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;
}
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 158,847评论 4 362
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,208评论 1 292
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 108,587评论 0 243
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,942评论 0 205
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,332评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,587评论 1 218
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,853评论 2 312
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,568评论 0 198
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,273评论 1 242
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,542评论 2 246
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,033评论 1 260
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,373评论 2 253
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,031评论 3 236
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,073评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,830评论 0 195
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,628评论 2 274
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,537评论 2 269

推荐阅读更多精彩内容