多项式的高级运算

基础：FFT与NTT

任意模数NTT（MTT）

有的时候质数不能表示为 (a * 2^k + 1) 的形式，无法进行NTT。如果质数比较小（比如 (10^4) 左右），就可以直接用 FFT 搞过去；但是如果质数到达 (10^9) 左右，就需要一些技巧了。

我们可以把所有系数化为 (ax+b) 的形式，其中 (x) 为我们指定的一个数（通常是 (2^{15})），这样的话 (a, b) 就都不超过 (2^{15}) 了，可以暴力求解。答案为 (aa'x ^2+ (ab'+ba')x + bb')，根据需要进行乘法即可。

一种减小常数的方法是：我们设复多项式 (P = (a, a'),Q=(b,b'),R=(a,-a'))，然后设 (F=PQ=(ab-a'b',ab'+a'b),F'=RQ=(ab+a'b',ab'-a'b))，再加加减减：(F+F'=(2ab,2ab'),F-F'=(-2a'b',2a'b))，这样，我们需要的东西就都有了，并且只用做 5 次FFT。

关键代码：

	for (register int i = 0; i <= n; ++i) {
		int a; read(a);
		P[i].x = a >> 15;
		P[i].y = a & mask;
		R[i].x = P[i].x, R[i].y = -P[i].y;
	}
	for (register int i = 0; i <= m; ++i) {
		int a; read(a);
		Q[i].x = a >> 15;
		Q[i].y = a & mask;
	}
    
	fft(P, 1), fft(Q, 1), fft(R, 1);
	for (register int i = 0; i < limi; ++i) {
		Q[i].x /= limi, Q[i].y /= limi;//提前除以 limi，fft就不用除了
		P[i] = P[i] * Q[i], R[i] = R[i] * Q[i];
	}
	fft(P, -1), fft(R, -1);
	for (register int i = 0; i <= n + m; ++i) {
		ll a1b1 = (ll)((P[i].x + R[i].x) / 2.0 + 0.5);
		ll a1b2 = (ll)((P[i].y + R[i].y) / 2.0 + 0.5);
		ll a2b1 = (ll)((P[i].y - R[i].y) / 2.0 + 0.5);
		ll a2b2 = (ll)((R[i].x - P[i].x) / 2.0 + 0.5);
		ll ans = a1b1 % Mod * ((1ll << 30) % Mod) % Mod;
		ans += (a2b1 + a1b2) % Mod * ((1ll << 15) % Mod) % Mod;
		ans += a2b2 % Mod;
		ans = (ans % Mod + Mod) % Mod;
		printf("%lld ", ans);
	}
	puts("");

多项式求乘法逆元

【模板】多项式乘法逆

牛顿迭代法：

[AB=1 ]

[AB - 1 = 0 ]

[B' = B - frac{AB-1}{A} ]

[=B-B(AB-1) ]

[=B(2-AB) ]

(Code：)

int n;
ll A[N], B[N], C[N], r[N];
ll limi, l;
inline ll quickpow(ll x, ll k)...
inline void ntt(ll *a, int type) {...//此处已经让type = 1的乘inv了
void sol(int deg, ll *a, ll *b) {//b is 逆元
	if (deg == 1) {
		b[0] = quickpow(a[0], P - 2);
		return ;
	}
	sol((deg + 1) >> 1, a, b);
	
	limi = 1, l = 0;
	while (limi <= (deg << 1))	limi <<= 1, ++l;
	for (register int i = 1; i <= limi; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
	for (register int i = 0; i < deg; ++i)	C[i] = a[i];//转移到C防变化 
	for (register int i = deg; i < limi; ++i)	C[i] = 0;//多次清空更保险 
	ntt(b, 1); ntt(C, 1);
	for (register int i = 0; i < limi; ++i)//B = 2B' - AB'B' = B'(2 - AB')
		b[i] = ((2ll - b[i] * C[i]) % P + P) % P * b[i] % P;
	ntt(b, -1);
	for (register int i = deg; i <= limi; ++i)//多次清空更保险 
		b[i] = 0;
}
int main() {
	read(n);
	for (register int i = 0; i < n; ++i)	read(A[i]);
	sol(n, A, B);
	for (register int i = 0; i < n; ++i)
		printf("%lld ", B[i]);
	return 0;
}

这里给一个简化的板子，供复习：

7.28 Update : 拿了个更简化的倍增版本

inline void get_inv(ll *a, ll *b, int d) {//d : 项数 
	b[0] = quickpow(a[0], P - 2);
	L = 0; 
	for (register int len = 1; len < (d << 1); len <<= 1) {//一定求完长度为 len/2 的 b 数组，想求长度为 len 的 b 数组
		limi = len << 1, ++L;
		for (register int i = 1; i < limi; ++i)
			r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
		
		for (register int i = 0; i < len; ++i)	C[i] = a[i], D[i] = b[i];
		ntt(C, 1), ntt(D, 1);
		for (register int i = 0; i < limi; ++i)
			b[i] = D[i] * (2 - C[i] * D[i] % P) % P, C[i] = D[i] = 0;
		ntt(b, -1);
		
		for (register int i = len; i < limi; ++i)	b[i] = 0;
	}
}

多项式开根号

• P5205 【模板】多项式开根

题意

• 给A(x)，其中a0 = 1，求B(x)，使得B(x)^2 = A(x) (mod x^n)

思路简析

与多项式求逆相同，由于一平方就mod x^m -> mod x^(2m)，我们考虑递归求解。

表达式

同样假设我们已经求出B(x)的一半b(x)，那么：

A * b = 1(mod x^m)
A * B = 1(mod x^m)
∴B - b = 0（mod x^m)
两边平方:
B^2 - 2 * B * b + b^2 = 0(mod x^(2*m))
据B^2 = A(mod x^(2*m))：
A - 2 * B * b + b^2 = 0(mod x^(2*m))
于是：
B = (A/b + b) / 2

配合多项式求逆解出B(x)。

（7.28 Update:）

这里给一个更方便，无需更多的聪明与技巧的方法：牛顿迭代:

[sqrt{A(x)}=B(x) ]

[A(x)=B^2(x) ]

[F(B(x)) = B^2(x)-A(x)=0 ]

[B_{2n}=B_n-frac{B^2-A}{2B} ]

[=frac{B^2 + A}{2B} ]

[=1/2 * (B + frac{A}{B}) ]

我们可以一遍 NTT 算出 (A, B, frac{1}{B}) 的点值表示，然后在运算并 NTT 搞回去，但是这需要 NTT 三次，因为第一个 (B) 与分式是加法关系，直接系数相加即可，这样就只用 NTT 两次。

递归边界

模板题非常友善，告诉我们 a0 = 1，于是递归边界为

if (deg == 1) {b[0] = 1; return ;}

如果题目没有那么友善，那么我们或许可以多random几个数 我们需要用二次剩余之类的麻烦的东西，或者考虑换一种算法。

浅谈二次剩余：

Code:

void get_sqrt(ll *a, ll *b, int deg) {
	if (deg == 1) {b[0] = 1; return ;}
	get_sqrt(a, b, (deg + 1) >> 1);
  //get_len
	limi = 1, len = 0;
	while (limi <= (deg << 1))	limi <<= 1, len++;
	for (register int i = 0; i <= limi; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
        
  //copy and multiply
	for (register int i = 0; i <= limi; ++i)	bn[i] = 0;//attention
	get_inv(b, bn, deg);
	for (register int i = 0; i < deg; ++i)	C[i] = a[i];
	for (register int i = deg; i <= limi; ++i)	C[i] = 0;
	ntt(C, 1); ntt(bn, 1);
	for (register int i = 0; i <= limi; ++i)
		C[i] = C[i] * bn[i] % P;
	ntt(C, -1);
	for (register int i = 0; i < deg; ++i)	b[i] = (C[i] + b[i]) * inv2 % P;
	for (register int i = deg; i <= limi; ++i)	b[i] = 0;
}

(7.28 Update)

倍增版本：（附新求逆模板）（我怕两个板子之间的 (limi) 互相串出错，因此每个函数单独写了个 (limi) 和 (L)。NTT 的时候直接传入 (limi)）

inline void get_inv(ll *a, ll *b, int d) {
	b[0] = quickpow(a[0], P - 2);
	int L = 0, limi = 1;
	for (register int len = 1; len < (d << 1); len <<= 1) {
		limi = (len << 1), ++L;
		for (register int i = 1; i < limi; ++i)
			r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
		
		for (register int i = 0; i < len; ++i)	E[i] = a[i], F[i] = b[i];
		ntt(E, 1, limi), ntt(F, 1, limi);
		for (register int i = 0; i < limi; ++i)
			b[i] = F[i] * (2 - E[i] * F[i] % P) % P, E[i] = F[i] = 0;
		ntt(b, -1, limi);
		
		for (register int i = len; i < limi; ++i)	b[i] = 0;
	}
}
inline void Sqrt(ll *a, ll *b, int d) {
	b[0] = 1;
	int limi = 1, L = 0;
	for (register int len = 1; len < (d << 1); len <<= 1) {
		get_inv(b, C, len);
		
		limi = len << 1, ++L;
		for (register int i = 1; i < limi; ++i)
			r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
		
		for (register int i = 0; i < len; ++i)	D[i] = a[i];
		ntt(D, 1, limi); ntt(C, 1, limi);
		for (register int i = 0; i < limi; ++i)	C[i] = D[i] * C[i] % P;
		ntt(C, -1, limi);
		for (register int i = 0; i < len; ++i)
			b[i] = INV2 * (b[i] + C[i]) % P, C[i] = D[i] = 0;
		for (register int i = len; i < limi; ++i)	b[i] = 0;
	}
}

注意：

1. 用数组前一定注意清空。我也不知道为什么，反正不清空就会出错。估计是NTT的祸吧。

多项式除法

• P4512 【模板】多项式除法

• A(x) * B(x) + C(x) = D(x)，给出D(x), A(x)，求B(x)，C(x)。（类似高精除）

• m = deg(A) < 100,000, n = deg(D) <= 100,000，m <= n

思路简析

我们发现神奇的事情：把A,B,C,D都翻转过来，（把C加到D后面），等式仍然成立。并且还有一个好处：C转过来后D的0~n-m项都不受C的影响，而B又肯定超不过n-m+1项。因此我们可以借助反转后的A,D数组算出B数组，然后什么都好搞了。

Code:

int main() {
	read(n); read(m); n++; m++;//n,m变成项数
	for (register int i = 0; i < n; ++i)	read(D[i]), Dbp[i] = D[i];
	for (register int i = 0; i < m; ++i)	read(A[i]), Abp[i] = A[i];//backup
	Reverse(A, m);
	Reverse(D, n);
	for (register int i = n - m + 1; i < n; ++i)
		A[i] = D[i] = 0;
	get_inv(A, An, n - m + 1);
	mul(D, An, B, n - m + 1, n - m + 1);
    //D(n-m+1项) * An(n-m+1项) -> B
	Reverse(B, n - m + 1);
	mul(Abp, B, AB, m, n - m + 1);
	for (register int i = 0; i < m - 1; ++i)
		C[i] = (Dbp[i] - AB[i] + P) % P;
    ...
}

注意

• 在算D*inv(A)时，保险起见，只保留D和A的0~n-m项，且对A,D做备份，算C时用。

• 什么时候用n-m，什么时候用n-m+1，要分清楚！

• 此时多项式变量名逐渐增多，注意区分，不要把An写成A!

分治FFT

（其实我觉得对于我这个几乎只写NTT的人来说，叫分治NTT比较好）

【模板】分治 FFT

简单说一下，分治FFT用到了CDQ分治的思想，但不用非得学CDQ分治，毕竟这个思想还是比较好理解的，之前也经常用到。简而言之，这里的CDQ分治思想就是：每次只考虑左半段对右半段的贡献，先递归解决左半段，然后让右半段加上左半段的贡献，再递归解决右半段。这样，一次次贡献的加和就组成了每个位置的值。

将题意稍稍转化，f[i] = g[0 ~ i]与f[0 ~ i]的卷积（多项式乘法）。剩下的推导部分可以通过手玩来推导。

最关键的两条语句：

    for (int i = l;i <= mid; i++) a[i-l] = f[i];
    for (int i = 0;i < len; i++) b[i] = g[i];

(Code:) my code

（主要篇幅是NTT，其实重点就在于以上两条语句，NTT只是工具）

多项式求ln

• P4725 【模板】多项式对数函数（多项式 ln）

在学习微积分后，我再学ln，感觉舒适了很多。

求ln很简单，两边求个导，用一下多项式求逆，再积分即可。

7.28 Update:

给一下推导：

[ln(A(x))=B(x) ]

[frac{A'(x)}{A(x)} = B(x) ]

[B(x) = intfrac{A'(x)}{A(x)} ]

思维难度低，代码量大。

好想好写，算是比较小清新的板子了。

（更新封装风格的代码）

inline void dao(ll *a, ll *b, int d) {
	for (register int i = 0; i < d; ++i)	b[i] = a[i + 1] * (i + 1) % P;
}
inline void ji(ll *a, ll *b, int d) {//利用线性推逆元，积分可以做到 O(n)
	for (register int i = 1; i < d; ++i)
		b[i] = a[i - 1] * quickpow(i, P - 2) % P;
}
inline void get_ln(ll *a, ll *b, int d) {
	get_inv(a, C, d);
	dao(a, D, d);
	int limi = 1, L = 0;
	while (limi < (d << 1))	limi <<= 1, ++L;
	for (register int i = 1; i < limi; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
	ntt(C, 1, limi), ntt(D, 1, limi);
	for (register int i = 0; i < limi; ++i)	C[i] = C[i] * D[i] % P;
	ntt(C, -1, limi);
	for (register int i = d; i < limi; ++i)	C[i] = 0;
	ji(C, b, d);
}

多项式求exp

继续牛顿迭代：

[e^{A(x)} = B(x) ]

[A(x) = ln(B(x)) ]

[A(x) - ln(B(x)) = 0 ]

[B_{2n}(x) = B_n(x) - frac{A(x) - ln(B(x))}{-frac{1}{B(x)}} ]

[B_{2n} = B(x)(1 - ln(B(x)) + A(x)) ]

然后就终极套娃即可。

(1 - ln(B(x)) + A(x)) 可以用系数表示法直接运算，但是记住这是系数表示法，不是点值表示，不是每一项都 + 1，而是只有常数项 +1

多项式快速幂

自然可以倍增快速幂，不过那是 (n~logn~logk) 的。既然多项式可以快速取 (ln,exp)，我们可以利用这一点来做到 (n~logn)

[A^k(x) = B(x) ]

[klnA(x)=lnB(x) ]

[B(x) = e^{klnA(x)} ]

套 ln 和 exp 即可。

下面给一个 快速幂 + exp + ln + 求逆求导积分 的代码：

ll inv_A[N], inv_B[N];
inline void get_inv(ll *a, ll *b, int d) {
	b[0] = quickpow(a[0], P - 2);
	int L = 0, limi = 1;
	for (register int len = 1; len < (d << 1); len <<= 1) {
		limi = (len << 1); ++L;
		for (register int i = 1; i < limi; ++i)
			r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
		
		for (register int i = 0; i < len; ++i)	inv_A[i] = a[i], inv_B[i] = b[i];
		ntt(inv_A, 1, limi), ntt(inv_B, 1, limi);
		for (register int i = 0; i < limi; ++i)
			b[i] = inv_B[i] * (2ll - inv_A[i] * inv_B[i] % P) % P, inv_A[i] = inv_B[i] = 0;
		ntt(b, -1, limi);
		
		for (register int i = len; i < limi; ++i)	b[i] = 0;
	}
}
inline void Dao(ll *a, ll *b, int d) {
	for (register int i = 0; i < d; ++i)	b[i] = a[i + 1] * (i + 1) % P;
	b[d - 1] = 0;
}
inline void Ji(ll *a, ll *b, int d) {
	for (register int i = 1; i < d; ++i)	b[i] = a[i - 1] * inv[i] % P;
	b[0] = 0;
}
ll Ln_A[N], Ln_B[N];
inline void get_Ln(ll *a, ll *b, int d) {
	b[0] = 0;
	Dao(a, Ln_A, d);
	get_inv(a, Ln_B, d);
	int limi = 1, L = 0;
	while (limi < (d << 1)) limi <<= 1, ++L;
	for (register int i = 1; i < limi; ++i)
		r[i] = (r[i>> 1] >> 1) | ((i & 1) << (L - 1));
	ntt(Ln_A, 1, limi); ntt(Ln_B, 1, limi);
	for (register int i = 0; i < limi; ++i)	Ln_A[i] = Ln_A[i] * Ln_B[i] % P;
	ntt(Ln_A, -1, limi);
	Ji(Ln_A, b, d);
	for (register int i = 0; i < limi; ++i)	Ln_A[i] = Ln_B[i] = 0;//Attention!!
}
ll Exp_A[N];
inline void get_Exp(ll *a, ll *b, int d) {
	b[0] = 1;
	int L = 0, limi = 1;
	for (register int len = 1; len < (d << 1); len <<= 1) {
		get_Ln(b, Exp_A, len);
		Exp_A[0] = (1 - Exp_A[0] + a[0]) % P;
		for (register int i = 1; i < len; ++i)	Exp_A[i] = (a[i] - Exp_A[i]) % P;
		limi = len << 1; ++L;
		for (register int i = 1; i < limi; ++i)
			r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
		ntt(b, 1, limi), ntt(Exp_A, 1, limi);
		for (register int i = 0; i < limi; ++i)	b[i] = (b[i] * Exp_A[i]) % P;
		ntt(b, -1, limi);
		
		for (register int i = len; i < limi; ++i)	b[i] = 0;
		for (register int i = 0; i < limi; ++i)	Exp_A[i] = 0;
	}
}
ll Pow_A[N];
inline void get_Pow(ll *a, ll *b, int d, ll k) {
	get_Ln(a, Pow_A, d);
	for (register int i = 0; i < d; ++i)	Pow_A[i] = Pow_A[i] * k % P;
	get_Exp(Pow_A, b, d);
}