haraduka's diary

やる気が欲しい

ATC001 C問題

今回は高速フーリエ変換の問題。
まさかこんな問題がフーリエ変換で解けるとは…

FFTはいくつかアルゴリズムがある。再帰のものは理解したけど、Cooley-Tukeyは理解するのが面倒だからやめた。サイトのプログラム大体そのまま使った。

typedef complex<double> Complex;
const double PI = acos(-1.0);

vector<Complex> FFT(double theta, const vector<Complex>& a){
    int n = (int)a.size();
    vector<Complex> ret = a;
    for(int m=n; m>=2; m>>=1){
        int mh = m>>1;
        for(int i=0; i<mh; i++){
            Complex w = exp(i*theta*Complex(0, 1));
            for(int j=i; j<n; j+=m){
                int k = j+mh;
                Complex x = ret[j] - ret[k];
                ret[j] += ret[k];
                ret[k] = w*x;
            }
        }
        theta *= 2;
    }
    int i=0;
    for(int j=1; j<n-1; j++){
        for(int k=n>>1; k>(i^=k); k>>= 1){}
        if(j < i) swap(ret[i], ret[j]);
    }
    return ret;
}

//畳み込み
template<class T>
vector<T> Convolution(const vector<T> &lhs, const vector<T> &rhs){
    int n = 1; //nはlhs,rhsのsizeより大きい2の累乗
    while(n < (int)max(lhs.size(), rhs.size())*2){n <<= 1;}
    vector<Complex> temp1(n);
    vector<Complex> temp2(n);
    //データ全てをinput
    for(int i=0; i<n/2; i++){
        if(i < (int)lhs.size()) temp1[i] = Complex(lhs[i], 0);
        if(i < (int)rhs.size()) temp2[i] = Complex(rhs[i], 0);
    }
    //FFTにかけて周波数分解
    temp1 = FFT(2.0*PI/n, temp1);
    temp2 = FFT(2.0*PI/n, temp2);
    //周波数同士で掛け算
    for(int i=0; i<n; i++) temp1[i] *= temp2[i];
    //逆FFTをかけて元に戻す
    temp1 = FFT(-2.0*PI/n, temp1);
    vector<T> ret(n);
    for(int i=0; i<n; i++){
        ret[i] = (T)(temp1[i].real()/n + 0.5); //T=intの時用
    }
    return ret;
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);

    int N; cin >> N;
    vector<int> a(100000), b(100000);
    for(int i=0; i<N; i++) cin >> a[i] >> b[i];
    vector<int> res = Convolution(a, b);
    cout << 0 << endl;
    for(int i=0; i<2*N-1; i++) cout << res[i] << endl;
}