bzoj_2179 FFT快速傅立叶

题目:

  给出两个n位10进制整数x和y,你需要计算x*y。

思路:

  • 关于快速傅立叶变换:
      设 。如果想要得到这样两个多项式的乘积 。最朴素的办法就是 的。但是如果把 看成一个 次的函数 ,用 个点对 去表示 ,这就是 的点值表示(显然它们之间是唯一确定的)。也用同样的方法表示 ,那么这样就能 地得到 的点值表示:。而快速傅立叶变换就是为了快速的在多项式和点值表示之间转换。
      首先假设 ,(如果不是就在后面补 )。快速傅立叶变换的特殊点就是,它代入的 个均匀分布在单位圆上的复数,当 时就是下图:
      
      根据复数的运算,如果从最左边的点开始编号为 ,那么 。设代入 得到的一组值为 。再把 代入得到一组系数为 ,那么 。这样就能把点值表示变回多项式了,推导就参考网上的文章吧(推荐小学生都能看懂的FFT!!!)。
      但是到这里计算的复杂度还是 的,要考虑优化,对于一个多项式 ,如果设 。那么 这样就能通过分治 变换了。具体实现的时候一般会预处理出最后被换到的位置 ,找规律可以发现 (也就是二进制反过来),这样就能迭代实现了。
  • 关于此题:
      把一个数字看成 。这样就能用上面的变换了。

代码:

#include <bits/stdc++.h>
#define mset(a,b) memset(a,b,sizeof a)
#define _(...) (void)(__VA_ARGS__)
#define rep(i,a,b) for(int i(a),i##_END_(b);i<=i##_END_;i++)
#define drep(i,a,b) for(int i(a),i##_END_(b);i>=i##_END_;i--)
template<class T> inline bool tomax(T &a,T b){return a<b?a=b,1:0;}
template<class T> inline bool tomin(T &a,T b){return b<a?a=b,1:0;}
const double PI=acos(-1);
const int N=131072+5;

struct Complex{
    double r,i;
    inline Complex operator + (const Complex _) const {
        return (Complex){r+_.r,i+_.i};
    }
    inline Complex operator - (const Complex _) const {
        return (Complex){r-_.r,i-_.i};
    }
    inline Complex operator * (const Complex _) const {
        return (Complex){r*_.r-i*_.i,r*_.i+i*_.r};
    }
};

int rev[N];
void fft(Complex a[],Complex w[],int n){
    rep(i,0,n-1)
        if(i<rev[i])
            std::swap(a[i],a[rev[i]]);
    for(int l=2;l<=n;l<<=1){
        int m=l/2;
        for(Complex *p=a;p!=a+n;p+=l)
            rep(i,0,m-1){
                Complex t=w[n/l*i]*p[i+m];
                p[i+m]=p[i]-t;
                p[i]=p[i]+t;
            }
    }
}

Complex a[N],b[N],inv[N],omg[N];
void init(int &n){
    int lim=0;
    while((1<<lim)<n)lim++;
    n=(1<<lim);
    rep(i,0,n-1)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1));
    rep(i,0,n-1)
        omg[i]=(Complex){cos(PI*2/n*i),sin(PI*2/n*i)};
    rep(i,0,n-1)
        inv[i]=omg[(n-i)%n];
}

#define readChar(c) do (c)=getchar();while(c=='\n'||c=='\r'||c==' ')
int res[N];
int main(){
    int n;char c;
    scanf("%d",&n);
    drep(i,n-1,0){
        readChar(c);
        a[i]=(Complex){c-'0',0};
    }
    drep(i,n-1,0){
        readChar(c);
        b[i]=(Complex){c-'0',0};
    }

    init(n<<=1);
    fft(a,omg,n);fft(b,omg,n);
    rep(i,0,n-1)
        a[i]=a[i]*b[i];
    fft(a,inv,n);
    rep(i,0,n-1)
        a[i].r/=n;

    rep(i,0,n-1)
        res[i]=a[i].r+0.5;
    rep(i,0,n-2){
        res[i+1]+=res[i]/10;
        res[i]%=10;
    }
    int p=n-1;
    while(res[p]>=10){
        res[p+1]+=res[p]/10;
        res[p++]%=10;
    }
    while(res[p]==0&&p)p--;
    drep(i,p,0)
        putchar(res[i]+'0');
    puts("");
    return 0;
}