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

概述

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

n次单位复数根

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

定义

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

wkn=e2πik/n，其中wkn是n个n次单位复数根中的第k个，i是虚数单位，i=−1−−−√

eiu=cos(u)+isin(u)，其中三角函数是弧度制。

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

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

性质

性质1

w1n是w的主n此单位根，由wkn=e2πik/n可以的得出win=wi−1n∗w1n，整个定义都是指数形式的，是等比数列。

性质0

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

性质2（群的性质）

wkn=wk%nn∗wn∗⌊k/n⌋n，其中%是取模的意思。

=wk%nn∗(wnn)⌊k/n⌋

=wk%nn∗1⌊k/n⌋（性质0）

=wk%nn

性质3（消去引理）

wdkdn=wkn，由wkn=e2πik/n可以得出，不就是上下约分吗。

性质4（折半引理）

(wk+n/2n)2=(wkn)2∗(wn/2n)2

=(wkn)2∗(wnn)2

=(wkn)2∗12（性质0）

=(wkn)2

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

性质5（求和引理）

∑n−1j=0(wkn)j=((wkn)n−1)/(wkn−1)，等比数列求和

=((wnn)k−1)/(wkn−1)

=(1k−1)/(wkn−1)（性质0）

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

点值

我们要求y(wn)=a0∗(wn)0+a1∗(wn)1……an∗(wn)n

令y0(wn)=a0∗(wn)0+a2∗(wn)2……

y1(wn)=a1∗(wn)1+a3∗(wn)3……

y(wn)=y0((wn)2)+w∗y1((wn)2)

=y0(wn/2)+w∗y1(wn/2)(性质3)

也就是说我们要求wn对应的n个点值，可以把系数分组后由wn/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列是wkjn，就是在多项式里自变量为wkn。这里我们有一个结论，将y矩阵变成x矩阵的矩阵v−1的第j行第k列是w−kjn/n

首先a∗v∗v−1=a，那么v∗v−1是一个值都为1的矩阵。

[v∗v−1]j,j′=∑n−1k=0wjknw−kj′n/n

=∑n−1k=0wk(j−j′)n/n

=∑n−1k=0(wj−j′n)k/n

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

当j=j′时，显然

∑n−1k=0(wj−j′n)k/n=∑n−1k=0w0n/n

=∑n−1k=01/n

=n/n

=1

有了这个结论后，我们只需要在插值的时候将wkn变成w−kn就可以直接套用点值时的做法，点值时的性质显然现在还符合。

这就结束了？

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

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

我们用gk%p代替w1n，其中p是质数，p=kn+1。g是原根。

首先要(gk)n%p=(w1n)n=1

由于费马小定理ap−1=1(%p)

(gk)n%p=gkn%p=gp−1%p=1

有了w1n我们就可以求出所有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;}