[c++] FFT 実装 AtCoder Tyipical Contest C 高速フーリエ変換 

##問題
https://atcoder.jp/contests/atc001/tasks/fft_c
##FFTについて
多項式$f(x)=c_0 +c_1 x+..+c_{n-1} x^{n-1}$の係数はわからないけど、$n+1$個の点$f(x_0)=f_0 ,..,f(x_n)=f_n$がわかっている状況を想定していて、
この時、連立方程式

```math
\begin{cases}
    f(x_0)=c_0 +x_0 c_1 +..+(x_0 )^{n-1}c_{n-1} \\
    .\\
    .\\
    .\\
    f(x_n)=c_0 +x_n c_1 +..+(x_n )^{n-1}c_{n-1} \\
\end{cases}
```
が解ければ、${c_i}$がわかるので、元の$f(x)$の形がわかる。この作業を補完っていうらしい。この操作の総称をDFT離散フーリエ変換って言うみたい。
この時、$x_i$と$n$をうまく選べば、すごく計算すべき量を減らせるよっていう手法がFFTで、結果的には$n$を自然数$r$のべき乗$n=r^m$にして、$x_i$を1のn乗根$\zeta_n$の$i$乗にすればうまく計算できるみたい。
$\zeta_n$は基本的に複素数になることに注意。
詳しいことは、[ここ](https://qiita.com/ageprocpp/items/0d63d4ed80de4a35fe79)を見てね。
##実装例
自前の複素数libraryを使うと早くなるよ！！と書いてあったので、@ageprocpp さんの助言も併せて実装すると、標準ライブラリを使った時より、最高200ms早くなりました！！

```c++
struct FFT{
    struct mycomplex{   
        double real;
        double imag;
        mycomplex(double a=0.0,double b=0.0):real(a),imag(b){}
        mycomplex operator-() const {return mycomplex(-real,-imag);}
        mycomplex &operator+=(const mycomplex &a){real+=a.real;imag+=a.imag;return *this;}
        mycomplex &operator-=(const mycomplex &a){real-=a.real;imag-=a.imag;return *this;}
        mycomplex &operator*=(const mycomplex &a){double nr=real*a.real-imag*a.imag;double ni=real*a.imag+imag*a.real;real=nr;imag=ni;return *this;}
        mycomplex &operator/=(const mycomplex &a){double nr=(real*a.real+imag*a.imag)/(a.real*a.real+a.imag*a.imag);double ni=(imag*a.real-real*a.imag)/(a.real*a.real+a.imag*a.imag);real=nr;imag=ni;return *this;}
        mycomplex &operator*=(const double &a){real*=a;imag*=a;return *this;}
        mycomplex &operator/=(const double &a){real/=a;imag/=a;return *this;}
        mycomplex operator+(const mycomplex &a) const {return mycomplex(*this)+=a;}
        mycomplex operator-(const mycomplex &a) const {return mycomplex(*this)-=a;}
        mycomplex operator*(const mycomplex &a) const {return mycomplex(*this)*=a;}
        mycomplex operator/(const mycomplex &a) const {return mycomplex(*this)/=a;}
        mycomplex operator*(const double &a) const {return mycomplex(*this)*=a;}
        mycomplex operator/(const double &a) const {return mycomplex(*this)/=a;}
        mycomplex conj() const {return mycomplex(real,-imag);}
        mycomplex polar(double r,double q) const {return mycomplex(r*cos(q),r*sin(q));}
    };
    typedef mycomplex clx;
    const double eps=1e-3;

    vector<clx> dft(vector<clx> f,int n,bool inverse=false){
        if(n==1)return f;
        vector<clx> a0,a1;
        rep(i,n)i%2==0?a0.push_back(f[i]):a1.push_back(f[i]);
        a0=dft(a0,n/2,inverse);
        a1=dft(a1,n/2,inverse);
        clx zeta(cos(2*pi/n),sin(2*pi/n));
        if(inverse)zeta=zeta.conj();
        clx x=clx(1);
        rep(i,n){
            f[i]=a0[i%(n/2)]+x*a1[i%(n/2)];
            x*=zeta;
        }
        return f;
    }
    vector<ll> operator()(vector<double> a,vector<double> b){
        int ma=a.size(),mb=b.size();
        ll N=1;
        for(N;N<=ma+mb+1;N*=2);
        //0を詰める
        vector<clx> aa(N,clx()),bb(N,clx());
        for(int i=0;i<ma;i++)aa[i]=clx(a[i]);
        for(int i=0;i<mb;i++)bb[i]=clx(b[i]);

        aa=dft(aa,N);
        bb=dft(bb,N);

        vector<clx> res;
        rep(i,N+1)res.push_back(aa[i]*bb[i]);

        res=dft(res,N,true);

        vector<ll> ans;//整数値に直す（本質的ではない）
        rep(i,N+1)ans.push_back(eps+res[i].real/N);//Nで割る
        return ans;
    }
};
```