写在前面
看懂这篇博客,需要掌握
- 一般程度的代数知识
- 入门程度的复数知识
多项式乘法与卷积
对于给定的两个多项式
我们称
为(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)的意义下的开根结果。