FFT和NTT学习

MCAdam · 2020-04-06 17:25:13 · 个人记录

原来的博客 非常trival

复数基础

• 欧拉公式：

首先有：

把ix带入到e^{x}中得到：\displaystyle e^{ix}=\sum_{n=0}^{\infty}\frac{(ix)^n}{n!}=\sum_{n=0}^{\infty}\frac{(-1)^nx^{2n}}{(2n)!}+\sum_{n=0}^{\infty}\frac{(-1)^nix^{2n+1}}{(2n+1)!}

可以发现欧拉公式：

令x=\pi得到：

• 极坐标

对于任何复数我们都可以用z=a+bi来表示，也就是复平面上的一个点

极坐标：用到原点的距离与角度（和x轴形成的角度）来表示(r,\theta)

• 因为e^{ix}=\cos x+i\sin x，那么在复平面上就是(\cos x,\sin x)，也就是一个圆！

也就是说我们可以用e^{ix}来表示复平面上单位圆的任意一个点，为了更直观的几何意义可以写作：e^{i\theta}（其中\theta是弧度制下的角）

那么对于复平面上任意一个点：(a,b)就可以用极坐标来表示为：re^{i\theta}

这样我们就能非常直观的理解复数乘法：r_1e^{i\theta_1}\times r_2e^{i\theta_2}=(r1\times r2)e^{i(\theta_1+\theta_2)}，即模长相乘俯角相加（除法同理）

单位根

若w^n=1，那么称w为一个n次单位根。根据代数基本定理（不会证），n次单位根有且仅有n个

如果直接考虑w的取值，那么就是在复数域上考虑；当然我们也可以把它放在同余系下考虑，那么就涉及到原根的知识

非常地不严谨：单位根就是分布在某个域上所谓的“单位圆”，并且将整个“单位圆”平均分的位置（应该是分圆多项式？小蓝书上看到过）

为了方便表述，记w_n^k表示n次单位根的第k个，其中k\in [0,n-1]。那么单位根满足以下几个性质：

\displaystyle w_n^i\times w_n^j=w_n^{i+j} \displaystyle w_n^{i}=w_n^{i\bmod n} \displaystyle w_{kn}^{ki}=w_n^i

（最后一个性质，直接用上面的几何意义来理解，即把一个圆分成n份取其中的\dfrac{k}{n}，等价于把一个圆分成pn份，取其中的\dfrac{pk}{pn}）

上面说到：e^{i\theta}=(\sin \theta,\cos \theta)，那么e^{i\times 2\pi}=1

那么w_n^k=e^{i\times 2\pi\times \frac{k}{n}}就是合法的n次单位根，也就是复平面上的单位圆

DFT与IDFT

用来解决多项式乘法：\displaystyle H(x)[p]=\sum_{i=0}^{n}F(x)[i]G(x)[p-i]，直接做是O(n^2)

显然有H(m)=F(m)G(m)，所以我们可以考虑：

（1）先把F,G由系数表达式转化成若干个点值\to {\rm DFT}

（2）相乘得到G的若干个点值

（3）然后在把点值转换成系数\to {\rm IDFT}

• DFT

考虑一个n次多项式，我们需要n+1个点值才能表示出它的系数，不妨直接设n为多项式的项数

也就是按照次数奇偶分组，为了让项数相等方便处理，可以让n为偶数（如果是奇数就增大1，F[n+1]为0没有贡献）

设F1[k]=F[2k],F2[k]=F[2k+1]，那么{\rm DFT}(n,F)[p]={\rm DFT}(n,F_1)[2p]+x{\rm DFT}(n,F_2)[2p]

因为w_{n}^{2p}=w_{n/2}^{p}，所以\displaystyle {\rm DFT}(n,F)[p]={\rm DFT}(n/2,F_1)[p]+w_n^p{\rm DFT}(n/2,F_2)[p]

此时可以发现问题规模已经减半了，所以可以直接分治计算。

注意单位根实际上只需要0\sim n-1就够了，因为当p\geq n/2时，w_{n}^{2p}=-w_{n}^{2p-n/2}，所以在计算前面一半的时候可以顺便计算后面一半（系数取-即可）

注意到每一层都会有项数不相等的问题，所以一开始可以把n调整到2的若干次幂，这样就保证了每层奇偶都是一样的

• IDFT

设F做完{\rm DFT}后得到G，现在需要对G做{\rm IDFT}得到F

设\rm IDFT的变换系数为c(p,k)，即：\displaystyle {\rm IDFT}(G,n)[p]=\sum_{k=0}^{n-1}c(p,k)G[k]

即\displaystyle [j=p]=\sum_{k=0}^{n-1}c(p,k)w_n^{jk}

由单位根反演可知：\displaystyle [j|n]=\frac{1}{n}\sum_{k=0}^{n-1}w_n^{jk }

因为j,p<n，所以[j=p]=[j\equiv p\pmod n]=[(j-p)|n]

所以\displaystyle [(j-p)|n]=\frac{1}{n}\sum_{k=0}^{n-1}w_n^{(j-p)k}=\sum_{k=0}^{n-1}w_n^{jk}\times \frac{w_n^{-pk}}{n}

观察上面两个式子，可以得到：c(p,k)=\dfrac{w_n^{-pk}}{n}

整理以下即：\displaystyle {\rm IDFT}(F)[p]=\frac{1}{n}\sum_{k=0}^{n-1}(w_n^{-p})^kF[k]

所以计算\rm IDFT的过程只需要在\rm DFT的基础上，把求出来的单位根取逆就行了，最后再除掉n

优化（在原来的博客里有）：

1、预处理单位根

2、递归变迭代，蝴蝶变换

3、三次变两次（谨慎使用）

• 蝴蝶变换

以长度为8的序列为例子，初始为(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)，那么每经过一次奇偶分组就会得到：

将下标的二进制形式写出来，那么就能够发现初始在p位置的，最后就到了p二进制翻转后的位置

这样就可以直接从底层迭代上来了，不需要递归

NTT

建议先阅读同余

这里所有的原根都是在质数p=k\times 2^t+1的同余系下定义的（其中t足够大）

定义同余系下的单位根：

保证t足够大，是因为要保证能够整除n

现在要证明，这里的单位根和上面FFT的虚数单位根有同样的性质

1、\omega_n^k=\omega_n^{k\bmod n}

证明：\omega_{n}^{k}\equiv g^{k(p-1)/n}\equiv g^{k(p-1)/n\bmod (p-1)}

g^{k(p-1)/n\bmod (p-1)}\equiv g^{(p-1)((k/n)\bmod 1)}\equiv g^{(p-1)(k\bmod n)/n}\equiv (\omega_n^1)^{k\bmod n}=\omega_n^{k\bmod n}

2、\omega_{2n}^{2k}=\omega_n^{k}

证明：\omega_{2n}^{2k}\equiv g^{2k(p-1)/2n}\equiv g^{k(p-1)/n}\equiv w_n^k

3、\omega_{n}^{k+n/2}=-w_n^k

\omega_{n}^{k+n/2}\equiv \omega_n^k\times \omega_n^{n/2} \because (\omega_{n}^{n/2})^2\equiv \omega_n^n\equiv 1 \therefore w_n^{n/2}=1/-1

又w_n^k互不相同，所以\omega_{n}^{k+n/2}=-w_n^k

4、\sum_{i=1}^{n-1}\omega_n^i\equiv 0

注意之前FFT就说过要把多项式弄成2的整数次幂，这就是为什么上面取的质数p的t要足够大

放个\rm NTT的板子

inline void NTT(vint &A,int lim,int tag)
{
    for(int i=0;i<lim;i++) if(i<cir[i]) swap(A[i],A[cir[i]]);
    int w,buf,inv0=power(lim,mod-2);
    for(int l=2;l<=lim;l<<=1)
    {
        int r=l>>1; w=1,buf=power(tag?g:invg,(mod-1)/l);
        for(int i=0;i<lim;i+=l,w=1)
        for(int j=i;j<i+r;j++,w=1ll*w*buf%mod)
        {
            int tmp=1ll*w*A[j+r]%mod;
            A[j+r]=(A[j]-tmp+mod)%mod;
            A[j]=(A[j]+tmp)%mod;
        }
    }
    if(tag^1) for(int i=0;i<lim;i++) A[i]=1ll*A[i]*inv0%mod;
}