【快速傅立叶变换fft&数论变换ntt学习小记】

概述

fft（快速傅立叶变换）是用来解决多项式乘法的nlog（n）算法，它的主要思想是先把多项式的多项式表达法转化成若干个二维点对（x，y）（点值），把相同x的y乘起来（计算），最后利用这些点对计算出多项式的多项式表达法的系数（插值）。这中间使用了n次单位复数根的一些特殊性质，采用分治的思想快速地完成点值和插值。

n次单位复数根

要想理解fft，首先要知道n次单位复数根有哪些神奇的性质。

定义

n次单位复数根是满足$w^n=1$的复数w。

$w_n^k=e^{2πik/n}$，其中$w_n^k$是n个n次单位复数根中的第k个，i是虚数单位，$i=\sqrt{-1}$

$e^{iu}=cos(u)+isin(u)$，其中三角函数是弧度制。

*以下内容是我自己脑补的

我们可以感性地认为w由两部分构成，复数的实数部和虚数部，分别设为x，y。w就是复（数）平面上的点（x，y），u代表的就是和x正轴夹角的弧度，通过sin，cos就可以求出坐标，表示出复数w。

性质

性质1

$w_n^1$是w的主n此单位根，由$w_n^k=e^{2πik/n}$可以的得出$w_n^i=w_n^{i-1}*w_n^1$，整个定义都是指数形式的，是等比数列。

性质0

还记得刚刚那个图吗，由图可以感性得出$w_n^n=w_n^0=1$，为什么等于1，cos（0）=1。

性质2（群的性质）

$w_n^k=w_n^{k\%n}*w_n^{n*\lfloor k/n\rfloor}$，其中%是取模的意思。

$=w_n^{k\%n}*(w_n^n)^{\lfloor k/n\rfloor}$

$=w_n^{k\%n}*1^{\lfloor k/n\rfloor}$（性质0）

$=w_n^{k\%n}$

性质3（消去引理）

$w_{dn}^{dk}=w_n^k$，由$w_n^k=e^{2πik/n}$可以得出，不就是上下约分吗。

性质4（折半引理）

$(w_n^{k+n/2})^2=(w_n^k)^2*(w_n^{n/2})^2$

$=(w_n^k)^2*(w_n^n)^2$

$=(w_n^k)^2*1^2$（性质0）

$=(w_n^k)^2$

还记得刚刚那个图吗，这两个复数对应的点关于原点对称，其实相当于相反数关系。

性质5（求和引理）

$\sum_{j=0}^{n-1}(w_n^k)^j=((w_n^k)^n-1)/(w_n^k-1)$，等比数列求和

$=((w_n^n)^k-1)/(w_n^k-1)$

$=(1^k-1)/(w_n^k-1)$（性质0）

$=0$，当k!=0和k!=n时

点值

我们要求$y(w_n)=a_0*(w_n)^0+a_1*(w_n)^1……a_n*(w_n)^n$

令$y_0(w_n)=a_0*(w_n)^0+a_2*(w_n)^2……$

$y_1(w_n)=a_1*(w_n)^1+a_3*(w_n)^3……$

$y(w_n)=y_0((w_n)^2)+w*y_1((w_n)^2)$

$=y_0(w_{n/2})+w*y_1(w_{n/2})$(性质3)

也就是说我们要求$w_n$对应的n个点值，可以把系数分组后由$w_{n/2}$对应的n/2个点值求出。

首先我们考虑怎么分组，我们把模2为0的放左边，把模2为1的放右边，相当于0，1分组，按低位分组，倒着按高位编号，可以用下面的图理解一下。

分解到剩一个时y值就是a。考虑合并答案，从小到大枚举块的大小，由性质4可以把当前块分成左右两部分，只用考虑一部分，两部分对应的两个u’,v’转移是相似的，因为他们的值是相反数，平方之后就没有影响，只是外面的系数w成相反数，所以可以同时处理。$y(u')=y(u)+w*y(v),y(v')=y(u)-w*y(v)$

所以整个点值的流程是，把每个位置的数按标号二进制反过来重新放置，从小到大枚举块的大小，每个块找出对应的两个位置，位置相差块的一半，然后转移一下就可以了。

插值

我们可以由定义得出将a矩阵变成y矩阵的矩阵v的第j行第k列是$w_n^{kj}$，就是在多项式里自变量为$w_n^k$。这里我们有一个结论，将y矩阵变成x矩阵的矩阵$v^{-1}$的第j行第k列是$w_n^{-kj}/n$

首先$a*v*v^{-1}=a$，那么$v*v^{-1}$是一个值都为1的矩阵。

$[v*v^{-1}]_{j,j'}=\sum_{k=0}^{n-1} w_n^{jk}w_n^{-kj'}/n$

$=\sum_{k=0}^{n-1}w_n^{k(j-j')}/n$

$=\sum_{k=0}^{n-1}(w_n^{j-j'})^k/n$

$=0$(性质5)，当$j!=j'$时

当$j=j'$时，显然

$\sum_{k=0}^{n-1}(w_n^{j-j'})^k/n=\sum_{k=0}^{n-1}w_n^0/n$

$=\sum_{k=0}^{n-1}1/n$

$=n/n$

$=1$

有了这个结论后，我们只需要在插值的时候将$w_n^k$变成$w_n^{-k}$就可以直接套用点值时的做法，点值时的性质显然现在还符合。

这就结束了？

使用模运算的fft（原来这叫ntt）（有原根才能用）

由于我们使用了浮点数，当数据变大时我们的误差就变得不可接受，所以我们考虑化浮点数为整数。

我们用$g^k\%p$代替$w_n^1$，其中p是质数，p=kn+1。g是原根。

首先要$(g^k)^n\%p=(w_n^1)^n=1$

由于费马小定理$a^{p-1}=1(\%p)$

$(g^k)^n\%p=g^{kn}\%p=g^{p-1}\%p=1$

有了$w_n^1$我们就可以求出所有w（性质1）。

由于复数和普通的数本质上没有什么不同，所以之前的性质照样适用。

在这里提供一对比较好的g，p。g=3，p=1004535809。这里有一个问题，这样求出来的k可能不是整数（又产生了误差），可是我们发现n是2的指数次幂，而p-1的质因子有很多个2，多到n足够大都不会除不尽（在可接受复杂度内的n）。

贴上高精度乘法的代码

code（without mod）

#include<cmath>
#include<cstdio>
#include<algorithm>
#define LD double
#define LL long long
#define min(a,b) ((a<b)?a:b)
#define max(a,b) ((a>b)?a:b)
#define fo(i,j,k) for(int i=j;i<=k;i++)
#define fd(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
int const maxn=1e5;
int n,m,A,B;LL c[maxn*4+10];
LD pi=acos(-1);
struct rec{
    LD x,y;
    rec(LD X=0,LD Y=0){x=X;y=Y;}
};
rec operator+(rec x,rec y){return rec(x.x+y.x,x.y+y.y);}
rec operator-(rec x,rec y){return rec(x.x-y.x,x.y-y.y);}
rec operator*(rec x,rec y){return rec(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
rec a[maxn*4+10],b[maxn*4+10],t[maxn*4+10];
void read(int &num,rec *a){
    num=0;int v=0;char ch=getchar();
    for(;(ch<'0')||(ch>'9');ch=getchar());
    for(;(ch>='0')&&(ch<='9');a[num++]=ch-'0',ch=getchar());
    int mx=num/2-1;fo(i,0,mx){rec tmp=a[i];a[i]=a[num-i-1];a[num-i-1]=tmp;}
}
int up(LD x){return int(x)+((int(x)==x)?0:1);}
void DFT(rec *a,LD tag){
    fo(i,0,n-1){
        int pos=0;
        for(int j=0,ii=i;j<m;pos=(pos<<1)+(ii&1),ii=ii>>1,j++);
        t[pos]=a[i];
    }
    for(int i=2;i<=n;i=i<<1){
        int half=i>>1;
        fo(j,0,half-1){
            rec w(cos(tag*pi*j/half),sin(tag*pi*j/half));
            for(int k=j;k<n;k+=i){
                rec x=t[k],y=w*t[k+half];
                t[k]=x+y;
                t[k+half]=x-y;
            }
        }
    }
    fo(i,0,n-1)a[i]=t[i];
}
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    read(A,a);read(B,b);
    m=up(log(max(A,B)<<1)/log(2));n=1<<m;
    DFT(a,1);DFT(b,1);
    fo(i,0,n-1)a[i]=a[i]*b[i];
    DFT(a,-1);
    fo(i,0,n-1)c[i]=(a[i].x+1e-6)/n;
    fo(i,0,n-1){
        c[i+1]+=c[i]/10;
        c[i]%=10;
    }
    for(;c[n]==0;n--);
    fd(i,n,0)putchar(c[i]+'0');
    return 0;
}

code（with mod）

#include<cmath>
#include<cstdio>
#include<algorithm>
#define LD double
#define LL long long
#define min(a,b) ((a<b)?a:b)
#define max(a,b) ((a>b)?a:b)
#define fo(i,j,k) for(int i=j;i<=k;i++)
#define fd(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
int const maxn=1e5,g=3,mo=1004535809;
int n,m,A,B,w[maxn*4+10],a[maxn*4+10],b[maxn*4+10],t[maxn*4+10];
void read(int &num,int *a){
    num=0;int v=0;char ch=getchar();
    for(;(ch<'0')||(ch>'9');ch=getchar());
    for(;(ch>='0')&&(ch<='9');a[num++]=ch-'0',ch=getchar());
    int mx=num/2-1;fo(i,0,mx)swap(a[i],a[num-i-1]);
}
int up(LD x){return int(x)+((int(x)==x)?0:1);}
void DFT(int *a,int tag){
    fo(i,0,n-1){
        int pos=0;
        for(int j=0,ii=i;j<m;pos=(pos<<1)+(ii&1),ii=ii>>1,j++);
        t[pos]=a[i];
    }
    for(int i=2;i<=n;i=i<<1){
        int half=i>>1;
        fo(j,0,half-1){
            int wi=(tag>0)?w[n/i*j]:w[n-n/i*j];
            for(int k=j;k<n;k+=i){
                int x=t[k],y=1ll*wi*t[k+half]%mo;
                t[k]=(x+y)%mo;
                t[k+half]=(x-y+mo)%mo;
            }
        }
    }
    fo(i,0,n-1)a[i]=t[i];
}
int Pow(int x,int y){
    int z=1;
    while(y){
        if(y&1)z=1ll*z*x%mo;
        x=1ll*x*x%mo;
        y=y>>1;
    }
    return z;
}
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    read(A,a);read(B,b);
    m=up(log(max(A,B)<<1)/log(2));n=1<<m;
    w[0]=1;w[1]=Pow(g,(mo-1)/n);
    fo(i,2,n)w[i]=1ll*w[i-1]*w[1]%mo;
    DFT(a,1);DFT(b,1);
    fo(i,0,n-1)a[i]=1ll*a[i]*b[i]%mo;
    DFT(a,-1);
    int ni=Pow(n,mo-2);
    fo(i,0,n-1)a[i]=1ll*a[i]*ni%mo;
    fo(i,0,n-1){
        a[i+1]+=a[i]/10;
        a[i]%=10;
    }
    for(;a[n]==0;n--);
    fd(i,n,0)putchar(a[i]+'0');
    return 0;
}