快速傅立葉變換FFT模板
阿新 • • 發佈:2018-05-28
return namespace double names http ++ main swap pre
遞歸版
UOJ34多項式乘法
//容易暴棧,但是很好理解 #include <cmath> #include <iostream> #include <cstdio> #include <cstdlib> #include <cstring> #include <algorithm> #include <queue> #include <map> const int maxlongint=2147483647; const int mo=1e9+7; const int N=400005; const double pi=acos(-1); using namespace std; struct arr { double x,y; arr() {x=y=0;} arr(double x,double y):x(x),y(y) {} }a[N],b[N],c[N]; arr operator +(arr x,arr y) {return arr(x.x+y.x,x.y+y.y);} arr operator -(arr x,arr y) {return arr(x.x-y.x,x.y-y.y);} arr operator *(arr x,arr y) {return arr(x.x*y.x-x.y*y.y,x.x*y.y+y.x*x.y);} int n,m,fn; void FFT(arr *y,int n,int t) { if(n==1) return; arr a0[n>>1],a1[n>>1]; for(int i=0;i<n;i+=2) a0[i>>1]=y[i],a1[i>>1]=y[i+1]; FFT(a0,n>>1,t),FFT(a1,n>>1,t); arr w1(cos(2*pi/n),t*sin(2*pi/n)),w0(1,0); for(int i=0;i<n>>1;i++,w0=w0*w1) y[i]=a0[i]+w0*a1[i],y[i+(n>>1)]=a0[i]-w0*a1[i]; } int main() scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&a[i].x); for(int i=0;i<=m;i++) scanf("%lf",&b[i].x); fn=1; while(fn<=n+m) fn<<=1; FFT(a,fn,1),FFT(b,fn,1); for(int i=0;i<fn;i++) c[i]=a[i]*b[i]; FFT(c,fn,-1); for(int i=0;i<=n+m;i++) printf("%.0lf ",abs(c[i].x/fn)); }
非遞歸版
BZOJ3527[Zjoi2014]力
#include <cmath> #include <iostream> #include <cstdio> #include <cstdlib> #include <cstring> #include <algorithm> #include <queue> #include <map> const int maxlongint=2147483647; const int mo=1e9+7; const int N=400005; const double pi=acos(-1); using namespace std; struct arr { double x,y; arr() {x=y=0;} arr(double x1,double y1) {x=x1,y=y1;}; }q[N],r[N],f[N],f1[N]; int n,fn; double qq[N]; arr operator + (arr x,arr y) {return arr(x.x+y.x,x.y+y.y);} arr operator - (arr x,arr y) {return arr(x.x-y.x,x.y-y.y);} arr operator * (arr x,arr y) {return arr(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);} void FFT(arr *a,int n,int t) { for(int i=0,p=0;i<n;i++) { if(i<p) swap(a[i],a[p]); for(int j=n>>1;(p^=j)<j;j>>=1); } for(int m=2;m<=n;m<<=1) { int half=m>>1; for(int i=0;i<half;i++) { arr w0(cos(i*pi*t/half),sin(i*pi*t/half)),aj; for(int j=i;j<n;j+=m) aj=a[j],a[j]=aj+w0*a[j+half],a[j+half]=aj-w0*a[j+half]; } } } int main() { scanf("%d",&n); for(int i=1;i<=n;i++) scanf("%lf",&qq[i]),r[i].x=1.0/i/i,q[i].x=qq[i]; for(fn=1;fn<n*2+1;fn<<=1); FFT(q,fn,1),FFT(r,fn,1); for(int i=0;i<fn;i++) f[i]=q[i]*r[i]; FFT(f,fn,-1); memset(q,0,sizeof(q)); memset(r,0,sizeof(r)); for(int i=1;i<=n;i++) r[i].x=1.0/i/i,q[i].x=qq[n-i+1]; FFT(q,fn,1),FFT(r,fn,1); for(int i=0;i<fn;i++) f1[i]=q[i]*r[i]; FFT(f1,fn,-1); for(int i=1;i<=n;i++) printf("%.3lf\n",(f[i].x-f1[n-i+1].x)/fn); }
快速傅立葉變換FFT模板