FFT Fast-Fourier-Transform（白话粗口18X） .

由于受到良心谴责，(用模版AC)，所以今日全日搞FFT OTZ



终于搞掂！！！！！！首先自我发泄一下！！！！！！！！！！

连续AC5次先！！！！！！！！！！



以下内容涉及暴力等18X情节，需家长陪同，纯粹9UP，有七雷同咩{= =}|||



Fast-Fourier-Transform 又称为FFT，又叫快速傅里叶转换｛0 0｝



系离散傅里叶正逆转换（（i）DFT）ge一种算法，系由两位牛人coolkey，turkey（冰火两鸡）

提出ge一种基于分治思想ge一种算法。转换算法时间复杂度为O（nlogn）呢度

log指2ge对数，因为系点乘，所以乘法过程系O（n）。

所以成个算法时间复杂度系lim (n~>infi) O(nlogn)+O(n)+O(nlogn)=O(nlogn)

普通高精度乘法大家都知道系O（n^2)

所以快左好西多。



想深入了解FFT，首先有几样野系要先了解ge:

一，多项式ge两种表达式：

1、系数表达式（平时ge十进制数其实系一个次数界为n，x=10 ge一个表达式f(x)=a0+a1x^1+……a2x^n-1）。

instance:123==1*10^2+2*10^1+3*10^0==(1,2,3)<-呢个就系系数表达式{!__!}。



2、点值表达式(姐系将f(x)变成直角坐标系ge一个图像，取出n个点做表达式,至于点解要n个点？你可以类比下解多元方程，要解n个元，自然要n个等式咯，系咪先~)。

instance：f(x)=1*x^2+2*x^1+3*x^0

                  x=1，f(1)=6；x=2，f(2)=11；x=10,f(10)=123(注意到10进制数都系一点)

                  {(1,6),(2,11),(10,123)}<-呢个就系点值表达式咯。



二,矩阵里面ge：

1、向量（呢度ge向量……）。

2、矩阵乘法。（百度啦）

3、范德蒙德矩阵。(自己百度睇下)

4、逆矩阵*正矩阵=n*n单位矩阵

其实矩阵系同逆转换有关噶，特别系3同4，简直就系神器{>0<}Y



三，单位虚根（呢个就系FFTge关键）：

定义：满足w^n=1ge复数。所以n次ge单位复根就有n个喇，距地就系e^(2*pi*ik/n),其中i为虚部果个野，k=0~n-1。呢旧野仲有周期性。

用三脚函数表示就系：e^(u*i)=cos(u)+i*sin(u) (就系用呢个表示复根咯)

虚根三神技（我系签名度提过）：

1、相消引理:对任何w^(d*k)==w^(k)(姐系有周期咯)

instance：n=4，k=0~4,w^0=1,w^2=i,w^3=-i,w^4=1



2、折半引理:对任何n>0,w^(n/2)=-1（呢个就系快ge原因）。

instance：w^[(k+n/2)]*2=w^(2*k+n)=w^(2*k)*w^n=w^(2*k)=w^(k)*2（神变，甘样ge话，求出w^k就知道w^(k+n/2)喇，所以FFT就系2分法咯，hoho{^__^}）



3、求和引理:对任意整数n>=1同埋吾可以比n整出ge非零整数k，有∑(j=0~n-1) w^(k*j)=0。(呢个要黎推IDFT公式ge)



四，二进制数反转置换（instance：1100->0011）(呢个系递归转递推ge关键，自己捻啦呢，可以当一条小题做，不过要做出O（logn）噶窝，唔好直接O(n^2)窝，O(n^2)你转毛啊)

计埋呢个总共都系：O(nlogn)



大家唔好觉得好深奥.其实同高数一9样……



跟住讲DFT（其实DFT同iDFT一9样ge，所以我剩系用做一个函数就得喇。）：

Discrete-Fourier-Tranform（离散傅里叶变换）

姐系系数表达式同点值表达式ge互相转换。

其实都唔洗我讲，好简单：

instance：对于系数表达式（1，2，3）转 {(1,6),(2,11),(10,123)}，其实就系搵左三个数代入去多项式度，就搞定喇，姐系：

（a0，a1……an-1）<->{(x1,y1),(x2,y2)……(xn-1,yn-1)}

直接拿点黎代ge话，要O（n^2）ge时间（因为n个x同n个系数相乘）

所以直接转ge话你不如返屋企拾屎好过。（OTZ）



所以为冰火两鸡就提出左FFT：

其实就系利用左虚根ge周期性同折半引理，计到一半就知道另一半，所以同快速排序一样用二分法，姐系都系O(nlogn)。

原理我理解左好耐，终于知道距搞咩……

下面比较抽象，睇唔明，系我ge错{T___T}

首先大家捻下：



将一条n界ge多项式：y（x）=a0+a1*x+a2*x^2..........+an-1*x^(n-1)，按下标奇偶性分成两组：

奇数组：y1(x)=a1+a3*x+a5*x^2+.............+an-1*x^(n/2-1)

偶数组：y2(x)=a0+a2*x+a4*x^2+.............+an-2*x^(n/2-1)

所以嘛：y(x)=y1(x^2)+x*y2(x^2)  .............................................................式a

系咪先，唔明自己推推咯。



跟住大家捻翻相消引理：w^k=w^(k+n/2)

所以将单位虚根w（w^n=1）代入去ge后果就系~~~~

系y(x)求n个点（w^k，y）（k=0~n-1）ge问题转化为系y1(x)，y2(x)求n/2个点（w^k，y）（k=0~n/2-1）ge问题{(.)(.)}！！

如果继续分落去ge话，后果就系~~~~~~~



求n个点变成求n/2又变成求n/4个点……………………最后分到叶子节点个阵（二叉树），y(x)就系系数本身。

某SB问：但系有两个多项式窝？求两个多项式n/2个点同求一个多项式n个点唔系一9样咩？



大家请睇下推导~~~~~

式一：y(w^k)=y1(w^(k*2))+w^k*y2(w^(k*2)) （睇下式a）

式二：y(w^(k+n/2))=y1(w^(k*2+n))+w^(k+n/2)*y2(w^(k*2+n))

                           =y1(w^(k*2))+[-w^(k)]*y2(w^(k*2))      (睇下折半引理同相消引理)

                           =y1(w^(k*2))-w^(k)*y2(w^(k*2))

留意式一和二，共同拥有三样野：

y1(w^(k*2))，w^k，y2(w^(k*2))     (w^k有个名叫螺旋因子，可能卷积就好似螺旋挂)

所以我地只要求出呢三旧野，自然就知道y(w^k)同y(w^(k+n/2))喇所以真系每层只需一般计算量，有n个数，有logn层，所以FFT就系O（nlogn）喇！！！！



所以对大数a，b进行左DFT之后，就可以进行O(n)ge点乘喇~~~



点乘之后，要再进行逆转换，变翻系数表达式.



跟住讲IDFT，又叫插值。

如果大家有百度过范德蒙德矩阵(Vn)ge话，你地就会知道系数表达式同点值表达式有一个关系：

{(x1,y1),(x2,y2)……(xn-1,yn-1)}=Vn*(a0,a1,a2....an)

所以我地只需要猥琐甘将求出ge点值表达式乘以一个Vn 逆矩阵就得喇。

y*Vn^(-1)=a*Vn*Vn^(-1)=a*In=a

其中Vn^(-1)为列矩阵，In为单位矩阵.a*In其实就系行列向量ge变换。



讲左甘多废话之后，比两条式比你地：

设aj属于(a0,a1,a2....an)，yk属于{(x1,y1),(x2,y2)……(xn-1,yn-1)}

正转换就系：(k=0~n-1) yk=∑(j=0~n-1) aj*w(k*j)｛w^n=1)｝

点乘就系：sum=∑(k=0~n-1) yak*ybk

逆转换就系：(j=0~n-1) aj=(1/n)∑(k=0~n-1) yk*w(-k*j)｛w^n=1)｝



大家注意到其实正逆转换就系a y位置互换同埋w系变左负幂，所以正逆转换其实可以用一个函数搞掂。（其实系我唔想写两个）



终于讲完鸟，下面送上我呕心沥血ge通俗易明代码(无位运算无复杂ge函数无复杂ge下标操作)

hdu 1402        AC           359ms          7136k        2277B       C++          10SGetETernal{=___=}|||

#include<cmath>
#include<string>
#include<iostream>
using namespace std;
#define pi acos(-1。0) （求圆周率）

struct complex （定义复数结构）
{
    double r,i;
    complex (){r=i=0.0;}
    complex (double real,double image)
 {
  r=real;
  i=image;
 }
    complex operator + (complex o) （定义三种虚数运算）
 {
  return complex(r+o.r,i+o.i);
 }
    complex operator - (complex o)
 {
  return complex(r-o.r,i-o.i);
 }
    complex operator * (complex o)
 {
  return complex(r*o.r-i*o.i,r*o.i+i*o.r);
 }
}op1[200001],op2[200001];    

char a[100001],b[100001];              
int sum[200001];                     //保存结果sum

void brc(complex *y,int l)        //二进制平摊反转置换O（logn）
{
    int i,j,k;
    j=l/2;
    for (i=1;i<l-1;i++)
    {
        if (i<j) swap(y[i],y[j]);      //交换互为下标反转ge两个元素，（i<j）保证只交换一次
        k=l/2;
        while (j>=k)                    //由最高位检索,遇1变0（姐系加n/2）遇到0就变1，跳出，
        {
            j-=k;
            k/=2;
        }
        if (j<k) j+=k;                 
    }
}

void fft(complex *y,int l,double on)  //FFT O(nlogn)   其中on=1时为DFT，on=-1时为IDFT
{
    int i,j,k,m;
    complex u,t;
    brc(y,l);                                         //调用反转置换
    for (m=2;m<=l;m*=2)                 //控制层数
    {
        complex wn(cos(on*2*pi/m),sin(on*2*pi/m));  //初始化单位复根,每层n吾同，所以要更新咯
        for (j=0;j<l;j+=m)                     //控制起始下表
        {
            complex w(1,0);                //初始化螺旋因子
            for (k=j;k<j+m/2;k++)        //配对过程
            {
                u=y[k];                            //u=y1(w^k) 结合上面我讲ge
                t=w*y[k+m/2];                //t=w^k*y2(w^k)
                y[k]=u+t;                        //y(w^k)=u+t
                y[k+m/2]=u-t;                //y(w^(k+n/2))=u-t
                w=w*wn;                       //更新螺旋因子
            }                                         //上面操作又叫做蝴蝶操作
        }
    } 
    if (on==-1) for (i=0;i<l;i++) y[i].r/=l;      //IDFT时公式ge 1/n
}

int main()
{
    int l1,l2,i,l;
    while (scanf("%s %s",a,b)!=EOF)
    {
        l1=strlen(a);
        l2=strlen(b);
        l=1;
        while (l<l1*2 || l<l2*2) l*=2;            //将次数界变成2^n,配合二分法同反转置换
        for (i=0;i<l1;i++)                             //倒置存入，配合我后面ge进位
        {
            op1[i].r=a[l1-i-1]-'0';
            op1[i].i=0.0;
        }
        for (;i<l;i++) op1[i].r=op1[i].i=0.0;//将多余次数界初始化为0
        for (i=0;i<l2;i++) 
        {
            op2[i].r=b[l2-i-1]-'0';
            op2[i].i=0.0;
        }
        for (;i<l;i++) op2[i].r=op2[i].i=0.0;
        fft(op1,l,1);                                     //DFT(a)
        fft(op2,l,1);                                     //DFT(b)
        for (i=0;i<l;i++) op1[i]=op1[i]*op2[i];  //点乘并将结果存入a
        fft(op1,l,-1);                                    //IDFT(a*b)
        for (i=0;i<l;i++) sum[i]=op1[i].r+0.5;   //由于FFT有精度问题，所以呢度要四舍五入
        for (i=0;i<l;i++)                               //对sum做进位操作
        {
            sum[i+1]+=sum[i]/10;
            sum[i]%=10;
        }
        l=l1+l2-1;          
        while (sum[l]<=0 && l>0) l--;       //检索最高非零位
        for (i=l;i>=0;i--) printf("%d",sum[i]);    //倒序输出
        printf("/n");
    }
    return 0;
}