## 题目描述

Given two sequences $a$ and $b$ of equal length $n$, find the

$$\sum_{x=1}^n \sum_{y=1}^n \sum_{z=1}^n \sum_{w=1}^n (a_x+a_y+a_z+a_w)^{(b_x \oplus b_y \oplus b_z \oplus b_w)}$$

## 算法分析

$a_x+a_y+a_z+a_w$的取值范围为$[4,2000]$，$b_x \oplus b_y \oplus b_z \oplus b_w$的取值范围为$[0,511]$，可以分别计算每种情况出现的次数再求和。

$$f_{k,i,j}=\sum_{i_1+i_2=i} \sum_{j_1 \oplus j_2=j} f_{k-1,i_1,j_1}f_{1,i_2,j_2}$$

$f_{4,i,j}$即要求的每种情况出现的次数。转移方程在一维上是乘法卷积，另一维上是异或卷积，可以分别用FFT和FWT进行处理。总时间复杂度为$O(a_{\max}b_{\max}\log(a_{\max}b_{\max}))$。

## 代码

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

int const N = 100005, M = 512, MOD = 998244353, G = 3, INV2 = 499122177;

int a[N], b[N], rev[M << 2];

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

int wn[M << 2], A[M << 2], f[M << 2][M];

void init(int n) {
int len = 1, p = 0;
for (; len < n; len <<= 1, ++p) ;
for (int i = 1; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1 | (i & 1) << p - 1;
}
wn[0] = 1, wn[1] = power(G, (MOD - 1) / len);
for (int i = 2; i < len >> 1; ++i) {
wn[i] = 1ll * wn[i - 1] * wn[1] % MOD;
}
}

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

void fwt(int *a, int len, bool inv = 0) {
for (int i = 1; i < len; i <<= 1) {
for (int j = 0; j < len; j += i << 1) {
for (int k = 0; k < i; ++k) {
int x = a[j + k], y = a[j + i + k];
a[j + k] = (x + y) % MOD;
a[j + i + k] = (MOD + x - y) % MOD;
if (inv) {
a[j + k] = 1ll * a[j + k] * INV2 % MOD;
a[j + i + k] = 1ll * a[j + i + k] * INV2 % MOD;
}
}
}
}
}

int main() {
int n;
scanf("%d", &n);
for (int i = 1; i <= n; ++i) {
scanf("%d", &a[i]);
}
for (int i = 1; i <= n; ++i) {
scanf("%d", &b[i]);
f[a[i]][b[i]] = f[a[i]][b[i]] + 1;
}
init(M << 2);
for (int i = 0; i < M << 2; ++i) {
for (int j = 0; j < M; ++j) {
A[j] = f[i][j];
}
fwt(A, M);
for (int j = 0; j < M; ++j) {
f[i][j] = A[j];
}
}
for (int i = 0; i < M; ++i) {
for (int j = 0; j < M << 2; ++j) {
A[j] = f[j][i];
}
fft(A, M << 2);
for (int j = 0; j < M << 2; ++j) {
f[j][i] = A[j];
}
}
for (int i = 0; i < M << 2; ++i) {
for (int j = 0; j < M; ++j) {
f[i][j] = power(f[i][j], 4);
}
}
for (int i = 0; i < M << 2; ++i) {
for (int j = 0; j < M; ++j) {
A[j] = f[i][j];
}
fwt(A, M, 1);
for (int j = 0; j < M; ++j) {
f[i][j] = A[j];
}
}
for (int i = 0; i < M; ++i) {
for (int j = 0; j < M << 2; ++j) {
A[j] = f[j][i];
}
fft(A, M << 2, 1);
for (int j = 0; j < M << 2; ++j) {
f[j][i] = A[j];
}
}
int ans = 0;
for (int i = 0; i < M << 2; ++i) {
for (int j = 0; j < M; ++j) {
if (f[i][j]) {
ans = (ans + 1ll * f[i][j] * power(i, j)) % MOD;
}
}
}
printf("%d\n", ans);
return 0;
}

## 题目描述

You have an unbalanced coin, the state of which (head or tail) after tossing depends on the previous state. If the head side is up now, the probability that the coin is still head up after tossing is $p$, and the probability that the tail side comes up is $1-p$. Similarly, if the tail side is up now, the probability that the coin is still tail up after tossing is $p$, and the probability that the head side comes up is $1-p$.

Assume the coin is head up initially. You will then toss it for $n$ times, during these $n$ times, every time the head side is up, you will get one score. What is the probability that you will get exactly $k$ score for $k=0,1,2,\cdots,n$?

## 算法分析

\begin{align} f_{i,0,j} &= f_{i-1,0,j-1} \times p+f_{i-1,1,j-1} \times (1-p) \\ f_{i,1,j} &= f_{i-1,0,j} \times (1-p)+f_{i-1,1,j} \times p \end{align}

\begin{align} F_{i,0} &= F_{i-1,0} \times px+F_{i-1,1} \times (1-p)x \\ F_{i,1} &= F_{i-1,0} \times (1-p)+F_{i-1,1} \times p \end{align}

$$$$\left[\begin{matrix} px & (1-p)x \\ 1-p & p\end{matrix}\right] \left[\begin{matrix} F_{i-1,0} \\ F_{i-1,1} \end{matrix}\right]= \left[\begin{matrix} F_{i,0} \\ F_{i,1} \end{matrix}\right]$$$$

## 代码

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

int const N = 200005;
double const PI = acos(-1);

int rev[N];

struct Complex {
double a, b;
Complex(double _a = 0, double _b = 0) : a(_a), b(_b) {}
friend Complex operator+(Complex const &x, Complex const &y) {
return Complex(x.a + y.a, x.b + y.b);
}
friend Complex operator-(Complex const &x, Complex const &y) {
return Complex(x.a - y.a, x.b - y.b);
}
friend Complex operator*(Complex const &x, Complex const &y) {
return Complex(x.a * y.a - x.b * y.b, x.a * y.b + x.b * y.a);
}
friend Complex operator/(Complex const &x, double const &y) {
return Complex(x.a / y, x.b / y);
}
} wn[N], A[2][2][N], B[2][2][N], C[2][2][N];

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

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

struct Matrix {
int n[2][2];
double *a[2][2];

friend Matrix operator*(Matrix const &x, Matrix const &y) {
int len = init(x.n[0][0] + y.n[0][0] - 1);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
for (int k = 0; k < len; ++k) {
A[i][j][k] = k < x.n[i][j] ? x.a[i][j][k] : 0;
B[i][j][k] = k < y.n[i][j] ? y.a[i][j][k] : 0;
C[i][j][k] = 0;
}
fft(A[i][j], len);
fft(B[i][j], len);
}
}
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < len; ++l) {
C[i][j][l] = C[i][j][l] + A[i][k][l] * B[k][j][l];
}
}
fft(C[i][j], len, 1);
}
}
Matrix ret;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
ret.n[i][j] = 0;
for (int k = 0; k < 2; ++k) {
ret.n[i][j] = std::max(ret.n[i][j], x.n[i][k] + y.n[k][j] - 1);
}
ret.a[i][j] = new double[ret.n[i][j]];
for (int k = 0; k < ret.n[i][j]; ++k) {
ret.a[i][j][k] = C[i][j][k].a;
}
}
}
return ret;
}
} base, ans;

int main() {
int n;
double p;
scanf("%d%lf", &n, &p);
base.n[0][0] = base.n[0][1] = 2;
base.n[1][0] = base.n[1][1] = 1;
base.a[0][0] = new double[2];
base.a[0][0][0] = 0;
base.a[0][0][1] = p;
base.a[0][1] = new double[2];
base.a[0][1][0] = 0;
base.a[0][1][1] = 1 - p;
base.a[1][0] = new double;
base.a[1][0][0] = 1 - p;
base.a[1][1] = new double;
base.a[1][1][0] = p;
ans.n[0][0] = ans.n[1][1] = 1;
ans.a[0][0] = new double;
ans.a[0][0][0] = 1;
ans.a[1][1] = new double;
ans.a[1][1][0] = 1;
int b = n;
for (; b;) {
if (b & 1) {
ans = ans * base;
}
if (b >>= 1) {
base = base * base;
}
}
for (int i = 0; i <= n; ++i) {
printf("%.10f\n", ans.a[0][0][i] + ans.a[1][0][i]);
}
return 0;
}

## 题目描述

Tsr is a cute boy with handsome moustache.

You are given a sequence with length $n$. Tsr wants you to calculate the sum of variance of each successive subsequence. Note: The variance in this problem should’t be divided by length.

Recall $\overline{a}_{l, r}=\frac{1}{r-l+1} \sum_{i=l}^r a_i$. Then you are supposed to calculate $\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2$.

## 算法分析

\begin{align} \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2 &= \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k^2-2a_k \overline{a}_{i, j}+\overline{a}_{i, j}^2) \\ &= \sum_{k=1}^n \sum_{i=1}^k \sum_{j=k}^n a_k^2-2\sum_{i=1}^n \sum_{j=i}^n \overline{a}_{i, j} \sum_{k=i}^j a_k+ \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 \\ &= \sum_{k=1}^n k(n-k+1)a_k^2-\sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 \end{align}

\begin{align} \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 &= \sum_{j=1}^n \sum_{i=1}^j \frac{(s_j-s_{i-1})^2}{j-i+1} \\ &= \sum_{j=1}^n \sum_{i=1}^j \frac{s_j^2-2s_js_{i-1}+s_{i-1}^2}{j-i+1} \\ &= \sum_{j=1}^n s_j^2 \sum_{i=1}^j \frac{1}{i}-2\sum_{j=1}^n s_j \sum_{i=1}^j \frac{s_{i-1}}{j-(i-1)}+\sum_{j=1}^n \sum_{i=1}^j \frac{s_{i-1}^2}{j-(i-1)} \end{align}

## 代码

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

int const N = 10005, T = 400005;
double const PI = acos(-1);

int a[N], rev[T];

struct Point {
double x, y;
Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
friend Point const operator + (Point const &a, Point const &b) {
return Point(a.x + b.x, a.y + b.y);
}
friend Point const operator - (Point const &a, Point const &b) {
return Point(a.x - b.x, a.y - b.y);
}
friend Point const 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);
}
friend Point const operator / (Point const &a, double const &n) {
return Point(a.x / n, a.y / n);
}
} wn[T], A[T], B[T], C[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 main() {
int T;
scanf("%d", &T);
for (; T--;) {
int n;
scanf("%d", &n);
double ans = 0;
for (int i = 1; i <= n; ++i) {
scanf("%d", &a[i]);
ans += 1. * a[i] * a[i] * i * (n - i + 1);
a[i] += a[i - 1];
}
double sum = 0;
for (int i = 1; i <= n; ++i) {
sum += 1. / i;
ans -= sum * a[i] * a[i];
}
int len = init(n << 1);
for (int i = 0; i < n; ++i) {
A[i] = Point(a[i]);
B[i] = Point(1. * a[i] * a[i]);
C[i] = Point(1. / (i + 1));
}
for (int i = n; i < len; ++i) {
A[i] = B[i] = C[i] = Point(0);
}
fft(A, len);
fft(B, len);
fft(C, len);
for (int i = 0; i < len; ++i) {
A[i] = A[i] * C[i];
B[i] = B[i] * C[i];
}
fft(A, len, 1);
fft(B, len, 1);
for (int i = 0; i < n; ++i) {
ans += 2 * a[i + 1] * A[i].x - B[i].x;
}
printf("%.10f\n", ans);
}
return 0;
}


## 题目描述

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?

## 算法分析

$$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}$$

$$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)$$

## 代码

/*
*/

#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;
}


## 题目描述

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).

## 算法分析

$$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)}}$$

• 多项式求逆：
求$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;
}


## 算法分析1

$$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$$

## 代码1

#include <cstdio>

#define int long long

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 --) {
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}$$

$$\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)!}}$$

## 代码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;
}