FFT代码详解
关于FFT原理部分的介绍,在网上已经有很多了,所以在此只讲代码实现部分的内容。
原理可以参考https://www.cnblogs.com/RabbitHu/p/FFT.html
推荐看完它的原理解释再来看这里的代码解释
废话不多说,上代码(多项式乘法)
#include <iostream> #include <cstdio> #include <cmath> #define N 4000001 using namespace std; struct cp//手写复数类可以卡常 { double real,imag; }; cp operator +(cp a,cp b) { return (cp){a.real+b.real,a.imag+b.imag}; } cp operator -(cp a,cp b) { return (cp){a.real-b.real,a.imag-b.imag}; }
复数乘法:设$R_{a}$表示$a$的实部系数,$I_{a}$表示$a$的虚部系数
则$a*b$
$=(R_{a}+I_{a})*(R_{b}+I_{b})$
$=R_{a}*R_{b}+R_{a}*I_{b}+R_{b}*I_{a}+I_{a}*I_{b}$
因为$i^2=-1$
所以结果的实部为$R_{a}*R_{b}-I_{a}*I_{b}$
虚部为$R_{a}*I_{b}+R_{b}*I_{a}$
cp operator *(cp a,cp b) { return (cp){a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real}; } double pi=acos(-1.0); int lim,rev[N],len; cp w[N],inv[N],a[N],b[N]; void get_w() { for(int i=0;i<=lim;i++) { double angle=(double)i*2*pi/lim; w[i].imag=sin(angle); w[i].real=cos(angle); inv[i]=(cp){w[i].real,-w[i].imag}; } }
fft参数解释
$arr:$系数数组,在$fft$后变为点值数组,$arr_{i}$表示将$w^i_n$带入多项式后求得的值
$w:$预处理好的w单位根,在$fft$的时候正常带入即可,在$idft$的时候带入单位根的倒数(具体参见$idft$)void fft(cp *arr,cp *w)
{ for(int i=0;i<lim;i++) { //处理每一个系数在分治过程中的实际位置; //if是因为只需交换一次,所以选择由小的一方来执行 if(i<rev[i]) swap(arr[i],arr[rev[i]]); } for(int i=2;i<=lim;i*=2)//枚举区间长度 { int l=i/2; for(int j=0;j<lim;j+=i)//枚举区间位置,这些区间是互不相交的 { //枚举带入的单位根w(k,l),k>=l的单位根也可以在这里一并求出 for(int k=0;k<l;k++) {
意义变更
在这里$arr$的意义从系数变为$w^k_i$的点值,$a_{j,j+i-1}$分别表示将$w^{0,i-1}_i$的点值
下面的的t相当于文首博客中的$w^k_n * A_2(w^k_{n \over 2})$
cp t=arr[j+k+l]*w[lim/i*k];//w(k,i)=w(k/i,1)=w(n*k/i,n) arr[j+k+l]=arr[j+k]-t; arr[j+k]=arr[j+k]+t; } } } } int main() { int n,m; cin>>n>>m; for(int i=0;i<=n;i++) scanf("%lf",&a[i].real); for(int i=0;i<=m;i++) scanf("%lf",&b[i].real); lim=1; while(lim<=n+m) len++,lim<<=1;//这样会比用cmath的log要快? for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1)); get_w(); fft(a,w); fft(b,w); for(int i=0;i<=lim;i++) a[i]=a[i]*b[i]; fft(a,inv); for(int i=0;i<=n+m;i++) printf("%d ", (int)(a[i].real/lim+0.5)); //除以lim的原因具体参见idft,0.5是为了四舍五入 }
看都看了,顺手点个推荐呗 :)