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

概述

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

n次单位复数根

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

定义

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

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

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

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

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

性质

性质1

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

性质0

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

性质2（群的性质）

wkn=wk%nn∗wn∗⌊k/n⌋n w_n^k=w_n^{k\%n}w_n^{n\lfloor k/n\rfloor}，其中%是取模的意思。

=wk%nn∗(wnn)⌊k/n⌋ =w_n^{k\%n}*(w_n^n)^{\lfloor k/n\rfloor}

=wk%nn∗1⌊k/n⌋ =w_n^{k\%n}*1^{\lfloor k/n\rfloor}（性质0）

=wk%nn =w_n^{k\%n}

性质3（消去引理）

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

性质4（折半引理）

(wk+n/2n)2=(wkn)2∗(wn/2n)2 (w_n^{k+n/2})^2=(w_n^k)^2*(w_n^{n/2})^2

=(wkn)2∗(wnn)2 =(w_n^k)^2*(w_n^n)^2

=(wkn)2∗12 =(w_n^k)^2*1^2（性质0）

=(wkn)2 =(w_n^k)^2

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

性质5（求和引理）

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

=((wnn)k−1)/(wkn−1) =((w_n^n)^k-1)/(w_n^k-1)

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

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

点值

我们要求 y(wn)=a0∗(wn)0+a1∗(wn)1……an∗(wn)n y(w_n)=a_0(w_n)^0+a_1(w_n)^1……a_n*(w_n)^n

令 y0(wn)=a0∗(wn)0+a2∗(wn)2…… y_0(w_n)=a_0(w_n)^0+a_2(w_n)^2……

y1(wn)=a1∗(wn)1+a3∗(wn)3…… y_1(w_n)=a_1(w_n)^1+a_3(w_n)^3……

y(wn)=y0((wn)2)+w∗y1((wn)2) y(w_n)=y_0((w_n)^2)+w*y_1((w_n)^2)

=y0(wn/2)+w∗y1(wn/2) =y_0(w_{n/2})+w*y_1(w_{n/2})(性质3)

也就是说我们要求 wn w_n对应的n个点值，可以把系数分组后由 wn/2 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) y(u')=y(u)+wy(v),y(v')=y(u)-wy(v)

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

插值

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

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

[v∗v−1]j,j′=∑n−1k=0wjknw−kj′n/n [v*v^{-1}]_{j,j'}=\sum_{k=0}^{n-1} w_n^{jk}w_n^{-kj'}/n

=∑n−1k=0wk(j−j′)n/n =\sum_{k=0}^{n-1}w_n^{k(j-j')}/n

=∑n−1k=0(wj−j′n)k/n =\sum_{k=0}^{n-1}(w_n^{j-j'})^k/n

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

当 j=j′ j=j'时，显然

∑n−1k=0(wj−j′n)k/n=∑n−1k=0w0n/n \sum_{k=0}^{n-1}(w_n^{j-j'})^k/n=\sum_{k=0}^{n-1}w_n^0/n

=∑n−1k=01/n =\sum_{k=0}^{n-1}1/n

=n/n =n/n

=1 =1

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

这就结束了？

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

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

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

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

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

(gk)n%p=gkn%p=gp−1%p=1 (g^k)^n\%p=g^{kn}\%p=g^{p-1}\%p=1

有了 w1n 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;
}

【快速傅立叶变换fft数论变换ntt学习小记】相关推荐

BZOJ_2179_FFT快速傅立叶_(FFT)
描述 http://www.lydsy.com/JudgeOnline/problem.php?id=2179 超大整数乘法分析 FFT模板题. 把数字看成是多项式,x是10.然后用FFT做多项式乘 ...
【学习笔记】超简单的快速数论变换（NTT）（FFT的优化）（含全套证明）
整理的算法模板合集: ACM模板目录一.前置知识二.快速数论变换(NTT) 三.NTT证明(和FFT的关系) 四.NTT模板数组形式的实现 vector形式的实现点我看多项式全家桶(●^◡_ ...
算法学习FFT系列(2):快速数论变换NTT bzoj3992: [SDOI2015]序列统计例题详解
bzoj3992: [SDOI2015]序列统计 Description 小C有一个集合S,里面的元素都是小于M的非负整数.他用程序编写了一个数列生成器,可以生成一个长度为N的数列,数列中的每个数都属 ...
比FFT还容易明白的NTT（快速数论变换）
NTT相关一种快速数论变换算法,这种算法是以数论为基础,对样本点为的数论变换,按时间抽取的方法,得到一组等价的迭代方程,有效高速简化了方程中的计算公式·与直接计算相比,大大减少了运算次数.(见快速傅 ...
FFT 快速傅里叶变换 NTT 快速数论变换
SDNU 1531 a*b III (FFT模板) Description 计算a乘b,多组输入(50组以内). Input 输入a b,数据范围0 <= a,b <= 10^100000 ...
解题报告（二）C、（darkBZOJ 2194）快速傅立叶之二（FFT、卷积的概念、常用变换）
繁凡出品的全新系列:解题报告系列 -- 超高质量算法题单,配套我写的超高质量题解和代码,题目难度不一定按照题号排序,我会在每道题后面加上题目难度指数(1∼51 \sim 51∼5),以模板题难度 11 ...
快速傅立叶变换(FFT)算法(原来这就是蝶形变换)
快速傅立叶变换(FFT)算法(原来这就是蝶形变换) 为了实现FFT的海面模拟,不得不先撸个FFT算法实现. 离散傅立叶变换(DFT) 学习FFT之前,首先要先了解什么是DFT,我们都知道傅立叶变换是将 ...
神经网络中快速傅立叶变换(FFT)的梯度传递
最近需要在神经网络中构造复数酉矩阵做权重系数,使用了快速傅立叶变换和反变换. 但是FFT不是theano的现成操作模块(有人写过对应的代码,所以应该会很快加进去了),所以想自己去写梯度传递来彻底搞清楚 ...
多项式算法2：NTT（快速数论变换）
多项式算法2:NTT(快速数论变换) 前言前置知识正文前言算法简介上次的FFT大量使用浮点数,有精度不好控制的问题,使用可以保证精度的方法又容易超时.这里NTT是在取模的一个前提下找到一个新 ...

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

概述

n次单位复数根

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

定义

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

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

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

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

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

性质

性质1

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

性质0

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

性质2（群的性质）

wkn=wk%nn∗wn∗⌊k/n⌋n w_n^k=w_n^{k\%n}*w_n^{n*\lfloor k/n\rfloor}，其中%是取模的意思。

=wk%nn∗(wnn)⌊k/n⌋ =w_n^{k\%n}*(w_n^n)^{\lfloor k/n\rfloor}

=wk%nn∗1⌊k/n⌋ =w_n^{k\%n}*1^{\lfloor k/n\rfloor}（性质0）

=wk%nn =w_n^{k\%n}

性质3（消去引理）

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

性质4（折半引理）

(wk+n/2n)2=(wkn)2∗(wn/2n)2 (w_n^{k+n/2})^2=(w_n^k)^2*(w_n^{n/2})^2

=(wkn)2∗(wnn)2 =(w_n^k)^2*(w_n^n)^2

=(wkn)2∗12 =(w_n^k)^2*1^2（性质0）

=(wkn)2 =(w_n^k)^2

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

性质5（求和引理）

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

=((wnn)k−1)/(wkn−1) =((w_n^n)^k-1)/(w_n^k-1)

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

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

点值

我们要求 y(wn)=a0∗(wn)0+a1∗(wn)1……an∗(wn)n y(w_n)=a_0*(w_n)^0+a_1*(w_n)^1……a_n*(w_n)^n

令 y0(wn)=a0∗(wn)0+a2∗(wn)2…… y_0(w_n)=a_0*(w_n)^0+a_2*(w_n)^2……

y1(wn)=a1∗(wn)1+a3∗(wn)3…… y_1(w_n)=a_1*(w_n)^1+a_3*(w_n)^3……

y(wn)=y0((wn)2)+w∗y1((wn)2) y(w_n)=y_0((w_n)^2)+w*y_1((w_n)^2)

=y0(wn/2)+w∗y1(wn/2) =y_0(w_{n/2})+w*y_1(w_{n/2})(性质3)

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

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

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

插值

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

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

[v∗v−1]j,j′=∑n−1k=0wjknw−kj′n/n [v*v^{-1}]_{j,j'}=\sum_{k=0}^{n-1} w_n^{jk}w_n^{-kj'}/n

=∑n−1k=0wk(j−j′)n/n =\sum_{k=0}^{n-1}w_n^{k(j-j')}/n

=∑n−1k=0(wj−j′n)k/n =\sum_{k=0}^{n-1}(w_n^{j-j'})^k/n

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

当 j=j′ j=j'时，显然

∑n−1k=0(wj−j′n)k/n=∑n−1k=0w0n/n \sum_{k=0}^{n-1}(w_n^{j-j'})^k/n=\sum_{k=0}^{n-1}w_n^0/n

=∑n−1k=01/n =\sum_{k=0}^{n-1}1/n

=n/n =n/n

=1 =1

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

这就结束了？

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

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

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

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

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

(gk)n%p=gkn%p=gp−1%p=1 (g^k)^n\%p=g^{kn}\%p=g^{p-1}\%p=1

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

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

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

贴上高精度乘法的代码

code（without mod）

code（with mod）

【快速傅立叶变换fft数论变换ntt学习小记】相关推荐

最新文章

热门文章

wkn=wk%nn∗wn∗⌊k/n⌋n w_n^k=w_n^{k\%n}w_n^{n\lfloor k/n\rfloor}，其中%是取模的意思。

我们要求 y(wn)=a0∗(wn)0+a1∗(wn)1……an∗(wn)n y(w_n)=a_0(w_n)^0+a_1(w_n)^1……a_n*(w_n)^n

令 y0(wn)=a0∗(wn)0+a2∗(wn)2…… y_0(w_n)=a_0(w_n)^0+a_2(w_n)^2……

y1(wn)=a1∗(wn)1+a3∗(wn)3…… y_1(w_n)=a_1(w_n)^1+a_3(w_n)^3……

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