FFT快速傅里叶变换

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
typedef complex<double> cp;                     //complex库
const double PI = acos(-1);

void fft(cp *a,int n,int flag)                  //作用：求出a的点值表达，存进a
{
    int i;
    cp a0[n/2+1],a1[n/2+1];

    if(n==1) return;
    cp w_n(cos(2*PI/n),sin(flag*2*PI/n));   //flag=1:求值  flag=2:插值
    cp w(1,0);

    for(i=0; i<n/2; i++) a0[i]=a[i*2],a1[i]=a[i*2+1];   //分治

    fft(a0,n/2,flag);
    fft(a1,n/2,flag);

    for(i=0; i<n/2; i++)
    {
        a[i]=a0[i]+w*a1[i];
        a[i+n/2]=a0[i]-w*a1[i];
        w=w*w_n;                                //递推单位复数根
    }

}

cp x[300005]= {0},y[300005]= {0};

char arr1[100005], arr2[100005];
int ans[200005];

int main()
{
    while(~scanf("%s%s", arr1, arr2))
    {
        for(int i=0; i<300005; i++) x[i]=cp(0, 0);
        for(int i=0; i<300005; i++) y[i]=cp(0, 0);

        int n, m, i;
        n=strlen(arr1)-1;
        m=strlen(arr2)-1;

        for(i=0; i<=n; i++) x[i]=cp(arr1[n-i]-'0', 0);
        for(i=0; i<=m; i++) y[i]=cp(arr2[m-i]-'0', 0);
        m+=n;
        for(n=1; n<=m; n=n*2);

        fft(x,n,1);                                 //求值
        fft(y,n,1);                                 //求值

        for(i=0; i<n; i++) x[i]=x[i]*y[i];          //点值乘法
        fft(x,n,-1);                                //插值

        memset(ans, 0, sizeof(ans));

        for(i=0; i<n || ans[i]; i++)
        {
            ans[i]+=(x[i].real()/n+0.5);
            ans[i+1]+=(ans[i]/10);
            ans[i]%=10;
        }

        while(i>0 && !ans[i])i--;
        for(; i>=0; i--) printf("%d", ans[i]);
        printf("\n");
    }
    return 0;
}

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
const double PI = acos(-1);

struct Complex
{
    double real, image;
    Complex(double real=0.0, double image=0.0): real(real), image(image) {}
    Complex operator +(const Complex &b)
    {
        return Complex(real+b.real,image+b.image);
    }
    Complex operator -(const Complex &b)
    {
        return Complex(real-b.real,image-b.image);
    }
    Complex operator *(const Complex &b)
    {
        return Complex(real*b.real-image*b.image,real*b.image+image*b.real);
    }
};

Complex x[300005]= {0},y[300005]= {0};
char arr1[100005], arr2[100005];
int ans[200005];
int n, m;

void reverse(Complex a[],int n) //对长度为的复数序列a进行DFT和IDFT之前必要的反转变换操作(比特反转),n一定是2的幂次
{
    for(int i=1,j=n/2,k; i<n-1; i++)
    {
        if(i<j) swap(a[i],a[j]);
        k=n/2;
        while(j>=k)
        {
            j-=k;
            k/=2;
        }
        if(j<k) j+=k;
    }
}

void fft(Complex *a,int n,int flag) //flag=1是求值，flag=-1是插值
{
    reverse(a,n);
    for(int i=1; i<n; i<<=1) //i=该层的相邻结点之间编号之差，logn次迭代模拟，自底向上
    {
        Complex wn=Complex(cos(PI/i),flag*sin(PI/i));
        for(int j=0; j<n; j+=(i<<1)) //编号为j的结点
        {
            Complex w=Complex(1,0);
            for(int k=0; k<i; k++,w=w*wn)
            {
                Complex x=a[j+k],y=w*a[j+k+i]; //x和y是当前节点左右递归下去后得到的fft值，当然，y(右儿子，奇数项)的fft值要乘上一个w
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
    if(flag==-1) //如果是插值，还得在前面加个1/n的系数
        for(int i=0; i<n; i++)
            a[i].real/=n;
}

void Fast_Fourier_Transform(char *arr1, char *arr2)
{
    int i;
    n=strlen(arr1)-1;
    m=strlen(arr2)-1;

    for(i=0; i<=n; i++) x[i]=Complex(arr1[n-i]-'0', 0);
    for(i=0; i<=m; i++) y[i]=Complex(arr2[m-i]-'0', 0);
    m+=n;
    for(n=1; n<=m; n=n*2);

    fft(x,n,1);                                 //求值
    fft(y,n,1);                                 //求值

    for(i=0; i<n; i++) x[i]=x[i]*y[i];          //点值乘法
    fft(x,n,-1);                                //插值
}

int main()
{
    while(~scanf("%s%s", arr1, arr2))
    {
        int i;
        for(i=0; i<300005; i++) x[i]=Complex(0, 0);
        for(i=0; i<300005; i++) y[i]=Complex(0, 0);

        Fast_Fourier_Transform(arr1, arr2);

        memset(ans, 0, sizeof(ans));

        for(i=0; i<n || ans[i]; i++)
        {
            ans[i]+=(x[i].real+0.5);
            ans[i+1]+=(ans[i]/10);
            ans[i]%=10;
        }

        while(i>0 && !ans[i])i--;
        for(; i>=0; i--) printf("%d", ans[i]);
        printf("\n");
    }
    return 0;
}