#include "TrueBASIC.h"
int main();
void FFT(double g[][2], int p);

void FFT(double g[][2], int p)
{
    int L, N, N2, del_i, i, ip, j, jmax, k, Itmp1_;
    double Wp[2], angle, factor[2], temp[2];

    // 複素数の入力データの組 g の高速フーリエ変換
    // 変換値は g に戻される
    N = (1 << p);                           // データ点の数
    N2 = N/2;
    // ビット順反転の方法による入力データの並べ換え
    j = 0;
    for(i = 0; i < N-1; ++i)
    {
        // g(i,1) は f((i-1)*del_t) の実数部分
        // g(i,2) は f((i-1)*del_t) の虚数部分
        // データが実数ならば，g(i,2) = 0
        if(i < j)
        {
            temp[0] = g[j][0];
            temp[1] = g[j][1];
            g[j][0] = g[i][0];
            g[j][1] = g[i][1];
            g[i][0] = temp[0];
            g[i][1] = temp[1];
        }
        k = N2;
        while(k <= j)
        {
            j -= k;
            k /= 2;
        }
        j += k;
    }
    // ダニエルソン-ランチョスのアルゴリズムの開始
    jmax = 1;
    for(L = 1; L <= p; ++L)
    {
        del_i = 2*jmax;
        Wp[0] = 1;                          // Wp を W^0 で初期化
        Wp[1] = 0;
        angle = pi/jmax;
        factor[0] = cos(angle);             // W^p の新しい値と古い値の比    
        factor[1] = -sin(angle);
        for(j = 0; j < jmax; ++j)
        {
            for(i = j; i < N; i += del_i)
            {
                // 長さ 2^L の変換を計算
                ip = i + jmax;
                temp[0] = g[ip][0]*Wp[0] - g[ip][1]*Wp[1];
                temp[1] = g[ip][0]*Wp[1] + g[ip][1]*Wp[0];
                g[ip][0] = g[i][0] - temp[0];
                g[ip][1] = g[i][1] - temp[1];
                g[i][0] = g[i][0] + temp[0];
                g[i][1] = g[i][1] + temp[1];
            }
            // 新しい W^p を計算
            temp[0] = Wp[0]*factor[0] - Wp[1]*factor[1];
            temp[1] = Wp[0]*factor[1] + Wp[1]*factor[0];
            for(Itmp1_ = 0; Itmp1_ < 2; ++Itmp1_)
                Wp[Itmp1_] = temp[Itmp1_];
        }
        jmax = del_i;
    }
}