快速傅里叶变换 FFT

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<vector>
#include<algorithm>
using namespace std;
#define MAXN 550000+100
#define PI 3.14159265358
#define EPS 1e-2
char a[MAXN],b[MAXN];
int n,m,top=0,label=0;
int ans[MAXN];
struct complex
{
    double R,I;
    complex(){}
    complex(double r0,double i0): R(r0),I(i0) {}
};
complex operator *(const complex& A,const complex &B)
{
    return complex(A.R*B.R-A.I*B.I,A.R*B.I+B.R*A.I);
}
complex operator +(const complex& A,const complex &B)
{
    return complex(A.R+B.R,A.I+B.I);
}
complex operator -(const complex& A,const complex &B)
{
    return complex(A.R-B.R,A.I-B.I);
}
vector<complex> C[3*MAXN];
complex p[MAXN];
int FFT(int first,int floor,int multi)
{
    int i;
    if(n==floor) 
    {
        C[++top].push_back(p[first]);
        return top;
    }
    complex wn=complex(cos(2*PI/n*floor),sin(2*PI/n*floor)*multi);
    complex w=complex(1,0);
    label++;
    int top0=FFT(first,floor*2,multi);
    int top1=FFT(first+floor,floor*2,multi);
    top++;
    for(i=0;i<n/floor/2;i++,w=w*wn)
        C[top].push_back(C[top0][i]+w*C[top1][i]);

    w=complex(1,0);
    for(i=0;i<n/floor/2;i++,w=w*wn)
        C[top].push_back(C[top0][i]-w*C[top1][i]);
    return top;
}
int make_int(double x)
{
    int A=x;
    if (fabs(x-A)<EPS)
        return A;
    return A+1;
}
void solve()
{
    int i=1;
    for(i=0;i<n;i++)
        p[i].R=a[n-i-1]-48;
    while(i<n) i*=2;
    n=i*2;
    
    i=1;
    while(i<m) i*=2;
    n=max(i*2,n);

    int top0=FFT(0,1,1);

    memset(p,0,sizeof(p));
    for(i=0;i<m;i++)
        p[i].R=b[m-i-1]-48;
    int top1=FFT(0,1,1);

    memset(p,0,sizeof(p));
    for(i=0;i<n;i++)
        p[i]=C[top0][i]*C[top1][i];
    int top2=FFT(0,1,-1);

    for(i=0;i<n;i++) 
    {
        ans[i]+=make_int(C[top2][i].R/n);
        ans[i+1]=ans[i]/10;
        ans[i]=ans[i]%10;
    }

    for(i=n-1;i>=0;i--) if(ans[i]!=0) break;
    for(i;i>=0;i--) printf("%d",ans[i]);
    printf("\n");
}

int main()
{
    int i;
    scanf("%s%s",a,b);
    n=strlen(a); m=strlen(b);
    solve();
}