你不可不读的多项式教程
写在前面
看懂这篇博客,需要掌握
- 一般程度的代数知识
- 入门程度的复数知识
多项式乘法与卷积
对于给定的两个多项式
我们称
为\(F,G\)的卷积。不难证明,这就是多项式乘法,即
(简单地说,\(F,G\)的乘积中次数为\(k\)的项系数为
计算\(F,G\)的卷积显然有\(O(n^2)\)的朴素思路,但是对于较高的数据范围这样的时间复杂度是无法接受的。
多项式的点值表达
除了常见的系数表达以外,多项式还可以用点值表达来表示。
对于已知次数为\(n\)的多项式\(F(x)\),我们可以用
来唯一确定其系数表达式。证明可以由代数基本定理推出。
多项式的点值表达相对于系数表达有一个优点:在计算卷积时时间复杂度缩减为\(O(n)\)。更准确地说,有
但是将系数表达转换为点值表达的过程仍然是\(O(n^2)\)的,于是我们想到对这个过程进行优化。
单位根
快速傅里叶变换能够在\(O(nlogn)\)的时间内解决系数表达到点值表达的过程,是因为引入了单位根的概念。称方程
的复数根为\(n\)次单位根。这样的单位根一共有\(n\)个,记为
两条基本性质如下:
\((1)\)的证明是显然的:复数根在复平面上的表示即为单位圆的半径。
\((2)\)可以通过
计算得出。本条性质还可以推出另外两条常用推论:
此外,单位根还有两个重要的性质
\((3)\)和\((4)\)的证明都可以通过直接套用\((*)\)完成。
快速傅里叶变换(FFT)
快速傅里叶变换最核心的环节(之一)就是多项式的系数表达与点值表达的转换,可以采用分治的方式解完成。(此处先介绍递归写法)
对于多项式(为什么为最高次为\(n-1\)次下面解释)
有
不妨设
那么有
(此时可以看出为什么此前将\(F\)最高次设为\(n-1\):美观)
现在我们要求\(F\)的点值表达,于是代入单位根,也就是
同理,可以得到
也就是说,如果我们求出了\(L,R\)在\(\omega_{n/2}^0,...,\omega_{n/2}^{n/2-1}\)的点值表达,我们就可以得到\(F\)在\(\omega_{n}^0,...,\omega_{n}^n-1\)的点值表达。
由于每一次长度减半,所以可以在\(O(nlogn)\)的时间内求出多项式\(F\)的点值表示。
对于多项式的点值表达到系数表达的转换,可以推出类似的公式。
\(FFT\)的递归实现如下:
const double pi=acos(-1);
typedef complex<double> cx;
void fft (vector<cx> x,int n,int type) {
if (n==1) return;
vector<cx> l,r;
for (register int i=0; i<n; i+=2)
l.push_back(x[i]),r.push_back(x[i+1]);
fft(l,n>>1,type),fft(r,n>>1,type);
cx wn(cos(2*pi/n),type*sin(2*pi/n)),w(1,0);
for (register int i=0; i<(n>>1); ++i,w*=wn) {
x[i]=l[i]+w*r[i],x[i+(n>>1)]=l[i]-w*r[i];
}
}
在得到了递归实现以后,我们发现其常数较大,因此考虑是否存在非递归实现。观察递归过程中的数组下标,其变换过程如下
观察其二进制表示
不难发现每一项递归之后的最终位置为其二进制表示翻转的结果。
那么我们可以通过预处理一个\(R\)数组来直接得到最终位置。不难发现,每次递归实质上是按照对应位进行排序,于是我们可以得到实现:
for (register int i=0; i<n; ++i) {
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
}
那么根据\(R\),我们可以写出\(FFT\)的非递归实现:
void fft (vector<cx> x,int n,int type) {
for (register int i=0; i<n; ++i) // 求出递归后的数组
if (i<R[i]) swap(x[i],x[R[i]]); // 如果不判断就会交换两次,等于没交换
for (register int i=1; i<n; i<<=1) { // 枚举区间长度的一半(当然也可以枚举区间长度)
cx wn(cos(pi/i),type*sin(pi/i)); // 此处消掉了公式里的2
for (register int j=0; j<n; j+=(i<<1)) { // 枚举区间,对每一个区间做修改
cx w(1,0);
for (register int k=0; k<i; ++k,w*=wn) { // 枚举区间左半侧,同时修改右半侧
cx t=x[j+k],s=x[j+k+i]*w;
x[j+k]=t+s,x[j+k+i]=t-s;
}
}
}
}
快速数论变换(NTT)
\(NTT\)与\(FFT\)的唯一区别在于将单位根改为了原根,从而有效避免了复数取模运算中的精度丢失。原根满足此前提到的单位根的性质,其细节不做赘述。实现如下:
void ntt (vector<int> x,int n,int type) {
for (register int i=0; i<n; ++i)
if (i<R[i]) swap(x[i],x[R[i]]);
for (register int i=1; i<n; i<<=1) {
int wn=ksm(type==1?G:invG,(MOD-1)/(i<<1));
for (register int j=0; j<n; j+=(i<<1)) {
int w=1;
for (register int k=0; k<i; ++k,w*=wn) {
int t=x[j+k],s=(x[j+k+i]*w)%MOD;
x[j+k]=(t+s)%MOD,x[j+k+i]=(t-s)%MOD;
}
}
}
}
\(G\)为模数的原根,\(invG\)为该原根在模数意义下的逆元。对于很大一部分题目(准确地说,对于满足\(MOD=a*2^b+1\)的题目),我们能够直接由暴力求出其原根,对于另一类任意模数题目,则需要利用中国剩余定理进行求解。暴力实现如下:
void getG () {
int tmp=MOD;
for (register int i=2; i*i<=MOD; ++i) {
if (tmp%i==0) {
fac[++cnt]=i;
while (tmp%i==0) tmp/=i;
}
}
if (tmp!=1) fac[++cnt]=tmp;
for (register int i=2; i<MOD; ++i) {
int flag=1;
for (register int j=1; j<=cnt; ++j) {
if (ksm(i,(MOD-1)/fac[j])==1) {
flag=0;
break;
}
}
if (flag) {
write(i);
break;
}
}
}
任意模数NTT暂时先坑着,回头写完多项式求逆和开根再填。
多项式求逆
对于\(n-1\)次多项式\(F\),其逆元为满足
的多项式\(G\)。
假设我们已经求出满足
的\(G^{'}\),由于显然有
那么有
如果\(F(x)\equiv0(mod~x^{\frac{n}{2}})\),结果是显然的。那么考虑
等式两侧同时平方,得到
两侧同乘\(A(x)\),有
所以
于是说明了我们可以用模\(x^{\frac{n}{2}}\)意义下逆元推出模\(x^n\)的意义下的逆元。
非递归的多项式求逆实现如下:
void inv (vector<int> x,vector<int> y,int n) {
b[0]=ksm(a[0],MOD-2);
int len;
for (len=1; len<(n<<1); len<<=1) {
for (register int i=0; i<len; ++i) A[i]=x[i],B[i]=y[i];
for (register int i=0; i<len<<1; ++i) R[i]=(R[i>>1]>>1)|((i&1)?len:0);
ntt(A,len<<1,1),ntt(B,len<<1,1);
for (register int i=0; i<len<<1; ++i) y[i]=((2-A[i]*B[i]%MOD)*B[i]%MOD+MOD)%MOD;
ntt(y,len<<1,1);
for (register int i=len; i<len<<1; ++i) y[i]=0;
}
for (register int i=0; i<len; ++i) A[i]=B[i]=0;
for (register int i=len; i<len<<1; ++i) y[i]=0;
}
多项式开根
对于\(n-1\)次多项式\(F\),其开根结果为满足
的多项式\(G\)。
假设我们已经求出满足
的\(G^{'}\),由于显然有
那么有
可以化简得到
所以
两侧同时平方,得到
化简得
所以
于是说明了我们可以用模\(x^{\frac{n}{2}}\)意义下开根结果推出模\(x^n\)的意义下的开根结果。