离散傅立叶变换与快速傅立叶变换（DFT与FFT）

#include <iostream>
#include <cmath>
#include <string.h>
using namespace std;
#define N 200005
#define PI acos(-1.0)
struct complex  
{  
    double r,i;  
    complex(double real=0.0,double image=0.0)  
    {  
        r=real;  
        i=image;  
    }  
    //以下为三种虚数运算的定义   
    complex operator+(const complex o)  
    {  
        return complex(r+o.r,i+o.i);  
    }  
    complex operator-(const complex o)  
    {  
        return complex(r-o.r,i-o.i);  
    }  
    complex operator*(const complex o)  
    {  
        return complex(r*o.r-i*o.i,r*o.i+i*o.r);  
    }  
}x1[N],x2[N];  
char a[N/2],b[N/2];  
int sum[N]; //结果存在sum里   


void bitrev(complex *y,int l) //二进制平摊反转置换 O(logn)   
{  
    register int i,j,k;  
    for(i=1,j=l/2;i<l-1;i++)  
    {  
        if(i<j)  swap(y[i],y[j]); //交换互为下标反转的元素    
                                 //i<j保证只交换一次   
        k=l/2;  
        while(j>=k) //由最高位检索，遇1变0，遇0变1，跳出   
        {  
            j-=k;  
            k/=2;  
        }  
        if(j<k)  j+=k;  
    }  
}  
void fft(complex *in,int n)
{
    int i,j,k;
    complex u,t;
    bitrev(in,n);
    for(int i=2;i<=n;i=i*2)
    {
        complex wn(cos((2*PI)/i),sin((2*PI)/i));//初始化单位复根
        for(j=0;j<n;j=j+i)
        {
            complex w(1,0);
            for(k=j;k<j+i/2;k++)
            {
                u=in[k];
                t=w*in[k+i/2];
                in[k]=u+t;
                in[k+i/2]=u-t;
                w=w*wn;
            }
        }
    }
}
void antifft(complex *in,int n)
{
    complex x,y;
    int i,j,k;
    bitrev(in,n);
    for(int i=2;i<=n;i=i*2)
    {
        complex init(cos((2*PI*-1)/i),sin((2*PI*-1)/i));
        for(j=0;j<n;j=j+i)
        {
            complex w(1,0);
            for(k=j;k<j+i/2;k++)
            {
                x=in[k];
                y=w*in[k+i/2];
                in[k]=x+y;in[k+i/2]=x-y;
                w=w*init;
            }
        }
    }
    for(int i=0;i<n;i++)
        in[i]=in[i].r/n;
}


int main() {

    while(~scanf("%s%s",a,b))
    {
        //提取系数并逆序
        int na=strlen(a);
        int nb=strlen(b);
    
        int n=1;
        //expand
        
        while(n<na*2 || n<nb*2)   n=n*2; //将次数界变成2^n
        for(int i=na-1;i>=0;i--)
        {
            x1[na-1-i].r=a[i]-'0';x1[na-1-i].i=0;
        }
        
        for(int i=nb-1;i>=0;i--)
        {
            x2[nb-1-i].r=b[i]-'0';x2[nb-1-i].i=0;
        }
        
        for(int i=na;i<n;i++)
        {
            x1[i].i=0;x1[i].r=0;
        }
        for(int i=nb;i<n;i++)
        {
            x2[i].i=0;x2[i].r=0;
        }
        
        //分别傅里叶变换
        fft(x1,n);
        fft(x2,n);
        //点值相乘得到新的傅里叶变换数值
        for(int i=0;i<n;i++)
            x1[i]=x1[i]*x2[i];
        //逆傅里叶变换得到系数，*10求和
        antifft(x1,n);
        memset(sum,0,sizeof(sum));
        //四舍五入
        for(int i=0;i<n;i++)
            sum[i]=x1[i].r+0.5;
        for(int i=0;i<n;i++) //进位    
        {  
            sum[i+1]+=sum[i]/10;  
            sum[i]%=10;
        }  
        while(sum[n]<=0 && n>0)   n--; //检索最高位   
        
        for(int i=n;i>=0;i--)  
            putchar(sum[i]+'0'); //倒序输出    
        printf("
");  
    }
    return 0;

}