洛谷p3803 FFT入门
ps:花了我一天的时间弄懂fft的原理,感觉fft的折半很神奇!
大致谈一谈FFT的基本原理:
对于两个多项式的卷积,可以O(n^2)求出来(妥妥的暴力)
显然一个多项式可以用a0+a1X+a2X2+a3X3+a4X^4……表示。
也可以用(x1,y1),(x2,y2),(x3,y3),(x4,y4)的点集来表示。
用点值表示有一个好处:两个多项式的卷积可以直接取相同的x值,y值相乘得到。
那么,怎么转化为点值表示呢?
直接代进去?显然也是O(n^2),没用……
设\(A(X)=a0+a1X+a2X^2+a3X^3+a4X^4……\)
\(A[0](X)=a0+a2X+a4X^2……\)
\(A[1](X)=a1+a3X+a5X^2……\)
那么A(X)=A0(X2)+X*A1(X2);
然而,X^2仍会有N种不同的取值。
引入一个复数根可以很好的解决这个问题:设i=\(\sqrt{-1}\)
\(cosθ+i∗sinθ\)有很多很好的性质:
例如\((cosθ+i∗sinθ)^k=coskθ+i*sinkθ\)
记\(Wn=cos\frac{2π}{n}+i∗sin\frac{2π}{n}\)
\(A(W{n}^{k})=A[0](((W{n}^{k})^2)+W{n}^{k}*A[1]((W{n}^{k})^2)\)
因为$$(W{n}{k})2=W{n}{2k}=cos\frac{2π*2k}{n}+i∗sin\frac{2π*2k}{n}=cos\frac{2πk}{\frac{n}{2}}+i∗sin\frac{2πk}{\frac{n}{2}}=(W{\frac{n}{2}{k}})$$
所以$$A(W{n}^{k})=A0+W{n}{k}*A[1]((W{n}^{k})^2)=A[0]W{\frac{n}{2}{k}}+W{n}{k}*A[1]W{\frac{n}{2}}{k}$$
根据三角函数的周期性:2π为一个周期。
\(W{n}^{k}=W{n}^{k mod n}\)
那么\(A(W{n}^{k})=A[0]W{\frac{n}{2}}^{k}+W{n}^{k}*A[1]W{\frac{n}{2}^{k}}\)可以拆成两部分,设\(t<\frac{n}{2}\)
\(A(W{n}^{t})=A[0]W{\frac{n}{2}}^{t}+W{n}^{t}*A[1]W{\frac{n}{2}^{t}}\)
\(A(W{n}^{t+\frac{n}{2}})=A[0]W{\frac{n}{2}}^{t}-W{n}^{t}*A[1]W{\frac{n}{2}^{t}}\)
因为\(W{n}^{t+\frac{n}{2}}=W{n}^{t}*W{n}^{\frac{n}{2}}\)
\(W{n}^{\frac{n}{2}}=W{2}=cos\frac{2π}{2}+i*sin\frac{2π}{2}=cosπ+i*sinπ=-1\)
证毕,啦啦啦!
一直折半下去算就可以将时间优化到O(nlogn)
上图摘自算法导论。算法导论下载,提取码qcfs
一题洛谷的裸题p3803
【代码】
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
const int N=1e6+10;
const double pi=acos(-1.0);
int n,m,nn,k;
struct Complex {
double real,image; //real+i*image
Complex(){}
Complex(double _real,double _image)
{
real=_real; image=_image;
}
friend Complex operator + (Complex A,Complex B) { return (Complex(A.real+B.real,A.image+B.image)); }
friend Complex operator - (Complex A,Complex B) { return (Complex(A.real-B.real,A.image-B.image)); }
friend Complex operator * (Complex A,Complex B) { return (Complex(A.real*B.real-A.image*B.image,A.real*B.image+A.image*B.real)); }
} A[N<<2],B[N<<2],C[N<<2];
int rev[N<<2];
void FFT(Complex *A,int n,int DFT)
{
for(int i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]);
for(int s=1;(1<<s)<=n;s++)
{
int mi=(1<<s);
Complex wn=Complex(cos(2*pi/mi),sin(2*pi/mi)*DFT);
for(int t=0;t<n;t+=mi)
{
Complex w=Complex(1,0);
for(int j=0;j<(mi>>1);j++)
{
Complex u=A[t+j],v=w*A[t+j+(mi>>1)];
A[t+j]=u+v; A[t+j+(mi>>1)]=u-v;
w=w*wn;
}
}
}
if(DFT==-1) for(int i=0;i<n;i++) A[i].real/=n,A[i].image/=n;
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++) { int xx; scanf("%d",&xx); A[i]=Complex(xx,0); }
for(int i=0;i<=m;i++) { int xx; scanf("%d",&xx); B[i]=Complex(xx,0); }
m+=n; k=0;
for(nn=1;nn<=m;nn<<=1) k++;
for(int i=0;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
FFT(A,nn,1); FFT(B,nn,1);
for(int i=0;i<=nn;i++) C[i]=A[i]*B[i];
FFT(C,nn,-1);
for(int i=0;i<=m;i++) printf("%d ",(int)(C[i].real+0.5));
printf("\n");
return 0;
}
ps:这篇总结写得真心累……