fftのソース

現状のfftソースを晒してみる。

#include <stdio.h> #include <stdlib.h> #include <math.h> typedef struct { double r ; double i ; } ad_t ; /* ===== fft ======== */ void *fft(int nx ,ad_t a[]) { int i, k, m, n; double tmpr, tmpi, r, cosw, sinw; n = 1<<nx; r = -M_PI / (double)n; for(m = n/2; m > 0 ; m /= 2 ) { r *= 2; for(k = 0 ; k < m; k++) { cosw = cos(r*(double)k); sinw = sin(r*(double)k); for(i= k ; i < n; i += m*2) { tmpr = a[i].r - a[m+i].r; tmpi = a[i].i - a[m+i].i; a[i].r = a[i].r + a[m+i].r; a[i].i = a[i].i + a[m+i].i; a[m+i].r = tmpr * cosw - tmpi * sinw; a[m+i].i = tmpi * cosw + tmpr * sinw; } } } return; } /* ===== ifft ======== */ void *ifft(int nx,ad_t a[]) { int i, k, m,n; double tmpr, tmpi, r, cosw, sinw; n = 1<< nx ; r = M_PI ; for(m = 1; m < n ; m *= 2 ) { for(k = 0; k < m; k++) { cosw = cos(r*(double)k); sinw = sin(r*(double)k); for(i= k ; i < n; i += m*2) { tmpr = cosw * a[m+i].r - sinw * a[m+i].i; tmpi = cosw * a[m+i].i + sinw * a[m+i].r; a[m+i].r = a[i].r - tmpr ; a[m+i].i = a[i].i - tmpi; a[i].r = a[i].r + tmpr ; a[i].i = a[i].i + tmpi; } } r *= 0.5; } return; } /* a * b -> a */ void *ri_mul(int nx,ad_t a[],ad_t b[]) { int m ; double tmpa , tmpb ; a[0].r = a[0].r*b[0].r; a[0].i = a[0].i*b[0].i; for(m=1;m<(1<<nx);m++){ tmpa = a[m].r; tmpb = b[m].r; a[m].r = tmpa*tmpb - a[m].i*b[m].i; a[m].i = tmpa*b[m].i + a[m].i*tmpb; } } /* スケーリングと繰越 n-2 -> 0 方向 */ void *ri_carry(int nx,ad_t a[],int redex) { double pn,tmp; int n , m; long car , t; n = 1<<nx ; pn = 1.0/(double)n; car = 0 ; for(m=n-2;m>=0;m--){ tmp = a[m].r*pn + (double)car ; if( tmp >= 0.0 ){ t = (long)(tmp+0.5) ; car = t / redex ; a[m].r =(double)( t % redex ); }else{ t = (long)(tmp-0.5) ; car = t / redex - 1; a[m].r =(double)( t % redex + redex ); } } a[n-1].r=(double)car; } /* スケーリングと繰越 0 -> n-2 方向 */ void *ri_carry_r(int nx,ad_t a[],int redex) { double pn,tmp; int n , m; long car , t; n = 1<<nx ; pn = 1.0/(double)n; car = 0 ; for(m=0;m<n-1;m++){ tmp = a[m].r*pn + (double)car ; if( tmp >= 0.0 ){ t = (long)(tmp+0.5) ; car = t / redex ; a[m].r =(double)( t % redex ); }else{ t = (long)(tmp-0.5) ; car = t / redex - 1; a[m].r =(double)( t % redex + redex ); } } a[n-1].r=(double)car; } /* 任意の値のビット反転 */ inline unsigned int bitrv_v(int nx, unsigned int a ) { a = ((a & 0x55555555)<< 1 ) | ((a >> 1) & 0x55555555 ); a = ((a & 0x33333333)<< 2 ) | ((a >> 2) & 0x33333333 ); a = ((a & 0x0f0f0f0f)<< 4 ) | ((a >> 4) & 0x0f0f0f0f ); a = ((a & 0x00ff00ff)<< 8 ) | ((a >> 8) & 0x00ff00ff ); a = ((a & 0x0000ffff)<<16 ) | ((a >>16) & 0x0000ffff ); return (a >> (32-nx)) ; } /* ビット反転 インクリメンタル */ inline unsigned int inc_rvbit(int nx , unsigned int a ) { unsigned int rb ; rb = 1<<nx; do { rb >>= 1; a ^= rb; /* xor は、繰越なしで1bitの足し算したものと同じ */ } while ( rb > a ); /* xor の結果が、1の時終わり */ return a; } /* ビット反転位置のデータ入れ替え */ void *bitrv(int nx , ad_t a[]) { int i ,j; int nm1 , nh ; ad_t t; nh = 1<<(nx-1); nm1 = (1<<nx) -1 ; j = 0 ; for (i = 1; i < nh ; i++) { j = inc_rvbit(nx , j ); if (j > i) { t = a[j]; a[j] = a[i]; a[i] = t; if( j < nh ){ t = a[nm1-j]; a[nm1-j]=a[nm1-i]; a[nm1-i]=t; } } } }

以下動作確認

input data n = 16 = 1<<4 0 a[ 0].r = 2.000000 a[ 0].i = 0.000000 b[ 0].r = 9.000000 b[ 0].i = 0.000000 1 a[ 1].r = 1.000000 a[ 1].i = 0.000000 b[ 1].r = 8.000000 b[ 1].i = 0.000000 2 a[ 2].r = 2.000000 a[ 2].i = 0.000000 b[ 2].r = 7.000000 b[ 2].i = 0.000000 3 a[ 3].r = 3.000000 a[ 3].i = 0.000000 b[ 3].r = 6.000000 b[ 3].i = 0.000000 4 a[ 4].r = 4.000000 a[ 4].i = 0.000000 b[ 4].r = 5.000000 b[ 4].i = 0.000000 5 a[ 5].r = 5.000000 a[ 5].i = 0.000000 b[ 5].r = 4.000000 b[ 5].i = 0.000000 6 a[ 6].r = 6.000000 a[ 6].i = 0.000000 b[ 6].r = 2.000000 b[ 6].i = 0.000000 7 a[ 7].r = 7.000000 a[ 7].i = 0.000000 b[ 7].r = 3.000000 b[ 7].i = 0.000000 8 a[ 8].r = 0.000000 a[ 8].i = 0.000000 b[ 8].r = 0.000000 b[ 8].i = 0.000000 9 a[ 9].r = 0.000000 a[ 9].i = 0.000000 b[ 9].r = 0.000000 b[ 9].i = 0.000000 10 a[10].r = 0.000000 a[10].i = 0.000000 b[10].r = 0.000000 b[10].i = 0.000000 11 a[11].r = 0.000000 a[11].i = 0.000000 b[11].r = 0.000000 b[11].i = 0.000000 12 a[12].r = 0.000000 a[12].i = 0.000000 b[12].r = 0.000000 b[12].i = 0.000000 13 a[13].r = 0.000000 a[13].i = 0.000000 b[13].r = 0.000000 b[13].i = 0.000000 14 a[14].r = 0.000000 a[14].i = 0.000000 b[14].r = 0.000000 b[14].i = 0.000000 15 a[15].r = 0.000000 a[15].i = 0.000000 b[15].r = 0.000000 b[15].i = 0.000000 fft(nx,a) fft(nx,b) 0 a[ 0].r = 30.000000 a[ 0].i = 0.000000 b[ 0].r = 44.000000 b[ 0].i = 0.000000 1 a[ 8].r = -7.137071 a[ 8].i = -20.109358 b[ 8].r = 17.920298 b[ 8].i = -24.812274 2 a[ 4].r = -2.000000 a[ 4].i = 9.656854 b[ 4].r = 4.707107 b[ 4].i = -9.949747 3 a[12].r = 4.380086 a[12].i = -5.986423 b[12].r = 5.530124 b[12].i = -7.699802 4 a[ 2].r = -2.000000 a[ 2].i = 4.000000 b[ 2].r = 5.000000 b[ 2].i = -3.000000 5 a[10].r = 5.276769 a[10].i = -2.672715 b[10].r = 5.398808 b[10].i = -4.971880 6 a[ 6].r = -2.000000 a[ 6].i = 1.656854 b[ 6].r = 3.292893 b[ 6].i = 0.050253 7 a[14].r = 5.480217 a[14].i = -0.795649 b[14].r = 7.150769 b[14].i = -2.084352 8 a[ 1].r = -2.000000 a[ 1].i = 0.000000 b[ 1].r = 2.000000 b[ 1].i = 0.000000 9 a[ 9].r = 5.480217 a[ 9].i = 0.795649 b[ 9].r = 7.150769 b[ 9].i = 2.084352 10 a[ 5].r = -2.000000 a[ 5].i = -1.656854 b[ 5].r = 3.292893 b[ 5].i = -0.050253 11 a[13].r = 5.276769 a[13].i = 2.672715 b[13].r = 5.398808 b[13].i = 4.971880 12 a[ 3].r = -2.000000 a[ 3].i = -4.000000 b[ 3].r = 5.000000 b[ 3].i = 3.000000 13 a[11].r = 4.380086 a[11].i = 5.986423 b[11].r = 5.530124 b[11].i = 7.699802 14 a[ 7].r = -2.000000 a[ 7].i = -9.656854 b[ 7].r = 4.707107 b[ 7].i = 9.949747 15 a[15].r = -7.137071 a[15].i = 20.109358 b[15].r = 17.920298 b[15].i = 24.812274 ri_mul(nx,a,b) 0 a[ 0].r = 1320.000000 a[ 0].i = 0.000000 1 a[ 8].r = -626.857348 a[ 8].i = -183.278730 2 a[ 4].r = 86.669048 a[ 4].i = 65.355339 3 a[12].r = -21.871852 a[12].i = -66.831453 4 a[ 2].r = 2.000000 a[ 2].i = 26.000000 5 a[10].r = 15.199846 a[10].i = -40.664931 6 a[ 6].r = -6.669048 a[ 6].i = 5.355339 7 a[14].r = 37.529354 a[14].i = -17.112207 8 a[ 1].r = -4.000000 a[ 1].i = 0.000000 9 a[ 9].r = 37.529354 a[ 9].i = 17.112207 10 a[ 5].r = -6.669048 a[ 5].i = -5.355339 11 a[13].r = 15.199846 a[13].i = 40.664931 12 a[ 3].r = 2.000000 a[ 3].i = -26.000000 13 a[11].r = -21.871852 a[11].i = 66.831453 14 a[ 7].r = 86.669048 a[ 7].i = -65.355339 15 a[15].r = -626.857348 a[15].i = 183.278730 ifft(nx,a) 0 a[ 0].r = 288.000000 a[ 0].i = -0.000000 1 a[ 1].r = 400.000000 a[ 1].i = -0.000000 2 a[ 2].r = 640.000000 a[ 2].i = -0.000000 3 a[ 3].r = 992.000000 a[ 3].i = 0.000000 4 a[ 4].r = 1440.000000 a[ 4].i = -0.000000 5 a[ 5].r = 1968.000000 a[ 5].i = 0.000000 6 a[ 6].r = 2528.000000 a[ 6].i = -0.000000 7 a[ 7].r = 3216.000000 a[ 7].i = -0.000000 8 a[ 8].r = 2672.000000 a[ 8].i = 0.000000 9 a[ 9].r = 2208.000000 a[ 9].i = 0.000000 10 a[10].r = 1744.000000 a[10].i = 0.000000 11 a[11].r = 1296.000000 a[11].i = -0.000000 12 a[12].r = 880.000000 a[12].i = -0.000000 13 a[13].r = 512.000000 a[13].i = -0.000000 14 a[14].r = 336.000000 a[14].i = -0.000000 15 a[15].r = -0.000000 a[15].i = 0.000000 copy a -> b ri_carry(nx,a,10) ri_carry_r(nx,b,10) 0 a[15].r = 2.000000 a[15].i = 0.000000 b[15].r = 2.000000 b[15].i = 0.000000 1 a[ 0].r = 0.000000 a[ 0].i = -0.000000 b[14].r = 4.000000 b[14].i = -0.000000 2 a[ 1].r = 9.000000 a[ 1].i = -0.000000 b[13].r = 8.000000 b[13].i = -0.000000 3 a[ 2].r = 7.000000 a[ 2].i = -0.000000 b[12].r = 4.000000 b[12].i = -0.000000 4 a[ 3].r = 2.000000 a[ 3].i = 0.000000 b[11].r = 3.000000 b[11].i = -0.000000 5 a[ 4].r = 4.000000 a[ 4].i = -0.000000 b[10].r = 4.000000 b[10].i = 0.000000 6 a[ 5].r = 0.000000 a[ 5].i = 0.000000 b[ 9].r = 6.000000 b[ 9].i = 0.000000 7 a[ 6].r = 9.000000 a[ 6].i = -0.000000 b[ 8].r = 8.000000 b[ 8].i = 0.000000 8 a[ 7].r = 9.000000 a[ 7].i = -0.000000 b[ 7].r = 8.000000 b[ 7].i = -0.000000 9 a[ 8].r = 1.000000 a[ 8].i = 0.000000 b[ 6].r = 1.000000 b[ 6].i = -0.000000 10 a[ 9].r = 9.000000 a[ 9].i = 0.000000 b[ 5].r = 2.000000 b[ 5].i = 0.000000 11 a[10].r = 7.000000 a[10].i = 0.000000 b[ 4].r = 6.000000 b[ 4].i = -0.000000 12 a[11].r = 6.000000 a[11].i = -0.000000 b[ 3].r = 6.000000 b[ 3].i = 0.000000 13 a[12].r = 8.000000 a[12].i = -0.000000 b[ 2].r = 2.000000 b[ 2].i = -0.000000 14 a[13].r = 4.000000 a[13].i = -0.000000 b[ 1].r = 6.000000 b[ 1].i = -0.000000 15 a[14].r = 1.000000 a[14].i = -0.000000 b[ 0].r = 8.000000 b[ 0].i = -0.000000 a 21234567 * 98765423 = 2097240991976841 b 76543212 * 32456789 = 2484346881266268

fft一本で逆変換(ifft)も可能だが、ビット反転位置のデータ入れ替えを省略するため、fft,ifftを使っている。
bitrv は、調査以外では使わないのでまともに動くかちと不安。
ri_carry と ri_carry_r は、繰越の方向が異なるだけ。どちらを使うかは状況による。
現在、実数fftについて検証中。メモリー使用量が半分ですむならそれに越したことは無い。

円周率計算のプログラムを作っているのだが・・

まだこのアカウント生きてたｗｗｗ 掲題の件、世界記録を目指すとかでは無いです。
チューニングとスレッドの勉強のため負荷のかかるものを探したらそうなった。
目標　　　　　１００億桁以上
アルゴリズム　アークタンジェントのテーラ展開ガウスの式。
環境
　CPU　i7-3820　4core
　メモリ　64GB
　ホストOS　Win7
　VM　VirtualBox
　VOS　Fedora 17
　VCPU　４～６
　VMEM　54GB
コンパイラ　gcc
言語　　　C
環境依存の最適化を行うけど、アセンブラは使わない予定
今のところ　GPU も使う予定なし。
主にfftのチューニングになるのかな？