ケロシのポリ - 多项式学习笔记

ケロシのポリ - 多项式学习笔记

持续更新中……

目录

• ケロシのポリ - 多项式学习笔记
  • 多项式乘法
  • 多项式乘法逆
  • 多项式开根
  • 多项式在 $\ln$ 和 $\exp$ 上的意义
  • 多项式对数函数 $\ln$
  • 多项式牛顿迭代
  • 多项式指数函数 $\exp$
  • 多项式快速幂
  • 多项式连续点值平移
  • Chirp Z-Transform
  • 任意模数多项式乘法
  • 快速阶乘
  • 第二类斯特林数行
  • 第二类斯特林数列
  • 第一类斯特林数行
  • 第一类斯特林数列

多项式乘法

给出 $F(x)$ 和 $G(x)$，求出 $H(x)=F(x)G(x)$。

考虑用 NTT 求出点值，而点值可以直接相乘。

然后对乘出来的点值做一遍逆变换即可。

时间复杂度 $O(n \log n)$。

ケロシの代码

void solve() {
	cin >> n >> m;
	FOR(i, 0, n) cin >> a[i];
	FOR(i, 0, m) cin >> b[i];
	int lim = 1, l = 0;
	while(lim <= (n + m)) lim <<= 1, l++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	NTT(a, lim, 1);
	NTT(b, lim, 1);
	REP(i, lim) Mul(a[i], b[i]);
	NTT(a, lim, -1);
	FOR(i, 0, n + m) cout << a[i] << " ";
	cout << endl;
}

多项式乘法逆

给出 $F(x)$，求 $G(x)$ 使得 $F(x)G(x) \equiv 1 (\bmod x^n)$。

首先 $G(x)_0=\frac{1}{F(x)_0}$，然后考虑倍增，用 $\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }$ 的答案 $G_0(x)$ 推 $\bmod x^n$ 的答案：

$$\begin{aligned} & \because F(x)G_0(x) \equiv 1 (\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }) \\ & \because F(x)G(x) \equiv 1 (\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }) \\ & \therefore G(x)-G_0(x) \equiv 0 (\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }) \\ & \therefore (G(x)-G_0(x))^2 \equiv 0 (\bmod x^{n}) \\ & \therefore G(x)^2-2G(x)G_0(x)+G_0^2(x) \equiv 0 (\bmod x^{n}) \\ & \therefore G(x)-2G_0(x)+F(x)G_0^2(x) \equiv 0 (\bmod x^{n}) \\ & \therefore G(x) \equiv 2G_0(x) - F(x)G_0^2(x) (\bmod x^{n}) \\ \end{aligned}$$

需要用到 $O(n \log n)$ 的多项式乘法，计算一下复杂度：

$$T(n)=T(\left \lceil \frac{n}{2} \right \rceil )+O(n \log n)$$

所以复杂度 $T(n)=O(n \log n)$。

ケロシの代码

int c[N];
void invs(int * a, int * b, int d) {
	if(d == 1) {
		b[0] = fp(a[0], P - 2);
		return;
	}
	invs(a, b, (d + 1) >> 1);
	int lim = 1, l = 0;
	while(lim < (d << 1)) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	REP(i, d) c[i] = a[i];
	FOR(i, d, lim - 1) c[i] = 0;
	NTT(b, lim, 1);
	NTT(c, lim, 1);
	REP(i, lim) Mul(b[i], sub(2, mul(b[i], c[i])));
	NTT(b, lim, -1);
	FOR(i, d, lim - 1) b[i] = 0;
}

多项式开根

给出 $F(x)$，求 $G(x)$ 使得 $G^2(x) \equiv F(x) (\bmod x^n)$。

先考虑常数项，那么 $G_0(x)$ 就是 $F_0(x)$ 的一个二次剩余。

然后还是考虑倍增，用 $\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }$ 的答案 $G_0(x)$ 推 $\bmod x^n$ 的答案：

$$\begin{aligned} & \because G_0^2(x) \equiv F(x) (\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }) \\ & \therefore G_0^2(x)-F(x) \equiv 0 (\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }) \\ & \therefore (G_0^2(x)-F(x))^2 \equiv 0 (\bmod x^{n}) \\ & \therefore (G_0^2(x)+F(x))^2 \equiv 4G_0^2(x)F(x) (\bmod x^{n}) \\ & \therefore \left ( \frac{G_0^2(x)+F(x)}{2G_0(x)} \right )^2 \equiv F(x) (\bmod x^{n}) \\ & \therefore \frac{G_0^2(x)+F(x)}{2G_0(x)}\equiv G(x) (\bmod x^{n}) \\ & \therefore \frac{G_0(x)}{2}+\frac{F(x)}{2G_0(x)}\equiv G(x) (\bmod x^{n}) \\ \end{aligned}$$

时间复杂度 $O(n \log n)$。

多项式在 $\ln$ 和 $\exp$ 上的意义

一般来说，$\ln$ 和 $\exp$ 都是对于实数的，那么要计算多项式等形式的 $\ln$ 和 $\exp$ 需要用到泰勒公式进行展开：

$$\exp(x)=\frac{x^0}{0!}+\frac{x^1}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots$$

$$\ln(x+1)=\frac{x^1}{1}-\frac{x^2}{2}+\frac{x^3}{3}-\frac{x^4}{4}+\cdots$$

我们把式子中的 $x$ 替换成多项式即可使运算有意义。

为了使其不无限循环下去，式子中的多项式的 $x_0$ 常数项必须为 $0$，使得 $x^n$ 及以后 $\bmod x^n$ 都等于 $0$。

所以当多项式 $F(x)$ 满足 $F(x)_0=1$ 时 $\ln(F(x))$ 有意义，
当多项式 $F(x)$ 满足 $F(x)_0=0$ 时 $\exp(F(x))$ 有意义。

多项式对数函数 $\ln$

给出 $F(x)$，求 $G(x)$ 使得 $G(x) \equiv \ln(F(x)) (\bmod x^n)$。

对于 $\ln(x)$，注意到其求导后即为 $\frac{1}{x}$，比较好算，所以考虑两边求导。

对于复合函数 $f(g(x))$，其求导需要遵循链式法则：$f(g(x))'=f'(g(x))g'(x)$。

$$\begin{aligned} & \because G(x) \equiv \ln(F(x)) (\bmod x^n) \\ & \therefore G'(x) \equiv \ln(F(x))' (\bmod x^n) \\ & \therefore G'(x) \equiv \ln'(F(x))F'(x) (\bmod x^n) \\ & \therefore G'(x) \equiv \frac{F'(x)}{F(x)} (\bmod x^n) \\ & \therefore G(x) \equiv \int \frac{F'(x)}{F(x)} \mathrm{d} x(\bmod x^n) \\ \end{aligned}$$

需要多项式求逆和多项式乘法，时间复杂度 $O(n \log n)$。

ケロシの代码

int p[N], q[N];
void lns(int * a, int * b, int d) {
	int lim = 1, l = 0;
	while(lim < (d << 1)) lim <<= 1, l++;
	REP(i, lim) q[i] = 0;
	invs(a, q, d);
	FOR(i, 1, d - 1) p[i - 1] = mul(i, a[i]);
	FOR(i, d - 1, lim - 1) p[i] = 0;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	NTT(p, lim, 1);
	NTT(q, lim, 1);
	REP(i, lim) Mul(p[i], q[i]);
	NTT(p, lim, - 1);
	FOR(i, 1, d) b[i] = mul(p[i - 1], inv[i]);
	b[0] = 0;
}

多项式牛顿迭代

给出多项式 $F(x)$，求 $F(x)$ 的零点 $G(x)$，即满足 $F(G(x))\equiv 0(\bmod x^n)$。

考虑使用牛顿迭代来进行倍增。

先单独考虑常数项 $F(G(x)_0)_0 = 0$。

接下来使用牛顿迭代，用 $\bmod x^{ \left \lceil \frac{n}{2} \right \rceil }$ 的答案 $G_0(x)$ 推 $\bmod x^n$ 的答案。

考虑泰勒展开公式：

$$f(x)=\frac{f(x_0)(x-x_0)^0}{0!}+\frac{f'(x_0)(x-x_0)^1}{1!}+\frac{f''(x_0)(x-x_0)^2}{2!}+\cdots$$

考虑将泰勒展开公式中的 $x$ 和 $x_0$ 分别代入 $G(x)$ 和 $G_0(x)$，那么 $(G(x)-G_0(x))^2$ 及以后的项 $\bmod x^n$ 都等于 $0$，所以直接代入泰勒展开公式：

$$F(G(x)) \equiv F(G_0(x))+F'(G_0(x))(G(x)-G_0(x))(\bmod x^n)$$

因为 $F(G(x)) \equiv 0 (\bmod x^n)$

所以可以推出 $G(x)$：

$$\begin{aligned} & \because F(G(x)) \equiv F(G_0(x))+F'(G_0(x))(G(x)-G_0(x)) (\bmod x^n) \\ & \therefore 0 \equiv F(G_0(x))+F'(G_0(x))(G(x)-G_0(x)) (\bmod x^n) \\ & \therefore -F(G_0(x)) \equiv F'(G_0(x))(G(x)-G_0(x)) (\bmod x^n) \\ & \therefore -\frac{F(G_0(x))}{F'(G_0(x))} \equiv G(x)-G_0(x) (\bmod x^n) \\ & \therefore G(x) \equiv G_0(x)-\frac{F(G_0(x))}{F'(G_0(x))} (\bmod x^n) \\ \end{aligned}$$

多项式指数函数 $\exp$

给出 $F(x)$，求 $G(x)$ 使得 $G(x) \equiv \exp(F(x)) (\bmod x^n)$。

考虑将求 $exp$ 转换成关于 $\ln$ 的零点，

$$\begin{aligned} & \because G(x) \equiv \exp(F(x)) (\bmod x^n) \\ & \therefore \ln(G(x)) \equiv F(x) (\bmod x^n) \\ & \therefore \ln(G(x)) - F(x) \equiv 0 (\bmod x^n) \\ \end{aligned}$$

那么 $G(x)$ 就是函数 $H(x)=\ln(x)+F$ 的零点。

接下来使用牛顿迭代，代入公式。

但是需要求得 $H'(x)$，因为自变量 $x$ 是多项式，$F$ 也是多项式，所以对于 $x$ 来说 $F$ 就是常数项，所以 $H'(x)=\ln'(x)=\frac{1}{x}$

然后直接代入公式即可：

$$\begin{aligned} & \because G(x) \equiv G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))} (\bmod x^n) \\ & \therefore G(x) \equiv G_0(x)-(\ln(G_0(x))-F(x))G_0(x) (\bmod x^n) \\ & \therefore G(x) \equiv G_0(x)(1-\ln(G_0(x))+F(x)) (\bmod x^n) \\ \end{aligned}$$

需要用到多项式乘法，多项式乘法逆和多项式对数函数 $\ln$。

时间复杂度 $O(n \log n)$。

ケロシの代码

int h[N];
void exps(int * a, int * b, int d) {
	if(d == 1) {
		b[0] = 1;
		return;
	}
	exps(a, b, (d + 1) >> 1);
	REP(i, d) h[i] = 0;
	lns(b, h, d);
	REP(i, d) h[i] = sub(a[i], h[i]);
	Add(h[0], 1);
	int lim = 1, l = 0;
	while(lim < (d << 1)) lim <<= 1, l++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	FOR(i, d, lim - 1) h[i] = 0;
	NTT(b, lim, 1);
	NTT(h, lim, 1);
	REP(i, lim) Mul(b[i], h[i]);
	NTT(b, lim, - 1);
	FOR(i, d, lim - 1) b[i] = 0;
}

多项式快速幂

给出 $F(x)$ 和 $k$，求 $G(x)$ 使得 $G(x) \equiv F(x)^k (\bmod x^n)$，保证 $F(x)_0=1$。

考虑两边同时取对数下放指数：

$$\begin{aligned} & \because G(x) \equiv F(x)^k (\bmod x^n) \\ & \therefore \ln (G(x)) \equiv \ln (F(x)^k) (\bmod x^n) \\ & \therefore \ln (G(x)) \equiv k \ln (F(x)) (\bmod x^n) \\ & \therefore G(x) \equiv \exp ( k \ln (F(x)) ) (\bmod x^n) \\ \end{aligned}$$

时间复杂度 $O(n \log n)$。

多项式连续点值平移

给定不超过 $n$ 次多项式的 $n+1$ 个点值：$f(0),f(1),f(2),\cdots,f(n+1)$，并且给定 $m$，满足 $m>n$，求出 $f(m),f(m+1),f(m+2),\cdots,f(m+n)$。

考虑拉格朗日插值求出 $x>n$ 的点值：

$$\begin{aligned} f(x)&=\sum_{i=0}^{n} f(i) \prod_{j \ne i} \frac{x-j}{i-j} \\ &=\sum_{i=0}^{n} f(i) \frac{x!}{(x-i)(x-n-1)!i!(n-i)!(-1)^{n-i}} \end{aligned}$$

考虑卷积优化，把式子拆成关于 $i$，$x-i$ 和 $x$ 的部分：

$$\begin{aligned} f(x)=&\sum_{i=0}^{n} f(i) \frac{x!}{(x-i)(x-n-1)!i!(n-i)!(-1)^{n-i}} \\ =&\frac{x!}{(x-n-1)!} \sum_{i=0}^{n} \frac{f(i)}{i!(n-i)!(-1)^{n-i}} \frac{1}{x-i} \\ \end{aligned}$$

为了方便卷积，对齐下标：

$$f(m+x)=\frac{(m+x)!}{(m+x-n-1)!} \sum_{i=0}^{n} \frac{f(i)}{i!(n-i)!(-1)^{n-i}} \frac{1}{m+x-i} \\ $$

拆成两个函数，注意下标必须要正的，所以将 $x-i$ 改成 $x+n-i$ 即可：

$$F(i)=\frac{f(i)}{i!(n-i)!(-1)^{n-i}}$$

$$G(x+n-i)=\frac{1}{m+x-i}$$

$$G(i)=\frac{1}{m-n+i}$$

然后卷起来即可：

$$f(m+x)=\prod_{i=0}^{n}(m+x-i) \sum_{i=0}^{n}F(i)G(x+n-i)$$

时间复杂度 $O(n \log n)$。

ケロシの代码

int n, m, a[N], F[N], G[N];
int fac[N], finv[N];
void solve() {
	cin >> n >> m;
	FOR(i, 0, n) cin >> a[i];
	fac[0] = 1;
	FOR(i, 1, n) fac[i] = mul(fac[i - 1], i);
	finv[n] = fp(fac[n], P - 2);
	ROF(i, n - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
	FOR(i, 0, n) F[i] = mul({a[i], finv[i], finv[n - i], sgn(n - i)});
	FOR(i, 0, n << 1) G[i] = fp(m - n + i, P - 2);
	int lim = 1, l = 0;
	while(lim <= (n << 1)) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	NTT(F, lim, 1);
	NTT(G, lim, 1);
	REP(i, lim) Mul(F[i], G[i]);
	NTT(F, lim, - 1);
	int res = 1;
	FOR(i, 0, n) Mul(res, m - i);
	FOR(i, 0, n) {
		cout << mul(F[n + i], res) << " ";
		Mul(res, fp(m + i - n, P - 2));
		Mul(res, m + i + 1);
	}
	cout << endl;
}

Chirp Z-Transform

给定一个 $n$ 次多项式 $F(x)$，给定 $c$ 和 $m$，求出 $F(c^0),F(c^1),F(c^2),\cdots,F(c^{m-1})$。

考虑：

$$F(c^i)=\sum_{j=0}^{n-1}F_jc^{ij}$$

我们希望用卷积加速这个式子，但是有有关 $ij$ 的项，不好处理。

注意到 ${ij}$ 是在指数上，拆成加减后会变成乘法，比较好做，所以考虑拆解 $ij$：

$$\begin{aligned} 2ij&=(i+j)^2-i^2-j^2 \\ ij&=\frac{(i+j)^2}{2}-\frac{i^2}{2}-\frac{j^2}{2} \end{aligned}$$

但是这个式子的指数有除以二，需要二次剩余，不好做。

考虑另外一种拆法：

$$\begin{aligned} 2ij&=(i+j)^2-i^2-j^2 \\ &=(i+j)^2-i^2-j^2-(i+j)+i+j \\ &=(i+j)(i+j-1)-i(i-1)-j(j-1) \\ ij&=\frac{(i+j)(i+j-1)}{2}-\frac{i(i-1)}{2}-\frac{j(j-1)}{2} \\ &=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2} \end{aligned}$$

这样就都是整数了，代入式子：

$$\begin{aligned} F(c^i)&=\sum_{j=0}^{n-1}F_jc^{ij} \\ &=\sum_{j=0}^{n-1}F_jc^{\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2} } \\ &=c^{-\binom{i}{2}}\sum_{j=0}^{n-1}F_jc^{-\binom{j}{2}}c^{\binom{i+j}{2}} \\ \end{aligned}$$

现在就是一部分只和 $j$ 有关，另一部分之和 $i+j$ 有关，变成卷积形式了，直接卷积即可。

假设 $n$ 与 $m$ 同阶，时间复杂度 $O(n \log n)$。

ケロシの代码

int n, c, m;
int a[N], pw1[B + 1], pw2[B + 1];
int F[N], G[N];
inline int pw(int x) {
	return mul(pw1[x % B], pw2[x / B]);
}
inline int C2(int x) {
	return (1ll * x * (x - 1) / 2) % (P - 1);
}
void solve() {
	cin >> n >> c >> m;
	REP(i, n) cin >> a[i];
	pw1[0] = pw2[0] = 1; 
	FOR(i, 1, B) pw1[i] = mul(pw1[i - 1], c);
	FOR(i, 1, B) pw2[i] = mul(pw2[i - 1], pw1[B]);
	REP(i, n + m) F[i] = pw(C2(i));
	REP(i, n) G[i] = mul(a[i], pw(P - 1 - C2(i)));
	reverse(G, G + n);
	int lim = 1, l = 0;
	while(lim <= (n + m)) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	NTT(F, lim, 1);
	NTT(G, lim, 1);
	REP(i, lim) Mul(F[i], G[i]);
	NTT(F, lim, - 1);
	REP(i, m) cout << mul(F[i + n - 1], pw(P - 1 - C2(i))) << " "; cout << endl;
}

任意模数多项式乘法

能做 NTT 的模数需要大量单位根，所以不是所有的数都能做 NTT 的模数。

考虑使用三个原根都是 $3$ 的模数：$998244353,469762049,1004535809$。

对于每个模数都做多项式乘法，然后使用 CRT 把答案解出来即可。

时间复杂度 $O(n \log n)$。

ケロシの代码

const int N = 1e6 + 5;
const int P1 = 998244353;
const int P2 = 469762049;
const int P3 = 1004535809;
const int I1 = 208783132;
const int I2 = 395249030;
inline int add1(int x, int y) { return (x + y < P1 ? x + y : x + y - P1); }
inline void Add1(int & x, int y) { x = (x + y < P1 ? x + y : x + y - P1); }
inline int sub1(int x, int y) { return (x < y ? x - y + P1 : x - y); }
inline void Sub1(int & x, int y) { x = (x < y ? x - y + P1 : x - y); }
inline int mul1(int x, int y) { return (1ll * x * y) % P1; }
inline void Mul1(int & x, int y) { x = (1ll * x * y) % P1; }
inline int add2(int x, int y) { return (x + y < P2 ? x + y : x + y - P2); }
inline void Add2(int & x, int y) { x = (x + y < P2 ? x + y : x + y - P2); }
inline int sub2(int x, int y) { return (x < y ? x - y + P2 : x - y); }
inline void Sub2(int & x, int y) { x = (x < y ? x - y + P2 : x - y); }
inline int mul2(int x, int y) { return (1ll * x * y) % P2; }
inline void Mul2(int & x, int y) { x = (1ll * x * y) % P2; }
inline int add3(int x, int y) { return (x + y < P3 ? x + y : x + y - P3); }
inline void Add3(int & x, int y) { x = (x + y < P3 ? x + y : x + y - P3); }
inline int sub3(int x, int y) { return (x < y ? x - y + P3 : x - y); }
inline void Sub3(int & x, int y) { x = (x < y ? x - y + P3 : x - y); }
inline int mul3(int x, int y) { return (1ll * x * y) % P3; }
inline void Mul3(int & x, int y) { x = (1ll * x * y) % P3; }
int fp1(int x, int y) {
	int res = 1;
	for(; y; y >>= 1) {
		if(y & 1) Mul1(res, x);
		Mul1(x, x);
	}
	return res;
}
int fp2(int x, int y) {
	int res = 1;
	for(; y; y >>= 1) {
		if(y & 1) Mul2(res, x);
		Mul2(x, x);
	}
	return res;
}
int fp3(int x, int y) {
	int res = 1;
	for(; y; y >>= 1) {
		if(y & 1) Mul3(res, x);
		Mul3(x, x);
	}
	return res;
}
struct cint {
	int x1, x2, x3;
	cint operator + (const cint & A) {
		return {add1(x1, A.x1), add2(x2, A.x2), add3(x3, A.x3)};
	}
	cint operator - (const cint & A) {
		return {sub1(x1, A.x1), sub2(x2, A.x2), sub3(x3, A.x3)};
	}
	cint operator * (const cint & A) {
		return {mul1(x1, A.x1), mul2(x2, A.x2), mul3(x3, A.x3)};
	}
	cint & operator += (const cint & A) {
		Add1(x1, A.x1); Add2(x2, A.x2); Add3(x3, A.x3);
		return * this;
	}
	cint & operator -= (const cint & A) {
		Sub1(x1, A.x1); Sub2(x2, A.x2); Sub3(x3, A.x3);
		return * this;
	}
	cint & operator *= (const cint & A) {
		Mul1(x1, A.x1); Mul2(x2, A.x2); Mul3(x3, A.x3);
		return * this;
	}
	int crt(int P) {
		ll x4 = 1ll * P1 * mul2(sub2(x2, x1 % P2), I1) + x1;
		int k4 = mul3(sub3(x3, x4 % P3), I2);
		return (x4 + 1ll * k4 * (1ll * P1 * P2 % P)) % P;
	}
};
int rev[N];
void NTT(cint * a, int lim, int o) {
	REP(i, lim) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1; i < lim; i <<= 1) {
		cint wn = {
			fp1(o == 1 ? 3 : (P1 + 1) / 3, (P1 - 1) / (i << 1)),
			fp2(o == 1 ? 3 : (P2 + 1) / 3, (P2 - 1) / (i << 1)),
			fp3(o == 1 ? 3 : (P3 + 1) / 3, (P3 - 1) / (i << 1))
		};
		for(int j = 0; j < lim; j += (i << 1)) {
			cint w = {1, 1, 1};
			for(int k = 0; k < i; k ++, w *= wn) {
				cint x = a[j + k];
				cint y = a[j + k + i] * w;
				a[j + k] = x + y;
				a[j + k + i] = x - y;
			}
		}
	}
	if(o == 1) return;
	cint inv = {fp1(lim, P1 - 2), fp2(lim, P2 - 2), fp3(lim, P3 - 2)};
	REP(i, lim) a[i] *= inv;
}
int n, m, P; cint F[N], G[N];
void solve() {
	cin >> n >> m >> P;
	FOR(i, 0, n) {
		int x; cin >> x; x %= P;
		F[i] = {x, x, x};
	}
	FOR(i, 0, m) {
		int x; cin >> x; x %= P;
		G[i] = {x, x, x};
	}
	int lim = 1, l = 0;
	while(lim <= (n + m)) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	NTT(F, lim, 1);
	NTT(G, lim, 1);
	REP(i, lim) F[i] *= G[i];
	NTT(F, lim, - 1);
	FOR(i, 0, n + m) cout << F[i].crt(P) << " "; cout << endl;
}

快速阶乘

给定 $n$ 和 $P$，求 $n! \bmod P$。

考虑分块，设快长 $B=\sqrt{n}$，并考虑求解一段数连乘的值，考虑设多项式 $F_d(x)=\prod_{i=1}^{d}(x+i)$。

那么 $n!=\prod_{i=0}^{B-1}F_B(iB)\prod_{i=B^2+1}^{n}i$，也就是求出 $F_B(x)$ 即可。

考虑求解 $F_B(x)$，但是在 $F_B(x)$ 需要的下标都是 $iB$，这启发直接计算点值。

考虑倍增求解点值，考虑 $F_d(iB)$ 如何转移到 $F_{2d}(iB)$，发现 $F_{2d}(x)=F_d(x)F_d(x+d)$。

所以直接将 $F_d(iB)$ 进行连续点值平移到 $F_d(iB+d)$ 也就是 $F_d((i+\frac{d}{B})B)$ 即可。

但这样是需要两只 $\log$。

考虑缩减范围，$F_d(x)$ 就是一个 $d+1$ 次多项式，只保留 $F_d(0),F_d(B),F_d(2B),\cdots,F_d(dB)$ 一共 $d+1$ 即可。

倍增到 $F_{2d}(x)$ 的时候再做一次点值平移得出 $F_d((d+1)B),F_d((d+1)B),\cdots,F_d(2dB)$ 即可。

这样时间复杂度就是 $O(\sqrt{n} \log n)$。

ケロシの代码

int n, P, B;
int a[N], b[N];
cint F[N], G[N];
int fac[N], finv[N];
void moves(int * a, int * b, int n, int m) {
	int lim = 1, l = 0;
	while(lim <= (n << 1)) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	REP(i, lim) F[i] = G[i] = Cint(0);
	FOR(i, 0, n) F[i] = Cint(mul({a[i], finv[i], finv[n - i], sgn(n - i)}));
	FOR(i, 0, n << 1) G[i] = Cint(inv(sub(m, n - i)));
	NTT(F, lim, 1);
	NTT(G, lim, 1);
	REP(i, lim) F[i] *= G[i];
	NTT(F, lim, - 1);
	int res = 1;
	FOR(i, 0, n) Mul(res, sub(m, i));
	FOR(i, 0, n) {
		b[i] = mul(F[n + i].crt(P), res);
		Mul(res, inv(sub(m, n - i)));
		Mul(res, add(m, i + 1));
	}
}
void facs(int d) {
	if(d == 1) {
		a[0] = 1; a[1] = B + 1;
		return;
	}
	facs(d >> 1);
	moves(a, b, d >> 1, (d >> 1) + 1);
	REP(i, d) a[(d >> 1) + i + 1] = b[i];
	moves(a, b, (d >> 1) << 1, mul(d >> 1, inv(B)));
	FOR(i, 0, (d >> 1) << 1) Mul(a[i], b[i]);
	if(d & 1) {
		REP(i, d) Mul(a[i], i * B + d);
		a[d] = 1;
		FOR(i, 1, d) Mul(a[d], d * B + i);
	}
}
void solve() {
	cin >> n >> P;
	B = sqrt(n);
	fac[0] = 1;
	FOR(i, 1, B) fac[i] = mul(fac[i - 1], i);
	finv[B] = fp(fac[B], P - 2);
	ROF(i, B - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
	facs(B);
	int ans = 1;
	REP(i, B) Mul(ans, a[i]);
	FOR(i, B * B + 1, n) Mul(ans, i);
	cout << ans << endl;
}

第二类斯特林数行

给出 $n$，对于 $i\in[0,n]$，求出 $\begin{Bmatrix} n \\ i \end{Bmatrix}$。

考虑第二类斯特林数的通项公式：

$$\begin{aligned} \begin{Bmatrix} n \\ i \end{Bmatrix} &= \sum_{j=0}^{i} \frac{(-1)^{i-j}j^n}{j!(i-j)!} \\ &=\sum_{j=0}^{i} \frac{j^n}{j!}\frac{(-1)^{i-j}}{(i-j)!} \end{aligned}$$

直接设多项式 $F_i=\frac{i^n}{i!}$，$G_i=\frac{(-1)^{i}}{i!}$，直接卷积求解即可。

时间复杂度 $O(n \log n)$。

ケロシの代码

int n;
int fac[N], finv[N];
int F[N], G[N];
void init(int n) {
	fac[0] = 1;
	FOR(i, 1, n) fac[i] = mul(fac[i - 1], i);
	finv[n] = fp(fac[n], P - 2);
	ROF(i, n - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
}
void solve() {
	cin >> n;
	init(n);
	int lim = 1, l = 0;
	while(lim <= (n << 1)) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	FOR(i, 0, n) F[i] = mul(fp(i, n), finv[i]);
	FOR(i, 0, n) G[i] = mul(sgn(i), finv[i]);
	NTT(F, lim, 1);
	NTT(G, lim, 1);
	REP(i, lim) Mul(F[i], G[i]);
	NTT(F, lim, - 1);
	FOR(i, 0, n) cout << F[i] << " "; cout << endl;
}

第二类斯特林数列

给出 $n$ 和 $k$，对于 $i\in[0,n]$，求出 $\begin{Bmatrix} i \\ k \end{Bmatrix}$。

考虑组合意义，先乘上 $k!$ 变成 $i$ 个有标号球放进 $k$ 个有标号非空盒子里，考虑每次加一个盒子，相当于选择放多少个球再分配标号，不难发现 EGF 的乘法可以做这个。

所以设只放一个盒子的 EGF，除了放 $0$ 个都有 $1$ 的方案，所以设 EGF $F=\sum_{i=1}\frac{x^i}{i!}=e^x-1$。

所以答案的 EGF 即为 $\frac{F^k}{k!}$。

使用多项式快速幂即可。

时间复杂度 $O(n \log n)$。

ケロシの代码

int n, k;
int F[N], G[N];
void solve() {
	cin >> n >> k;
	if(n < k) {
		REP(_, n + 1) cout << 0 << " ";
		return;
	}
	init(n);
	int len = n - k + 1;
	REP(i, len) F[i] = finv[i + 1];
	lns(F, G, len);
	REP(i, len) Mul(G[i], k);
	exps(G, F, len);
	REP(_, k) cout << 0 << " ";
	REP(i, len) cout << mul({F[i], fac[i + k], finv[k]}) << " "; cout << endl;
}

第一类斯特林数行

给出 $n$，对于 $i\in[0,n]$，求出 $\begin{bmatrix} n \\ i \end{bmatrix}$。

$$x^\overline{n} =\sum_{i=0}^{n} \begin{bmatrix} n \\ i \end{bmatrix} x_i$$

所以只要求 $x^\overline{n}$ 即可，考虑倍增：

$$x^\overline{2n}=x^\overline{n}(x+n)^\overline{n}$$

然后考虑 $(x+n)^\overline{n}$ 怎么求：

$$\begin{aligned} (x+n)^\overline{n}&=\sum_{i=0}^{n} \begin{bmatrix} n \\ i \end{bmatrix} (x+n)^i \\ &=\sum_{i=0}^{n} \begin{bmatrix} n \\ i \end{bmatrix} \sum_{j=0}^{i} \binom{i}{j} x^{j} n^{i-j} \\ &=\sum_{i=0}^{n} \begin{bmatrix} n \\ i \end{bmatrix} \sum_{j=0}^{i} \frac{i!}{j!(i-j)!} x^{j} n^{i-j} \\ &=\sum_{j=0}^{n} \frac{x^{j}}{j!} \sum_{i=j}^{n} \begin{bmatrix} n \\ i \end{bmatrix} \frac{i!}{(i-j)!} n^{i-j} \\ &=\sum_{j=0}^{n} \frac{x^{j}}{j!} \sum_{i=j}^{n} i! \begin{bmatrix} n \\ i \end{bmatrix} \frac{n^{i-j}}{(i-j)!}\\ \end{aligned}$$

右边很明显是一个卷积形式，设 $F_i=i! \begin{bmatrix} n \\ i\end{bmatrix}$ 和 $G_{n-i}=\frac{n^i}{i!}$，直接卷即可。

时间复杂度 $O(n \log n)$。

ケロシの代码

int n, a[N], b[N];
int fac[N], finv[N];
int F[N], G[N];
void init(int n) {
	fac[0] = 1;
	FOR(i, 1, n) fac[i] = mul(fac[i - 1], i);
	finv[n] = fp(fac[n], P - 2);
	ROF(i, n - 1, 0) finv[i] = mul(finv[i + 1], i + 1);
}
void slv(int d) {
	if(d == 1) {
		a[1] = 1;
		return;
	}
	if(d & 1) {
		slv(d - 1);
		FOR(i, 0, d) F[i] = 0;
		REP(i, d) Add(F[i + 1], a[i]);
		REP(i, d) Add(F[i], mul(a[i], d - 1));
		FOR(i, 0, d) a[i] = F[i];
		return;
	}
	slv(d >> 1);
	int len = d >> 1;
	int lim = 1, l = 0;
	while(lim <= d) lim <<= 1, l ++;
	REP(i, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
	REP(i, lim) F[i] = G[i] = 0;
	FOR(i, 0, len) F[i] = mul(fac[i], a[i]);
	FOR(i, 0, len) G[len - i] = mul(finv[i], fp(len, i));
	NTT(F, lim, 1);
	NTT(G, lim, 1);
	REP(i, lim) Mul(F[i], G[i]);
	NTT(F, lim, - 1);
	REP(i, lim) b[i] = 0;
	FOR(i, 0, len) b[i] = mul(finv[i], F[len + i]);
	NTT(a, lim, 1);
	NTT(b, lim, 1);
	REP(i, lim) Mul(a[i], b[i]);
	NTT(a, lim, - 1);
}
void solve() {
	cin >> n;
	init(n);
	slv(n);
	FOR(i, 0, n) cout << a[i] << " "; cout << endl;
}

第一类斯特林数列

给出 $n$ 和 $k$，对于 $i\in[0,n]$，求出 $\begin{bmatrix} i \\ k \end{bmatrix}$。

考虑使用与第二类斯特林数列相同的做法，但是注意单个环的方案是 $(i-1)!$，所以 EGF $F=\sum_{i=1}\frac{(i-1)!x^i}{i!}=\sum_{i=1}\frac{x^i}{i}=\ln(\frac{1}{1-x})$。

时间复杂度 $O(n \log n)$。

ケロシの代码

int n, k;
int F[N], G[N];
void solve() {
	cin >> n >> k;
	if(n < k) {
		REP(_, n + 1) cout << 0 << " ";
		return;
	}
	init(n);
	int len = n - k + 1;
	REP(i, len) F[i] = inv[i + 1];
	lns(F, G, len);
	REP(i, len) Mul(G[i], k);
	exps(G, F, len);
	REP(_, k) cout << 0 << " ";
	REP(i, len) cout << mul({F[i], fac[i + k], finv[k]}) << " "; cout << endl;
}