ATC 001 C 高速フーリエ変換

2015/06/20 (Sat) ATC AtCoder 数学

問題

日本語

方針

スライドの内容を理解して実装する．とても勉強になりました． $e^n$ が全ての元を網羅し， $e^1$ に戻ってくるような集合ならなんでもできるらしい．(例えば素数の冪乗) これなら誤差死しなさそう．

実装

typedef double Real;
typedef complex<Real> Complex;
 
Real const PI = acos(-1);
 
namespace FastFourierTransform {
    Complex z_[1<<19];
    Complex * z = z_+(1<<18);
    int len;
    void FastFourierTransform(vector<Complex> & f, int inv = 0){
        int n = f.size();
        if(n==1) return;
        vector<Complex> f0(n>>1), f1(n>>1);
        rep(i,n>>1) f0[i] = f[i<<1], f1[i] = f[i<<1|1];
        FastFourierTransform(f0,inv);
        FastFourierTransform(f1,inv);
        int k = len/n, m = n/2-1;
        rep(i,n) f[i] = f0[i&m] + z[inv ? -k*i : k*i] * f1[i&m];
        return;
    }
 
    vector<Complex> Convolution(vector<Complex> g, vector<Complex> h){
        int n_ = max(g.size(), h.size());
        int n = 1;
        while(n < n_) n<<=1;
        g.resize(n);
        h.resize(n);
        len = n+1;
        rep(i,len) z[i] = polar<Real>(1,2*PI/n*i), z[-i] = conj(z[i]);
        FastFourierTransform(g);
        FastFourierTransform(h);
        vector<Complex> f(n);
        rep(i,n) f[i] = g[i]*h[i];
        FastFourierTransform(f,1);
        rep(i,n) f[i] /= n;
        return f;
    }
}
 
signed main(){
    int n;
    scanf("%d",&n);
    vector<Complex> g(n*2+1), h(n*2+1);
    rep(i,n){
        int a,b;
        scanf("%d%d",&a,&b);
        g[i+1] = a, h[i+1] = b;
    }
    vector<Complex> f = FastFourierTransform::Convolution(g,h);
    for(int i=1;i<=n*2;i++) printf("%d\n",(int)(f[i].real()+0.5));
}