快速傅里叶变换 FFT 模板【bzoj2179】 FFT快速傅立叶

来源:互联网 发布:淘宝网屏风富贵牡丹 编辑:程序博客网 时间:2024/05/17 09:23

你眼中看似落叶纷飞变化无常的世界,实际只是躺在上帝怀中一份早已谱好的乐章。                                                            ——韩昊

时域是什么呢?正常人在观察这个世界都是已时间为参照,就比如声音,最简单的声音在物理中的理解就是:
这里写图片描述
这样一条以时间为横轴的美美哒的正弦函数,这就是声音在时域中的表现方式。

那什么是频域呢?那就是忽略掉时间来观察这个世界咯。
要如何忽略掉时间呢?最简单的方法就是你站在时间轴上,你不就看不到它了吗o(∩▽∩)o~
就比如说在频域上看这个正弦函数,站在时间轴正无穷处,沿着时间轴的方向来看它,它是什么样的呢(很好想象吧)?
这里写图片描述

诶,它差不多就这样(我就说很好想象)……

——然后我们回归正题——

对于一个高次函数:
比如: f(x)=2x²+x+3 (这个就是时域上的这个函数啦)
那频域上怎么表示这个函数呢?
我们知道三个点可以求出一条抛物线,那么这个函数就可以用三个点来表示:
(1 , 6)(2 ,13)( 3, 24)
这就是这个函数在频域的表示方法啦。
对于其他高次函数,几次函数就用几个点表示,就是点值表示法。

在时域中的一个奇妙无比的表达式,在频域中也就是几个点……

然而恰好在两个高次函数在时域中相乘和在频域中相乘是等效的(如果不理解找几个二次函数,乘一乘,代一代就明白了)。

然而你用表达式算乘积的话,要花O(n^2) 的时间,而用点值表示法相乘只需要话费O(n)的时间,所以在频域中计算,明显效率更高啊!

可是需要的实在时域中的表达式,所以就涉及到了两个域中的相互转换。这个理解起来简单,给你一个表达式,然后求上面的随便n个,或者给你n个点,求一个n次表达式。宝宝初中就会了!

但是在计算机中按照数学方式的代入法求点,和设未知数然后解出所有的系数显然是不现实的,因为太慢了……这样转完了比原来n^2的算法还慢,就得不偿失了。

所以选择点的时候就不能随便选了,要选择有代表性的点。
这时引入了单位根的概念。
对于 x^n=1 成立的所有x都为n次单位根,当然x是在复数范围内的。
就比如说当n等于4的时候x就可以取: 1,-1,i,-i。
下图是n等于8时的两个单位根。
这里写图片描述
因为单位圆上的每一个点到原点的距离都为1,所以把这个圆n等分,就是我们要求的n单位根了,这n个单位根就是我们要选择的n个点。

因为这些点具有一些奇怪的性质,可以在nlogn的时间内实现高次函数在时域与频域内的变换。
具体的实现方法可以看代码,但是至于为什么这样做和怎么做,笨笨的我只能感性的理解,不能够解释清楚QAQ。
如果想要透彻的领悟,可以去picks的博客学习。
(果然这篇博客最后还是个大坑,隔靴搔痒,重要的东西全都没有说QAQ)

———END————

题目大意:
给两个巨大的数,求乘积。

题目分析:(看题目猜题意:FFT)

把两个巨大的数看成两个x=10的高次函数,然后快速傅里叶变换求乘积输出系数就好了(注意最后要进位哦)。

此代码借鉴自某非常厉害的神犇的博客,戳进去可了解更多。

代码如下:

#include<cstdio>#include<algorithm>#include<cmath>#define N 1<<17using namespace std;const double pi=acos(-1);const double DFT=2.0,IDFT=-2.0;struct complex{    double a,b;    complex(const double &a=0.0,const double &b=0.0):a(a),b(b){}    complex operator + (const complex &c) const {return complex(a+c.a,b+c.b); }    complex operator - (const complex &c) const {return complex(a-c.a,b-c.b); }    complex operator * (const complex &c) const {return complex(a*c.a-b*c.b,a*c.b+b*c.a); }}factor1[N],factor2[N],product[N];int pos[N];void initialization(int len){    for(int i=0;i<len;i++)    {        pos[i]=pos[i>>1]>>1;        if(i&1) pos[i]|=len>>1;    }    return;}void Fast_Fourier_Transform(complex x[],int len,int mode){    for(int i=0;i<len;i++)    if(i<pos[i]) swap(x[i],x[pos[i]]);    for(int i=2;i<=len;i<<=1)    {        int step=i>>1;        complex wm(cos(2*pi/i),sin(mode*pi/i));        for(int j=0;j<len;j+=i)        {            int lim=j+step;            complex w(1,0);            for(int k=j;k<lim;k++)            {                complex a=x[k],b=w*x[k+step];                x[k]=a+b,x[k+step]=a-b;                w=w*wm;            }        }    }    if(mode==IDFT)    for(int i=0;i<len;i++) x[i].a/=len;    return;}#define FFT Fast_Fourier_Transformint n,ans[N],top=0;char s[N];int main(){    scanf("%d",&n);    int len=1;    while(len<(n<<1)) len<<=1;    initialization(len);    scanf("%s",s);    for(int i=0;i<n;i++) factor1[i].a=s[i]-'0';    scanf("%s",s);    for(int i=0;i<n;i++) factor2[i].a=s[i]-'0';    FFT(factor1,len,DFT);    FFT(factor2,len,DFT);    for(int i=0;i<len;i++) product[i]=factor1[i]*factor2[i];    FFT(product,len,IDFT);    for(int i=(n<<1)-2;~i;i--) ans[top++]=floor(product[i].a+0.5);    for(int i=0;i<top;i++) ans[i+1]+=ans[i]/10,ans[i]%=10;    while(!ans[top]) top--;    for(int i=top;~i;i--) printf("%d",ans[i]);    printf("\n");    return 0;}
0 0
原创粉丝点击