FFT的迭代程序实现——hdu1402

    《快速傅里叶变换FFT的迭代实现》描述了最简单的FFT的迭代实现，在此基础上可以用它进行大整数乘法或者多项式乘法。不过，还需要考虑IDFT的快速实现。IDFT有2种实现方式。第一种仿照FFT，观察IDFT的定义式，和DFT的定义本质上没有区别，利用单位复根的性质可以写出IFFT。第二种方法则利用共轭的性质：

      即：逆变换等于共轭的变换的共轭，所以可以复用FFT的实现代码。关键算法确定后，还要实现诸如确定补齐长度、位倒置、复数运算、输入输出转换等具体问题。hdu1402，用FFT做大数乘法。

     当然，这个程序还有很多可以优化的地方，在循环控制等方面，优化以后的效果也相当显著。大家可以去UOJ的模板题看看跑的最快的FFT是怎么写的。

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

#define EPS 1E-6
#define is0(x) ( -EPS < (x) && (x) < EPS )

double const PI = acos(-1.0);

//定义复数及其运算
struct complex_t{
    double real;
    double imag;
    complex_t(double a=0.0,double b=0.0):real(a),imag(b){}
    complex_t(complex_t const&rhs):real(rhs.real),imag(rhs.imag){}
};

complex_t operator + (complex_t const&lhs,complex_t const&rhs){
    return complex_t( lhs.real + rhs.real, lhs.imag + rhs.imag );
}

complex_t operator - (complex_t const&lhs,complex_t const&rhs){
    return complex_t( lhs.real - rhs.real, lhs.imag - rhs.imag );
}

complex_t operator * (complex_t const&lhs,complex_t const&rhs){
    return complex_t(
            lhs.real * rhs.real - lhs.imag * rhs.imag,
            lhs.real * rhs.imag + lhs.imag * rhs.real
        );
}

//定义全局数组
#define SIZE 65536
int Size = 0;
char A[SIZE];
char B[SIZE];
char C[SIZE+SIZE] = {'\0'};
complex_t Xa[SIZE+SIZE],Ya[SIZE+SIZE];
complex_t Xb[SIZE+SIZE],Yb[SIZE+SIZE];
complex_t Xc[SIZE+SIZE],Yc[SIZE+SIZE];

//将字符串转为复数序列，n为序列长度，非字符串长度
void string2Complex(char const s[],int n,complex_t c[]){
    fill(c,c+n,complex_t());
    int len = strlen(s);
    for(int i=0;s[i];++i){
        c[i].real = (double)(s[len-1-i]-'0');
    }
    return;
}

//雷德算法，调整系数位置，n为数组长度，从0开始
void Rader(complex_t const a[],int n,complex_t b[]){
    copy(a,a+n,b);
    for(int i=1,j=n>>1,k;i<n-1;++i){
        if ( i < j ) swap(b[i],b[j]);
        k = n >> 1;
        while( j >= k ) j -= k, k >>= 1;
        if ( j < k ) j += k;
    }
}

//FFT，n为序列长度
void FFT(complex_t const x[],int n,complex_t y[]){
    Rader(x,n,y);
    //最外层循环
    for(int i=1;(1<<i)<=n;++i){
        int m = 1 << i;
        complex_t wm(cos(2.0*PI/(double)m),-sin(2.0*PI/(double)m));

        //中间循环n/m次
        for(int j=0;j<n;j+=m){
            complex_t w(1.0);
            for(int k=0;k<(m>>1);++k){
                complex_t t( w * y[j+k+m/2] );
                complex_t u(y[j+k]);
                y[j+k] = u + t;
                y[j+k+m/2] = u - t;
                w = w * wm;
            }
        }
    }
}


int main(){
    while( EOF != scanf("%s%s",A,B) ){
        //确定序列长度
        int la = strlen(A);
        int lb = strlen(B);
        Size = la + lb;
        for(int i=0;;++i){
            if( Size <= (1<<i) ){
                Size = 1<<i;
                break;
            }
        }
        //转换
        string2Complex(A,Size,Xa);
        string2Complex(B,Size,Xb);
        //FFT
        FFT(Xa,Size,Ya);
        FFT(Xb,Size,Yb);
        //逐项乘，顺便共轭
        for(int i=0;i<Size;++i){
            Yc[i] = Ya[i] * Yb[i];
            Yc[i].imag = - Yc[i].imag;
        }
        //FFT
        FFT(Yc,Size,Xc);
        Xc[Size].real = 0.0;

        //再对Xc共轭除以N得到结果，但此处只关心实部
        for(int i=0;i<Size;++i) Xc[i].real /= (double)Size;
        //将复数序列转化为合理的输出
        int carry = 0;
        int head = SIZE + SIZE - 2;
        for(int i=0;i<Size;++i){
            int t = (int)(Xc[i].real + 0.5) + carry;
            C[SIZE+SIZE-2-i] = (char)(t % 10) + '0';
            if ( t % 10 ) head = SIZE + SIZE - 2 - i;
            carry = t / 10;
        }
        if ( carry ) head = SIZE + SIZE - 1 - Size;
        while( carry ){
            --head;
            C[head] = '0' + (char)(carry%10);
            carry /= 10;
        }
        printf("%s\n",C+head);
    }
    return 0;
}