多项式学习笔记

Posted on 2023-07-13 In XCPC , Note

之前学习多项式相关算法的学习笔记, 从原博客搬运而来.

(以下为原内容)

这篇博客会讲一下我学多项式的板子（可能不太完善, 还请谅解）.

板子参考了这里.

学习资料:

多项式与生成函数学习笔记

Karry5307's Blog

写在前面的话（板子使用指南）

基本所有的运算都可以大胆使用, 不用要求结果数组为空.

把用来计算的数组叫做操作数组, 存储结果的数组叫做结果数组. 则这个板子操作数组一般可以相同, 操作数组与结果数组之间一般不可以相同.

除多项式除法外, 所有的n意为多项式的项数（或者说是界）, 而非度数（应为度数+1）.

所有的代码都附了一个检验方法, 用来判断是否打对, 建议每完成1~2个模块就检验一次.

FFT & NTT

多项式乘法(FFT/NTT)

这篇博客里有详细代码. 这里给一个下面都会用到的"封装版"多项式乘法（预处理了单位根和逆元（方便积分））:

检验方法:

找一个多项式平方. 我喜欢找 \(1+x+x^2+\cdots\), 它的平方相当于求前缀和.

#define clr(f, s, t) memset(f + (s), 0x00, sizeof(int) * ((t) - (s)))
#define cpy(f, g, n) memcpy(g, f, sizeof(int) * (n))
const int MAXN = (1 << 18) + 5, bas = 1 << 18, P = 998244353, G = 3, invG = 332748118;
int pls(int a, int b) {return a + b < P ? a + b : a + b - P;}
int mns(int a, int b) {return a < b ? a + P - b : a - b;}
int mul(int a, int b) {return 1ll * a * b % P;}
int qpow(int a, int n) {int ret = 1; for(; n; n >>= 1, a = mul(a, a)) if(n & 1) ret = mul(ret, a); return ret;}
int _g[2][MAXN], tf, tr[MAXN], inv[MAXN];
void init() {
	for(int i = 0; i < bas; i++) {
		_g[1][i] = qpow(G, (P-1) / bas * i);
		_g[0][i] = qpow(invG, (P-1) / bas * i);
	}
	inv[1] = 1; for(int i = 2; i < MAXN; i++) inv[i] = mul(P - P/i, inv[P % i]);
}
int getlim(int n) {
	int lim = 1; for(; lim < n + n; lim <<= 1);
	return lim;
}
void tpre(int lim) {
	if(tf == lim) return ;
	tf = lim; for(int i = 0; i < lim; i++) tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
}
void NTT(int* f, int lim, int fl) {
	tpre(lim); for(int i = 0; i < lim; i++) if(i < tr[i]) swap(f[i], f[tr[i]]);
	for(int l = 2, k = 1; l <= lim; l <<= 1, k <<= 1) 
		for(int i = 0; i < lim; i += l) 
			for(int j = i; j < i+k; j++) {
				int tt = mul(f[j+k], _g[fl][bas / l * (j-i)]);
				f[j+k] = mns(f[j], tt);
				f[j] = pls(f[j], tt);
			}
	if(!fl) for(int i = 0; i < lim; i++) f[i] = mul(f[i], inv[lim]);
}
void Mul(int* f, int* g, int* h, int n) {
	static int a[MAXN], b[MAXN];
	int lim = getlim(n);
	cpy(f, a, n); clr(a, n, lim);
	cpy(g, b, n); clr(b, n, lim);
	NTT(a, lim, 1); NTT(b, lim, 1);
	for(int i = 0; i < lim; i++) h[i] = mul(a[i], b[i]);
	NTT(h, lim, 0);
}
/*
int n, f[MAXN], g[MAXN];
void test() {
	n = read();
	for(int i = 0; i < n; i++) f[i] = read();
	Mul(f, f, g, n);
	for(int i = 0; i < n; i++) printf("%d ", g[i]);
}
Input:
5
1 1 1 1 1
Output:
1 2 3 4 5
*/

任意模数NTT（MTT）

P4245 【模板】任意模数多项式乘法

给 2 个多项式 \(F(x),G(x)\), 求 \(F(x)G(x)\). 系数对 \(p\) 取模, 不保证 \(p\) 是 NTT 模数.

也就是MTT, 使用 4 次 FFT 完成任意模数的多项式乘法.

设 \(K=2^{15}\), 我们把多项式每项系数分为两部分（高低位）.
\[ F(x)=K\cdot F_1(x)+F_0(x) \\ G(x)=K\cdot G_1(x)+G_0(x) \\ \therefore F(x)G(x)=K^2\cdot F_1(x)G_1(x)+K\cdot [F_1(x)G_0(x)+F_0(x)G_1(x)]+F_0(x)G_0(x) \] 如何快速得到这四个多项式的点值表示？

构造 \[ P(x)=F_0(x)+iG_0(x) \\ Q(x)=F_0(x)-iG_0(x) \] 我们惊奇地发现: \[ \mathrm{DFT}(P)[j]=P(w_n^j)=F_0(w_n^j)+iG_0(w_n^j) \\ =\sum_{k=0}^{n-1}F_0[k]w_n^{kj}+i\sum_{k=0}^{n-1}G_0[k]w_n^{kj} \\ =\sum_{k=0}^{n-1}(F_0[k]+iG_0[k])(\cos(\dfrac {2\pi kj}{n})+i\sin(\dfrac{2\pi kj}n)) \\ =\sum_{k=0}^{n-1}(F_0[k]\cos(\dfrac{2\pi kj}n)-G_0[k]\sin(\dfrac {2\pi kj}n))+\\ i\sum_{k=0}^{n-1}(F_0[k]\sin(\dfrac {2\pi kj}n)+G_0[k]\sin(\dfrac {2\pi kj}n)) \]

同理 \[ \mathrm{DFT}(Q)[n-j]=P(w_n^{-j})=F_0(w_n^{-j})-iG_0(w_n^{-j}) \\ =\sum_{k=0}^{n-1}F_0[k]w_n^{-kj}-i\sum_{k=0}^{n-1}G_0[k]w_n^{-kj} \\ =\sum_{k=0}^{n-1}(F_0[k]-iG_0[k])(\cos(\dfrac {2\pi kj}{n})-i\sin(\dfrac{2\pi kj}n)) \\ =\sum_{k=0}^{n-1}(F_0[k]\cos(\dfrac{2\pi kj}n)-G_0[k]\sin(\dfrac {2\pi kj}n))+\\ i\sum_{k=0}^{n-1}(F_0[k]\sin(\dfrac {2\pi kj}n)+G_0[k]\sin(\dfrac {2\pi kj}n)) \] 故 \(P\) 的第 \(j\) 项点值与 \(Q\) 的第 \(n-j\) 项点值共轭.

于是我们可以使用 1 次 FFT 得到 \(P(x)\) 和 \(Q(x)\) 的点值, 再解方程就可得到 \(F_0(x)\) 和 \(G_0(x)\) 的点值.

同样地可得到 \(F_1(x),G_1(x)\), 使用了 2 次FFT.

然后考虑怎么求解回系数.

构造 \[ P(x)=F_1(x)G_1(x)+i(F_1(x)G_0(x)+F_0(x)G_1(x)) \\ Q(x)=F_0(x)G_0(x) \] 做两次 IDFT 即可.

下面是与上面相对的另一个板子. 如果要把其它多项式运算运用到任意模数下, 要求把所有的乘法用 Mul 代替.

检验方法: 同上

#define clr(f, s, t) memset(f + (s), 0x00, sizeof(int) * ((t) - (s)))
#define cpy(f, g, n) memcpy(g, f, sizeof(int) * (n))
const int MAXN = (1 << 19) + 5, bas = 1 << 19;
const db PI = acos(-1.0);
int P;
int pls(int a, int b) {return a + b < P ? a + b : a + b - P;}
int mns(int a, int b) {return a < b ? a + P - b : a - b;}
int mul(int a, int b) {return 1ll * a * b % P;}
int qpow(int a, int n) {int ret = 1; for(; n; n >>= 1, a = mul(a, a)) if(n & 1) ret = mul(ret, a); return ret;}
struct cp {db x, y;};
cp operator + (const cp& a, const cp& b) {return (cp){a.x + b.x, a.y + b.y};}
cp operator - (const cp& a, const cp& b) {return (cp){a.x - b.x, a.y - b.y};}
cp operator * (const cp& a, const cp& b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
cp operator * (const cp& a, const db& k) {return (cp){a.x * k, a.y * k};}
const cp I = (cp){0, 1};
cp _g[2][MAXN];
int tr[MAXN], tf;
void init() {
	for(int i = 0; i < bas; i++) {
		db a = cos(2 * PI * i / bas), b = sin(2 * PI * i / bas);
		_g[1][i] = (cp){a, b};
		_g[0][i] = (cp){a, -b};
	}
}
int getlim(int n) {
	int lim = 1; for(; lim < n + n; lim <<= 1);
	return lim;
}
void tpre(int lim) {
	if(tf == lim) return;
	tf = lim; for(int i = 0; i < lim; i++) tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
}
ll tran(db x) {return ((ll)(x > 0 ? x + .5 : x - .5) % P + P) % P;}
void FFT(cp* f, int lim, int fl) {
	tpre(lim); for(int i = 0; i < lim; i++) if(i < tr[i]) swap(f[i], f[tr[i]]);
	for(int l = 2, k = 1; l <= lim; l <<= 1, k <<= 1)
		for(int i = 0; i < lim; i += l)
			for(int j = i; j < i+k; j++) {
				cp tt = f[j+k] * _g[fl][(j-i) * (bas / l)];
				f[j+k] = f[j] - tt;
				f[j] = f[j] + tt;
			}
	if(!fl) for(int i = 0; i < lim; i++) f[i].x /= lim, f[i].y /= lim;
}
void Mul(int* f, int* g, int* h, int n) {
	static cp f0[MAXN], f1[MAXN], g0[MAXN], g1[MAXN];
	int lim = getlim(n);
	for(int i = 0; i < n; i++) f0[i].x = f[i] >> 15, f0[i].y = f[i] & 32767;
	for(int i = 0; i < n; i++) g0[i].x = g[i] >> 15, g0[i].y = g[i] & 32767;
	for(int i = n; i < lim; i++) f0[i] = (cp){0, 0};
	for(int i = n; i < lim; i++) g0[i] = (cp){0, 0};
	FFT(f0, lim, 1); FFT(g0, lim, 1);
	for(int i = 0; i < lim; i++) {
		f1[i] = f0[i ? lim - i : 0], f1[i].y *= -1;
		g1[i] = g0[i ? lim - i : 0], g1[i].y *= -1;
	}
	for(int i = 0; i < lim; i++) {
		cp a = (f0[i] + f1[i]) * 0.5;		//f0
		cp b = (f1[i] - f0[i]) * 0.5 * I;	//f1
		cp c = (g0[i] + g1[i]) * 0.5;		//g0
		cp d = (g1[i] - g0[i]) * 0.5 * I;	//g1
		f0[i] = a * c + I * (a * d + b * c);
		g0[i] = b * d;
	}
	FFT(f0, lim, 0); FFT(g0, lim, 0);
	for(int i = 0; i < n; i++)
		h[i] = (1ll * tran(f0[i].x) * (1 << 30) + 1ll * tran(f0[i].y) * (1 << 15) % P + tran(g0[i].x)) % P;
	clr(h, n, lim);
}

多项式乘法逆

P4238 【模板】多项式乘法逆

已知 \(F(x)\), 要求 \(G(x)\), 使得 \(F(x)G(x)\equiv 1\pmod {x^n}\)

保证 \(F[0]\ne 0\). 系数对 998244353取模.

我们采用倍增的思想.

假设我们已知多项式 \(G_0(x)\) 满足:

\[ G_0(x)F(x)\equiv 1\pmod{x^{n/2}} \] 又 \[ G(x)F(x)\equiv 1\pmod {x^{n/2}} \] 则由于 \(F(x)\not\equiv 0\pmod{x^{n/2}}\), 作差得

\[ G(x)-G_0(x)\equiv 0\pmod{x^{n/2}} \] 平方得

\[ G^2(x)-2G_0(x)G(x)+G_0^2(x)\equiv 0\pmod{x^n} \] 则 \[ G^2(x)\equiv 2G_0(x)G(x)-G_0^2(x)\pmod {x^n} \] 两边同乘 \(F(x)\) ,由于 \(G(x)F(x)\equiv 1\pmod {x^n}\), 得

\[ G(x)\equiv 2G_0(x)-F(x)G_0^2(x)\equiv G_0(x)[2-F(x)G_0(x)]\pmod {x^n} \] 这就是递归计算式了. 使用NTT转成点值进行计算.

最后的终止条件为 \(F[0]\equiv (G[0])^{-1}\).

多项式有逆的充要条件是常数项有乘法逆.

时间复杂度:

\[ T(n)=T(n/2)+O(n\log n)\\ T(n)=O(n\log n) \]

检验方法:

求 \(1-x\) 的逆为 \(\dfrac 1 {1-x}=1+x+x^2+\cdots\)

void Inv(int* f, int* g, int n) {
	static int a[MAXN];
	if(n == 1) {g[0] = qpow(f[0], P-2); return ;}
	Inv(f, g, (n+1)>>1); int lim = getlim(n);
	clr(g, (n+1)>>1, lim);
	cpy(f, a, n); clr(a, n, lim);
	NTT(a, lim, 1); NTT(g, lim, 1);
	for(int i = 0; i < lim; i++) g[i] = mul(g[i], mns(2, mul(a[i], g[i])));
	NTT(g, lim, 0); 
}
/*
int n, f[MAXN], g[MAXN];
void test() {
	n = read();
	for(int i = 0; i < n; i++) f[i] = read();
	Inv(f, g, n);
	for(int i = 0; i < n; i++) printf("%d ", g[i]);
}
Input:
5
1 998244352 0 0 0
Output:
1 1 1 1 1
*/

多项式除法 & 取模

P4512 【模板】多项式除法

已知 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次多项式 \(G(x)\) , 求 \(Q(x)\), \(R(x)\) 满足

\(Q(x)\) 的度数为 \(n-m\), \(R(x)\) 的度数 \(<m\)

\(F(x)=Q(x)G(x)+R(x)\)

对 998244353 取模

不妨就认为 \(\mathrm{deg}\ R(x)=m-1\), 最高位用 \(0\) 补齐.

定义 \(A_T(x)=x^{\mathrm{deg}\ A(x)}A(\dfrac 1 x)\), 可以发现 \(A_T(x)\) 就是把 \(A(x)\) 的系数反过来的多项式.

故有 \[ F(\dfrac 1x)=Q(\dfrac 1x)G(\dfrac 1 x)+R(\dfrac 1x) \\ x^nF(\dfrac 1x)=x^{n-m}Q(\dfrac 1x)\cdot x^mG(\dfrac 1x)+x^{n-m+1}\cdot x^{m-1}R(\dfrac 1x) \\ F_T(x)=Q_T(x)G_T(x)+x^{n-m+1}R_T(x) \\ F_T(X)=Q_T(x)G_T(x) \pmod{x^{n-m+1}} \\ Q_T(x)=\dfrac {F_T(x)}{G_T(x)}\pmod {x^{n-m+1}} \] 又由于 \(\mathrm{deg}\ Q(x)=n-m<n-m+1\), 所以得到的结果是正确的.

然后 \(R(x)=F(x)-Q(x)G(x)\).

（特别注意: 代码中的 \(n,m\) 指的是度数而不是界！！！）

检验方法: 除一除就好了嘛!

void Div(int* f, int* g, int* q, int* r, int n, int m) {
	static int a[MAXN], b[MAXN];
	int lim = getlim(n+1);
	cpy(g, a, m+1); reverse(a, a+m+1); clr(a, m+1, lim);
	Inv(a, b, n-m+1);
	cpy(f, a, n+1); reverse(a, a+n+1); clr(a, n+1, lim);
	Mul(a, b, q, n-m+1); reverse(q, q+n-m+1); clr(q, n-m+1, lim);
	Mul(g, q, a, m);
	for(int i = 0; i < m; i++) r[i] = mns(f[i], a[i]);
	clr(r, m, lim);
}
/*
int n, m, f[MAXN], g[MAXN], q[MAXN], r[MAXN];
void test() {
	n = read(); m = read();
	for(int i = 0; i <= n; i++) f[i] = read();
	for(int i = 0; i <= m; i++) g[i] = read();
	Div(f, g, q, r, n, m);
	for(int i = 0; i <= n-m; i++) printf("%d ", q[i]);
	printf("\n");
	for(int i = 0; i <= m-1; i++) printf("%d ", r[i]);
	printf("\n");

}
Input:
2 1
2 2 1
1 1
Output:
1 1
1
*/

多项式ln

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

已知 \(F(x)\), 求 \(G(x)\equiv\ln F(x)\pmod{x^n}\).

保证 \(F[0]=1\). 对 998244353 取模

这怎么求呢？我们可以先求微, 得 \[ dG(x)\equiv d\ln F(x)\equiv \dfrac{F'(x)dx}{F(x)}\pmod{x^{n}} \] 再两边积分得 \[ G(x)\equiv\ln F(x)\equiv\int_0^x\dfrac{F'(x)dx}{F(x)}\pmod {x^n} \] 而 \(G[0]\equiv \ln F[0]\).

当 \(F[0]=1\) 时, \(G[0]\equiv0\);

当 \(F[0]\ne 1\) 时, \(G[0]\equiv \ln a(a\ne 1\wedge a\in\mathbb N_+)\not\in \mathbb Q\), 在模意义下没有意义.

这也是题目要保证 \(F[0]=1\) 的原因.

时间复杂度为 \(O(n\log n)\)

检验方法:

求 \(1-x\) 的对数为 \(-\sum_{k\ge 1}\dfrac {x^k} k\)

void Deriv(int* f, int* g, int n) {
	for(int i = 1; i < n; i++) g[i-1] = mul(f[i], i);
	g[n-1] = 0;
}
void Integ(int* f, int* g, int n) {
	for(int i = 1; i < n; i++) g[i] = mul(f[i-1], inv[i]);
	g[0] = 0;
}
void Ln(int* f, int* g, int n) {
	static int a[MAXN], b[MAXN];
	int lim = getlim(n);
	Deriv(f, a, n); clr(a, n, lim);
	Inv(f, b, n); clr(b, n, lim);
	NTT(a, lim, 1); NTT(b, lim, 1);
	for(int i = 0; i < lim; i++) a[i] = mul(a[i], b[i]);
	NTT(a, lim, 0); Integ(a, g, n);
}
/*
int n, f[MAXN], g[MAXN];
void test() {
	n = read();
	for(int i = 0; i < n; i++) f[i] = read();
	Ln(f, g, n);
	for(int i = 0; i < n; i++) printf("%lld ", qpow(g[i], P-2));
}
Input:
5
1 998244352 0 0 0
Output:
0 998244352 998244351 998244350 998244349 
*/

泰勒展开

在进入 exp 前, 我们先来看看泰勒展开 &牛顿迭代

对于一个函数 \(f(x)\), 我们可以利用其在 \(x=x_0\) 处的多阶导函数的值来多项式逼近 \(f(x)\). \[ f(x)=f(x_0)+\dfrac{f'(x_0)}{1!}(x-x_0)+\dfrac{f''(x_0)}{2!}(x-x_0)^2+\cdots+\dfrac{f^{(n)}(x_0)}{n!}(x-x_0)^n+\cdots \] 怎么理解？我们不妨让右边那串为 \(g(x)\), 验证可知: \[ f(x_0)=g(x_0)\\ f'(x_0)=g'(x_0)\\ \cdots\\ f^{(n)}(x_0)=g^{(n)}(x_0)\\ \cdots \] 这就叫 "逼近".

牛顿迭代（多项式复合求根）

已知 \(G(x)\), 求 \(F(x)\)使 \(G(F(x))\equiv 0\pmod {x^n}\)

考虑使用倍增法. 假如我们已经知道了 \(G(F_0(x))\equiv 0\pmod {x^{n/2}}\).

将 \(G(F(x))\) 在 \(F(x)=F_0(x)\) 处泰勒展开, 有

\[ G(F(x))=G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))+\dfrac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+\cdots \] 又因为

\[ G(F(x))\equiv 0\pmod {x^n})\\ \mathrm{and}\\ F(x)-F_0(x)\equiv0\pmod {x^{n/2}}\\\Rightarrow(F(x)-F_0(x))^2\equiv 0\pmod {x^n} \] 代入得 \[ 0\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\pmod {x^n} \] 移项得 \[ F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n} \] 这便是递推公式.

多项式exp

已知 \(F(x)\), 求 \(G(x)\equiv\exp F(x)\pmod {x^n}\)

保证 \(F[0]=0\).

直接求并不好求, 我们不妨先取对数.

\[ \ln G(x)-F(x)\equiv0\pmod{x^n} \] 直接套用牛顿迭代, 令 \(H(X)=\ln X-F(x)\).

则求解 \[ H(G(x))\equiv 0\pmod {x^n} \] 则有

\[ G(x)\equiv G_0(x)-\dfrac{H(G_0(x))}{H'(G_0(x))}\\ \equiv G_0(x)-G_0(x)[\ln G_0(x)-F(x)]\\ \equiv G_0(x)[1-\ln G_0(x)+F(x)] \] 而又有 \(G[0]\equiv \exp{F[0]}\).

可以递归求解. 看起来很麻烦, 但时间复杂度还是 \(O(n\log n)\) （虽然比较慢）.

类似地, 由于 \(G[0]\equiv \exp F[0]\), 当 \(F[0]\ne 0\) 时在模意义下无意义, 所以题目保证了 \(F[0]=0\).

检验方法:

求 \(x\) 的指数为 \(\sum_{k\ge 0}\dfrac{x^k}{k!}\)

void Exp(int* f, int* g, int n) {
	static int a[MAXN];
	if(n == 1) {g[0] = 1; return ;}
	Exp(f, g, (n+1)>>1); int lim = getlim(n);
	clr(g, (n+1)>>1, lim);
	Ln(g, a, n); 
	for(int i = 0; i < n; i++) a[i] = mns(f[i], a[i]);
	a[0] = pls(a[0], 1); clr(a, n, lim);
	NTT(a, lim, 1); NTT(g, lim, 1);
	for(int i = 0; i < lim; i++) g[i] = mul(a[i], g[i]);
	NTT(g, lim, 0); 
}
/*
int n, f[MAXN], g[MAXN];
void test() {
	n = read();
	for(int i = 0; i < n; i++) f[i] = read();
	Exp(f, g, n);
	for(int i = 0; i < n; i++) printf("%d ", qpow(g[i], P-2));
}
Input:
5
0 1 0 0 0
Output:
1 1 2 6 24
*/

多项式快速幂 1

P5245 【模板】多项式快速幂

已知 \(F(x)\), 求 \(G(x)\equiv F^k(x)\pmod {x^n}\)

保证 \(F[0]=1\).

答案对 998244353 取模

我们根据对数恒等式, \[ F^k(x)\equiv\exp(k\ln F(x)) \] 先求\(\ln\), 再做数乘, 再求\(\exp\) 即可. 时间复杂度 \(O(n\log n)\). 常数巨大.

从这里我们也知道了, 多项式快速幂不符合费马小定理, 指数应当 \(\mod P\) 才对.

多项式快速幂 2

P5273 【模板】多项式幂函数 (加强版)

已知 \(F(x)\), 求 \(G(x)\equiv F^k(x)\), 对 998244353 取模.

如果 \(F[0]\equiv 1\), 那么直接仿照上面计算即可. 如果\(F[0]\not\equiv 1\), 则需要进行一些转化.

找到 \(F(x)\) 的第一项非零项 \(a_tx^t\) , 把其提出, 得到

\(F(x)\equiv a_tx^tH(x)\), 其中 \(H[n]\equiv\dfrac {F[n+t]}{a_k}\).

则 \(G(x)\equiv F^k(x)\equiv a_t^kx^{tk}H^k(x)\)

\(H^k(x)\) 再如上计算即可. 值得注意的是, \(H^k(x)\) 的 \(k\) 应 \(\bmod P\), 而 \(a_t^k\) 中的 \(k\) 应 \(\bmod \varphi(P)\)

代码极丑:

检验方法:

找一个多项式随便平方一下吧.

void Qpow(int* f, int* g, int n, int k1, int k2, bool out) {
    //k1 表示 k mod P, k2 表示 k mod \phi(P), out 表示 k 是否取过模. 
	static int a[MAXN];
	int shift = 0;
	if(k1 == 0) {
		clr(g, 0, n); g[0] = qpow(f[0], k2);
		return;
	}
	while(shift < n && !f[shift]) shift++;
	if((shift && out) || 1ll * shift * k1 >= n) {
		clr(g, 0, n);
		return;
	}
	int in = qpow(f[shift], P-2), t = qpow(f[shift], k2);
	for(int i = 0; i < n; i++) f[i] = i + shift < n ? mul(f[i + shift], in) : 0;
	Ln(f, a, n);
	for(int i = 0; i < n; i++) a[i] = mul(a[i], k1);
	Exp(a, f, n);
	shift *= k1;
	for(int i = 0; i < n; i++) g[i] = i - shift >= 0 ? mul(f[i - shift], t) : 0;
}
/*
char sk[MAXL];
int n, f[MAXN], g[MAXN], len, out;
ll k1, k2;
void test() {
	n = read(); scanf("%s", sk+1); len = strlen(sk+1);
	for(int i = 0; i < n; i++) f[i] = read();
	for(int i = 1; i <= len; i++) {
		k1 = k1 * 10 + int(sk[i] - '0');
		if(k1 >= P) k1 %= P, out = 1;
	}
	for(int i = 1; i <= len; i++) {
		k2 = k2 * 10 + int(sk[i] - '0');
		if(k2 >= P-1) k2 %= P-1, out = 1;
	}
	Qpow(f, g, n, k1, k2, out);
	for(int i = 0; i < n; i++) printf("%d ", g[i]);
}
Input:
10 3
0 1 1 1 0 0 0 0 0 0
Output:
0 0 0 1 3 6 7 6 3 1
*/

多项式开根

P5205 【模板】多项式开根

P5277 【模板】多项式开根（加强版）

已知 \(F(x)\), 求 \(G(x)\equiv \sqrt{F(x)}\pmod {x^n}\).

保证 \(F[0]\) 是模998244353意义下的二次剩余.

答案对 998244353 取模

你当然可以向下面讲到的多项式快速幂一样计算, 但是那样太慢.

原式化为\(G^2(x)-F(x)\equiv 0\pmod {x^n}\)

牛顿迭代可得 \[ G(x)\equiv G_0(x)-\dfrac{G^2_0(x)-F(x)}{2G_0(x)}\equiv \dfrac{G_0(x)}2+\dfrac{F(x)}{2G_0(x)} \] 递归终点是 \[ G[0]\equiv \sqrt{F[0]}\pmod p \] 这里给出加强版的代码（其中Sqrt_P::solve(n) 是求 \(n\) 的二次剩余）.

时间复杂度为 \(O(n\log n)\)

检验方法:

求一个多项式的平方根.

namespace Sqrt_P {
	int a, I2;
	struct F2 {int x, y;};
	F2 operator * (const F2& a, const F2& b) {return (F2){pls(mul(a.x, b.x), mul(mul(a.y, b.y), I2)), pls(mul(a.x, b.y), mul(a.y, b.x))};}
	F2 qpow2(F2 a, int n) {F2 ret = (F2){1, 0}; for(; n; n >>= 1, a = a * a) if(n & 1) ret = a * ret;  return ret;}
	int solve(int n) {
		n %= P;
		if(!n) return 0;
		if(qpow(n, (P-1) / 2) == P-1) return -1;
		while(1) {
			a = rand() % P;
			I2 = mns(mul(a, a), n);
			if(qpow(I2, (P-1) / 2) == P-1) break;
		}
		int ret = qpow2((F2){a, 1}, (P+1) / 2).x;
		return min(ret, P-ret);
	}
}

void Sqrt(int* f, int* g, int n) {
	static int a[MAXN], b[MAXN];
	if(n == 1) {g[0] = Sqrt_P::solve(f[0]); return ;}
	Sqrt(f, g, (n+1) >> 1); int lim = getlim(n);
	clr(g, (n+1) >> 1, lim);
	cpy(f, a, n); clr(a, n, lim);
	Inv(g, b, n); clr(b, n, lim);
	Mul(a, b, a, n);
	for(int i = 0; i < n; i++) g[i] = mul(inv[2], pls(g[i], a[i]));
	clr(g, n, lim);
}
/*
int n, f[MAXN], g[MAXN];
void test() {
	n = read();
	for(int i = 0; i < n; i++) f[i] = read();
	Sqrt(f, g, n);
	for(int i = 0; i < n; i++) printf("%d ", g[i]);
}
Input:
3
4 8 4
Output:
2 2 0
*/