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

---

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

---

时域是什么呢？正常人在观察这个世界都是已时间为参照，就比如声音，最简单的声音在物理中的理解就是：

这样一条以时间为横轴的美美哒的正弦函数，这就是声音在时域中的表现方式。

那什么是频域呢？那就是忽略掉时间来观察这个世界咯。
要如何忽略掉时间呢？最简单的方法就是你站在时间轴上，你不就看不到它了吗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的时间内实现高次函数在时域与频域内的变换。
具体的实现方法可以看代码，但是至于为什么这样做和怎么做，笨笨的我只能感性的理解，不能够解释清楚ＱＡＱ。
如果想要透彻的领悟，可以去picks的博客学习。
（果然这篇博客最后还是个大坑，隔靴搔痒，重要的东西全都没有说QAQ）

———END————

题目大意：
给两个巨大的数，求乘积。

题目分析：（看题目猜题意：FFT）

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

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

代码如下：

#include<cstdio>
#include<algorithm>
#include<cmath>
#define N 1<<17
using 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_Transform

int 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;
}