快速傅立叶变换FFT模板

Posted 2020-11-09 无尽的蓝黄

递归版

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);
}