猫史档案馆


复杂数据处理:傅里叶变换相关

用户:爵士OIer爵士OIer查看:121 回复:88 评论:121 创建时间:2021-11-13T20:20:11


内容同时发布于我的博客(洛谷uid:246979)。

非常感谢慕斯的帮助。

center_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_imagecenter_image

//重载运算符扩域至复数 
struct Complex
{
    Complex (double xx=0,double yy=0){x=xx,y=yy;}double x,y;
    Complex operator + (Complex const &B) const {return Complex(x+B.x,y+B.y);}
    Complex operator - (Complex const &B) const {return Complex(x-B.x,y-B.y);}
    Complex operator * (Complex const &B) const {return Complex(x*B.x-y*B.y,x*B.y+y*B.x);}
}f[2700000],g[2700000];

//FFT,flag控制求值或插值(正/逆)
void fft(Complex *f,int len,bool flag)
{
    if(len==1)return;
    Complex *fl=f,*fr=f+len/2;
    for(int k=0;k<len;k++)tmp[k]=f[k];
    for(int k=0;k<len/2;k++)
        fl[k]=tmp[k<<1],fr[k]=tmp[k<<1|1];
    //递归分治求解 
    fft(fl,len/2,flag),fft(fr,len/2,flag);
    //buf为k单位根,tG为起始1次单位根 
    Complex tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0);
    if(flag)tG.y*=-1;//DFT为负,ISFT为正
    //求出而序列的点值 
    for(int k=0;k<len/2;k++)
        tmp[k]=fl[k]+buf*fr[k], 
        tmp[k+len/2]=fl[k]-buf*fr[k],
        buf=buf*tG;
    for(int k=0;k<len;k++)f[k]=tmp[k];
}

center_image

//重载运算符,复数 
struct Complex
{
    Complex (double xx=0,double yy=0){x=xx,y=yy;}double x,y;
    Complex operator + (Complex const &B) const {return Complex(x+B.x,y+B.y);}
    Complex operator - (Complex const &B) const {return Complex(x-B.x,y-B.y);}
    Complex operator * (Complex const &B) const {return Complex(x*B.x-y*B.y,x*B.y+y*B.x);}
}f[2700000],g[2700000];
//FFT,flag控制求值或插值(正/逆)
void fft(Complex *f,bool flag)
{
    //蝴蝶变换优化
    for(int i=0;i<n;i++)
        if(i<tr[i])swap(f[i],f[tr[i]]);
    //迭代,枚举处理多项式长度 
    for(int p=2;p<=n;p<<=1)
    {
        int len=p>>1;
        Complex tG(cos(2*Pi/p),sin(2*Pi/p));
        if(flag)tG.y*=-1;
        //处理以p为段长的每一段 
        for(int k=0;k<n;k+=p)
        {
            Complex buf(1,0);
            //代入单位根逐次数求解
            for(int l=k;l<k+len;l++)
            {
                Complex tt=buf*f[len+l];
                f[len+l]=f[l]-tt,f[l]=f[l]+tt,buf=buf*tG;
            }
        }
    }
}

center_image

void mul(){
    for(m+=n,n=1;n<=m;n<<=1);
    for(int i=0;i<n;i++)
        tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
    fft(f,1),fft(g,1);
    for(int i=0;i<n;i++)f[i]=f[i]*g[i];
    fft(f,0);
    for(int i=0;i<=m;i++)
        ans[i]=(int)(f[i].x/n+0.5);
}

 

祝大家学习愉快(

可以在评论区回复。


回复

上一页1 页 / 共 3下一页
Asheep233Asheep233

喵d

点赞0


评论


冷鱼闲风冷鱼闲风

爵士出贴必属精品。顶!傅立叶大学生的恶梦imgsrc="https://static.codemao.cn/emoji/codemao/%E7%BC%96%E7%A8%8B%E7%8C%AB_%E7%B4%A7%E5%BC%A0.gif"alt="emotion_编程猫_紧张"。(hh

点赞2


评论


传闻中的三叶传闻中的三叶

写的好

点赞0


评论


嗯哼君 渠源嗯哼君 渠源

我看不懂,但我大受震撼

未来金牌得主的精品教程get√

点赞2


评论


爵士OIer爵士OIer

ddd

点赞0


评论


爵士OIer爵士OIer

ddd

点赞0


评论


爵士OIer爵士OIer

分治那里应该是这样center_image

之前写的时候太晚了,没注意检查

点赞0


评论


传闻中的三叶传闻中的三叶

牛牛牛

点赞0


评论


f(hxr)f(hxr)

666

点赞0


评论


NaughtNaught

666

点赞0


评论


爵士OIer爵士OIer

ddd

点赞0


评论


AlcalaAlcala

看不懂()

点赞0


评论


kill_godkill_god

你居然还在

点赞0


评论


QBZ_QBZ_

我看不懂但我大受震撼

点赞1


评论


Architect_煎饼Architect_煎饼

《天 书》

(或许考上高中就能看懂???)

无穷小是什么emotion_doge

能讲点初一看的懂得吗()()()

点赞0


评论


妙龙QWQ妙龙QWQ

我看不懂,但是大受震汗。

点赞0


评论


_韩佳霖_韩佳霖

我看不懂但我大受震撼(膜拜橙名dalao,,我也想橙名T_T

点赞0


评论


coding小霸王coding小霸王

hardcore

点赞0


评论


不歪刻晴真君_不歪刻晴真君_

赶紧收藏起来

点赞0


评论


umbrallawwumbrallaww

我居然看懂了一半)

点赞1


评论


YOF冇人GGYOF冇人GG

what’ up,虽然我看不懂,但是我大受震撼

点赞0


评论


coding小霸王coding小霸王

爵士几年级了啊

点赞0


评论


芜湖咸鱼芜湖咸鱼

我看不懂但我大受震撼

点赞3


评论


NomandNomand

我看不懂,但我大受震撼

点赞3


评论


_韩佳霖_韩佳霖

爵士OIer_Chtholly 2021-11-27 10:02:39 授予  自由发言 权限
满一年。   啊这,,,

点赞0


评论


_韩佳霖_韩佳霖

能对一个5年级蒟蒻友善一点吗()

点赞0


评论


史莱姆awa史莱姆awa

我看不懂,但我也看不懂emotion_doge

点赞0


评论


史莱姆awa史莱姆awa

就代码能看懂些emotion_doge

点赞0


评论


离开的十五分之六老咸鱼离开的十五分之六老咸鱼

我14年白活了()

点赞0


评论


不知道该叫啥的一只萌新不知道该叫啥的一只萌新

我看不懂但我大受震撼

点赞0


评论