FFT模板(多项式乘法)
标签: FFT
扯淡
一晚上都用来捣鼓这个东西了......
这里贴一位神犇的博客,我认为讲的比较清楚了。(刚好适合我这种复数都没学的)
http://blog.csdn.net/leo_h1104/article/details/51615710
题解
不写点什么也不好,我就简单的说一下吧。
我们首先得知道DFT(离散傅里叶变换)和IDFT(逆离散傅里叶变换)。
一个多项式有很两种表示方法:
法一:(f(x)=sum_{i=0}^n A_i*x^i)
法二:图像上的任(n+1)个点,如(f(x)=x+1)就可以用(0,1),(1,2),(2,3)来表示。
法二其实是很适合两个函数的相乘,只需要对应横坐标点的纵坐标相乘。
DFT其实就是将表示法 法一转换成法二,IDFT则相反。
假如DFT使用的点坐标基于实数,那么复杂度为(O(n^2)),相比较与基于复数的FFT,效率十分底下。
FFT就使用了单位根的(0~n-1)次方作为点的横坐标(这里n需要补成2的次幂),再利用单位根的某些性质,把规模减小一半。可以实现(O(nlogn))的计算。
Code
#include<complex>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
#define ll long long
#define REP(i,a,b) for(int i=(a),_end_=(b);i<=_end_;i++)
#define DREP(i,a,b) for(int i=(a),_end_=(b);i>=_end_;i--)
#define EREP(i,a) for(int i=start[(a)];i;i=e[i].next)
inline int read()
{
int sum=0,p=1;char ch=getchar();
while(!(('0'<=ch && ch<='9') || ch=='-'))ch=getchar();
if(ch=='-')p=-1,ch=getchar();
while('0'<=ch && ch<='9')sum=sum*10+ch-48,ch=getchar();
return sum*p;
}
const int maxn=3e6+20;
int n,m,l,rev[maxn];
complex <double> a[maxn],b[maxn];
void init()
{
n=read();m=read();
REP(i,0,n)a[i]=read();
REP(i,0,m)b[i]=read();
m+=n;
for(n=1;n<=m;n<<=1)l++;
REP(i,0,n-1)rev[i]=(rev[i>>1]>>1) | ((i&1)<<(l-1));
}
const double Pi=acos(-1);
void FFT(complex <double> *p,int opt)
{
REP(i,0,n-1)if(i<rev[i])swap(p[i],p[rev[i]]);
for(int i=1;i<n;i<<=1)
{
complex <double> W(cos(Pi/i),opt*sin(Pi/i));
for(int P=i<<1,j=0;j<n;j+=P)
{
complex <double> w(1,0);
for(int k=j;k<i+j;k++,w*=W)
{
complex <double> x=p[k],y=w*p[k+i];
p[k]=x+y;
p[k+i]=x-y;
}
}
}
if(opt==-1)REP(i,0,n)p[i]/=n;
}
void doing()
{
FFT(a,1);
FFT(b,1);
REP(i,0,n)a[i]=a[i]*b[i];
FFT(a,-1);
REP(i,0,m)printf("%d ",(int)(a[i].real()+0.5));
}
int main()
{
freopen("FFT.in","r",stdin);
freopen("FFT.out","w",stdout);
init();
doing();
return 0;
}