高速フーリエ変換

クラスで実装されているが、実行するときはconvolutionを呼び出せばいい

多項式に見立てたvectorを\(2\)つ与えると、それらに対し高速フーリエ変換を用いた畳み込み演算を行い一つの多項式にし、vectorとして返す

\(\{1,2,3\}*\{1,3,5\}=\{1,5,14,19,15\}\)といった調子

実装上doubleを経由しているので誤差死が怖い

計算量は\(N=max(\vert A \vert,\vert B\vert)\)として、\(Ο(NlogN)\)

github

class fast_Fourier_transform{
    // Copyright (c) 2023 0214sh7
    // https://github.com/0214sh7/library/
    private:
    const double PI=3.14159265358979; 
    
    std::vector<std::complex<double>> DFT(std::vector<std::complex<double>> A){
        int N=A.size();
        if(N==1)return A;
        std::vector<std::complex<double>> A0(N/2),A1(N/2),iA0,iA1,iA(N);
        for(int i=0;i<N;i++){
            if(i%2==0){
                A0[i/2]=A[i];
            }else{
                A1[i/2]=A[i];
            }
        }
        iA0=DFT(A0);
        iA1=DFT(A1);
        
        for(int i=0;i<N;i++){
            std::complex<double> ith_zeta = std::complex<double>(cos(2*PI*i/N),sin(2*PI*i/N));
            iA[i]=(iA0[i%(N/2)]+ ith_zeta*iA1[i%(N/2)]);
        }
        return iA;
    }
    
    std::vector<std::complex<double>> iDFT(std::vector<std::complex<double>> iA){
        int N=iA.size();
        std::vector<std::complex<double>> A,dA,rA;
        dA=DFT(iA);
        for(int i=0;i<N;i++){
            rA.push_back(dA[(N-i)%N]);
            A.push_back(rA[i]/std::complex<double>(N,0));
        }
        return A;
    }
    
    public:
    std::vector<long long> convolution(std::vector<long long> A,std::vector<long long> B){
        int deg = A.size() + B.size() -1;
        long long N=1;
        while(N<deg){N<<=1;}
        A.resize(N);B.resize(N);
        std::vector<std::complex<double>> dC(N),iC(N),dA,iA,dB,iB;
        std::vector<long long> C(N);
        for(int i=0;i<A.size();i++){
            dA.push_back(A[i]);
        }
        for(int i=0;i<B.size();i++){
            dB.push_back(B[i]);
        }
        iA=DFT(dA);iB=DFT(dB);
        for(int i=0;i<N;i++){
            iC[i]=iA[i]*iB[i];
        }
        dC=iDFT(iC);
        for(int i=0;i<dC.size();i++){
            C[i]=(0.1+dC[i].real());
        }
        C.resize(deg);
        return C;
    }
    
};

使用例

---

実行コード

#include <bits/stdc++.h>

class fast_Fourier_transform{/*省略*/};

int main(void){
    
    fast_Fourier_transform FFT;
    std::vector<long long> A={1,2,3},B={1,3,5},C;
    
    C = FFT.convolution(A,B);
    for(int i=0;i<C.size();i++){
        std::cout << C[i] << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

出力

1 5 14 19 15 

更新履歴

---

日時	内容
2023/06/29	ライセンスのコメントアウトを変更
2021/03/25	バグを修正
2021/03/25	使用例を追加
2020/04/04	高速フーリエ変換を追加