目的・方法

よく紹介されているFFT(高速フーリエ変換)のアルゴリズムでは、要素数が2の冪乗のCooley-Tukey法が使われている。これを使って要素数が2の冪乗でないデータをフーリエ変換しようと思うと、要素数より大きい最小の2の冪乗を求め、配列を生成し、データを移し、ゼロ埋めしてFFTしなければならない。また、ゼロ埋による副作用も考えなくてはならない。任意の要素数の配列をそのまま扱えたほうが、FFTを使う側としては楽である。

Cooley-Tukey法は、要素数が合成数ならば使える。サンプル数が素数のときは、レーダーのアルゴリズムがある。これを組み合わせて、任意要素数のFFTのプログラムを書く。

DFT(離散フーリエ変換)とは

サンプル数nのデータ列 $x[i]$ から以下の演算をして、$y[i]$を得る。

y_i = \sum_{j=0}^{n-1}x_j \omega_n^{-i j}

ここで $\omega_n = e^{2 \pi i/n}$ である。
この計算は $O(n^2)$ だが、FFTのアルゴリズムを用いると同じ計算を $O(n\log n)$ の計算で行うことができる。

Cooley-Tukey法とは

原理

$n=n_1 n_2$ のとき、$i=i_1+i_2 n_1$ , $j=j_1 n_2 + j_2$ と置くことで以下の計算で離散フーリエ変換ができる。
内側の丸括弧の中は要素数$n_1$の離散フーリエ変換になっているので、再帰的にこの計算を用いることができる。

\begin{align}
y_{i_1+i_2 n_1} &= \sum_{j_2=0}^{n_2-1} \sum_{j_1=0}^{n_1-1} x_{j_1 n_2+j_2} \omega_{n}^{-(i_1+i _2 n_1) (j_1 n_2 + j_2)} \\
&= \sum_{j_2=0}^{n_2-1} \left( \sum_{j_1=0}^{n_1-1} x_{j_1 n_2+j_2} \omega_{n}^{-i_1 j_1 n_2}\right) \omega_n^{-(i_1+i _2 n_1)j_2} \\
&= \sum_{j_2=0}^{n_2-1} \left( \sum_{j_1=0}^{n_1-1}x_{j_1 n_2+j_2} \omega_{n_1}^{-i_1 j_1}\right) \omega_n^{-(i_1+i _2 n_1) j_2}
\end{align}

実はこの式には「添字の並び変え」の操作が暗黙に含まれている。一番内側にある$x_{j_1 n_2 + j_2}$ は$j_1$を1ずつ増やすと、$n_2$ずつ飛び飛びに増えていく。なので、添字の並び変えをしなければ、離散フーリエ変換がそのまま使えない。

なぜ速くなるのか

計算の数が減るからである。
実は、DFTの式の $x_j \omega_n^{ij}$ は異なるiについて同じ値になることがあるため、そのままだと同じ計算を複数回行うことになってしまうのである。Cooley-Tukey法はこの重複を無くすことで計算を減らしている。

図解

$\omega_n^{ij}$は回転因子と呼ばれ、回転するもので表せる。下図は要素数6のDFTを回転する時計で表したものである。この行列のi行j列目にある時計が$\omega_n^{ij}$である。左図は元の行列、右図は添字の並び変えを行った後の行列である。

$x[j] \omega_n^{ij}$は同じ列に同じ時計があるときに同じ値になる。左図を見ると、いくつか重複している。
Cooley-Tukey法による添字の並び変えを行うと、この重複がうまく抽出されて、同じパターンが複数現れる。

右図の左側は要素数3のフーリエ変換の行列と全く同じである。右側はそうなってはいないが、実はこれも、要素数3のフーリエ変換の行列を行ごとに「ひねっていく」ことで作れるのである。これが式の最後の$\omega_n^{-(i_1+i_2 n_1) j_2}$である。

対して素数の要素数の場合はどうだろうか。下は要素数5の行列である。

周波数ゼロを除き、全ての角度が全ての列に重複なく分配されている。まったくスキがない。

Raderのアルゴリズムとは

原理

nが素数のとき、nには原始根が存在し、nを法とする数は全てこの原始根の冪乗で表せる。例えばn=5の場合、原始根は2である。2を5を法として累乗していくと、{2,4,3,1}となり、1~4までの数字全てが現れる。n=7の場合、原始根は3で{3,2,6,4,5,1}となり1~6全ての数字が現れる。

DFTの式再掲

y_i = \sum_{j=0}^{n-1}x_j\omega_n^{-i j}

nの原始根をgとし、添字を原始根を使って$i=g^q$ , $j=g^{-p}$と表すと

\begin{align}
y_0 &= \sum_{j=0}^{n-1}x_j \\
y_{g^p} &= x_0 + \sum_{q=0}^{n-2}x_{g^{-p}}\omega_n^{g^{q-p}}
\end{align}

これは実は長さn-1の巡回畳み込みという計算になっている。そして巡回畳み込みは畳み込み定理によりFFTで計算できるため、高速になる。畳み込みをFFTで行うとき、ある方法でゼロ埋めをすれば、配列のサイズを2の冪乗に変更しても結果が変わらないようにすることができる。

図解

このようになる。
ポイントとしては、行と列、両方の入れ替えが必要なことと、$x_0$だけ別に計算することである。
添字が1から始まっているので、{2,4,3,1}ではなく{3,5,4,2}になっている。

群論との関係

Raderのアルゴリズムは群論の言葉で "素数 p に対する有限体 $F_p=\mathbb Z/p\mathbb Z$の乗法群は位数p-1の巡回群と同型である：$F^\times \cong C_{p-1}$" ということを利用している。
どういうことかというと、pで割った余りの世界の「掛け算の表」はもしpが素数ならば、p-1で割った余りの世界の「足し算の表」と、「添字の並び替えと添字のはりかえれをすれば同じ」になるということである。

面白いことに、FFTのアルゴリズムは群論の定理と関係しているようである。

FFT アルゴリズム	群論の定理
Raderのアルゴリズム	$(\mathbb Z/p\mathbb Z)^\times \cong C_{p-1}$
PrimeFactor型FFT	$\mathbb{Z}/nm\mathbb{Z} \cong \mathbb{Z}/n\mathbb{Z} \times \mathbb{Z}/m\mathbb{Z}$ (中国の剰余定理)
おそらくCooley-Tukey法も群論の何らかの性質と関係していると思う。なぜ群論がフーリエ変換の高速化に使えるかというと、フーリエ変換で用いる回転因子の行列は$(\mathbb{Z}/n\mathbb{Z})^\times$ と同じ対称性を持っているので、フーリエ変換で行われる$n^2$個の計算にも、この対称性が反映されるのである。対称性があるということは、$n^2$個の計算は完全にバラバラではなく、重複する計算があったり、共通する「中間の計算」を持っていたりするので、これを利用すれば、計算量を$n \log n$個に減らせるということなのである。

参考ページ
余り(剰余)の性質をプログラムに活かす
 有理数とか有限体とかのはなし

実装

実装のポイント

予め三角関数の値を配列に格納する。
Raderのアルゴリズムは配列を生成する必要があるため、要素数がある程度大きくないと逆に遅くなる。このためある程度要素数が大きくなったら使う。

入力としてstd::vector<std::complex<double>>型、std::complex<double>*型などが使えるように、テンプレートを使っている。

三角関数を計算するルーチンの実装

DFTでは$2\pi$の等分した角度の整数倍の三角関数の値を使うので、値が正確に0,1,-1になるときがある。少しでも正確な結果にするために、そういうときにきっかりとした値を返す関数を作った。
この関数は最初に回転因子の配列を作るときに呼び出されるだけで、メインループの中で何度も呼び出すわけではないので、if文により少し遅くなっても問題ない。

三角関数計算ルーチン

double fft_cos(int n, int m){
	if(n == 0){return 1;}
	if(n == m){/*2 Pi*/return 1;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return -1;}
		if(m/n == 4){/*Pi/2*/return 0;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return 0;}
	}
	return cos(2*M_PI*n/m);
}
double fft_sin(int n, int m){
	if(n == 0){return 0;}
	if(n == m){/*2 Pi*/return 0;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return 0;}
		if(m/n == 4){/*Pi/2*/return 1;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return -1;}
	}
	return sin(2*M_PI*n/m);
}

### DFTの実装

DFTの実装

template<typename T>
void dft(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n
	){
	int m = size/n;
	for(int i=0;i<size;i++){
		temp[i] = 0;
	}
	for(int i=0;i<n;i++){
		for(int j=0;j<m;j++){
			for(int k=0;k<m;k++){
				//std::cout << std::floor(size*std::arg(omega[((j*k)%m)*n])/(2*M_PI)) << "\t";
				temp[i*m+j] += array[i*m+k]*omega[((j*k)%m)*n];
			}
			//std::cout << "\n";
		}
		//std::cout << "\n";
	}
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
}

このルーチンは、与えられた配列をn等分した部分配列にそれぞれDFTをほどこす。
tempは作業用配列、omegaは回転因子$\omega_n$が格納されている配列である。
このルーチンはFFTで呼び出す用途なので、FFTで生成したこれらの配列を受け取るようになっている。

Cooley-Tukeyの実装

template<typename T> void fft(
	T &array,
	int size,
	int sign
){
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,size));
		omega[i].imag(sign*fft_sin(i,size));
	}
	int primes[32];
	int np=0;
	int n=1;
	for(;;){
		int m = size/n;
		int p = 2 + (m&1);
		for(;p*p<=m;p+=2){
			if(m%p==0){
				/*m is divisible by p*/
				int q = m/p;
				/*permutation*/
				for(int i=0;i<n;i++){
					for(int j=0;j<p;j++){
						for(int k=0;k<q;k++){
							temp[i*m+q*j+k] = array[i*m+p*k+j];
						}
					}
				}
				for(int i=0;i<size;i++){
					array[i] = temp[i];
				}
				primes[np++]=p;
				n *= p;
				if(q <= FFT_HARDCORD_SIZE){break;}
				goto next;
			}
		}
		break;
		next:;
	}
	/*bottom of recursion*/
	/*perform dft on n sub-arrays*/
	{
		if(size/n <= FFT_HARDCORD_SIZE){
			fft_hardcord(array,size,n,sign);
		}else if(size/n <= FFT_PRIME_THRESHOLD){
			dft(array,temp,omega,size,n);
		}else{
			fft_prime(array,size,n,sign);
		}
	}
	/*sum up with twiddle factors*/
	for(int h=np-1;h>=0;h--){
		int m = size/n;
		int p = primes[h];
		int q = n/p;
		int r = size/q;
		
		for(int i=0;i<size;i++){temp[i] = 0;}
		for(int i=0;i<q;i++){
			for(int j=0;j<p;j++){
				for(int k=0;k<m;k++){
					for(int l=0;l<p;l++){
						temp[i*p*m+(j*m+k)] += array[i*p*m+(l*m+k)]*omega[(l*(j*m+k)%r)*q];
					}
				}
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
		n = q;
	}
}

まずはじめのforループでは、試し割りによって要素数を因数分解しながら、並べかえを行う。
見つけた素因数は配列primesに格納していく。要素数はint型なので素因数の数が３２を超えることはない。これで元のDFTを小さなサイズのDFTに分解していくわけだが、サイズが4くらいになると、分解を続けるよりも、一気にやってしまったほうが速いので、サイズがある程度（FFT_HARDCORD_SIZEによって決める）になったらそこで分解をやめる。
分解が終わったら、分割されたDFTのサイズは素数になっている。このサイズが十分小さかったら、ハードコードのルーチンを呼び出す。逆にある程度大きかったらRaderのアルゴリズムを使う。その中間では通常のDFTのルーチンを呼び出す。
最後に回転因子を掛けながら分割された部分のDFTの結果を足し合わせて、全体のDFTを得る。

Raderの実装

template<typename T> void fft_prime(T &array,int size,int n,int sign){
	int p = size/n;
	int len=p-1;
	bool power2=true;
	int i=0;
	for(;len>1;len>>=1,i++){
		if(len & 1){power2=false;}
	}
	len = power2 ? p-1 : (1<<(i+2));
	/*
	std::cout << "p-1 : " << p-1 << "\n";
	std::cout << "len : " << len << "\n";
	*/
	int pad = len - (p-1);
	int g = primePrimitiveRoot(p);
	int ig = powerMod(g,p-2,p);
	std::vector<Complex> data(len),temp(len),omega(len),omega2(len);
	for(int i=0;i<len;i++){
		omega[i].real(fft_cos(i,len));
		omega[i].imag(-fft_sin(i,len));
		omega2[i].real(fft_cos(powerMod(ig,i%(p-1),p),p));
		omega2[i].imag(sign*fft_sin(powerMod(ig,i%(p-1),p),p));
	}
	fft2(omega2,temp,omega);
	for(int i=0;i<n;i++){
		data[0] = array[i*p+1];
		for(int j=1;j<len;j++){
			if(j <= pad){
				data[j] = 0;
			}else{
				data[j] = array[i*p+powerMod(g,j-pad,p)];
			}
		}
		/*===Convolution theorem===*/
		fft2(data,temp,omega);
		for(int j=0;j<len;j++){temp[j] = data[j]*omega2[j];}
		for(int j=0;j<len;j++){
			data[j] = temp[j];
			omega[j].imag(-omega[j].imag());
		}
		fft2(data,temp,omega);
		/*add array[i*p] term*/
		for(int j=0;j<len;j++){
			data[j].real(data[j].real()/len + array[i*p].real());
			data[j].imag(data[j].imag()/len + array[i*p].imag());
			//data[j] = data[j]/(Complex(len)) + array[i*p];
			//omega[j].imag(-omega[j].imag());
		}
		/*===Convolution theorem end===*/
		
		/*DC term*/
		temp[0]=0;
		for(int j=0;j<p;j++){temp[0] += array[i*p+j];}
		array[i*p] = temp[0];
		
		/*non-DC term*/
		for(int j=0;j<p-1;j++){array[i*p+powerMod(ig,j,p)] = data[j];}
	}
}

このルーチンは、与えられた配列をn等分した部分配列にそれぞれDFTをほどこす。以下はざっくりした説明。

まずp-1が2の冪乗になっているか確認し、もしなっていなければ、2*(p-1)より大きな最小の2の冪乗を計算する。これはゼロ埋めしてFFTをしたとき、畳込みが正しく復元できるようにするために必要な長さである。
関数primePrimitiveRootによってpの単位元を取得し、変数gに入れる。
変数igにgの逆元を入れる。pを法としてgのp-2乗をするとgの逆元になる。例えばp=5なら、g=2は{2,4,3,1}という数列を生み出し、ig=3は{3,4,2,1}という数列を生み出す。要するに、igはgと逆方向に進む数列を生み出す。
回転因子$\omega_n^{g^{q-p}}$の配列omega2と、FFTで使うための回転因子の配列omegaを作る。
長さを2の冪乗にした配列にゼロ埋めした元の配列を格納する。またこのときに並べ替えもすます。ゼロ埋めのやり方だが、普通に元の配列の後ろをゼロで埋めるのではなく、元の配列の１番目と２番目の間をゼロで埋める。畳み込みを正しく復元するためにはこうするといいのである。(参考:Haskell FFT 9: Prime-Length FFT With Rader's Algorithm
畳み込み定理によって畳み込みを計算する。サイズが2の冪乗専用のFFTルーチンfft2_modを呼び出す。
最後に並べ替えを戻し、周波数ゼロの成分と合わせて完成

原始根を求めるルーチン、冪剰余の実装

Primitive Rootを参考にした。

原始根を求めるルーチン、冪剰余の実装

int primePrimitiveRoot(int p){
	int i,t;
	std::vector<int> list;
	
	for(i=2,t=p-1;i*i<=t;i++){
		if(t % i == 0){
			list.push_back(i);
			do{
				t/=i;
			}while(t%i==0);
		}
	}
	
	for(std::vector<int>::iterator it = list.begin();it != list.end();it++){
		 *it = (p-1) / (*it);
	}
	for(i=2;i<=p-1;i++){
		for(std::vector<int>::iterator it = list.begin();it != list.end();it++){
			if(powerMod(i,*it,p)==1){
				goto loopend;
			}
		}
		break;
	  loopend: ;
	}
	return i;
}

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;
}

例えばp=7の場合、原始根の候補は{1,2,3,4,5,6}である。これらの生成する数列をみると以下のようになる。

{1,1,1,1,1,1}
{2,4,1,2,4,1}
{3,2,6,4,5,1}
{4,2,1,4,2,1}
{5,4,6,2,3,1}
{6,1,6,1,6,1}

これみると、1から6までの全ての数字を生み出しているのは3と5なので、3と5が7の原始根であることが分かる。他の数列をみると、全て周期的になっていて、その周期は1,2,3,6のどれかになっている。
実はpの剰余が生成する数列は必ずp-1の約数の周期を持っている。これは群論のラグランジュの定理として知られる。なのでp-1の約数で冪乗したときに1になってるかどうかを調べることで、原始根かどうかが調べられる。

なぜ１番目と２番目の間をゼロ埋めするのか

巡回畳み込み積分は以下の式で表される。

(f \ast g)_i = \sum_{j=0}^{n-1}f_jg_{i-j}

ただし$g_{-j} = g_{n-j}$である。これを書き下してみると

\begin{align}
(f \ast g)_0 &= f_0 g_0 + f_1 g_{n-1} + f_2 g_{n-2} + f_3 g_{n-3} + \cdots \\
(f \ast g)_1 &= f_0 g_1 + f_1 g_0 + f_2 g_{n-1} + f_3 g_{n-2} + \cdots \\
(f \ast g)_2 &= f_0 g_2 + f_1 g_1 + f_2 g_0 + f_3 g_{n-1} + \cdots \\
&\cdots
\end{align}

これを見ると$g$の先頭と末尾はfの１番目と２番目と掛け合わされることがわかる。Raderのアルゴリズムの場合、$g_i$はある周期$n$を持っている。今$p$個ゼロ埋めをしたとすると、データ数はn+pになり、畳み込みの式は、

\begin{align}
(f \ast g)_0 &= f_0 g_0 + f_1 g_{n+p-1} +  f_2 g_{n+p-2} ... +  f_p g_{n} + f_{p+1} g_{n-1} + f_{p+2} g_{n-2} + \cdots \\
(f \ast g)_1 &= f_0 g_1 + f_1 g_0 +  f_2 g_{n+p-1} ... +  f_{p+1} g_{n} + f_{p+2} g_{n-1} + f_{p+3} g_{n-2} + \cdots \\
(f \ast g)_2 &= f_0 g_2 + f_1 g_1 +  f_2 g_0 + f_3 g_{n+p-1}  ... +  f_{p+2} g_{n} + f_{p+3} g_{n-1} + f_{p+4} g_{n-2} + \cdots \\
&\cdots
\end{align}

となる。$g_{n-1}$に掛かかっている$f$の添字に注目すると、元の畳込みと同じになるようにするには、{$f_1,f_2,..f_p$}までをゼロ埋めすればいいということがわかる。

サイズが2の冪乗専用のFFTルーチンの実装

template<typename T> void fft2(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size
	){
	int n=1;
	for(;n<size;n <<= 1){
		int m = size/n;
		int q = m>>1;
		/*permutation*/
		for(int i=0;i<n;i++){
			for(int k=0;k<q;k++){
				temp[i*m+k] = array[i*m+2*k];
				temp[i*m+q+k] = array[i*m+2*k+1];
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
	}
	for(;n>1;n >>= 1){
		int m = size/n;
		int q = n>>1;
		int r = size/q;
		/*adding up with twiddle factors*/
		for(int i=0;i<size;i++){temp[i] = 0;}
		for(int i=0;i<q;i++){
			for(int k=0;k<m;k++){
				temp[2*m*i+k] += array[2*m*i+k];
				temp[2*m*i+m+k] += array[2*m*i+k];
				temp[2*m*i+k] += array[2*m*i+m+k]*omega[(k%r)*q];
				temp[2*m*i+m+k] += array[2*m*i+m+k]*omega[((m+k)%r)*q];
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
	}
}

サイズが合成数のルーチンから、因数分解の計算を除き、ループも展開してある。

ハードコートのルーチンの実装

ハードコート部分を生成するコードfft_hardcorderを作ってコードを生成した。これによってハードコートを生成し、FFTのプログラムからはヘッダファイルとしてインクルードする。

ハードコートのルーチンの実装

fft_hardcorder.cpp

# include <iostream>
# include <cmath>
double fft_cos(int n, int m){
	if(n == 0){return 1;}
	if(n == m){/*2 Pi*/return 1;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return -1;}
		if(m/n == 4){/*Pi/2*/return 0;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return 0;}
	}
	return cos(2*M_PI*n/m);
}
double fft_sin(int n, int m){
	if(n == 0){return 0;}
	if(n == m){/*2 Pi*/return 0;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return 0;}
		if(m/n == 4){/*Pi/2*/return 1;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return -1;}
	}
	return sin(2*M_PI*n/m);
}

int main(int argc, char *argv[]){
	FILE *fp;
	char *fileout;
	int i,j,k;
	int n,p,q;
	int flag=0;
	if(argc < 3){
		printf("usage : %s [outputfile] [number]\n",argv[0]);
		return -1;
	}
	n = atoi(argv[2]);
	
	fp = fopen(argv[1],"w");
	fprintf(fp,"#ifndef FFT_HARDCORD_H\n");
	fprintf(fp,"#define FFT_HARDCORD_H\n");
	fprintf(fp,"#include <vector>\n");
	fprintf(fp,"using Complex = std::complex<double>;\n");

	fprintf(fp,"#define FFT_HARDCORD_SIZE %d\n",n);
	fprintf(fp,"inline void fft_hardcord(T &array,int size,int n,int sign){\n\n");
	fprintf(fp,"int m = size/n;");
	fprintf(fp,"switch(m){\n");
	fprintf(fp,"\t  case 1:\n");
	fprintf(fp,"\t\t break;\n");
	for(i=2;i<=n;i++){
		fprintf(fp,"\t  case %d:\n",i);
		fprintf(fp,"\t\tfor(int i=0;i<n;i++){\n");
		for(j=0;j<i;j++){
			fprintf(fp,"\t\t\tdouble r%d = ",j);
			flag = 0;
			for(k=0;k<i;k++){
				p=(k*j)%i; q=i;
				if(p==0){
					if(!flag){flag=1;}else{fprintf(fp," + ");}
					fprintf(fp,"array[i*%d+%d].real()",i,k);
					continue;
				}
				if(q%p==0){
					if(q/p==2){
						if(!flag){flag=1;}else{fprintf(fp," + ");}
						fprintf(fp,"(-array[i*%d+%d].real())",i,k);
						continue;
					}
					if(q/p==4){
						continue;
					}
				}
				if(q%(q-p)==0){
					if(q/(q-p)==4){
						continue;
					}
				}
				if(!flag){flag=1;}else{fprintf(fp," + ");}
				fprintf(fp,"(%.14e)*(array[i*%d+%d].real())",cos(2*M_PI*p/q),i,k);
			}
			fprintf(fp,"\n\t\t\t");
			flag=0;
			for(k=0;k<i;k++){
				p=(k*j)%i; q=i;
				if(p==0){
					continue;
				}
				if(q%p==0){
					if(q/p==2){
						continue;
					}
					if(q/p==4){
						if(!flag){flag=1;fprintf(fp," + ");}else{fprintf(fp," + ");}
						fprintf(fp,"(-sign)*(array[i*%d+%d].imag())",i,k);
						continue;
					}
				}
				if(q%(q-p)==0){
					if(q/(q-p)==4){
						if(!flag){flag=1;fprintf(fp," + ");}else{fprintf(fp," + ");}
						fprintf(fp,"sign*(array[i*%d+%d].imag())",i,k);
						continue;
					}
				}
				if(!flag){flag=1;fprintf(fp," + ");}else{fprintf(fp," + ");}
				fprintf(fp,"(-sign)*(%.14e)*(array[i*%d+%d].imag())",sin(2*M_PI*p/q),i,k);
			}
			fprintf(fp,";\n");
		}
		for(j=0;j<i;j++){
			fprintf(fp,"\t\t\tdouble i%d = ",j);
			flag=0;
			for(k=0;k<i;k++){
				p=(k*j)%i; q=i;
				if(p==0){
					continue;
				}
				if(q%p==0){
					if(q/p==2){
						continue;
					}
					if(q/p==4){
						if(!flag){flag=1;}else{fprintf(fp," + ");}
						fprintf(fp,"sign*(array[i*%d+%d].real())",i,k);
						continue;
					}
				}
				if(q%(q-p)==0){
					if(q/(q-p)==4){
						if(!flag){flag=1;}else{fprintf(fp," + ");}
						fprintf(fp,"(-sign)*(array[i*%d+%d].real())",i,k);
						continue;
					}
				}
				if(!flag){flag=1;}else{fprintf(fp," + ");}
				fprintf(fp,"sign*(%.14e)*(array[i*%d+%d].real())",sin(2*M_PI*p/q),i,k);
			}
			fprintf(fp,"\n\t\t\t");
			flag=0;
			for(k=0;k<i;k++){
				p=(k*j)%i; q=i;
				if(p==0){
					if(!flag){flag=1;fprintf(fp," + ");}else{fprintf(fp," + ");}
					fprintf(fp,"(array[i*%d+%d].imag())",i,k);
					continue;
				}
				if(q%p==0){
					if(q/p==2){
						if(!flag){flag=1;fprintf(fp," + ");}else{fprintf(fp," + ");}
						fprintf(fp,"(-array[i*%d+%d].imag())",i,k);
						continue;
					}
					if(q/p==4){
						continue;
					}
				}
				if(q%(q-p)==0){
					if(q/(q-p)==4){
						continue;
					}
				}
				if(!flag){flag=1;fprintf(fp," + ");}else{fprintf(fp," + ");}
				fprintf(fp,"(%.14e)*(array[i*%d+%d].imag())",cos(2*M_PI*p/q),i,k);
			}
			fprintf(fp,";\n");
		}
		for(int j=0;j<i;j++){
			fprintf(fp,"\t\t\tarray[i*%d+%d].real(r%d);\n",i,j,j);
			fprintf(fp,"\t\t\tarray[i*%d+%d].imag(i%d);\n",i,j,j);
		}
		fprintf(fp,"\t\t}\n");
		fprintf(fp,"\t\tbreak;\n");
	}
	fprintf(fp,"}\n");
	fprintf(fp,"\n};\n");
	fprintf(fp,"#endif\n");
	fclose(fp);
	return 0;
}

fft_hardcorderにより生成されたハードコード

fft_hardcord.h

# ifndef FFT_HARDCORD_H
# define FFT_HARDCORD_H
# include <vector>
using Complex = std::complex<double>;
# define FFT_HARDCORD_SIZE 4
inline void fft_hardcord(std::vector<Complex> &array,std::vector<Complex> &temp,int n,int sign){

int m = array.size()/n;
switch(m){
	  case 1:
		 break;
	  case 2:
		for(int i=0;i<n;i++){
			double r0 = array[i*2+0].real() + array[i*2+1].real()
			;
			double r1 = array[i*2+0].real() + (-array[i*2+1].real())
			;
			double i0 = 
			 + (array[i*2+0].imag()) + (array[i*2+1].imag());
			double i1 = 
			 + (array[i*2+0].imag()) + (-array[i*2+1].imag());
			array[i*2+0].real(r0);
			array[i*2+0].imag(i0);
			array[i*2+1].real(r1);
			array[i*2+1].imag(i1);
		}
		break;
	  case 3:
		for(int i=0;i<n;i++){
			double r0 = array[i*3+0].real() + array[i*3+1].real() + array[i*3+2].real()
			;
			double r1 = array[i*3+0].real() + (-5.00000000000000e-01)*(array[i*3+1].real()) + (-5.00000000000000e-01)*(array[i*3+2].real())
			 + (-sign)*(8.66025403784439e-01)*(array[i*3+1].imag()) + (-sign)*(-8.66025403784438e-01)*(array[i*3+2].imag());
			double r2 = array[i*3+0].real() + (-5.00000000000000e-01)*(array[i*3+1].real()) + (-5.00000000000000e-01)*(array[i*3+2].real())
			 + (-sign)*(-8.66025403784438e-01)*(array[i*3+1].imag()) + (-sign)*(8.66025403784439e-01)*(array[i*3+2].imag());
			double i0 = 
			 + (array[i*3+0].imag()) + (array[i*3+1].imag()) + (array[i*3+2].imag());
			double i1 = sign*(8.66025403784439e-01)*(array[i*3+1].real()) + sign*(-8.66025403784438e-01)*(array[i*3+2].real())
			 + (array[i*3+0].imag()) + (-5.00000000000000e-01)*(array[i*3+1].imag()) + (-5.00000000000000e-01)*(array[i*3+2].imag());
			double i2 = sign*(-8.66025403784438e-01)*(array[i*3+1].real()) + sign*(8.66025403784439e-01)*(array[i*3+2].real())
			 + (array[i*3+0].imag()) + (-5.00000000000000e-01)*(array[i*3+1].imag()) + (-5.00000000000000e-01)*(array[i*3+2].imag());
			array[i*3+0].real(r0);
			array[i*3+0].imag(i0);
			array[i*3+1].real(r1);
			array[i*3+1].imag(i1);
			array[i*3+2].real(r2);
			array[i*3+2].imag(i2);
		}
		break;
	  case 4:
		for(int i=0;i<n;i++){
			double r0 = array[i*4+0].real() + array[i*4+1].real() + array[i*4+2].real() + array[i*4+3].real()
			;
			double r1 = array[i*4+0].real() + (-array[i*4+2].real())
			 + (-sign)*(array[i*4+1].imag()) + sign*(array[i*4+3].imag());
			double r2 = array[i*4+0].real() + (-array[i*4+1].real()) + array[i*4+2].real() + (-array[i*4+3].real())
			;
			double r3 = array[i*4+0].real() + (-array[i*4+2].real())
			 + sign*(array[i*4+1].imag()) + (-sign)*(array[i*4+3].imag());
			double i0 = 
			 + (array[i*4+0].imag()) + (array[i*4+1].imag()) + (array[i*4+2].imag()) + (array[i*4+3].imag());
			double i1 = sign*(array[i*4+1].real()) + (-sign)*(array[i*4+3].real())
			 + (array[i*4+0].imag()) + (-array[i*4+2].imag());
			double i2 = 
			 + (array[i*4+0].imag()) + (-array[i*4+1].imag()) + (array[i*4+2].imag()) + (-array[i*4+3].imag());
			double i3 = (-sign)*(array[i*4+1].real()) + sign*(array[i*4+3].real())
			 + (array[i*4+0].imag()) + (-array[i*4+2].imag());
			array[i*4+0].real(r0);
			array[i*4+0].imag(i0);
			array[i*4+1].real(r1);
			array[i*4+1].imag(i1);
			array[i*4+2].real(r2);
			array[i*4+2].imag(i2);
			array[i*4+3].real(r3);
			array[i*4+3].imag(i3);
		}
		break;
}

};
# endif

コード全体

二次元FFTのルーチンfft2d、Prime-Factor型fftコードfft_pf(Prime-Factor FFT)、再帰関数を使うfft_recも入っている。

コード全体

fft.h

# ifndef FFT_H
# define FFT_H
# include <iostream>
# include <complex>
# include <vector>
# include <list>
# include <fft_hardcord.h>
# define FFT_PRIME_THRESHOLD 256
using Complex = std::complex<double>;

template<typename T> void fft(T &array, int size);
template<typename T> void ifft(T &array, int size);
template<typename T> void fft2d(T &array, int row, int column);
template<typename T> void ifft2d(T &array, int row, int column);
template<typename T> void dft(T &array, int size);
template<typename T> void idft(T &array, int size);

void fft(std::vector<Complex> &array);
void ifft(std::vector<Complex> &array);
void dft(std::vector<Complex> &array);
void idft(std::vector<Complex> &array);

/*internal functions*/
template<typename T>
void dft(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n
);

template<typename T> void fft2(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size
);

template<typename T> void fft_prime(T &array,int size,int n,int sign);
template<typename T> void fft(T &array, int size, int sign);
template<typename T> void fft2d(T &array, int row, int column, int sign);

double fft_cos(int n, int m);
double fft_sin(int n, int m);
int primePrimitiveRoot(int p);
int powerMod(int a, int b, int m);
int eulerPhi(int n);
int primitiveRoot(int n);
int primitiveRoot(int n, int phi);

std::vector<std::vector<int>> factorInteger(int size);
std::vector<std::vector<int>> tally(std::vector<int> list);
void print(std::vector<std::vector<int>> &list);

void fft_pf(std::vector<Complex> &array);
template<typename T> void fft_pf(T &array, int size);
template<typename T> void fft_pf(T &array, int size, int sign);
template<typename T> void fft_pf(
	T &array,
	std::vector<Complex> temp,
	std::vector<Complex> omega,
	int size,
	int sign,
	std::vector<std::vector<int>> list,
	int index
);

void fft_rec(std::vector<Complex> &array);
template<typename T> void fft_rec(T &array, int size);
template<typename T> void fft_rec(T &array, int size, int sign);
template<typename T> void fft_rec(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n, 
	int sign
);
# endif

fft.cpp

# include <fft.h>
# include <iostream>

void dft(std::vector<Complex> &array){dft(array,array.size());}
template void dft(std::vector<Complex> &array, int size);
template void dft(Complex* &array, int size);

template void dft(
	std::vector<Complex> &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n
);
	template void dft(
	Complex* &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n
);

template<typename T>
void dft(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n
	){
	int m = size/n;
	for(int i=0;i<size;i++){
		temp[i] = 0;
	}
	for(int i=0;i<n;i++){
		for(int j=0;j<m;j++){
			for(int k=0;k<m;k++){
				//std::cout << std::floor(size*std::arg(omega[((j*k)%m)*n])/(2*M_PI)) << "\t";
				temp[i*m+j] += array[i*m+k]*omega[((j*k)%m)*n];
			}
			//std::cout << "\n";
		}
		//std::cout << "\n";
	}
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
}
	
template<typename T> void dft(T &array, int size){
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,omega.size()));
		omega[i].imag(-fft_sin(i,omega.size()));
	}
	dft(array,temp,omega,size,1);
}
	
template<typename T> void idft(T &array, int size){
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,omega.size()));
		omega[i].imag(fft_sin(i,omega.size()));
	}
	dft(array,temp,omega,size,-1);
	for(int i=0;i<size;i++){array[i] /= size;}
}
	
template<typename T> void fft2(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size
	){
	int n=1;
	for(;n<size;n <<= 1){
		int m = size/n;
		int q = m>>1;
		/*permutation*/
		for(int i=0;i<n;i++){
			for(int k=0;k<q;k++){
				temp[i*m+k] = array[i*m+2*k];
				temp[i*m+q+k] = array[i*m+2*k+1];
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
	}
	for(;n>1;n >>= 1){
		int m = size/n;
		int q = n>>1;
		int r = size/q;
		/*adding up with twiddle factors*/
		for(int i=0;i<size;i++){temp[i] = 0;}
		for(int i=0;i<q;i++){
			for(int k=0;k<m;k++){
				temp[2*m*i+k] += array[2*m*i+k];
				temp[2*m*i+m+k] += array[2*m*i+k];
				temp[2*m*i+k] += array[2*m*i+m+k]*omega[(k%r)*q];
				temp[2*m*i+m+k] += array[2*m*i+m+k]*omega[((m+k)%r)*q];
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
	}
}
	
void fft2(
	std::vector<Complex> &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega
	){
	fft2(array,temp,omega,array.size());
}

template<typename T> void fft_prime(T &array,int size,int n,int sign){
	int p = size/n;
	int len=p-1;
	bool power2=true;
	int i=0;
	for(;len>1;len>>=1,i++){
		if(len & 1){power2=false;}
	}
	len = power2 ? p-1 : (1<<(i+2));
	/*
	std::cout << "p-1 : " << p-1 << "\n";
	std::cout << "len : " << len << "\n";
	*/
	int pad = len - (p-1);
	int g = primePrimitiveRoot(p);
	int ig = powerMod(g,p-2,p);
	std::vector<Complex> data(len),temp(len),omega(len),omega2(len);
	for(int i=0;i<len;i++){
		omega[i].real(fft_cos(i,len));
		omega[i].imag(-fft_sin(i,len));
		omega2[i].real(fft_cos(powerMod(ig,i%(p-1),p),p));
		omega2[i].imag(sign*fft_sin(powerMod(ig,i%(p-1),p),p));
	}
	fft2(omega2,temp,omega);
	for(int i=0;i<n;i++){
		data[0] = array[i*p+1];
		for(int j=1;j<len;j++){
			if(j <= pad){
				data[j] = 0;
			}else{
				data[j] = array[i*p+powerMod(g,j-pad,p)];
			}
		}
		/*===Convolution theorem===*/
		fft2(data,temp,omega);
		for(int j=0;j<len;j++){temp[j] = data[j]*omega2[j];}
		for(int j=0;j<len;j++){
			data[j] = temp[j];
			omega[j].imag(-omega[j].imag());
		}
		fft2(data,temp,omega);
		/*add array[i*p] term*/
		for(int j=0;j<len;j++){
			data[j].real(data[j].real()/len + array[i*p].real());
			data[j].imag(data[j].imag()/len + array[i*p].imag());
			//data[j] = data[j]/(Complex(len)) + array[i*p];
			//omega[j].imag(-omega[j].imag());
		}
		/*===Convolution theorem end===*/
		
		/*DC term*/
		temp[0]=0;
		for(int j=0;j<p;j++){temp[0] += array[i*p+j];}
		array[i*p] = temp[0];
		
		/*non-DC term*/
		for(int j=0;j<p-1;j++){array[i*p+powerMod(ig,j,p)] = data[j];}
	}
}

template<typename T> void fft(
	T &array,
	int size,
	int sign
){
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,size));
		omega[i].imag(sign*fft_sin(i,size));
	}
	int primes[32];
	int np=0;
	int n=1;
	for(;;){
		int m = size/n;
		int p = 2 + (m&1);
		for(;p*p<=m;p+=2){
			if(m%p==0){
				/*m is divisible by p*/
				int q = m/p;
				/*permutation*/
				for(int i=0;i<n;i++){
					for(int j=0;j<p;j++){
						for(int k=0;k<q;k++){
							temp[i*m+q*j+k] = array[i*m+p*k+j];
						}
					}
				}
				for(int i=0;i<size;i++){
					array[i] = temp[i];
				}
				primes[np++]=p;
				n *= p;
				if(q <= FFT_HARDCORD_SIZE){break;}
				goto next;
			}
		}
		break;
		next:;
	}
	/*bottom of recursion*/
	/*perform dft on n sub-arrays*/
	{
		if(size/n <= FFT_HARDCORD_SIZE){
			fft_hardcord(array,size,n,sign);
		}else if(size/n <= FFT_PRIME_THRESHOLD){
			dft(array,temp,omega,size,n);
		}else{
			fft_prime(array,size,n,sign);
		}
	}
	/*sum up with twiddle factors*/
	for(int h=np-1;h>=0;h--){
		int m = size/n;
		int p = primes[h];
		int q = n/p;
		int r = size/q;
		
		for(int i=0;i<size;i++){temp[i] = 0;}
		for(int i=0;i<q;i++){
			for(int j=0;j<p;j++){
				for(int k=0;k<m;k++){
					for(int l=0;l<p;l++){
						temp[i*p*m+(j*m+k)] += array[i*p*m+(l*m+k)]*omega[(l*(j*m+k)%r)*q];
					}
				}
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
		n = q;
	}
}

void fft(std::vector<Complex> &array){fft(array,array.size());}
template void fft(std::vector<Complex> &array, int size, int sign);
template void fft(std::vector<Complex> &array, int size);
template void fft(Complex* &array, int size, int sign);
template void fft(Complex* &array, int size);
template<typename T> void fft(T &array, int size){fft(array,size,-1);}

template<typename T> void ifft(T &array, int size){
	fft(array,size,1);
	for(int i=0;i<size;i++){array[i] /= size;}
}
	
template<typename T> void fft2d(T &array, int row, int column, int sign){
	int size = row*column;
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,size));
		omega[i].imag(sign*fft_sin(i,size));
	}

	int primes[32];
	int np=0;
	int n=row;
	for(;;){
		int m = size/n;
		int p = 2 + (m&1);
		for(;p*p<=m;p+=2){
			if(m%p==0){
				/*m is divisible by p*/
				int q = m/p;
				/*permutation*/
				for(int i=0;i<n;i++){
					for(int j=0;j<p;j++){
						for(int k=0;k<q;k++){
							temp[i*m+q*j+k] = array[i*m+p*k+j];
						}
					}
				}
				for(int i=0;i<size;i++){
					array[i] = temp[i];
				}
				primes[np++]=p;
				n *= p;
				if(q <= FFT_HARDCORD_SIZE){break;}
				goto next1;
			}
		}
		break;
		next1:;
	}
	/*bottom of recursion*/
	/*perform dft on n sub-arrays*/
	{
		if(size/n <= FFT_HARDCORD_SIZE){
			fft_hardcord(array,size,n,sign);
		}else if(size/n <= FFT_PRIME_THRESHOLD){
			dft(array,temp,omega,size,n);
		}else{
			fft_prime(array,size,n,sign);
		}
	}
	/*integrate array elements multiplied by twiddle factors*/
	for(int h=np-1;h>=0;h--){
		int m = size/n;
		int p = primes[h];
		int q = n/p;
		int r = size/q;
		
		for(int i=0;i<size;i++){temp[i] = 0;}
		for(int i=0;i<q;i++){
			for(int j=0;j<p;j++){
				for(int k=0;k<m;k++){
					for(int l=0;l<p;l++){
						temp[i*p*m+(j*m+k)] += array[i*p*m+(l*m+k)]*omega[(l*(j*m+k)%r)*q];
					}
				}
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
		n = q;
	}
	
	/*transpose*/
	for(int i=0;i<row;i++){
		for(int j=0;j<column;j++){
			temp[j*row+i] = array[i*column+j];
		}
	}
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
	
	np=0;
	n=column;
	for(;;){
		int m = size/n;
		int p = 2 + (m&1);
		for(;p*p<=m;p+=2){
			if(m%p==0){
				/*m is divisible by p*/
				int q = m/p;
				/*permutation*/
				for(int i=0;i<n;i++){
					for(int j=0;j<p;j++){
						for(int k=0;k<q;k++){
							temp[i*m+q*j+k] = array[i*m+p*k+j];
						}
					}
				}
				for(int i=0;i<size;i++){
					array[i] = temp[i];
				}
				primes[np++]=p;
				n *= p;
				if(q <= FFT_HARDCORD_SIZE){break;}
				goto next2;
			}
		}
		break;
		next2:;
	}
	/*bottom of recursion*/
	/*perform dft on n sub-arrays*/
	{
		if(size/n <= FFT_HARDCORD_SIZE){
			fft_hardcord(array,size,n,sign);
		}else if(size/n <= FFT_PRIME_THRESHOLD){
			dft(array,temp,omega,size,n);
		}else{
			fft_prime(array,size,n,sign);
		}
	}
	/*sum up with twiddle factors*/
	for(int h=np-1;h>=0;h--){
		int m = size/n;
		int p = primes[h];
		int q = n/p;
		int r = size/q;
		
		for(int i=0;i<size;i++){temp[i] = 0;}
		for(int i=0;i<q;i++){
			for(int j=0;j<p;j++){
				for(int k=0;k<m;k++){
					for(int l=0;l<p;l++){
						temp[i*p*m+(j*m+k)] += array[i*p*m+(l*m+k)]*omega[(l*(j*m+k)%r)*q];
					}
				}
			}
		}
		for(int i=0;i<size;i++){array[i] = temp[i];}
		n = q;
	}
	
	/*transpose*/
	for(int i=0;i<column;i++){
		for(int j=0;j<row;j++){
			temp[j*column+i] = array[i*row+j];
		}
	}
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
}
	
template<typename T> void fft2d(T &array, int row, int column){fft2d(array,row,column,-1);}
template<typename T> void ifft2d(T &array, int row, int column){
	fft2d(array,row,column,1);
	for(int i=0;i<row*column;i++){array[i] /= (row*column);}
}

void ifft(std::vector<Complex> &array){ifft(array,array.size());}
void idft(std::vector<Complex> &array){idft(array,array.size());}

double fft_cos(int n, int m){
	if(n == 0){return 1;}
	if(n == m){/*2 Pi*/return 1;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return -1;}
		if(m/n == 4){/*Pi/2*/return 0;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return 0;}
	}
	return cos(2*M_PI*n/m);
}
double fft_sin(int n, int m){
	if(n == 0){return 0;}
	if(n == m){/*2 Pi*/return 0;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return 0;}
		if(m/n == 4){/*Pi/2*/return 1;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return -1;}
	}
	return sin(2*M_PI*n/m);
}

int eulerPhi(int n){
	int res = n;
	for(int i=2;i*i<=n;i++){
		if(n%i==0){
			while(n%i==0){n/=i;}
			res -= res/i;
		}
	}
	if(n>1){
		res -= res/n;
	}
	return res;
}

int primitiveRoot(int n){
	return primitiveRoot(n,eulerPhi(n));
}

int primitiveRoot(int n, int phi){
	int i,t;
	std::vector<int> list;
	
	for(i=2,t=phi;i*i<=t;i++){
		if(t % i == 0){
			list.push_back(i);
			do{
				t/=i;
			}while(t%i==0);
		}
	}
	
	for(std::vector<int>::iterator it = list.begin();it != list.end();it++){
		 *it = phi / (*it);
	}
	
	for(i=2;i<=phi;i++){
		for(std::vector<int>::iterator it = list.begin();it != list.end();it++){
			if(powerMod(i,*it,phi)==1){
				goto loopend;
			}
		}
		break;
	  loopend:;
	}
	return i;
}

int primePrimitiveRoot(int p){
	int i,t;
	std::vector<int> list;
	
	for(i=2,t=p-1;i*i<=t;i++){
		if(t % i == 0){
			list.push_back(i);
			do{
				t/=i;
			}while(t%i==0);
		}
	}
	
	for(std::vector<int>::iterator it = list.begin();it != list.end();it++){
		 *it = (p-1) / (*it);
	}
	for(i=2;i<=p-1;i++){
		for(std::vector<int>::iterator it = list.begin();it != list.end();it++){
			if(powerMod(i,*it,p)==1){
				goto loopend;
			}
		}
		break;
	  loopend:;
	}
	return i;
}

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;
}

template void ifft(std::vector<Complex> &array, int size);
template void fft2d(std::vector<Complex> &array, int row, int column);
template void ifft2d(std::vector<Complex> &array, int row, int column);
template void idft(std::vector<Complex> &array, int size);

template void fft2(
	std::vector<Complex> &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size
);

template void fft_prime(std::vector<Complex> &array,int size,int n,int sign);
template void fft2d(std::vector<Complex> &array, int row, int column, int sign);
template void ifft(Complex* &array, int size);
template void fft2d(Complex* &array, int row, int column);
template void ifft2d(Complex* &array, int row, int column);
template void idft(Complex* &array, int size);
	
template void fft2(
	Complex* &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size
);
template void fft_prime(Complex* &array,int size,int n,int sign);
template void fft2d(Complex* &array, int row, int column, int sign);

std::vector<std::vector<int>> factorInteger(int size){
	std::vector<int> primes;
	int n=1;
	for(;;){
		int m = size/n;
		int p = 2 + (m&1);
		for(;p*p<=m;p+=2){
			if(m%p==0){
				/*m is divisible by p*/
				primes.push_back(p);
				n *= p;
				goto next;
			}
		}
		break;
		next:;
	}
	primes.push_back(size/n);
	return tally(primes);
}
	
std::vector<std::vector<int>> tally(std::vector<int> list){
	std::vector<std::vector<int>> res;
	for(std::vector<int>::iterator it = list.begin();it!=list.end();it++){
		for(std::vector<std::vector<int>>::iterator it2 = res.begin();it2!=res.end();it2++){
			if(it2->front() == *it){
				it2->back() += 1;
				goto loopend;
			}
		}
		res.push_back(std::vector<int>{*it,1});
		loopend:;
	}
	return res;
}

void print(std::vector<std::vector<int>> &list){
	std::cout << "{";
	for(std::vector<std::vector<int>>::iterator it = list.begin();it!=list.end();it++){
		std::cout << "{";
		for(std::vector<int>::iterator it2 = it->begin();it2!=it->end();it2++){
			std::cout << *it2;
			if(std::next(it2) != it->end()){std::cout << ",";}
		}
		std::cout << "}";
		if(std::next(it) != list.end()){std::cout << ",";}
	}
	std::cout << "}\n";
}

int power(int a, int b){
	int res = 1;
	for(int i=0;i<b;i++){res *= a;}
	return res;
}

void fft_pf(std::vector<Complex> &array){fft_pf(array,array.size());}
template<typename T> void fft_pf(T &array, int size){fft_pf(array,size,-1);}
template void fft_pf(std::vector<Complex> &array, int size);
template void fft_pf(Complex* &array, int size);
template void fft_pf(std::vector<Complex> &array, int size, int sign);
template void fft_pf(Complex* &array, int size, int sign);

template<typename T> void fft_pf(T &array, int size, int sign){
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,size));
		omega[i].imag(sign*fft_sin(i,size));
	}
	fft_pf(array,temp,omega,size,sign,factorInteger(size),0);
}

template<typename T> void fft_pf(
	T &array,
	std::vector<Complex> temp,
	std::vector<Complex> omega,
	int size,
	int sign,
	std::vector<std::vector<int>> list,
	int index
){
	//print(list);
	if(index==list.size()){return;}
	
	int nn = 1;
	for(int i=0;i<index;i++){
		nn *= power(list[i][0],list[i][1]);
	}
	int mm =size/nn;
	
	int n = power(list[index][0],list[index][1]);
	int m = mm/n;
		
	/*CRT mapping*/
	
	{
		for(int ii=0;ii<nn;ii++){
			int j=0,k=0;
			for(int i=0;i<mm;i++,j++,k++){
				if(j>=n){j=0;} if(k>=m){k=0;}
				temp[ii*mm+j+n*k] = array[ii*mm+i];
			}
		}
	}
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
	
	/*===FFT2D===*/
	
	{
		//dft(array,temp,omega,size,size/n);
		fft_rec(array,temp,omega,size,size/n,sign);
	}
	
	/*transpose*/
	for(int ii=0;ii<nn;ii++){
		for(int i=0;i<n;i++){
			for(int j=0;j<m;j++){
				temp[ii*mm+j+i*m] = array[ii*mm+i+j*n];
			}
		}
	}
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
	
	{
		fft_pf(array,temp,omega,size,sign,list,index+1);
	}
	
	/*===FFT2D end===*/
	
	/*good mapping*/
	
	{
		for(int ii=0;ii<nn;ii++){
			int i=0;
			for(int j=0;j<n;j++){
				for(int k=0;k<m;k++,i++){
					temp[ii*mm+(m*j+n*k)%mm] = array[ii*mm+i];
				}
			}
		}
	}
	
	for(int i=0;i<size;i++){
		array[i] = temp[i];
	}
}

void fft_rec(std::vector<Complex> &array){fft_rec(array,array.size());}
template<typename T> void fft_rec(T &array, int size){fft_rec(array,size,-1);}
template void fft_rec(std::vector<Complex> &array, int size);
template void fft_rec(Complex* &array, int size);
template void fft_rec(std::vector<Complex> &array, int size, int sign);
template void fft_rec(Complex* &array, int size, int sign);

template void fft_rec(std::vector<Complex> &array, int size, int sign);
template void fft_rec(Complex* &array, int size, int sign);
template void fft_rec(
	std::vector<Complex> &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n, 
	int sign
);
template void fft_rec(
	Complex* &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n, 
	int sign
);

template<typename T> void fft_rec(T &array, int size, int sign){
	std::vector<Complex> temp(size),omega(size);
	for(int i=0;i<size;i++){
		omega[i].real(fft_cos(i,size));
		omega[i].imag(sign*fft_sin(i,size));
	}
	fft_rec(array,temp,omega,size,1,sign);
}

template<typename T> void fft_rec(
	T &array,
	std::vector<Complex> &temp,
	std::vector<Complex> &omega,
	int size,
	int n, 
	int sign
){	
	int m = size/n;
	int p = 2 + (m&1);
	
	for(;p*p<=m;p+=2){
		if(m%p==0){
			/*m is divisible by p*/
			int q = m/p;
			/*permutation*/
			for(int i=0;i<n;i++){
				for(int j=0;j<p;j++){
					for(int k=0;k<q;k++){
						temp[i*m+q*j+k] = array[i*m+p*k+j];
					}
				}
			}
			for(int i=0;i<size;i++){
				array[i] = temp[i];
			}
			n *= p;
			m /= p;
			fft_rec(array,temp,omega,size,n,sign);
			q = n/p;
			int r = size/q; 
		
			for(int i=0;i<size;i++){temp[i] = 0;}
			for(int i=0;i<q;i++){
				for(int j=0;j<p;j++){
					for(int k=0;k<m;k++){
						for(int l=0;l<p;l++){
							temp[i*p*m+(j*m+k)] += array[i*p*m+(l*m+k)]*omega[(l*(j*m+k)%r)*q];
						}
					}
				}
			}
			for(int i=0;i<size;i++){array[i] = temp[i];}
			return;
		}
	}
	/*bottom of recursion*/
	/*perform dft on n sub-arrays*/
	if(size/n <= FFT_HARDCORD_SIZE){
		fft_hardcord(array,size,n,sign);
	}else if(size/n <= FFT_PRIME_THRESHOLD){
		dft(array,temp,omega,size,n);
	}else{
		fft_prime(array,size,n,sign);
	}
}

動作テストコードとその他

fft_test.cpp

# include <fft.h>
# include <Array2D.h>
# include <vector>
int main(int argc, char *argv[]){
	Array2D a,b;
	if(argc > 1){
		a.read(argv[1]);
		b.set(a.rows,2*a.columns);
		if(a.columns >= 2){
			std::vector<Complex> array(a.rows);
			for(int i=0;i<a.rows;i++){
				array[i].real(a[i][0]);
				array[i].imag(a[i][1]);
			}
			fft(array);
			for(int i=0;i<b.rows;i++){
				b[i][0] = array[i].real();
				b[i][1] = array[i].imag();
			}
			for(int i=0;i<a.rows;i++){
				array[i].real(a[i][0]);
				array[i].imag(a[i][1]);
			}
			fft_pf(array);
			for(int i=0;i<b.rows;i++){
				b[i][2] = array[i].real();
				b[i][3] = array[i].imag();
			}

			std::streambuf *buf;
			if(argc > 2){
				buf = std::ofstream(argv[2]).rdbuf();
			}else{
				buf = std::cout.rdbuf();
			}
			std::ostream os(buf);
			b.print(os);
		}
	}else{
		std::cout << "usage : " << argv[0] << " [input file] [output file]\n";
	}
}

Array2D.h

# ifndef ARRAY2D_H
# define ARRAY2D_H
# include <iostream>
# include <fstream>
# include <list>
# include <string>
# include <string.h>
# include <sys/stat.h>

struct Array2D{
	struct ReadOptions{
		enum {automatic=-1};
		int rows;
		int columns;
		int rowOffset;
		int columnOffset;
		std::list<std::string> lineSeparators;
		std::list<std::string> fieldSeparators;
		bool repeatedSeparators;
		bool ignoreEmptyLines;
		
		ReadOptions();
	};
	enum precision{p_double,p_float};
	double **data;
	int rows;
	int columns;
	
	Array2D();
	Array2D(int,int);
	Array2D(const Array2D&);
	Array2D(Array2D&&);
	Array2D& operator=(const Array2D &array);
	Array2D& operator=(Array2D &&array);
	virtual ~Array2D();
	
	void set(int rows, int columns);
	void clear();
	void print();
	void print(const char *fileName);
	void print(std::ostream &os);
	void exportBMAT(const char *fileName,precision p);
	Array2D& read(const char *fileName, ReadOptions &options);
	Array2D& read(const char *fileName);
	double interpolate(double x);
	double min();	
	double minRow(int row);
	double minColumn(int column);
	double*& operator[] (const int index);
  private:
	struct Matcher{
		Matcher(std::list<std::string>&);
		~Matcher();
		void clear();
		int count(char *str,bool repeatedSeparators);
		char *replace(char *p, char c);
		char *match(char *p);
		char *tokenize(char **str,bool repeatedSeparators);
		void tokenize(char *str,std::list<char*> &list,bool repeatedSeparators);
		size_t size;
		const char **list[2];
	};
};
# endif

Array2D.cpp

# include <Array2D.h>

Array2D::Array2D(){
	this->rows=0;
	this->columns=0;
	this->data=nullptr;
}

Array2D::Array2D(int rows, int columns){
	this->set(rows,columns);
}

Array2D::Array2D(const Array2D &array){
	this->set(array.rows,array.columns);
	for(int i=0;i<rows;i++){
		for(int j=0;j<columns;j++){
			data[i][j] = array.data[i][j];
		}
	}
}

Array2D::Array2D(Array2D &&array){
	rows = array.rows;
	columns = array.columns;
	data = array.data;
}

Array2D& Array2D::operator=(const Array2D &array){
	Array2D tmp(array);
	*this = std::move(tmp);
	tmp.data=nullptr;
	return *this;
}

Array2D& Array2D::operator=(Array2D &&array){
	if(this==&array){
		return *this;
	}
	rows = array.rows;
	columns = array.columns;
	data = array.data;
	return *this;
}

Array2D::~Array2D(){
	this->clear();
}

void Array2D::set(int rows, int columns){
	this->rows = rows;
	this->columns = columns;
	data = new double*[rows];
	data[0] = new double[rows*columns];
	for(int i=1;i<rows;i++){
		data[i] = &(data[0][i*columns]);
	}
}

void Array2D::clear(){
	if(data){
		delete[] data[0];
		delete[] data;
		data = nullptr;
	}
}

void Array2D::print(){
	std::ostream os(std::cout.rdbuf());
	print(os);
}

void Array2D::print(const char *fileName){
	std::ostream os(std::ofstream(fileName).rdbuf());
	print(os);
}

void Array2D::print(std::ostream &os){
	os.precision(17);
	for(int i=0;i<rows;i++){
		for(int j=0;j<columns;j++){
			os << std::fixed << data[i][j] << "\t";
		}
		os << "\n";
	}
}

void Array2D::exportBMAT(const char *fileName, Array2D::precision p){
	FILE *fp = fopen(fileName,"wb");
	if(!fp){std::cout << "exportBMAT : unable to open file " << fileName << "\n"; return;}
	fwrite(&p,sizeof(precision),1,fp);
	fwrite(&rows,sizeof(int),1,fp);
	fwrite(&columns,sizeof(int),1,fp);
	int w1;
	float *f;
	switch(p){
	  case p_double:
		if((w1=fwrite(data[0],sizeof(double),rows*columns,fp))!=rows*columns){
			std::cout << "exportBMAT : tried to write " <<
			rows*columns*sizeof(double) << " bytes, but " << 
			w1*sizeof(double) << " bytes were written\n";
		}
		break;
	  case p_float:
		f = new float[rows*columns];
		for(int i=0;i<rows;i++){
			for(int j=0;j<columns;j++){
				f[i*columns+j] = data[i][j];
			}
		}
		if((w1=fwrite(f,sizeof(float),rows*columns,fp))!=rows*columns){
			std::cout << "exportBMAT : tried to write " <<
			 rows*columns*sizeof(float) << " bytes, but " << 
			 w1*sizeof(float) << " bytes were written\n";
		}
		delete[] f;
		break;
	}
	fclose(fp);
}

Array2D::ReadOptions::ReadOptions() :
lineSeparators({"\r\n","\n","\r"}),
fieldSeparators({"\t"," "})
{
	rows = ReadOptions::automatic;
	columns = ReadOptions::automatic;
	rowOffset = 0;
	columnOffset = 0;
	repeatedSeparators = true;
	ignoreEmptyLines = false;
}

Array2D& Array2D::read(const char *fileName, ReadOptions &options){
	struct stat st;
	size_t size;
	if(stat(fileName,&st)==0){size = st.st_size+1;}
	FILE *fp = fopen(fileName,"r");
	if(!fp){
		std::cout << "unable to open " << fileName << "\n";
		return *this;
	}
	setbuf(fp,NULL);
	char *buffer = new char[size];
	fread(buffer,1,size,fp);
	buffer[size-1] = '\0';
	fclose(fp);
	
	Matcher lineMatcher(options.lineSeparators);
	Matcher fieldMatcher(options.fieldSeparators);
	
	std::list<char*> listRow;
	lineMatcher.tokenize(buffer,listRow,options.repeatedSeparators);

	for(std::list<char*>::iterator it=listRow.begin();it!=listRow.end();it++){
		if(fieldMatcher.count(*it,options.repeatedSeparators)==0){
			*it = NULL;
		}
	}
	listRow.remove(NULL);

	rows = listRow.size();
	int rowOffset = options.rowOffset < 0 ? 
					rows + options.rowOffset : 
					options.rowOffset;

	if(options.rows == ReadOptions::automatic){
		if(rows < rowOffset){
			std::cout << "row offset " << rowOffset << "exceeds the row size " << rows << "\n";
			return *this;
		}
		rows = rows - rowOffset;
	}else{
		if(rows < options.rows + rowOffset){
			std::cout << "row size "<< options.rows << 
			" + rowOffset " << rowOffset << " exceeds the row size of input " << rows << "\n";
			return *this;
		}
		rows = options.rows;
	}
	
	std::list<char*>::iterator it = listRow.begin();
	for(int i=0;it != listRow.end() && i<rowOffset;it++,i++){}
	
	columns = fieldMatcher.count(*it,options.repeatedSeparators);
	int columnOffset = options.columnOffset < 0 ?
					   columns + options.columnOffset :
					   options.columnOffset;

	if(options.columns == ReadOptions::automatic){
		if(columns < columnOffset){
			std::cout << "column offset " << columnOffset << "exceeds the column size " << rows << "\n";
			return *this;
		}
		columns = columns - columnOffset;
	}else{
		if(columns < options.columns + columnOffset){
			std::cout << "column size "<< options.columns << 
			" + columnOffset " << columnOffset << " exceeds the column size of input " << columns << "\n";
			return *this;
		}
		columns = options.columns;
	}
	//std::cout << "rows : " << rows << std::endl;
	//std::cout << "columns : " << columns << "\n";
	
	this->set(rows,columns);
	
	for(int i=0;it!=listRow.end() && i<rows;it++,i++){
		char *p = *it;
		for(int j=0;j<columnOffset;j++){
			fieldMatcher.clear();
			fieldMatcher.tokenize(&p,options.repeatedSeparators);
		}
		for(int j=0;j<columns;j++){
			fieldMatcher.clear();
			//fieldMatcher.tokenize(&p,options.repeatedSeparators);
			data[i][j] = atof(fieldMatcher.tokenize(&p,options.repeatedSeparators));
		}
	}
	return *this;
}

Array2D& Array2D::read(const char *fileName){
	ReadOptions options;
	return read(fileName,options);
}

double Array2D::interpolate(double x){
	// assumes x data are in 0th column, y data are in 1th column;
	// assumes x data are sorted in increasing order.
	
	if(x < data[0][0]){
		return (x - data[0][0])*(data[1][1]-data[0][1])/(data[1][0]-data[0][0]) + data[0][1];
	}
	if(x > data[rows-1][0]){
		return (x - data[rows-1][0])*(data[rows-1][1]-data[rows-2][1])/(data[rows-1][0]-data[rows-2][0]) + data[rows-1][1];
	}
	
	int i=0;
	for(;i<rows && x >= data[i][0];i++){}
	return (x-data[i][0])*(data[i][1]-data[i-1][1])/(data[i][0]-data[i-1][0]) + data[i][1];
}

double*& Array2D::operator[] (const int index){
	return data[index];
}

double Array2D::min(){
	double min = data[0][0];
	for(int i=0;i<rows*columns;i++){
		if(min > data[0][i]){min = data[0][i];}
	}
	return min;
}

double Array2D::minRow(int row){
	double min = data[row][0];
	for(int i=0;i<columns;i++){
		if(min > data[row][i]){min = data[row][i];}
	}
	return min;
}

double Array2D::minColumn(int column){
	double min = data[0][column];
	for(int i=0;i<rows;i++){
		if(min > data[i][column]){min = data[i][column];}
	}
	return min;
}

Array2D::Matcher::Matcher(std::list<std::string> &list){
	size = list.size();
	this->list[0] = new const char*[size];
	this->list[1] = new const char*[size];
	int i=0;
	for(std::list<std::string>::iterator it=list.begin();
		it != list.end();
		it++,i++
	){
		this->list[0][i] = this->list[1][i] = it->c_str();
	}
}
Array2D::Matcher::~Matcher(){
	delete[] list[0];
	delete[] list[1];
}

void Array2D::Matcher::clear(){
	for(int i=0;i<size;i++){list[1][i] = list[0][i];}
}

char *Array2D::Matcher::match(char *p){
	for(int i=0;i<size;i++){
		if(*list[1][i]==*p){
			list[1][i]++;
			if(*list[1][i]=='\0'){
				for(int j=0;j<size;j++){list[1][j] = list[0][j];}
				int len = strlen(list[0][i]);
				return p-len+1; 
			}
		}else{
			list[1][i] = list[0][i];
		}
	}
	return NULL;
}

char *Array2D::Matcher::replace(char *p, char c){
	for(int i=0;i<size;i++){
		if(*list[1][i]==*p){
			list[1][i]++;
			if(*list[1][i]=='\0'){
				for(int j=0;j<size;j++){list[1][j] = list[0][j];}
				int len = strlen(list[0][i]);
				char *pp = p;
				for(;p-pp < len;pp--){*pp=c;}
				return pp+1; 
			}
		}else{
			list[1][i] = list[0][i];
		}
	}
	return NULL;
}

int Array2D::Matcher::count(char *str,bool repeatedSeparators){
	int counter=0;
	char *p,*ps,*pe,*pm,*tail;
	tail = pe = strchr(str,'\0');
	for(ps=p=str;*p;p++){
		pm = this->match(p);
		if(pm){
			pe = pm;
			if(!repeatedSeparators || ps!=pe){++counter;}
			ps = p+1;
		}
	}
	if(*ps && ps < tail){counter++;}
	return counter;
}

char *Array2D::Matcher::tokenize(char **str,bool repeatedSeparators){
	char *p,*ps,*pe,*pm;
	pe = strchr(*str,'\0');
	for(ps=p=*str;*p;p++){
		pm = this->replace(p,'\0');
		if(pm){
			pe = pm;
			if(!repeatedSeparators || ps!=pe){
				*str = p+1; return ps;
			}
			ps = p+1;
		}
	}
	return ps;
}

void Array2D::Matcher::tokenize(char *str, std::list<char*> &list, bool repeatedSeparators){
	char *p,*ps,*pe,*pm,*tail;
	
	tail = pe = strchr(str,'\0');
	for(ps=p=str;*p;p++){
		pm = this->replace(p,'\0');
		if(pm){
			pe = pm;
			if(!repeatedSeparators || ps != pe){
				list.push_back(ps);
			}
			ps = p+1;
		}
	}
	if(*ps && ps < tail){list.push_back(ps);}
}

ベンチマーク

赤 : FFT
青 : 通常のDFT
緑 : 素数要素数のFFT(Raderが呼び出される場合)

FFTは通常のDFTよりも約100倍速い。
合成数の場合は約数の組み合わせによって速さにバラつきがでる。最速なのはサイズが小さな素数の冪乗のときである。
Raderのアルゴリズムは合成数のFFTと比べて遅い。ただ、それでも通常のDFTよりは10倍ほど速い。
Raderのアルゴリズムはゼロ埋めでサイズを2の冪乗にするため、計算時間は階段状に増えていく。
素数のサイズが256以下のときは通常のDFTのほうがRaderのアルゴリズムより早かった。これはサイズが小さいと内部で配列を生成するコストのほうが大きくなるためと思われる。