BZOJ2179: FFT快速傅立叶 FFT

题目大意：给两个60000位以内的整数，求乘积
快速傅立叶变换入门题。递归常数巨大，因此需要采取迭代写法。
考虑进行递归FFT的过程：
①：将偶数为划分到左边，奇数位划分到右边
②：递归左右两边（设为a,b），设当前长度为len，则c[i]=a[i]+wi * b[i]，c[i+(len/2)]=a[i]-wi * b[i]
观察这个过程，如果把b当成在a右边len/2个长度，可以发现每次操作所需要的两个数，操作后又正好落回了这两个位置。这说明操作可以只用一个数组完成。将递归过程想象成一棵树，根节点长度为n，有两个长度为n/2的儿子，如此一直到底。。。如果能预处理出每个元素在最底层的位置，就可以自下向上递推。
通过前人的观察与证明，对于长度为2^n的序列，初始第i项在最底层的位置是i的后n位二进制反转。可以用pos[i]=(pos[i>>1]>>1);if(i&1) pos[i]|=maxlen>>1;来O(n)递推。
迭代过程：第一层循环自下到上枚举每一层的节点长度i，最底层不用管，从2一直枚举到maxlen，每次乘以2，同时计算出该层的单位根wi。设step=i/2，则第x个元素与第x+step个元素配对。第二层循环枚举每个节点j，从0到 < maxlen，每次增加i，同时初始化w=1。第三层循环处理当前节点的所有元素k，从j到 < j+step（也就是这个节点的前一半），每次把k和k+step配对，使a=c[k],b=w*c[k+step]，则c[k]=a+b,c[k+step]=a-b，同时将w乘上单位根wi。
不要忘记IDFT后每一位要除以maxlen。
注意由于算的是高精度乘法，还要对每一位取整后进行十进制进位。

#include<cstdio>#include<cmath>#include<algorithm>#define gm 1<<17using namespace std;struct complex{    double a,b;    complex(const double &a=0.0,const double &b=0.0):a(a),b(b){}    complex operator + (const complex &y) const    {        return complex(a+y.a,b+y.b);    }    complex operator - (const complex &y) const    {        return complex(a-y.a,b-y.b);    }    complex operator * (const complex &y) const    {        return complex(a*y.a-b*y.b,a*y.b+b*y.a);    }}a[gm],b[gm],c[gm];const double pi=acos(-1);const double DFT=2.0,IDFT=-2.0;double transform_mode;size_t pos[gm];inline void fast_transform(complex x[],size_t len){    for(size_t i=0;i<len;++i)    if(i<pos[i]) swap(x[i],x[pos[i]]);    for(size_t i=2;i<=len;i<<=1)    {        size_t step=i>>1;        complex wm(cos(2*pi/i),sin(transform_mode*pi/i));        for(size_t j=0;j<len;j+=i)        {            size_t lim=j+step;            complex w(1,0);            for(size_t 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(transform_mode==IDFT)    for(size_t i=0;i<len;++i) x[i].a/=len;}size_t n;char __c;size_t ans[gm],top=0;int main(){    scanf("%u",&n);    #define c __c    do c=getchar();while(c<'0'||c>'9');    for(size_t i=0;i<n;++i)    a[i].a=floor(c-'0'),c=getchar();    do c=getchar();while(c<'0'||c>'9');    for(size_t i=0;i<n;++i)    b[i].a=floor(c-'0'),c=getchar();    #undef c    size_t len=1;    while(len<(n<<1)) len<<=1;    for(size_t i=0;i<len;++i)    {        pos[i]=pos[i>>1]>>1;        if(i&1) pos[i]|=len>>1;    }    transform_mode=DFT;    fast_transform(a,len);    fast_transform(b,len);    for(size_t i=0;i<len;++i) c[i]=a[i]*b[i];    transform_mode=IDFT;    fast_transform(c,len);    for(size_t i=(n<<1)-2;~i;--i) ans[top++]=floor(c[i].a+0.5);    for(size_t i=0;i<top;++i) ans[i+1]+=ans[i]/10,ans[i]%=10;    if(ans[top]) printf("%u",ans[top]);    for(size_t i=top-1;~i;--i) printf("%u",ans[i]);    putchar('\n');    return 0;}