高速フーリエ変換(FFT)


高速フーリエ変換は新しい変換ではなく、フーリエ変換の高速アルゴリズムです.このアルゴリズムは通常の離散Fourier変換(DFT)の時間的複雑度O(n*n)をOに下げることができる.(n log n)は、フーリエ変換の速度を大幅に向上させた.高速フーリエ変換アルゴリズムの提案により、フーリエ変換は通信分野において極めて大きな運用と発展を遂げた.ACMでは、高速フーリエ変換は通常、大数の乗算に用いられる.2つの数がJAVAのBIにも耐えられないほど大きい場合には、高速フーリエアルゴリズムを用いる.このアルゴリズムは、2つの数を、いずれも多項式行列として書き,次いで離散フーリエ変換を行った.通信分野には、時間領域のボリューム、周波数領域の乗算という古典的な言葉がある.従って,2つの数の多項式係数行列を離散フーリエ変換した後,直接乗算することができる.乗算した結果をフーリエ逆変換(DTFT)する.
########私はこれまでどうしても高速フーリエアルゴリズムを理解できませんでした.当時は複雑すぎて難しいと思っていたので、見ても読めませんでした.しかし、私の専攻は電子情報工学なので、各専攻の授業でフーリエ変換と会うのは避けられない.先週、私たちのデジタル信号処理の先生は、いくつかの授業の時間を費やして、高速フーリエアルゴリズムを説明してくれました.私はまじめに聞いて、自分で理解したと思っていました.しかし、自分の頭の中の高速フーリエアルゴリズムで大数乗算の問題を解きたいと思ったとき、手がつけられないことに気づいた.どうやって運用するか全然分かりません.その时、私はこのアルゴリズムをマスターしなければならないと思っていました.ACMに現れるからだけでなく、私の専攻でも極めて重要な地位を占めているからです.もし私がフーリエ変換さえできなかったら、これから自分が電子情報専門だと言う顔はありません.そこで、私は2日かけて、資料を調べて、やっとこの強大なアルゴリズムを理解しました.私はいくつかの大神が書いたFFTに関する文章が確かにとても良くて、はっきりしていて分かりやすくて、私自身もコレクションして、みんなに推薦して、みんなで一緒に勉強します.
超分かりやすいフーリエ変換
高速フーリエ変換(FFT)
###JAVAコード
import java.util.Arrays;
import java.util.Scanner;
import java.util.Vector;

/**
 * Created by jal on 2017/12/6 0006.
 */
public class DIT_FFT {
    private static int maxn;
    public static void main(String[] args) {
        Scanner scanner = new Scanner(System.in);
        String numstr1 = scanner.next();
        String numstr2 = scanner.next();
        int len =numstr1.length()+numstr2.length();
        maxn = 0;
        double temp = log2(len);
        double floortemp = Math.floor(temp);
        if(floortemp == temp){
            maxn = len;
        }else {
            maxn = (int)Math.pow(2,floortemp+1);
        }
        Complex []arra = createArray(maxn);
        Complex []arrb =createArray(maxn);
        Complex []arrc =createArray(maxn);
        Complex []arrA = createArray(maxn);
        Complex []arrB = createArray(maxn);
        Complex []arrC = createArray(maxn);

        char []a = numstr1.toCharArray();
        int j = 0;
        for(int i = a.length - 1; i >= 0; i--){
            arra[j++] = new Complex(a[i] - '0');
        }
        j = 0;
        char []b = numstr2.toCharArray();
        for(int i = b.length - 1; i >= 0; i--){
            arrb[j++] = new Complex(b[i] - '0');
        }
        //System.out.println(Arrays.toString(arra));

        //invert(arra);



        arrA = fft(arra);
        System.out.println(Arrays.toString(arrA));
        //invert(arrb);

        arrB = fft(arrb);
        //System.out.println(Arrays.toString(arrB));
        for(int i = 0; i < arrC.length; i++){
            arrC[i] = arrA[i].times(arrB[i]);
        }
        //System.out.println(Arrays.toString(arrC));
        arrc = ifft(arrC);
        Vector vector = new Vector<>();
        vector = toIntOfString(arrc);
        String str = "";
        char vectorCHar[] = new char[vector.size()];
        for(int i = 0; i < vectorCHar.length; i++){
            vectorCHar[i] = (char)(vector.get(i) + '0');
        }
        str =String.valueOf(vectorCHar);
        str = new StringBuffer(str).reverse().toString();
        //System.out.println("str:"+str);
        str = trim(str);
        System.out.println(str);

    }
    private static String trim(String value) {

        int len = value.length();
        int st = 0;
        char[] val = value.toCharArray();    /* avoid getfield opcode */

        while ((st < len) && (val[st] <= '0')) {
            st++;
        }
//            while ((st < len) && (val[len - 1] <= ' ')) {
//                len--;
//            }
        return value.substring(st, len);

    }
    private static Vector toIntOfString(Complex[] arrc) {
        Vector result = new Vector<>();
        int n = arrc.length;
        int arrayTemp [] = new int[n+1];
        for(int i = 0; i < arrc.length; i++){
            arrayTemp[i] = (int)arrc[i].re();
        }
        int i;
        arrayTemp[n] = 0;
        for(i = 0; i < n; i++){
            result.add(arrayTemp[i]%10);
            arrayTemp[i+1] += arrayTemp[i] / 10;
        }
        int t = arrayTemp[n];
        while(t > 0){
            result.add(t % 10);
            t/= 10;
        }
        return result;
    }

    private static Complex[] ifft(Complex[] arrC) {
        int n = arrC.length;
        Complex arrayResult[] = new Complex[n];
        for(int i = 0; i < arrC.length; i++){
            arrC[i] = arrC[i].conjugate();
        }
        arrayResult= fft(arrC);
        for(int i = 0; i < arrayResult.length; i++){
            arrayResult[i] = arrayResult[i].conjugate().divides(new Complex(n));
        }
        return arrayResult;
    }

    private static Complex[] fft(Complex[] arrA) {
        int N = arrA.length;
        if(N == 1){
            return arrA;
        }
        Complex [] arrayEven = new Complex[N/2];
        Complex [] arrayOdd = new Complex[N/2];
        for(int i = 0; i < arrayEven.length; i++){
            arrayEven[i] = arrA[2*i];
            arrayOdd[i] = arrA[2*i+1];
        }
        arrayEven = fft(arrayEven);
        arrayOdd = fft(arrayOdd);
        Complex[]arrayResult = new Complex[N];
        for(int i = 0; i < N/2; i++){
            Complex W = new Complex(0,-2*Math.PI*i/N).exp();
            arrayResult[i] = arrayEven[i].plus(arrayOdd[i].times(W));
            arrayResult[i+N/2] = arrayEven[i].minus(arrayOdd[i].times(W));
        }
        return arrayResult;
    }

    private static  void bit_reverse_swap(Complex [] a) {
        int n = a.length;
        for (int i = 1, j = n >> 1, k; i < n - 1; ++i) {
            if (i < j) swap(a,i,j);
            // tricky
            for (k = n >> 1; j >= k; j -= k, k >>= 1)  // inspect the highest "1"
                ;
            j += k;
        }
    }

    private static void invert(Complex[] arra) {
        int n = (int)log2(maxn);
        for(int i = 0; i < arra.length; i++){
            String temp = Integer.toBinaryString(i);
            temp = new StringBuffer(temp).reverse().toString();
            int len = n - temp.length();
            char [] zeros = new char[len];
            Arrays.fill(zeros,'0');
            temp+=String.valueOf(zeros);

            int j = Integer.valueOf(temp,2);
            System.out.println("i:" + i +"j:" + j);
            if(j>i){
                swap(arra,i,j);
            }
        }
    }

    private static void swap(Complex[] arra, int i, int j) {
        Complex tmp = arra[i];
        arra[i] = arra[j];
        arra[j] = tmp;
    }
    private static Complex[] createArray(int maxn) {
        Complex[] result =new Complex[maxn];
        for(int i=0;i
 *  For additional documentation, see Section 9.9 of
 *  Algorithms, 4th Edition by Robert Sedgewick and Kevin Wayne.
 *
 *  @author Robert Sedgewick
 *  @author Kevin Wayne
 */
class Complex {
    private double re;   // the real part
    private  double im;   // the imaginary part

    /**
     * Initializes a complex number from the specified real and imaginary parts.
     *
     * @param real the real part
     * @param imag the imaginary part
     */
    public Complex(double real, double imag) {
        re = real;
        im = imag;
    }

    public Complex(double re) {
        this.re = re;
        this.im = 0;
    }

    /**
     * Returns a string representation of this complex number.
     *
     * @return a string representation of this complex number,
     *         of the form 34 - 56i.
     */
    public String toString() {
        if (im == 0) return re + "";
        if (re == 0) return im + "i";
        if (im <  0) return re + " - " + (-im) + "i";
        return re + " + " + im + "i";
    }

    /**
     * Returns the absolute value of this complex number.
     * This quantity is also known as the modulus or magnitude.
     *
     * @return the absolute value of this complex number
     */
    public double abs() {
        return Math.hypot(re, im);
    }

    /**
     * Returns the phase of this complex number.
     * This quantity is also known as the angle or argument.
     *
     * @return the phase of this complex number, a real number between -pi and pi
     */
    public double phase() {
        return Math.atan2(im, re);
    }

    /**
     * Returns the sum of this complex number and the specified complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this + that)}
     */
    public Complex plus(Complex that) {
        double real = this.re + that.re;
        double imag = this.im + that.im;
        return new Complex(real, imag);
    }

    /**
     * Returns the result of subtracting the specified complex number from
     * this complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this - that)}
     */
    public Complex minus(Complex that) {
        double real = this.re - that.re;
        double imag = this.im - that.im;
        return new Complex(real, imag);
    }

    /**
     * Returns the product of this complex number and the specified complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this * that)}
     */
    public Complex times(Complex that) {
        double real = this.re * that.re - this.im * that.im;
        double imag = this.re * that.im + this.im * that.re;
        return new Complex(real, imag);
    }

    /**
     * Returns the product of this complex number and the specified scalar.
     *
     * @param  alpha the scalar
     * @return the complex number whose value is {@code (alpha * this)}
     */
    public Complex scale(double alpha) {
        return new Complex(alpha * re, alpha * im);
    }

    /**
     * Returns the product of this complex number and the specified scalar.
     *
     * @param  alpha the scalar
     * @return the complex number whose value is {@code (alpha * this)}
     * @deprecated Replaced by {@link #scale(double)}.
     */
    @Deprecated
    public Complex times(double alpha) {
        return new Complex(alpha * re, alpha * im);
    }

    /**
     * Returns the complex conjugate of this complex number.
     *
     * @return the complex conjugate of this complex number
     */
    public Complex conjugate() {
        return new Complex(re, -im);
    }

    /**
     * Returns the reciprocal of this complex number.
     *
     * @return the complex number whose value is {@code (1 / this)}
     */
    public Complex reciprocal() {
        double scale = re*re + im*im;
        return new Complex(re / scale, -im / scale);
    }

    /**
     * Returns the real part of this complex number.
     *
     * @return the real part of this complex number
     */
    public double re() {
        return re;
    }

    /**
     * Returns the imaginary part of this complex number.
     *
     * @return the imaginary part of this complex number
     */
    public double im() {
        return im;
    }

    /**
     * Returns the result of dividing the specified complex number into
     * this complex number.
     *
     * @param  that the other complex number
     * @return the complex number whose value is {@code (this / that)}
     */
    public Complex divides(Complex that) {
        return this.times(that.reciprocal());
    }

    /**
     * Returns the complex exponential of this complex number.
     *
     * @return the complex exponential of this complex number
     */
    public Complex exp() {
        return new Complex(Math.exp(re) * Math.cos(im), Math.exp(re) * Math.sin(im));
    }

    /**
     * Returns the complex sine of this complex number.
     *
     * @return the complex sine of this complex number
     */
    public Complex sin() {
        return new Complex(Math.sin(re) * Math.cosh(im), Math.cos(re) * Math.sinh(im));
    }

    /**
     * Returns the complex cosine of this complex number.
     *
     * @return the complex cosine of this complex number
     */
    public Complex cos() {
        return new Complex(Math.cos(re) * Math.cosh(im), -Math.sin(re) * Math.sinh(im));
    }

    /**
     * Returns the complex tangent of this complex number.
     *
     * @return the complex tangent of this complex number
     */
    public Complex tan() {
        return sin().divides(cos());
    }



}

###C++コード
#include 
#include 
#include 
#include 
#include 
#include 
#include 
using namespace std;
double log2(int n){
    return log(n)/log(2);
}
const int MAXN = 1 << 20;
complex arra[MAXN],arrb[MAXN],arrC[MAXN];
complex* arrA;
complex* arrB;//,arrC[MAXN];
complex* arrc;
const double PI = 3.141592653;
complex* DIT_FFT(complex* arr,int len){
    if(len == 1)return arr;
    complex *arrayEven = new complex[len/2];
    complex *arrayOdd = new complex[len/2];
    for(int i = 0; i < len/2; i++){
        arrayEven[i] = arr[2*i];
        arrayOdd[i] = arr[2*i+1];
    }
    arrayEven = DIT_FFT(arrayEven,len/2);
    arrayOdd = DIT_FFT(arrayOdd,len/2);
    complex *arrayResult = new complex[len];
    for(int i = 0; i < len/2; i++){
        //TODO
        //help me
        //help me
        // what is 'W' ?
        complex W = exp(complex(0,-2*PI *i / len));//ecomplex(cos(-2*PI *i / len),sin(-2*PI *i / len));
        cout<* IFFT(complex*arrC,int len){
    int n = len;
    complex* arrayResult = new complex[len];
    for(int i = 0; i < len; i++){
        arrC[i] = conj(arrC[i]);
    }
    arrayResult= DIT_FFT(arrC,len);
    for(int i = 0; i < len; i++){
        arrayResult[i] = conj(arrayResult[i]);
        arrayResult[i]/=n;
    }
    return arrayResult;

}
int main(){

    string str1,str2;
    cin>>str1>>str2;
    int len = str1.size() + str2.size();
    int maxn = 0;
    double floortemp =log2(len);
    if(floortemp == floor(floortemp)){
        maxn = len;
    }
    else{
        maxn = pow(2,(int)(floor(floortemp) + 1) );
    }

    for(int i = str1.size() - 1; i >= 0; i--){
        arra[i] = complex(str1[i]-'0',0);
    }
    for(int i = str1.size() - 1; i >= 0; i--){
        arrb[i] = complex(str2[i]-'0',0);
    }
    arrA = DIT_FFT(arra,maxn);
    arrB = DIT_FFT(arrb,maxn);
    for(int i= 0; i < maxn; i++){
        arrC[i] = arrA[i]*arrB[i];
    }
    arrc = IFFT(arrC,len);
    vectorv;
    for(int i = 0; i < maxn;i++){
        v.push_back((int)arrc[i].real());
    }

}