【配列へのアクセスを減らす最適化】

　セクション２ で完成版 Stockham のアルゴリズムを示しましたが、ちょっとした工夫で、 配列へのアクセスを減らすことができます。 まずは、リスト７：完成版 Stockham のアルゴリズムを再掲しておきましょう。

#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
// n  : 系列長
// s  : ストライド
// eo : eo == 0 か false なら x が出力、eo == 1 か true なら y が出力
// x  : フーリエ変換する入力系列(eo == 0 のとき出力)
// y  : 作業用配列(eo == 1 のとき出力)
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) { if (eo) for (int q = 0; q < s; q++) y[q] = x[q]; }
    else {
        for (int p = 0; p < m; p++) {
            const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft0(n/2, 2*s, !eo, y, x);
    }
}

void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
    complex_t* y = new complex_t[N];
    fft0(N, 1, 0, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] /= N;
}

void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
    for (int p = 0; p < N; p++) x[p] = conj(x[p]);
    complex_t* y = new complex_t[N];
    fft0(N, 1, 0, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}

　よく見ると fft0 は n == 1 の時、 何も計算せずに配列のコピーを行っているだけです。 つまり、このステップの前の計算でコピー先の配列へ計算結果を代入してやれば、 このコピーは消去できます。そのように書き換えると以下のようになります。

#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
// n  : 系列長
// s  : ストライド
// eo : eo == 0 か false なら x が出力、eo == 1 か true なら y が出力
// x  : フーリエ変換する入力系列(eo == 0 のとき出力)
// y  : 作業用配列(eo == 1 のとき出力)
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 2) {
        complex_t* z = eo ? y : x;
        for (int q = 0; q < s; q++) {
            const complex_t a = x[q + 0];
            const complex_t b = x[q + s];
            z[q + 0] = a + b;
            z[q + s] = a - b;
        }
    }
    else if (n >= 4) {
        for (int p = 0; p < m; p++) {
            const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft0(n/2, 2*s, !eo, y, x);
    }
}

void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
    complex_t* y = new complex_t[N];
    fft0(N, 1, 0, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] /= N;
}

void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
    for (int p = 0; p < N; p++) x[p] = conj(x[p]);
    complex_t* y = new complex_t[N];
    fft0(N, 1, 0, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}

　副作用として、n == 2 の時のバタフライ演算が簡単になるため、 1 を掛けるだけの無駄な乗算も消去できています。