多項式入門教程

\(\newcommand{\rd}{{\rm 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) \rd x\\ A'(x) = \frac{\rd}{\rd x}A(x) \]

以及如下定義：

\[B_{n}(x) = B(x) \bmod {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)^2\equiv 0 \pmod {x^{2n}}\\ A^2B_n^2-2AB_n+1\equiv 0 \pmod {x^{2n}}\\ 1\equiv 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(n\log n) = O(n\log 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)^ 2\frac{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)^2\frac{F''(A_n)}{2!}+\cdots \]

由於 \(A - A_n\) 的 \(0\cdots 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(n\log n)=O(n\log 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)} \rd x \]

時間復雜度 $$O(n\log 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(n\log 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^{k\ln(G(x))} \]

時間復雜度

\[O(n\log n) \]

注意點

若 \(k\) 極大（如洛谷的"多項式快速冪加強版"），則我們需要計算三個值：（為了表達清晰，這里用 \(p\) 表示模數）

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

那么

\[F^k(x)\equiv b^{k_1}x^{ak_3} e^{k_2\ln(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(n\log 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(n\log n) \]

示例代碼

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

多點求值與快速插值

多點求值

給定最高次項為 \(x^{m-1}\) 的函數 \(F(x)\) ，以及 \(n\) 個參數 \(x_{1\cdots n}\) ，求

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

做法

\[F(x_i)=(F(x)\bmod (x-x_i))[0] \]

即多項式 \(F(x)\bmod (x-x_i)\) 的常數項。
設

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

則對於 \(i\leq n/2\) ，\(F(x_i)=(F\bmod L)(x_i)\) ；

對於 \(i> n/2\) ，\(F(x_i)=(F\bmod R)(x_i)\) ；

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

時間復雜度

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

快速插值

給定 \(n\) 對 \((x_i,y_i)\) ，求最高次項為 \(x^{n-1}\) 的多項式 \(F(x)\) 滿足

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

做法

考慮拉格朗日插值法

\[F(x) = \sum_{i=1}^{n}\frac {\prod_{j\ne i}(x-x_j)} {\prod_{j\ne 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\neq i}(x_i-x_j)}=\frac{y_i}{M_i(x_i)} \]

根據洛必達法則，當 \(x\rightarrow x_i\) 時，

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

故

\[t_i=y_i/M'(x_i) \]

設

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

則

\[F(x) = \sum_{i=1}^{n}t_i\prod_{j\ne i}(x-x_j)\\ =\left(\sum_{i=1}^{\lfloor \frac n2 \rfloor}t_i\prod_{j\neq i}(x-x_j)\right)R(x)+\left(\sum_{i=\lfloor \frac n2 \rfloor+1}^{n}t_i\prod_{j\neq i}(x-x_j)\right)L(x) \]

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

時間復雜度

\[T(n)=2T(n/2)+O(n\log n)=O(n\log^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