多项式入门教程

( ewcommand{ d}{{ m d}})原文链接 https://www.cnblogs.com/zhouzhendong/p/polynomial.html

UPD(2020-08-27): 做了大量更新

多项式基础操作

目录

1. 多项式求逆
2. 牛顿迭代
3. 二次剩余
4. 多项式开根
5. 多项式对数函数
6. 多项式指数函数
7. 多项式快速幂
8. 多项式除法、取模
9. 多点求值与快速插值

前置技能

1. 快速傅里叶变换(FFT) 和 快速数论变换(NTT)
   https://www.cnblogs.com/zhouzhendong/p/Fast-Fourier-Transform.html

关于代码

下面这个链接是博主给出的示例代码。如果在各个部分的示例代码中有未定义的变量名或者函数名，请到代码里找。

https://loj.ac/submission/872547

为了方便，这里也给出的一些重要的变量名/函数名的含义

mod = 998244353;
Pow,Add,Del等自行感受;
Fac[i] = i!;
Inv[i] = 1 / Fac[i];
Iv[i] = 1 / i;
typedef vector <int> vi;//代码用来表示多项式类型
vi fix(vi a,int n);//将a的size调整为n（给a末尾补0或者删除末尾直至a的size为n）
vi operator + (vi a,vi b);//多项式加法
vi operator - (vi a,vi b);//多项式减法
vi operator * (vi a,int b);//a*b
vi operator * (vi a,vi b);//a*b
vi polyinv(vi a);//多项式求逆
vi Der(vi a);//多项式求导
vi Int(vi a);//多项式积分
vi polyln(vi a);//多项式ln
vi polyexp(vi a);//多项式exp
vi operator / (vi a,vi b);//多项式除法
vi operator % (vi a,vi b);//多项式取模
vi polypow(vi a,int k);//多项式快速幂
int Sqrt(int x);//求x在模mod意义下的二次剩余
vi polysqrt(vi a);//多项式开根
namespace eval_inter;//多点求值和快速插值

多项式求导和积分

vi Der(vi a){
	For(i,1,(int)a.size()-1)
		a[i-1]=(LL)a[i]*i%mod;
	a.pop_back();
	return a;
}
vi Int(vi a){
	a.pb(0);
	Fod(i,(int)a.size()-1,1)
		a[i]=(LL)a[i-1]*Iv[i]%mod;
	a[0]=0;
	return a;
}

一些记号

为了方便，我们先声明以下记号的含义：

多项式 (A(x)) 和其对应的 (i) 次项系数 (a_i)

[A(x) = sum_{i=0}^{n-1}a_ix^i ]

多项式积分和求导

[int A(x) d x\ A'(x) = frac{ d}{ d x}A(x) ]

以及如下定义：

[B_{n}(x) = B(x) mod {x^n} ]

为了方便，在后面的推导中可能会简化多项式的表达，将 (A(x)) 简写为 (A)，以及将 (A_n(x)) 简写为 (A_n)。同样地，为了方便，在采用这种简写方式时，我可能会将 对 (x^?) 取模的同余式 写成等式，但是一般来说意义仍是明确的。

多项式求逆

给定多项式 (A(x)) ，求出多项式 (B(x)) ，使得

[A(x)B(x) equiv 1 pmod{x^n} ]

其中多项式 (A(x))、(B(x)) 只需要保留前 (n) 项（也就是 (x^0) 到 (x^{n-1})）

做法

假设我们已经得到了 (B_n)。

[AB_n equiv 1 pmod{x^{n}} ]

则

[(AB_n-1)^2equiv 0 pmod {x^{2n}}\ A^2B_n^2-2AB_n+1equiv 0 pmod {x^{2n}}\ 1equiv 2AB_n-A^2B_n^2 pmod {x^{2n}} ]

左右同除以 (A) ，得

[B_{2n} = 2B_n - AB_n^2 ]

因为我们知道 (B_1(x) = a_0^{-1})，所以不难通过倍增法得到 (B_n)。

时间复杂度

[T(n) = T(n/2) + O(nlog n) = O(nlog n) ]

示例代码

vi polyinv(vi a){
	int n=(int)a.size();
	if (n==1)
		return (vi){Pow(a[0],mod-2)};
	vi b=polyinv(fix(a,(n+1)/2));
	return fix(b*2-b*b*a,n);
}

这里有一个卡常技巧：如果我们已经知道了 (B_n)，那么我们在做 b*b*a 时可以不关心前 (n) 位的系数，于是我们可以直接利用 dft 是循环卷积的性质来缩小 dft 长度。实现如下：

vi polyinv(vi a){
	if (a.size()==1)
		return vi(1,Pow(a[0],mod-2));
	int n=a.size();
	vi b=polyinv(fix(a,(n+1)/2));
	int len=1<<(int)(log(n+(int)b.size())/log(2)+1);
	vi c=fix(b,len),d=fix(a,len);
	fft::init(len),FFT(&c[0],len,1),FFT(&d[0],len,1);
	For(i,0,len-1)
		c[i]=(LL)c[i]*c[i]%mod*d[i]%mod;
	FFT(&c[0],len,-1);
	b.resize(n);
	For(i,(n+1)/2,n-1)
		Del(b[i],c[i]);
	return b;
}

多项式牛顿迭代

设有可导函数 (F(G)) ，其中 (G) 是一个多项式。我们想要快速求其零点。

将 (F(G)) 泰勒展开得到：

[F(G) = F(H) + (G-H)F'(H) + (G-H)^ 2frac{F''(H)}{2!}+cdots ]

假设 (F(A)=0)，那么我们来考虑怎么通过上式来由 (A_n) 导出 (A_{2n})。

[0 = F(A) = F(A_n) + (A-A_n)F'(A_n)+(A-A_n)^2frac{F''(A_n)}{2!}+cdots ]

由于 (A - A_n) 的 (0cdots n-1) 次项系数为 (0)，所以

[0 = F(A_n) + (A-A_n)F'(A_n) pmod {x ^ {2n}}\ F'(A_n)A = A_nF'(A_n) - F(A_n)pmod {x ^ {2n}}\ A_{2n} = A_n - cfrac{F(A_n)}{F'(A_n)} pmod {x ^ {2n}} ]

这启发我们使用倍增法来求一些关于多项式的函数。

多项式开根

P.S. 这里需要求模意义下二次剩余，可以采用 BSGS 等算法解决。

给定多项式 (A(x)) ，求出多项式 (B(x)) 使得

[B^2(x)equiv A(x) pmod{x^n} ]

做法

我们考虑套用牛顿迭代：

令 (F(B) = B^2 - A)，我们要求的便是 (F(B)) 的零点。则：

[B_{2n} = B_n - cfrac{F(B_n)}{F'(B_n)}\ = B_n - cfrac{B_n^2 - A}{2B_n}\ = frac 1 2(B_n + frac{A}{B_n}) ]

因为我们知道

[B_1 = sqrt{a_0} ]

所以这里我们可以使用任意求二次剩余的算法和倍增法解决多项式开根问题。

时间复杂度

[T(n)=T(n/2)+O(nlog n)=O(nlog n) ]

示例代码

vi polysqrt(vi a){
	if (a.size()==1)
		return vi(1,Sqrt(a[0]));
	int n=a.size();
	vi b=fix(polysqrt(fix(a,(n+1)/2)),n);
	return fix((b+a*polyinv(b))*((mod+1)/2),n);
}

多项式对数函数

给定多项式 (F(x))，求 (ln(F(x))pmod{x^n}) 。

做法

该问题的解法由下式直接得出：

[ln(F(x)) = int cfrac{F'(x)}{F(x)} d x ]

时间复杂度 $$O(nlog n)$$。

注意点

(F(x)) 的常数项必须是 (1) 。

否则设 (G(x)=kF(x)) ，则 (ln(G(x))=ln(F(x))+ln(k)) ，其中 (ln(k)) 难以用模意义下的数来表示。

容易得知 (ln(F(x))) 的常数项是 (0)，也就是说，(ln(F(x))) 相对于 (F(x)) 损失了常数项，而在高位上没有损失。

示例代码

vi polyln(vi a){
	return Int(fix(Der(a)*polyinv(a),(int)a.size()-1));
}

多项式指数函数

给定多项式 (F(x)) ，求 (e^{F(x)} pmod{x^n}) 。

做法

首先，设

[e^{F(x)}=G(x) pmod{x^n} ]

则

[ln(G(x))-F(x)=0 pmod{x^n} ]

设函数 (Q(x)=ln(x) - F) ，利用牛顿迭代得到：

[G_{2n} = G_n - cfrac{Q(G_n)}{Q'(G_n)}\ = G_n - cfrac{ln(G_n) - F}{frac{1}{G_n}}\ = G_n(1 - ln(G_n) + F) ]

时间复杂度仍然是

[O(nlog n) ]

示例代码

vi polyexp(vi a){
	if (a.size()==1)
		return vi(1,1);
	int n=a.size();
	vi b=polyexp(fix(a,(n+1)/2));
	return fix(b+b*(a-polyln(fix(b,n))),n);
}

多项式快速幂

给定多项式 (F(x)) ，求 (F^k(x)pmod{x^n})。

做法

令 (F(x)=bx^aG(x)) ，其中 (a,b) 为常数，(G(x)) 的常数项为 1 。则

[F^k(x) = b^kx^{ak}G^k(x) = b^kx^{ak}e^{kln(G(x))} ]

时间复杂度

[O(nlog n) ]

注意点

若 (k) 极大（如洛谷的"多项式快速幂加强版"），则我们需要计算三个值：（为了表达清晰，这里用 (p) 表示模数）

[k_1 = k mod varphi(p)\ k_2 = k mod p\ k_3 = min(k,{ m n}) ]

那么

[F^k(x)equiv b^{k_1}x^{ak_3} e^{k_2ln(G(x))} pmod {x^n} pmod {p} ]

示例代码

若保证了 (f_0 = 1)，则

vi polypow(vi a,int k){
	return polyexp(polyln(a)*k);
}

否则：

vi polypow(vi a,int k){
//	return polyexp(polyln(a)*k);
	int p=0;
	while (p<(int)a.size()&&!a[p])
		p++;
	if ((LL)p*k>=(int)a.size())
		return vi();
	vi b(a.begin()+p,a.end());
	int coef=b[0],inv=Pow(coef,mod-2),powcoef=Pow(coef,k);
	b=polyexp(polyln(b*inv)*k)*powcoef;
	vi c(p*k+(int)b.size());
	For(i,0,(int)b.size()-1)
		c[i+p*k]=b[i];
	return fix(c,a.size());
}

多项式除法

给定多项式 (F(x),G(x)) ，求多项式 (Q(x))，使得

[F(x)=G(x)Q(x)+R(x) ]

其中 (F(x),G(x)) 分别是 (n,m) 次多项式。

(Q(x),R(x)) 分别是 (n-m+1,m-1) 次多项式。

做法

定义 (F^R(x)) 表示多项式 (F(x)) 系数翻转之后得到的结果。设 (F(x)) 最高次项为 (x^{n-1}) ，则

[F^R(x)=x^{n-1}F(frac 1x) ]

于是可得

[F(x)=G(x)Q(x)+R(x)\ F(frac 1x)=G(frac 1x)Q(frac 1x)+R(frac 1x)\ x^{n-1}F(frac 1x)=x^{m-1}G(frac 1x)x^{n-m}Q(frac 1x)+x^{n-m+1}x^{m-2}R(frac 1x)\ F^R(x)=G^R(x)Q^R(x)+x^{n-m+1}R^R(x) ]

又因为 (Q^R(x)) 最高次项为 (x^{n-m})，所以

[F^R(x)=G^R(x)Q^R(x)pmod{x^{n-m+1}}\ frac{F^R(x)}{G^R(x)}=Q^R(x)pmod{x^{n-m+1}} ]

时间复杂度

[O(nlog n) ]

示例代码

这里需要注意的是 (m>n) 要特判。

vi operator / (vi a,vi b){
	int n=a.size(),m=b.size();
	if (n<m)
		return vi();
	reverse(a.begin(),a.end());
	reverse(b.begin(),b.end());
	a=fix(fix(a,n-m+1)*polyinv(fix(b,n-m+1)),n-m+1);
	reverse(a.begin(),a.end());
	return a;
}

多项式取模

给定多项式 (F(x),G(x)) ，求多项式 (R(x))，使得

[F(x)=G(x)Q(x)+R(x) ]

其中 (F(x),G(x)) 分别是 (n,m) 次多项式。

(Q(x),R(x)) 分别是 (n-m+1,m-1) 次多项式。

做法

[R(x)=F(x)-G(x)Q(x) ]

时间复杂度

[O(nlog n) ]

示例代码

vi operator % (vi a,vi b){
	return fix(a-a/b*b,(int)b.size()-1);
}

多点求值与快速插值

多点求值

给定最高次项为 (x^{m-1}) 的函数 (F(x)) ，以及 (n) 个参数 (x_{1cdots n}) ，求

[F(x_1),F(x_2),cdots,F(x_n) ]

做法

[F(x_i)=(F(x)mod (x-x_i))[0] ]

即多项式 (F(x)mod (x-x_i)) 的常数项。
设

[L(x) = prod_{1leq ileq lfloor frac n2 floor} (x-x_i)\ R(x) = prod_{lfloor frac n2 floor<ileq n} (x-x_i) ]

则对于 (ileq n/2) ，(F(x_i)=(Fmod L)(x_i)) ；

对于 (i> n/2) ，(F(x_i)=(Fmod R)(x_i)) ；

先预处理得到所有的 (L(x)) 和 (R(x)) ，然后再分治求出答案即可。

时间复杂度

[T(n)=2T(n/2)+O(nlog n)=O(nlog^2 n) ]

快速插值

给定 (n) 对 ((x_i,y_i)) ，求最高次项为 (x^{n-1}) 的多项式 (F(x)) 满足

[forall 1leq ileq n, f(x_i) = y_i ]

做法

考虑拉格朗日插值法

[F(x) = sum_{i=1}^{n}frac {prod_{j e i}(x-x_j)} {prod_{j e i}(x_i-x_j)}y_i ]

令 (M(x) = prod_{i=1}^{n}(x-x_i))

[M_i(x) = M(x)/(x-x_i) ]

令

[t_i=frac{y_i}{prod_{j eq i}(x_i-x_j)}=frac{y_i}{M_i(x_i)} ]

根据洛必达法则，当 (x ightarrow x_i) 时，

[M_i(x)=frac{M(x)}{x-x_i} ightarrow M'(x) ]

故

[t_i=y_i/M'(x_i) ]

设

[L(x) = prod_{1leq ileq lfloor frac n2 floor} (x-x_i)\ R(x) = prod_{lfloor frac n2 floor<ileq n} (x-x_i) ]

则

[F(x) = sum_{i=1}^{n}t_iprod_{j e i}(x-x_j)\ =left(sum_{i=1}^{lfloor frac n2 floor}t_iprod_{j eq i}(x-x_j) ight)R(x)+left(sum_{i=lfloor frac n2 floor+1}^{n}t_iprod_{j eq i}(x-x_j) ight)L(x) ]

先预处理得到所有的 (L(x)) 和 (R(x)) ，然后再分治求出答案即可。

时间复杂度

[T(n)=2T(n/2)+O(nlog n)=O(nlog^2 n) ]

示例代码

namespace eval_inter{
	int n;
	vi f,x,y,prod[N*4];
	void getprod(int rt,int L,int R){
		if (L==R){
			prod[rt]={Del(-x[L]),1};
			return;
		}
		int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;
		getprod(ls,L,mid);
		getprod(rs,mid+1,R);
		prod[rt]=prod[ls]*prod[rs];
	}
	void gety(int rt,int L,int R,vi f){
		f=fix(f%prod[rt],(int)prod[rt].size()-1);
		if (L==R)
			return (void)(y[L]=f[0]);
		int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;
		gety(ls,L,mid,f);
		gety(rs,mid+1,R,f);
	}
	vi eval(vi _f,vi _x){
		n=_x.size();
		if (!n)
			return vi();
		f=_f,x=_x,y.resize(n);
		getprod(1,0,n-1);
		gety(1,0,n-1,f);
		return y;
	}
	vi getf(int rt,int L,int R){
		if (L==R)
			return vi(1,y[L]);
		int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;
		return getf(ls,L,mid)*prod[rs]+getf(rs,mid+1,R)*prod[ls];
	}
	vi inter(vi _x,vi _y){
		n=_x.size();
		if (!n)
			return vi();
		x=_x,y.resize(n);
		getprod(1,0,n-1);
		vi M=Der(prod[1]);
		gety(1,0,n-1,M);
		For(i,0,n-1)
			y[i]=(LL)_y[i]*Pow(y[i],mod-2)%mod;
		return getf(1,0,n-1);
	}
}

模板题

LOJ150 挑战多项式
https://loj.ac/problem/150

NFLSOJ里的
https://acm.nflsoj.com/problems/template