## 题目描述

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

## 题目描述

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