快速傅里叶变换FFT

是lzh学长讲过以后，又看了小迪的博客，才学会的fft

小迪这个博客太推荐了，一学就会https://www.cnblogs.com/RabbitHu/p/FFT.html

模板

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<complex>
#define maxn 4000010
#define PI (acos(-1.0))
using namespace std;
complex<double>a[maxn],b[maxn];
int id[maxn];
void fft(complex<double> *p,int N,int f){
    for(int i=0;i<=N;i++)
        if(id[i]>i)swap(p[i],p[id[i]]);
    for(int k=1;k<N;k<<=1){
        complex<double>wn(cos(PI/k),f*sin(PI/k));
        for(int j=0;j<N;j+=k<<1){
            complex<double>w(1,0);
            for(int i=0;i<k;i++,w=w*wn){
                complex<double>x=p[i+j];
                complex<double>y=p[i+j+k]*w;
                p[i+j]=x+y;
                p[i+j+k]=x-y;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<N;i++)p[i]=p[i]/(double)N;
}
int main(){
    int N,M;
    scanf("%d%d",&N,&M);
    double x;
    for(int i=0;i<=N;i++){
        scanf("%lf",&x);
        a[i].real(x);
    }
    for(int i=0;i<=M;i++){
        scanf("%lf",&x);
        b[i].real(x);
    }
    M=N+M;N=1;
    int l=0;
    while(N<=M)N<<=1,l++;
    for(int i=0;i<=N;i++)id[i]=(id[i>>1]>>1)|((i&1)<<(l-1));
    fft(a,N,1);fft(b,N,1);
    for(int i=0;i<=N;i++)a[i]=a[i]*b[i];
    fft(a,N,-1);
    for(int i=0;i<=M;i++)
        printf("%d ",(int)(a[i].real()+0.5));
    return 0;
}

 

A - A * B Problem Plus

题意：求两个长度为50000位的大整数的乘积

解：fft之后进位，并特判为0的情况

注意fft计算时一定要用double类型，结果四舍五入取整数

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<complex>
#define maxn 400010
#define PI (acos(-1.0))
using namespace std;
complex<double>a[maxn],b[maxn];
char s1[maxn],s2[maxn];
int id[maxn],ans[maxn];
void fft(complex<double> *p,int N,int f){
    for(int i=0;i<=N;i++)
        if(id[i]>i)swap(p[i],p[id[i]]);
    for(int k=1;k<N;k<<=1){
        complex<double>wn(cos(PI/k),f*sin(PI/k));//点值化为系数表示，用共轭复数 
        for(int j=0;j<N;j+=k<<1){
            complex<double>w(1,0);
            for(int i=0;i<k;i++,w=w*wn){
                complex<double>x=p[i+j];
                complex<double>y=p[i+j+k]*w;
                p[i+j]=x+y;
                p[i+j+k]=x-y;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<N;i++)p[i]=p[i]/(double)N;
}

int main(){
    int N,M;
//    scanf("%d%d",&N,&M);
    int x;
    while(scanf("%s",s1)!=EOF){
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        scanf("%s",s2);
        N=strlen(s1);
        M=strlen(s2); 
        for(int i=0;i<N;i++){
            x=s1[i]-'0';
            a[i].real(x);
        }
        for(int i=0;i<M;i++){
            x=s2[i]-'0';
            b[i].real(x);
        }
        M=N+M;N=1;
        int l=0;
        while(N<=M)N<<=1,l++;
        for(int i=0;i<=N;i++)id[i]=(id[i>>1]>>1)|((i&1)<<(l-1));
        fft(a,N,1);fft(b,N,1);
        for(int i=0;i<=N;i++)a[i]=a[i]*b[i];
        fft(a,N,-1);
        for(int i=0;i<M-1;i++)ans[i]=(int)(a[i].real()+0.5);
        for(int i=M-2;i>0;i--){
            if(ans[i]>=10){
                ans[i-1]+=ans[i]/10;
                ans[i]%=10;
            }
        }
        bool flag=0;
        for(int i=0;i<M-1;i++){
            if(ans[i]!=0)flag=1;
            if(flag)printf("%d",ans[i]);
        }
        if(!flag)printf("0");
        puts("");
    } 
    
    return 0;
}

 

力

OJ数据出问题了，我也不知道代码能不能过

具体思路见这个博客：https://blog.csdn.net/qq_33929112/article/details/54590319

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<complex>
#define PI (acos(-1.0))
#define maxn 400010
using namespace std;
double q[maxn];
complex<double>f[maxn],g[maxn],h[maxn];
int id[maxn],n;

void fft(complex<double> *p,int N,int f){
    for(int i=0;i<=N;i++)
        if(id[i]>i)swap(p[i],p[id[i]]);
    for(int k=1;k<N;k<<=1){
        complex<double>wn(cos(PI/k),f*sin(PI/k));
        for(int j=0;j<N;j+=k<<1){
            complex<double>w(1,0);
            for(int i=0;i<k;i++,w=w*wn){
                complex<double>x=p[i+j];
                complex<double>y=p[i+j+k]*w;
                p[i+j]=x+y;
                p[i+j+k]=x-y;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<N;i++)p[i]=p[i]/(double)N;
}

int main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++){
        scanf("%lf",&q[i]);
        f[i]=q[i];
        h[n-i-1]=q[i];
    }
    for(int i=1;i<n;i++){
        g[i]=1.0/i/i;
    }
    int l=0,N=1,M=(n-1)*2;
    while(N<=M)N<<=1,l++;
    for(int i=0;i<N;i++)id[i]=(id[i>>1]>>1)|((i&1)<<(l-1));
    fft(f,N,1);fft(g,N,1);fft(h,N,1);
    for(int i=0;i<N;i++)f[i]=f[i]*g[i];
    for(int i=0;i<N;i++)h[i]=h[i]*g[i];
    fft(f,N,-1);fft(h,N,-1);
    for(int i=0;i<n;i++){
        double ans=f[i].real()-h[n-i-1].real();
        printf("%.10lf\n",ans);
    }
    return 0;
}