FFT模板

// luogu-judger-enable-o2
/*(
 FFT 模板
 luogu P3803
in:
 n , m
 0,...,a[n]
 0,...,b[m]

example:
 1 2
 1 2
 1 2 1
 1 4 5 2
*/
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
void in(long long &x){
    x=0;char c=getchar();
    long long y=1;
    while(c<'0'||c>'9'){if(c=='-')y=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();}
    x*=y;
}
void o(long long x){
    if(x<0){x=-x;putchar('-');}
    if(x>9)o(x/10);
    putchar(x%10+'0');
}
#define int long long
const int _ = 1e7+10;
const double Pi = acos(-1.0);
struct complex{
    double x,y;
    complex(double xx=0.0,double yy = 0.0):x(xx),y(yy){};
};
complex operator + (complex a,complex b){
    return complex(a.x+b.x,a.y+b.y);
}
complex operator - (complex a,complex b){
    return complex(a.x-b.x,a.y-b.y);
}
complex operator * (complex a,complex b){
    return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
int n,m,N;
int limit = 1 , bit;
complex a[_],b[_];
int r[_];
void fft(complex *A,int type)
{
    for(int i=0;i<limit;i++){
        if(i<r[i])
            swap(A[i],A[r[i]]);
    }
    for(int mid = 1;mid < limit;mid<<=1){//当前区间的一半
        complex wn( cos(Pi/mid),type*sin(Pi/mid));
        for(int R = mid << 1,j = 0;j < limit;j+=R){//回溯区间大小 最大不超过 limit
            complex w(1,0);//幂?
            for(int k = 0 ;k<mid;k++,w=w*wn){
                complex x = A[j+k],y = w*A[j+mid+k];
                A[j+k]=x+y;
                A[j+k+mid]=x-y;
            }
        }

    }
}
signed main() {
    in(n);in(m);
    for(int i = 0;i <= n;i++)cin>>(a[i].x);
    for(int j = 0;j <= m;j++)cin>>(b[j].x);
    while (limit<=n+m)limit<<=1,bit++;
    for(int i = 0 ;i<limit;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
    }
    fft(a,1);
    fft(b,1);
    for(int i=0;i<=limit;i++)a[i]=a[i]*b[i];//y;
    fft(a,-1);
    for(int i=0;i<=n+m;i++){
        printf("%lld ",(long long)(a[i].x/limit + 0.5));
    }
    return 0;
}
原文地址:https://www.cnblogs.com/yesuweiYYYY/p/15037543.html