快速傅里叶变换

原创

胡桃小孩儿 2022-10-11 23:00:55 博主文章分类：algorithm_数论 ©著作权

©著作权归作者所有：来自51CTO博客作者胡桃小孩儿的原创作品，请联系作者获取转载授权，否则将追究法律责任

快速傅里叶变换不能三言两语能解释清楚，自己看了一些资料，仍不敢说完全掌握了。
快速傅里叶变换 (FFT)的作用及解释：

http://blog.jobbole.com/70549/

编程实现：

http://wenku.baidu.com/link?url=Ntpg6kmz98PiLYCo4ymYCNDEQ2iziPIaq4sCqhbxLuZ16ONotUgwLvJ0q8AG-mpTlpoCFr-M7dkf0iZ7tnws45W8Z-JD2g_0VgoSblbU4jK###

大神们的博客：

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-15

卷积是两个变量在某范围内相乘后求和的结果。如果卷积的变量是序列x(n)和h(n)，则卷积的结果

快速傅里叶变换_#include

，其中星号*表示卷积。

多项式乘法就可被FFT高效解决。

快速傅里叶变换提供了一种新的视觉，从不同的角度看待问题，让复杂问题简单化。把数字统统放到复数域内，由单位复根W的三大特性，对称性，可约性，周期性实现普通离散变换效率上的飞跃。

复根：

快速傅里叶变换_#include_02

对称性奇偶分类，分治处理，此外

周期性：

快速傅里叶变换_fft_03

可约性：

快速傅里叶变换_#include_04

帮助减少了运算量。

快速傅里叶变换解决问题的过程可以说就是先将函数表达变成点值关系，进行运算，再转成系数表示。

A * B Problem Plus
http://acm.hdu.edu.cn/showproblem.php?pid=1402

Calculate A * B. 数据很大

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int N=5e5+10;
const double PI=acos(-1.0);
struct comx{ // 复数
    double r,i;
    comx(double r=0.0,double i=0.0){
        this->r=r;
        this->i=i;
    }
    comx operator +(const comx x){
        return comx(r+x.r,i+x.i);
    }
    comx operator -(const comx x){
        return comx(r-x.r,i-x.i);
    }
    comx operator *(const comx x){
        return comx(r*x.r-i*x.i,r*x.i+i*x.r);
    }
};
// 倒位序·雷德算法 100 --> 001
void Rader(comx F[],int len){
    int j=len/2;
    for(int i=1;i<len-1;i++){
        if(i<j) swap(F[i],F[j]);
        int k=len>>1;
        while(j>=k){
            j-=k;
            k>>=1;
        }
        if(j<k) j+=k;
    }
}
//  FFT 快速傅里叶变换
void FFT(comx F[],int len,int f){
    Rader(F,len); // 倒位序 奇偶分类
    for(int i=2;i<=len;i<<=1){  // 分治
        comx I(cos(-f*2*PI/i),sin(-f*2*PI/i));
        for(int j=0;j<len;j+=i){
            comx w(1,0); // 旋转因子
            for(int k=j;k<j+i/2;k++){
                comx u=F[k];
                comx t=w*F[k+i/2];
                F[k]=u+t;
                F[k+i/2]=u-t;   //蝴蝶
                w=w*I;
            }
        }
    }
    if(f==-1) {
        for(int i=0;i<len;i++) F[i].r/=len;
    }
}
// 卷积
void Conv(comx a[],comx b[],int len){
    FFT(a,len,1); // DFT 变成点值
    FFT(b,len,1);
    for(int i=0;i<len;i++) a[i]=a[i]*b[i];
    FFT(a,len,-1);  // IDFT 变回关系式
}

char str1[N],str2[N];
comx a[N],b[N];
int ans[N],len;

int main()
{
    while(~scanf("%s%s",str1,str2)){
        int len1=strlen(str1),len2=strlen(str2);
        len=1;
        while(len<2*len1 || len<2*len2) len<<=1;
        int i;
        for(i=0;i<len1;i++){  //大数存入两个复数数组中
            a[i].r=str1[len1-1-i]-'0';
            a[i].i=0;
        }
        while(i<len){
            a[i].r=a[i].i=0;
            i++;
        }
        for(i=0;i<len2;i++){
            b[i].r=str2[len2-1-i]-'0';
            b[i].i=0;
        }
        while(i<len){
            b[i].r=b[i].i=0;
            i++;
        }
        Conv(a,b,len);
        for(i=0;i<len;i++) ans[i]=a[i].r+0.5; //四舍五入
        for(i=0;i<len;i++){
            ans[i+1]+=ans[i]/10; //高精度处理
            ans[i]%=10;
        }
        int length=0;
        for(i=len-1;i>=0;i--){
            if(ans[i]) {
                length=i;
                break;
            }
        }
        for(i=length;i>=0;i--) printf("%d",ans[i]);
        printf("\n");
    }
    return 0;
}

拓展：
将卷积部分改成单纯的求和，即成为高精度求和

求和（前提是两个大数均是非负数）:

void Conv(comx a[],comx b[],int len){
    FFT(a,len,1); // DFT 变成点值
    FFT(b,len,1);
    for(int i=0;i<len;i++) a[i]=a[i]+b[i];
    FFT(a,len,-1);  // IDFT 变回关系式
}

但是不能拓展成为高精度减法：

void Conv(comx a[],comx b[],int len){
    FFT(a,len,1); // DFT 变成点值
    FFT(b,len,1);
    for(int i=0;i<len;i++) a[i]=a[i]-b[i];
    FFT(a,len,-1);  // IDFT 变回关系式
}