円周率.jp - 実装方式

実装方式

FFT を素直に実装すると計算量の導出に書いたような処理になる。ただしこれらのコードは、再帰呼び出しになっているため、並べ替えのために単純にコピーを行ったり、同じ $\omega$ の計算を繰り返し行ったりしていて効率がよくない。そのため一般的には再帰呼び出しを取り払い、そのような無駄を取り除いた方式の実装を FFT の実装として考えることが多い。ここでは広く知られている 2 種類の実装方式を紹介する。

プログラムを実装する際にどちらの方が良いかは細かな実装やユーザの需要次第となるので決定的な言及はできないが，大まかには小メモリ動作を目指すなら Cooley-Tukey 型を，速さを求めるなら Stockham 型を使用する印象である。

なお、FFT の高速演算を行うために様々な方法が提案されているが、それらは全て実装方法がどちらの方式であっても適用することができる。 in-place 型 FFT での効率化実装については[JW06]が詳しいので，本サイトでは Stockham 型 FFT を用いて説明している。

Stockham 型 FFT / self-sorting 型 FFT

冒頭部分で述べたような、再帰呼び出しや同じ $\omega$ の計算を素直に省いた実装を self-sorting 型 FFT、もしくは開発者の名をとって Stockham 型 FFT という。 self-sorting の名の通りデータが常に並んでいるため計算がしやすい。

時間間引き FFT

void InTimeFFTStCore(int n, int m, Complex* x, Complex* y) {
  const int l = n / m / 2;
  for (int k = 0; k < m; ++k) {
    double wt = -M_PI * k / m;
    Complex w(cos(wt), sin(wt));
    for (int j = 0; j < l; ++j) {
      int xi0 = j * m + k, xi1 = xi0 + l * m;
      int yi0 = j * 2 * m + k, yi1 = yi0 + m;
      Complex t = x[xi1] * w;
      y[yi1] = x[xi0] - t;
      y[yi0] = x[xi0] + t;
    }
  }
}

void InTimeFFTSt(int n, Complex* x, Complex* y) {
  for (int m = 1, i = 0; m < n; m *= 2, ++i) {
    if (i % 2) {
      InTimeFFTStCore(n, m, y, x);
    } else {
      InTimeFFTStCore(n, m, x, (m == n / 2) ? x : y);
    }
  }
}

周波数間引き FFT

void InFreqFFTStCore(int n, int m, Complex* x, Complex* y) {
  const int l = n / m / 2;
  for (int k = 0; k < m; ++k) {
    double wt = -M_PI * k / m;
    Complex w(cos(wt), sin(wt));
    for (int j = 0; j < l; ++j) {
      int xi0 = j * 2 * m + k, xi1 = xi0 + m;
      int yi0 = j * m + k, yi1 = yi0 + l * m;
      Complex t0 = x[xi0] + x[xi1], t1 = (x[xi0] - x[xi1]) * w;
      y[yi0] = t0;
      y[yi1] = t1;
    }
  }
}

void InFreqFFTSt(int n, Complex* x, Complex* y) {
  int m = n / 2;
  if (static_cast<int>(log2(n)) % 2) {
    InFreqFFTStCore(n, m, x, x);
    m /= 2;
  }
  for (int i = 0; m >= 1; m /= 2, ++i) {
    if (i % 2)
      InFreqFFTStCore(n, m, y, x);
    else
      InFreqFFTStCore(n, m, x, y);
  }
}

Cooley-Tukey 型 FFT / in-place 型 FFT

FFT は各ステップでの計算の際のデータ依存を考えるとそもそも元のデータ領域だけで計算を完結させることができることが分かる。そのようにしてメモリ領域の無駄を省いた実装を in-place 型 FFT、もしくは Cooley-tukey 型 FFT という。ただしこの方式はメモリの無駄を無くした一方で、データ並べ替えを別作業として行う必要が有るため、厳密には CPU 使用量が少し多くなるはずである。

void BitReverse(int n, Complex* x) {
  for (int i = 1, j = 0; i < n; ++i) {
    for (int b = n >> 1; b; b >>= 1) {
      j ^= b;
      if (j & b) break;
    }
    if (i < j)
      swap(x[i], x[j]);
  }
}

時間間引き FFT

void InTimeFFTCT(int n, Complex* x) {
  BitReverse(n, x);
  for (int m = 1; m < n; m *= 2) {
    for (int k = 0; k < m; ++k) {
      double wt = -M_PI * k / m;
      Complex w(cos(wt), sin(wt));
      for (int j = 0; j < n; j += 2 * m) {
        Complex &x0 = x[k + j], &x1 = x[k + j + m];
        Complex x1w = x1 * w;
        x1 = x0 - x1w;
        x0 = x0 + x1w;
      }
    }
  }
}

周波数間引き FFT

void InFreqFFTCT(int n, Complex* x) {
  for (int m = n / 2; m >= 1; m /= 2) {
    for (int k = 0; k < m; ++k) {
      double wt = -M_PI * k / m;
      Complex w(cos(wt), sin(wt));
      for (int j = 0; j < n; j += 2 * m) {
        Complex &x0 = x[k + j], &x1 = x[k + j + m];
        Complex x1w = (x0 - x1) * w;
        x0 = x0 + x1;
        x1 = x1w;
      }
    }
  }
  BitReverse(n, x);
}