Kyoya and Train

题目描述

Kyoya Ootori wants to take the train to get to school. There are $n$ train stations and $m$ one-way train lines going between various stations. Kyoya is currently at train station $1$, and the school is at station $n$. To take a train, he must pay for a ticket, and the train also takes a certain amount of time. However, the trains are not perfect and take random amounts of time to arrive at their destination. If Kyoya arrives at school strictly after $t$ time units, he will have to pay a fine of $x$.
Each train line is described by a ticket price, and a probability distribution on the time the train takes. More formally, train line $i$ has ticket cost $c_i$, and a probability distribution $p_{i, k}$ which denotes the probability that this train will take $k$ time units for all $1 \le k \le t$. Amounts of time that each of the trains used by Kyouya takes are mutually independent random values (moreover, if Kyoya travels along the same train more than once, it is possible for the train to take different amounts of time and those amounts are also independent one from another).
Kyoya wants to get to school by spending the least amount of money in expectation (for the ticket price plus possible fine for being late). Of course, Kyoya has an optimal plan for how to get to school, and every time he arrives at a train station, he may recalculate his plan based on how much time he has remaining. What is the expected cost that Kyoya will pay to get to school if he moves optimally?

题意概述

有$n$座城市和$m$条铁路,第$i$条铁路从$a_i$号城市出发,到达$b_i$号城市,票价为$c_i$。某人要在$t$个单位时间内从$1$号城市到$n$号城市,如果超过时间就会被罚款$x$。已知每次搭乘第$i$条铁路有$p_{i, k}$的概率需要$k \; (1 \le k \le t)$个单位时间。求在采取最优策略的情况下总花费的期望。
数据范围:$2 \le n \le 50, \; 1 \le m \le 100, \; 1 \le t \le 20000, \; 0 \le x \le 10^6$。

算法分析

首先考虑最暴力的DP。令$f_{i, j}$表示当前时刻为$j$,在最优策略下从$i$号城市到$n$号城市的总花费的期望。那么
$$
f_{i, j}=
\begin{cases}
\min(c_e+\sum_{k=1}^t p_{e, k} \cdot f_{b_e, j+k} \mid e \in [1, m] \land a_e=i), & i \lt n \land j \le t \\
d_{i, n}+x , & i \lt n \land j \gt t \\
0, & i=n \land j \le t \\
x, & i=n \land j \gt t
\end{cases}
$$
第二部分中的$d_{i, n}$表示在只考虑铁路票价的情况下从$i$号城市到$n$号城市的最小花费(因为已经超时了,所以不如走总票价最少的路)。
这样做的复杂度是$O(mt^2)$的,需要对它进行优化。令
$$S_{e, j}=\sum_{k=1}^t p_{e, k} \cdot f_{b_e, j+k}$$
转移方程的第一部分变为
$$f_{i, j}=\min(c_e+S_{e, j} \mid e \in [1, m] \land a_e=i)$$
对于$S$,它类似于卷积,可以将其中一部分翻转后用FFT求值;对于$f$,可以用CDQ分治,统计较大的$j$对较小的$j$的贡献。
具体来说,假设我们要计算$l \le j \le r$的$f_{i, j}$,令$mid=\lfloor {l+r \over 2} \rfloor$,那么先计算$mid \lt j \le r$的$f_{i, j}$,用这些来更新$l \le j \le mid$的$S_{e, j}$($m$条边分别用FFT更新,复杂度$O(m(r-l)\log (r-l))$),最后计算$l \le j \le mid$的$f_{i, j}$。当$l=r$时$f$可以直接由$S$得到。总复杂度降至$O(mt\log^2t)$。
实现时细节很多,要注意DP边界条件和FFT下标变换。

代码

/*
 * Your nature demands love and your happiness depends on it.
 */

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>

static int const N = 55;
static int const M = 105;
static int const T = 100000;
static double const PI = acos(-1);
int rev[T];

class Point {
public:
  double x, y;
  Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
  std::pair<double, double> pair() {
    return std::make_pair(x, y);
  }
  friend Point operator+(Point const &a, Point const &b) {
    return Point(a.x + b.x, a.y + b.y);
  }
  friend Point operator-(Point const &a, Point const &b) {
    return Point(a.x - b.x, a.y - b.y);
  }
  friend Point operator*(Point const &a, Point const &b) {
    return Point(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
  }
  Point operator/(double const &n) {
    return Point(x / n, y / n);
  }
} wn[T], A[T], B[T];

int init(int n) {
  int m = n, l = 0;
  for (n = 1; n <= m; n <<= 1, ++l)
    ;
  for (int i = 1; i < n; ++i)
    rev[i] = rev[i >> 1] >> 1 | (i & 1) << l - 1;
  for (int i = 0; i < n >> 1; ++i)
    wn[i] = Point(cos(2 * PI / n * i), sin(2 * PI / n * i));
  return n;
}

void fft(Point *a, int n, int inv = 0) {
  for (int i = 0; i < n; ++i)
    if (i < rev[i])
      std::swap(a[i], a[rev[i]]);
  for (int i = 1; i < n; i <<= 1)
    for (int j = 0; j < n; j += i << 1)
      for (int k = 0; k < i; ++k) {
        Point x = a[j + k], y = wn[n / (i << 1) * k] * a[j + k + i];
        a[j + k] = x + y, a[j + k + i] = x - y;
      }
  if (inv) {
    std::reverse(a + 1, a + n);
    for (int i = 0; i < n; ++i)
      a[i] = a[i] / n;
  }
}

int n, m, t, x, mp[N][N];
double p[M][T], s[M][T], f[N][T];
struct Line {
  int a, b, c;
} li[M];

void update(int l, int r) {
  int mid = l + r >> 1, len = init(r - l + r - mid - 2);
  for (int i = 1; i <= m; ++i) {
    for (int j = 0; j < len; ++j)
      A[j] = B[j] = 0;
    for (int j = mid + 1; j <= r; ++j)
      A[j - mid - 1] = f[li[i].b][r - j + mid + 1];
    for (int j = 1; j <= r - l; ++j)
      B[j - 1] = p[i][j];
    fft(A, len), fft(B, len);
    for (int j = 0; j < len; ++j)
      A[j] = A[j] * B[j];
    fft(A, len, 1);
    for (int j = l; j <= mid; ++j)
      s[i][j] += A[r - j - 1].x;
  }
}

void solve(int l, int r) {
  if (l == r) {
    for (int i = 1; i <= m; ++i)
      f[li[i].a][l] = std::min(f[li[i].a][l], s[i][l] + li[i].c);
    return;
  }
  int mid = l + r >> 1;
  solve(mid + 1, r), update(l, r), solve(l, mid);
}

int main() {
  scanf("%d%d%d%d", &n, &m, &t, &x);
  memset(mp, 0x3f, sizeof mp);
  for (int i = 1; i <= m; ++i) {
    scanf("%d%d%d", &li[i].a, &li[i].b, &li[i].c);
    mp[li[i].a][li[i].b] = std::min(mp[li[i].a][li[i].b], li[i].c);
    for (int j = 1; j <= t; ++j)
      scanf("%lf", &p[i][j]), p[i][j] /= 100000;
  }
  for (int i = 1; i <= n; ++i)
    for (int j = 1; j <= n; ++j)
      for (int k = 1; k <= n; ++k)
        mp[j][k] = std::min(mp[j][k], mp[j][i] + mp[i][k]);
  for (int i = 1; i <= n; ++i)
    for (int j = 0; j <= t << 1; ++j)
      if (i == n)
        if (j <= t)
          f[i][j] = 0;
        else
          f[i][j] = x;
      else
        if (j <= t)
          f[i][j] = 1e9;
        else
          f[i][j] = x + mp[i][n];
  update(0, t << 1), solve(0, t), printf("%.8lf\n", f[1][0]);
  return 0;
}

The Child and Binary Tree

题目描述

Our child likes computer science very much, especially he likes binary trees.
Consider the sequence of $n$ distinct positive integers: $c_1, c_2, \ldots, c_n$. The child calls a vertex-weighted rooted binary tree good if and only if for every vertex $v$, the weight of $v$ is in the set ${c_1, c_2, \ldots, c_n}$. Also our child thinks that the weight of a vertex-weighted tree is the sum of all vertices’ weights.
Given an integer $m$, can you for all $s \; (1 \le s \le m)$ calculate the number of good vertex-weighted rooted binary trees with weight $s$? Please, check the samples for better understanding what trees are considered different.
We only want to know the answer modulo $998244353$ ($7 \times 17 \times 2^{23}+1$, a prime number).

题意概述

有$n$种点,每种点有无限个,第$i$种点的权值为$c_i$。定义一棵二叉树的权值等于它所有点的权值之和。求对于所有$s \in [1, m]$,权值为$s$的二叉树有几棵。两棵二叉树不同当且仅当它们左子树或右子树不同,或者根节点权值不同。
数据范围:$1 \le n, m, c_i \le 10^5$。

算法分析

令$f(x)$表示权值为$x$的二叉树个数,$F(x)$为其生成函数($F(x)=\sum_{i \ge 0} f(i)x^i$)。
令$C(x)$为给定$c$的集合的生成函数($C(x)=\sum_{i=1}^n x^{c_i}$)。
根据DP转移方程,易知
$$
f(x)=\sum_{w \in {c_1, c_2, \ldots, c_n}} \sum_{i=0}^{x-w} f(i)f(x-w-i)
$$

$$
F(x)=C(x)F(x)^2+1
$$
解得
$$
F(x)={1 \pm \sqrt{1-4C(x)} \over 2C(x)}={2 \over 1 \pm \sqrt{1-4C(x)}}
$$
显然,若取减号,则当$x$趋近$0$时分母为$0$,因此只能取加号。接着就是多项式开根和多项式求逆了。

  • 多项式求逆:
    求$GF \equiv 1 \pmod {x^n}$。

    假设已知$G_0F \equiv 1 \pmod {x^{\lceil n/2 \rceil}}$
    $G-G_0 \equiv 0 \pmod {x^{\lceil n/2 \rceil}}$
    $G^2-2GG_0+G_0^2 \equiv 0 \pmod {x^n}$
    $G-2G_0+G_0^2F \equiv 0 \pmod {x^n}$
    $G \equiv 2G_0-G_0^2F \pmod {x^n}$

  • 多项式开根:
    求$G^2 \equiv F \pmod {x^n}$。

    假设已知$G_0^2 \equiv F \pmod {x^{\lceil n/2 \rceil}}$
    $(G_0^2-F)^2 \equiv 0 \pmod {x^n}$
    $(G_0^2+F)^2 \equiv 4G_0^2F \pmod {x^n}$
    $\left({G_0^2+F \over 2G_0}\right)^2 \equiv F \pmod {x^n}$
    $G \equiv {G_0+G_0^{-1}F \over 2} \pmod {x^n}$

代码

#include <algorithm>
#include <cstdio>
#include <cstring>

static const int N = 500000;
static const int MOD = 998244353;
static const int G = 3;
static const int INV2 = 499122177;
int n, m, c[N], C[N], rev[N], wn[N], tmp[N], tmp2[N], tmp3[N];

int power(int a, int b) {
  int ret = 1;
  for (a %= MOD, b %= MOD - 1; b; b >>= 1)
    b & 1 && (ret = 1ll * ret * a % MOD), a = 1ll * a * a % MOD;
  return ret;
}

void init(int &n) {
  int m = n << 1, l = 0;
  for (n = 1; n < m; n <<= 1, ++l)
    ;
  for (int i = 1; i < n; ++i)
    rev[i] = rev[i >> 1] >> 1 | (i & 1) << (l - 1);
}

void ntt(int *a, int n, bool inv) {
  for (int i = 0; i < n; ++i)
    if (i < rev[i])
      std::swap(a[i], a[rev[i]]);
  wn[0] = 1, wn[1] = power(G, (MOD - 1) / n);
  for (int i = 2; i < n >> 1; ++i)
    wn[i] = 1ll * wn[i - 1] * wn[1] % MOD;
  for (int i = 1; i < n; i <<= 1)
    for (int j = 0; j < n; j += i << 1)
      for (int k = 0; k < i; ++k) {
        int x = a[j + k], y = 1ll * wn[n / (i << 1) * k] * a[j + k + i] % MOD;
        a[j + k] = (x + y) % MOD, a[j + k + i] = (MOD + x - y) % MOD;
      }
  if (inv) {
    for (int i = 1; i < n >> 1; ++i)
      std::swap(a[i], a[n - i]);
    int rec = power(n, MOD - 2);
    for (int i = 0; i < n; ++i)
      a[i] = 1ll * a[i] * rec % MOD;
  }
}

void get_inv(int *f, int *g, int n) {
  if (n == 1)
    return void(g[0] = power(f[0], MOD - 2));
  int rec = n;
  get_inv(f, g, (n + 1) >> 1), init(n);
  for (int i = (rec + 1) >> 1; i < n; ++i)
    g[i] = 0;
  for (int i = 0; i < rec; ++i)
    tmp[i] = f[i];
  for (int i = rec; i < n; ++i)
    tmp[i] = 0;
  ntt(g, n, 0), ntt(tmp, n, 0);
  for (int i = 0; i < n; ++i)
    g[i] = 1ll * g[i] * (MOD + 2 - 1ll * g[i] * tmp[i] % MOD) % MOD;
  ntt(g, n, 1);
  for (int i = rec; i < n; ++i)
    g[i] = 0;
}

void get_sqrt(int *f, int *g, int n) {
  if (n == 1)
    return void(g[0] = 1);
  int rec = n;
  get_sqrt(f, g, (n + 1) >> 1);
  for (int i = 0; i<(n + 1)>> 1; ++i)
    tmp2[i] = g[i];
  for (int i = (n + 1) >> 1; i < n; ++i)
    tmp2[i] = 0;
  get_inv(tmp2, tmp3, n), init(n);
  for (int i = (rec + 1) >> 1; i < n; ++i)
    g[i] = 0;
  for (int i = 0; i < rec; ++i)
    tmp2[i] = f[i];
  for (int i = rec; i < n; ++i)
    tmp2[i] = 0;
  ntt(tmp2, n, 0), ntt(tmp3, n, 0);
  for (int i = 0; i < n; ++i)
    tmp3[i] = 1ll * tmp3[i] * tmp2[i] % MOD;
  ntt(tmp3, n, 1);
  for (int i = 0; i < rec; ++i)
    g[i] = 1ll * (g[i] + tmp3[i]) * INV2 % MOD;
  for (int i = rec; i < n; ++i)
    g[i] = 0;
}

int main() {
  scanf("%d%d", &n, &m);
  for (int i = 0; i < n; ++i)
    scanf("%d", &c[i]);
  for (int i = 0; i < n; ++i)
    if (c[i] <= m)
      ++C[c[i]];
  for (int i = 1; i <= m; ++i)
    C[i] = (MOD - (C[i] << 2)) % MOD;
  C[0] = 1;
  get_sqrt(C, c, m + 1), ++c[0], get_inv(c, C, m + 1);
  for (int i = 1; i <= m; ++i)
    printf("%d\n", (C[i] << 1) % MOD);
  return 0;
}

序列求和 V4

题意概述

给定$n$和$k$,求$\sum_{i=1}^n i^k \bmod 10^9+7$。有$T$组数据。
数据范围:$1 \le T \le 500, \; 1 \le n \le 10^{18}, \; 1 \le k \le 50000$。

算法分析1

设答案为$f(n)$。这是一个$(k+1)$次多项式,只要确定$(k+2)$个点就可以确定$f$。
我们可以令$x=1\ldots k+2$,在$O(k\log k)$时间内计算出它们对应的$y$。
利用拉格朗日插值法。假设有$n$个点$(x_i, y_i)$,那么
$$f(x)=\sum_{i=1}^n {\prod_{j \in [1, n] \land j \neq i} (x-x_j) \over \prod_{j \in [1, n] \land j \neq i} (x_i-x_j)}y_i$$
考虑分子部分。在给定$x$的情况下,这一部分可以$O(k)$预处理$O(1)$求值。
考虑分母部分。它形如$\prod_{j \in [i-k-2, i-1] \land j \neq 0} j$,可以在插值时$O(1)$维护。
总时间复杂度为$O(Tk\log k)$。

代码1

#include <cstdio>

#define int long long

void read(int &n) {
  char c; while ((c = getchar()) < '0' || c > '9') ; n = c - '0';
  while ((c = getchar()) >= '0' && c <= '9') (n *= 10) += c - '0';
}

static const int K = 50005;
static const int MOD = 1000000007;
int n, k, a[K], p[K], q[K], inv[K], base[K];

int power(int a, int b) {
  int ret = 1; a %= MOD, b %= MOD - 1;
  while (b) b & 1 && ((ret *= a) %= MOD), (a *= a) %= MOD, b >>= 1;
  return ret;
}

int calc(int n) {
  if (n <= k + 2) return a[n]; n %= MOD;
  int w = power(base[k + 2], MOD - 2), ans = 0;
  p[0] = q[k + 3] = 1;
  for (int i = 1; i <= k + 2; ++ i) p[i] = p[i - 1] * (n - i) % MOD;
  for (int i = k + 2; i; -- i) q[i] = q[i + 1] * (n - i) % MOD;
  for (int i = 1; i <= k + 2; ++ i) (ans += a[i] * w % MOD * p[i - 1] % MOD * q[i + 1]) %= MOD, w = w * (i - k - 2) % MOD * inv[i] % MOD;
  return (ans + MOD) % MOD;
}

signed main() {
  int T; read(T), base[1] = 1;
  for (int i = 1; i < K; ++ i) inv[i] = power(i, MOD - 2);
  for (int i = 2; i < K; ++ i) base[i] = base[i - 1] * (MOD + 1 - i) % MOD;
  while (T --) {
    read(n), read(k);
    for (int i = 1; i < k + 3; ++ i) a[i] = (a[i - 1] + power(i, k)) % MOD;
    printf("%lld\n", calc(n));
  }
  return 0;
}

算法分析2

此题也可以用伯努利数解决,因为
$$\sum_{i=1}^n i^k={1 \over k+1} \sum_{i=0}^k (-1)^i {k+1 \choose i} B_in^{k+1-i}$$
其中伯努利数由其指数型生成函数${x \over e^x-1}$定义。易知
$$\sum_{i=0}^{\infty} B_i{x^i \over i!}={x \over e^x-1}={x \over \sum_{i=1}^{\infty} {x^i \over i!}}={1 \over \sum_{i=0}^{\infty} {x^i \over (i+1)!}}$$
只要对分母求多项式逆元,即可得到伯努利数,时间复杂度$O(k\log k)$。其它部分均可$O(k)$预处理,答案也可以$O(k)$计算。因此总时间复杂度是$O(k\log k+Tk)$。

代码2

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>

#define int long long

static int const N = 200000;
static int const MOD = 1000000007;
static int const M = 31622;
static double const PI = acos(-1);
int b[N], rev[N], fac[N], inv[N];
class complex {
private:
  double x, y;

public:
  complex(double _x = 0, double _y = 0) : x(_x), y(_y) {}
  double real() { return x; }
  friend complex operator+(complex const &a, complex const &b) {
    return complex(a.x + b.x, a.y + b.y);
  }
  friend complex operator-(complex const &a, complex const &b) {
    return complex(a.x - b.x, a.y - b.y);
  }
  friend complex operator*(complex const &a, complex const &b) {
    return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
  }
} wn[N], A[N], B[N], C[N], D[N], E[N], F[N], G[N];

int power(int a, int b) {
  int ret = 1;
  for (; b; b >>= 1)
    b & 1 && ((ret *= a) %= MOD), (a *= a) %= MOD;
  return ret;
}

int init(int m) {
  int n = 1, l = 0;
  for (; n <= m; n <<= 1, ++l)
    ;
  for (int i = 1; i < n; ++i)
    rev[i] = rev[i >> 1] >> 1 | (i & 1) << l - 1;
  for (int i = 0; i < n >> 1; ++i)
    wn[i] = complex(cos(2 * PI / n * i), sin(2 * PI / n * i));
  return n;
}

void fft(complex *a, int n, bool inv = 0) {
  for (int i = 0; i < n; ++i)
    if (i < rev[i])
      std::swap(a[i], a[rev[i]]);
  for (int i = 1; i < n; i <<= 1)
    for (int j = 0; j < n; j += i << 1)
      for (int k = 0; k < i; ++k) {
        complex x = a[j + k], y = wn[n / (i << 1) * k] * a[j + k + i];
        a[j + k] = x + y, a[j + k + i] = x - y;
      }
  if (inv) {
    std::reverse(a + 1, a + n);
    for (int i = 0; i < n; ++i)
      a[i] = a[i].real() / n;
  }
}

int get_mod(double t) {
  return (int)(fmod(t, MOD) + MOD + 0.5) % MOD;
}

void get_inv(int *f, int *g, int n) {
  if (n == 1)
    return void(g[0] = power(f[0], MOD - 2));
  get_inv(f, g, n + 1 >> 1);
  int len = init(n << 1);
  for (int i = 0; i < n + 1 >> 1; ++i)
    A[i] = g[i] / M, B[i] = g[i] % M;
  for (int i = n + 1 >> 1; i < len; ++i)
    A[i] = B[i] = 0;
  for (int i = 0; i < n; ++i)
    C[i] = f[i] / M, D[i] = f[i] % M;
  for (int i = n; i < len; ++i)
    C[i] = D[i] = 0;
  fft(A, len), fft(B, len), fft(C, len), fft(D, len);
  for (int i = 0; i < len; ++i) {
    E[i] = 0 - A[i] * C[i];
    F[i] = 0 - A[i] * D[i] - B[i] * C[i];
    G[i] = 2 - B[i] * D[i];
  }
  fft(E, len, 1), fft(F, len, 1), fft(G, len, 1);
  for (int i = 0; i < n; ++i) {
    int x = get_mod(E[i].real()) * M % MOD * M % MOD;
    int y = get_mod(F[i].real()) * M % MOD;
    int z = get_mod(G[i].real());
    int w = (x + y + z) % MOD;
    C[i] = w / M, D[i] = w % M;
  }
  for (int i = n; i < len; ++i)
    C[i] = D[i] = 0;
  fft(C, len), fft(D, len);
  for (int i = 0; i < len; ++i) {
    E[i] = A[i] * C[i];
    F[i] = A[i] * D[i] + B[i] * C[i];
    G[i] = B[i] * D[i];
  }
  fft(E, len, 1), fft(F, len, 1), fft(G, len, 1);
  for (int i = 0; i < n; ++i) {
    int x = get_mod(E[i].real()) * M % MOD * M % MOD;
    int y = get_mod(F[i].real()) * M % MOD;
    int z = get_mod(G[i].real());
    g[i] = (x + y + z) % MOD;
  }
}

int get_c(int n, int m) {
  return fac[n] * inv[m] % MOD * inv[n - m] % MOD;
}

signed main() {
  int T;
  fac[0] = 1;
  for (int i = 1; i < N; ++i)
    fac[i] = fac[i - 1] * i % MOD;
  inv[N - 1] = power(fac[N - 1], MOD - 2);
  for (int i = N - 1; i; --i)
    inv[i - 1] = inv[i] * i % MOD;
  get_inv(inv + 1, b, 50001);
  for (int i = 0; i < 50001; ++i)
    (b[i] *= fac[i]) %= MOD;
  for (scanf("%lld", &T); T--;) {
    int n, k, ans = 0;
    scanf("%lld%lld", &n, &k), n %= MOD;
    for (int i = k, f = k & 1 ? -1 : 1, p = n; ~i; --i, f = -f, (p *= n) %= MOD)
      ans += (MOD + f) * get_c(k + 1, i) % MOD * b[i] % MOD * p % MOD;
    printf("%lld\n", ans % MOD * power(k + 1, MOD - 2) % MOD);
  }
  return 0;
}