レーダーのアルゴリズムとは、MITのリンカーン研究所のチャールズ・M・レーダーにより考案された高速フーリエ変換のアルゴリズムである 。このアルゴリズムではサイズが素数の離散フーリエ変換（DFT）を巡回畳み込みに置き換えることで計算コストを減らす（ のアルゴリズムもまたDFTを巡回畳み込みと置き換えるアルゴリズムである)。
レーダーのアルゴリズムはDFTカーネルの周期性 のみを利用しているので、このアルゴリズムは同様の特徴をもつ（素数サイズの） やに対しても適用することができる。
このアルゴリズムはデータの並び変えを少し変更し、実数データに対して更なる高速化が可能である 。また、実数データについてに対する別の最適化として、 によって示された離散ハートレイ変換を用いたアルゴリズムがある。
ウィノグラードはレーダーのアルゴリズムを拡張し、素数の冪乗のサイズのDFTに適用できるアルゴリズムを考案した (; ）。このため、レーダーのアルゴリズムは、ウィノグラードのアルゴリズムの特別なケースとみなされる。ウィノグラードのアルゴリズムは乗算的FFTアルゴリズム とも言われる。このアルゴリズムは素数の冪乗以外のサイズにも適用できるが、合成数のサイズのDFTについては、クーリー・ターキーのFFTアルゴリズムを用いるほうがより単純で実装も容易であるため、このアルゴリズムは通常、クーリー・ターキーのアルゴリズムの再帰分割により得られたDFTのうち、大きな素数サイズのDFTにのみ適用される 。
アルゴリズム
サムネイル|RaderのFFT algoritmの視覚的表現。色つきの時計はDFT行列の位相を示す。11の原始根である2が生成する数列(2,4,8,5,10,9,7,3,6,1)を用いて行と列を先頭を除いて入れ替えることで、DFT行列は巡回行列となる。巡回行列をデータ列にかけることは巡回畳み込みと同値である。
N が素数の場合、n =1,...,N–1の添字はNを法とする乗算について群をなす。数論によるとこのような群については生成元（原始根）が存在し、これを gとすると、すべての n は q = 0,...,N–2 を用いて n = gq (mod N) と表される。(qとnの間に全単射が存在する。) n と k の添字をこの置き換えを用いて新しい添字 q と p により表すと、定式は以下のように置き換えられる。
上式は以下に示される長さ N–1(q =0,...,N–2) の数列 aq と bq の巡回畳み込み積分となっている。
畳み込み積分の計算 畳み込み定理 により、FFTを用いて計算することができる。また、N−1 は合成数のため、より高速なクーリーターキー型のFFTアルゴリズムが適用できる。しかしながら、N−1 が大きな素数の素因数をもつ場合、レーダーのアルゴリズムを再帰的に使う必要がある。その代わりに、長さ N−1 のDFTをゼロ埋めによって2の冪乗にすることで、高速なクーリーターキーのアルゴリズムを用いることができる。ゼロ埋め後のデータサイズが最低 2(N−1)−1 以上あれば元のデータを完全に復元することができ、ゼロ埋めによる結果の変化は生じない。
畳み込み積分の計算にゼロ埋めを用いたFFTを用いることでのこのアルゴリズムは加算について O(N) の計算時間とそれに加えと畳み込み積分についてと O(N log N) の計算時間で実行できる。しかしながら、このアルゴリズムは付近の合成数のサイズのFFTに比べ、およそ3から10倍の時間がかかる。
畳み込み積分の計算をゼロ埋めなしのFFTで行う場合、計算時間は N の性質に強く依存する。最悪の場合、N−1 が素数 N2 により N−1= 2N2 と表され、また N2−1 が素数 N3 により N2−1 = 2N3 と表され、以下同様に続いていく場合である。このような場合、レーダーのアルゴリズムの再帰が連続することになり、O(N2) の計算時間がかかる可能性がある。このような性質をもつ N は ソフィー・ジェルマン素数と呼ばれ、上記の数列は一次の Cunningham(ビル-カニンガム)チェーンと呼ばれる。しかしながら、これまでの研究ではカニンガムチェーンの成長は log2(N) よりも遅いことが分かっているため、レーダーのアルゴリズムの再帰によりかかる計算時間は O(N2) よりかは速いと思われる。幸いにも、畳み込み計算にゼロ埋めを用いたFFTを使えば計算時間は O(N log N) のオーダーになることが保証されている。
プログラムの例 以下はC言語による実装例である。
#include <stdio.h> #include <stdlib.h> #include <math.h> int powerMod(int a, int b, int m); int primePrimitiveRoot(int p); void DFT(double **x, int n); void IDFT(double **x, int n);
void FFT_prime(double **x, int n){ int i,len,flag,pad; int g,ig; double X = {0,0}; double **b,**a,**c;
/*ゼロ埋めのサイズの計算*/ flag = 0; len=1; for(i=n-1;i>1;i>>=1){ if(i & 1){flag=1;} len <<= 1;
if(flag){ len <<= 2; /*2(n-1) - 1 以上の2の冪乗*/ }else{ len = n-1; /*n-1が2の冪乗の場合ゼロ埋めはしない*/ } pad = len - (n-1); /*ゼロ埋めの数*/ g = primePrimitiveRoot(n); /*nの原始根*/ ig = powerMod(g,n-2,n); /*nの原始根の逆数*/
/*数列 a_q の作成*/ a = malloc(len*sizeof(double*)); for(i=0;i<len;i++){ a[i] = malloc(2*sizeof(double)); /*1番目と2番目の要素の間をゼロ埋めする。*/ if(i == 0){ a = x; a = x; }else if(i <= pad){ a[i] = 0; a[i] = 0; }else{ a[i] = x[powerMod(g,i-pad,n)]; a[i] = x[powerMod(g,i-pad,n)];
/*数列 b_q の作成*/ b = malloc(len*sizeof(double*)); for(i=0;i<len;i++){ b[i] = malloc(2*sizeof(double)); b[i] = cos(2*M_PI*powerMod(ig,i,n)/n); b[i] = -sin(2*M_PI*powerMod(ig,i,n)/n);
/*離散フーリエ変換*/ DFT(a,len); DFT(b,len);
/*離散フーリエ変換されたaとbの要素ごとの積を計算し、配列cに格納する。*/ c = malloc(len*sizeof(double*)); for(i=0;i<len;i++){ c[i] = malloc(2*sizeof(double)); /*複素数の積の計算*/ c[i] = a[i]*b[i] - a[i]*b[i]; c[i] = a[i]*b[i] + a[i]*b[i];
/*逆離散フーリエ変換*/ IDFT(c,len);
/*x0をXg^-qに足す*/ for(i=0;i<len;i++){ c[i] += x; c[i] += x;
/*X0の計算*/ for(i=0;i<n;i++){ X += x[i]; X += x[i];
/*X0,Xg^-q,を元の配列xに格納する。*/ x = X; x = X;
/*添字の並び変え*/ for(i=0;i<n-1;i++){ x[powerMod(ig,i,n)] = c[i]; x[powerMod(ig,i,n)] = c[i];
for(i=0;i<len;i++){free(a[i]);} free(a); for(i=0;i<len;i++){free(b[i]);} free(b); for(i=0;i<len;i++){free(c[i]);} free(c); }
/*冪剰余を求める関数*/ int powerMod(int a, int b, int m){ int i; for(i=1;b;b>>=1){ if(b & 1){ i = (i*a)%m;
a = (a*a)%m; } return i; }
/*素数の最小の原始根を見つける関数*/ int primePrimitiveRoot(int p){ int i,j,size=0; int *list;
list = malloc((p-1)*sizeof(int)); /*p-1の素因数を見つけ、listに格納する。*/ for(i=2,j=p-1;i*i<=j;i++){ if(j % i == 0){ list[size] = i; size++; do{ j /= i; }while(j%i==0);
for(i=0;i<size;i++){ list[i] = (p-1)/list[i]; } for(i=2;i<=p-1;i++){ for(j=0;j<size;j++){ if(powerMod(i,list[j],p)==1){ goto next;
break; next: ; } free(list); return i; }
/*離散フーリエ変換の関数*/ void DFT(double **x, int n){ int i,j; double **c;
c = malloc(n*sizeof(double*)); for(i=0;i<n;i++){ c[i] = malloc(sizeof(double)); c[i] = 0; c[i] = 0;
for(i=0;i<n;i++){ for(j=0;j<n;j++){ c[i] += x[j]*cos(2*M_PI*i*j/n) + x[j]*sin(2*M_PI*i*j/n); c[i] += -x[j]*sin(2*M_PI*i*j/n) + x[j]*cos(2*M_PI*i*j/n);
for(i=0;i<n;i++){ x[i] = c[i]; x[i] = c[i];
for(i=0;i<n;i++){free(c[i]);} free(c); }
/*逆離散フーリエ変換の関数*/ void IDFT(double **x, int n){ int i,j; double **c;
c = malloc(n*sizeof(double*)); for(i=0;i<n;i++){ c[i] = malloc(sizeof(double)); c[i] = 0; c[i] = 0;
for(i=0;i<n;i++){ for(j=0;j<n;j++){ c[i] += x[j]*cos(2*M_PI*i*j/n) - x[j]*sin(2*M_PI*i*j/n); c[i] += x[j]*sin(2*M_PI*i*j/n) + x[j]*cos(2*M_PI*i*j/n);
for(i=0;i<n;i++){ x[i] = c[i]/n; x[i] = c[i]/n;
for(i=0;i<n;i++){free(c[i]);} free(c); } 複素数データを表現するため、一列目を実数部、二列目を虚数部とした二次元配列を用いている。 レーダーのアルゴリズムではgの逆数が生成する数列も必要であるため、フェルマーの小定理からgの逆数を求めている。 畳み込み積分を不変に保つゼロ埋めの手法として、aについて先頭要素とその次の要素の間にゼロを埋め、bについてはgの生成する数列がn-1の周期をもつことを自然に用いて元の数列を巡回させている。ゼロ埋めの手法はこの方法以外にも考えられる。 畳込み積分の計算に通常の離散フーリエ変換を用いているが、これを高速フーリエ変換のルーチンで置き換えることができる。 原始根を求めるアルゴリズムは原始根を参照。
参考文献