多项式相乘快速算法原理及相应C代码实现---用到fft

最近认真研究了一下算法导论里面的多项式乘法的快速计算问题，主要是用到了FFT，自己也实现了一下，总结如下。

1.多项式乘法

两个多项式相乘即为多项式乘法，例如：3*x^7+4*x^5+1*x^2+5与8*x^6+7*x^4+6*x^3+9两个式子相乘，会得到一个最高次数项为13的多项式。一般来说，普通的计算方法是：把A多项式中的每一项与B中多项式中的每一项相乘，得到n个多项式，再把每个多项式相加到一起，得到最终的结果，不妨假设A，B的最高次项都为n-1，长度都为n，那么计算最终的结果需要o(n^2)时间复杂度。而使用快速傅里叶变换（FFT），则可以将时间复杂度降低到o(nlog n)。这是因为，对一个复数序列做正/反快速傅里叶变换的时间复杂度都是o(nlog n)，而变换后的序列逐项相乘即为原序列做多项式乘法的结果（多项式乘法相当于卷积）。所以，FFT可以降低多项式相乘运算的时间复杂度，具体的解释和证明在《算法导论》或者其他任何一本相关的算法书中都有详细描述，在此不再赘述。另外，需要注意的是，一些其他的运算也可以转化成多项式乘法，进而利用FFT来加快运算。例如：1.数字乘法运算，和多项式乘法类似，A*B的操作就是用A每一位上的数字乘以B每一位上的数字。尤其是在大数乘法中，FFT可以大幅度加快运算。2.给定A到B的不同长度路径详细数据，B到C的不同路径长度详细数据，求A到C不同长度路径的数量。可以把A到B和B到C不同长度的路径看成不同次数的项，例如：A到B有3条长度为4,2条长度为5的路径，B到C有1条长度为2,4条长度为3的路径，那么A到C不同长度路径的数量等于(3*x^4+2*x^5)*(4*x^3+1*x^2)得到的各项的系数，转化成多项式相乘问题之后，就可以利用FFT来加快运算速度了。

2.FFT

大多数人应该是只需要会用FFT即可，但是这个算法比较基础，因此我自己编程实现了一下，总的代码只有150行左右，其实不算长，当然，对输入序列长度不是2的整数次幂这种情况我没有相应的预处理，算是偷懒了。其实只要对长度取一下对数即可，例如：输入长度如果是37，首先把37*2，拓展成74（这个是FFT必须的），然后对74取log₂上取整即可，得到2⁷=128，因此在74后面再添加54个0。

另外，需要注意的一点是，reverse函数在FFT和IFFT中是必须的，不过鉴于多项式乘法需要成对进行FFT和IFFT，所以在做多项式乘法的时候，reverse应该是可以省略的（当然，这个函数的耗时很小）。w的值应该提前计算出来，这样在FFT和IFFT中蝶形计算每一项的时候，就不用重复计算w了，可以节省很多时间。

具体FFT的原理和解释，可以查维基、信息论、数字信号处理、随机过程等任一领域的教科书。

3.代码

程序主要包括FFT，IFFT函数，以及一些复数运算

3.1复数的定义和相关运算定义

//复数
struct Complex{
    double real;
    double image;
};
Complex a1[MAX_SIZE],a2[MAX_SIZE],result[MAX_SIZE],w[MAX_SIZE];
//复数相乘计算
Complex operator*(Complex a,Complex b){
    Complex r;
    r.real=a.real*b.real-a.image*b.image;
    r.image=a.real*b.image+a.image*b.real;
    return r;
}
//复数相加计算
Complex operator+(Complex a,Complex b){
    Complex r;
    r.real=a.real+b.real;
    r.image=a.image+b.image;
    return r;
}
//复数相减计算
Complex operator-(Complex a,Complex b){
    Complex r;
    r.real=a.real-b.real;
    r.image=a.image-b.image;
    return r;
}
//复数除法计算
Complex operator/(Complex a,double b){
    Complex r;
    r.real=a.real/b;
    r.image=a.image/b;
    return r;
}
//复数虚部反计算
Complex operator~(Complex a){
    Complex r;
    r.real=a.real;
    r.image=0-a.image;
    return r;
}

3.2FFT和IFFT函数及相关函数

其实FFT和IFFT的原理一样，只是IFFT多了一个除法步骤，也可以把两个合并成一个函数。Reverse用于重新排列输入数组的元素下标，例如输入数组长度为8，则0,1,2,3,4,5,6,7下标的元素经过重新排列后变为0,4,2,6,1,5,3,7下标的元素。Compute_W用于预先计算FFT中需要的w值。

//重新排列方法2，效率较高
void Reverse(int* id,int size,int m){
    for(int i=0;i<size;i++){
        for(int j=0;j<(m+1)/2;j++){
            int v1=(1<<(j)&i)<<(m-2*j-1);
            int v2=(1<<(m-j-1)&i)>>(m-2*j-1);
            id[i]|=(v1|v2);
        }
    }
};

//重新排列方法1，该方法是用pow函数效率比较低
void Reverse(int* id,int size,int m){
    for(int i=0;i<size;i++){
        for(int j=0;j<m;j++){
            int exp=(i>>j)&1;
            id[i]+=exp*(int)pow((double)2,(double)(m-j-1));
        }
    }
};

//计算并存储需要乘的w值
void Compute_W(Complex w[],int size){
    for(int i=0;i<size/2;i++){
        w[i].real=cos(2*PI*i/size);
        w[i].image=sin(2*PI*i/size);
        w[i+size/2].real=0-w[i].real;
        w[i+size/2].image=0-w[i].image;
    }
};
//快速傅里叶
void FFT(Complex in[],int size){
    int* id=new int[size];
    memset(id,0,sizeof(int)*size);
    int m=log((double)size)/log((double)2);
    Reverse(id,size,m);    //将输入重新排列，符合输出
    Complex *resort= new Complex[size];
    memset(resort,0,sizeof(Complex)*size);
    int i,j,k,s;
    for(i=0;i<size;i++)
        resort[i]=in[id[i]];
    for(i=1;i<=m;i++){
        s=(int)pow((double)2,(double)i);
        for(j=0;j<size/s;j++){
            for(k=j*s;k<j*s+s/2;k++){        
                Complex k1=   resort[k]+w[size/s*(k-j*s)]*resort[k+s/2];
                resort[k+s/2]=resort[k]-w[size/s*(k-j*s)]*resort[k+s/2];
                resort[k]=k1;
            }
        }
    }
    for(i=0;i<size;i++)
        in[i]=resort[i];
    delete[] id;
    delete[] resort;
};
//快速逆傅里叶
void IFFT(Complex in[],int size){
    int* id=new int[size];
    memset(id,0,sizeof(int)*size);
    int m=log((double)size)/log((double)2);
    Reverse(id,size,m);    //将输入重新排列，符合输出
    Complex *resort= new Complex[size];
    memset(resort,0,sizeof(Complex)*size);
    int i,j,k,s;
    for(i=0;i<size;i++)
        resort[i]=in[id[i]];
    for(i=1;i<=m;i++){
        s=(int)pow((double)2,(double)i);
        for(j=0;j<size/s;j++){
            for(k=j*s;k<j*s+s/2;k++){
                Complex k1=(resort[k]+(~w[size/s*(k-j*s)])*resort[k+s/2]);
                resort[k+s/2]=(resort[k]-(~w[size/s*(k-j*s)])*resort[k+s/2]);
                resort[k]=k1;
            }
        }
    }
    for(i=0;i<size;i++)
        in[i]=resort[i]/size;
    delete[] id;
    delete[] resort;
};

3.3主函数

输入两个多项式的系数（长度必须都是2的整数次幂），输出两个多项式相乘的结果

int main(){
    //输入两个多项式数列
    int size,size1,size2,i;
    memset(a1,0,sizeof(a1));
    memset(a2,0,sizeof(a2));
    memset(w,0,sizeof(w));
    memset(result,0,sizeof(result));
    scanf("%d%d",&size1,&size2);
    for(i=0;i<size1;i++)
        scanf("%lf",&a1[i].real);
    for(i=0;i<size2;i++)
        scanf("%lf",&a2[i].real);
    size=size1>size2?size1*2:size2*2;
    Compute_W(w,size);
    FFT(a1,size);
    FFT(a2,size);
    for(i=0;i<size;i++)
        result[i]=a1[i]*a2[i];
    IFFT(result,size);
    for(i=0;i<size1+size2-1;i++)
        printf("%.2lf ",result[i].real);
    printf("\n");
    return 0;
}

下面是完整的代码

CPP文件下载