/*
Copyright(C) 1995-98 Ryuji Matumoto (matumoto@pluto.ai.kyutech.ac.jp)

ライセンス: GNU General Public Licence(GPL)
 
このプログラムについて:
 このプログラムは，ソースコードが公開されている，作業用テーブルを使わ
ないタイプのFFTのプログラムとしては，私の調べた限り最速です．
(とおもったが世の中にはまだ速いソースがありました。)

FFTの構成:
  fftb (Complex * x, int depth, int f)
       周波数間引き型 FFT 再帰+ループ版．計算速度が速いが，計算精度悪
       い  

  fftb2 (Complex * x, int depth, int f)
       京都大学の大浦托哉氏の公開しているFFTのコードを元に三角関数の計
       算法を改良．精度改善．fftb()より速度が僅かに遅くなる．
         (http://momonga.t.u-tokyo.ac.jp/~ooura/index-j.html)

更新:
   97/7/1	関数名を統一
   98/3/2	公開用に作業関数を一つのファイルに取り込む．
                コメントを追加．
   98/3/4       三角関数の計算法を改善． 精度が大幅に改善．

使用環境:
gcc-2.7.2 以降が必要． gccの組み込み複素数型を使ってます． 

    typedef __complex__ double Complex;

と定義して使用してます．最適化フラグを付けてコンパイルしてください． 

動作確認機種: 
  Sun OS 4.1.4
  Linux(Intel) 2.0.33 + Slackware3.1

使用法:

  void fftb (Complex * x, int depth, int f)
  void fftb2 (Complex * x, int depth, int f)

   Complex *x : 入力データ．計算結果は，上書き出力される．
               配列の大きさは， 2^depth．
   int f : +1 の時は逆フーリエ変換
           -1 の時は順フーリエ変換


プログラムの解説:

キャッシュが内蔵されているCPUでは，キャッシュに入っているデータは，
高速にアクセス可能．

数値計算の書籍に上がっているFFTのプログラムは大部分がループで構成
されている．しかしこれらのプログラムはキャッシュミスが多発する．

FFTのバタフライ計算をよく調べてみると，メモリアクセスの依存性が
２進ツリーのようになっている．

このプログラムでは，キャッシュに入りきるまで，２進ツリーのように
FFTを分解し，そしてループで計算している．

このプログラムは周波数間引き型であるが，時間間引き型 FFTでも同様の
構成が可能．


今後の予定:

データが2次キャッシュより大きくなったら，David H.Bailey氏のfour-step
FFT アルゴリズムを使えば良いのではないかと思うけど，それはまた今度．
(http://science.nas.nasa.gov/Pubs/TechReports/RNRreports/dbailey/RNR-89-004/RNR-89-004.o.html)
*/

#include <stdio.h>
#include <math.h>

/* このプログラムは gcc の組み込み複素数型を使用してます．
   gcc-2.6より前では，複素数型に虫がいます．
   gcc-2.7以降でないと正常にコンパイルできません．
*/

#if !defined(__GNUC__) || __GNUC__ < 2
#error "Need GCC Ver 2.7 or later"
#endif


/* 
CACHE_SIZE_Lの決め方について．

このFFTのプログラムは，扱うデータの大きさが一次データキャッシュに入る
ようになったら別戦略を取るように設計されてます．パラメータCACHE_SIZE_L
でその分岐点を決めます．


最近の多くのマシンは一次キャッシュは，データキャッシュと命令キャッ
シュに分かれます．

   一次キャッシュ = 命令キャッシュ + データキャッシュ

例: Pentiumは 一次キャッシュは16kb
   (データキャッシュは8kb, 命令キャッシュは8kb)

一次データキャッシュが 8kbのマシン は CACHE_SIZE_L = 9
(Pentium, PentiumPro, Sparcなど)

一次データキャッシュが 16kbのマシン は CACHE_SIZE_L = 10
(PowerPCなど)

計算法: 次の一次方程式を解いてください．

   一次データキャッシュ = sizeof(double) * 2 * 2^CACHE_SIZE_L

   # sizeof(double) * 2 == 複素数型の大きさ．
   # 2^CACHE_SIZE_L == 配列の大きさ

一次データキャッシュの大きさでパラメータを計算できるはずですが，，，
世の中計算通りにはいかないようです．実測してください．

私の試した限りでは，Pentium*とSparcは CACHE_SIZE_L = 8 を推薦．
PowerPCは CACHE_SIZE_L = 9 を推薦．

*/

#define CACHE_SIZE_L 8 /* or 9 */


#define ipow2(X) (1 << (X))

#define Real double
typedef __complex__ Real Complex;

#define RE(X) __real__(X)
#define IM(X) __imag__(X)
#define I (1.0i)


static Complex
c_arg (Real arg)
{
  extern Real sin(Real);
  extern Real cos(Real);
  Complex ret;

  RE (ret) = cos (arg);
  IM (ret) = sin (arg);

  return ret;
}

static void
fftutl_div_n (Complex * x, int size)
{
  int i;

  for (i = 0; i < size; i++)
    x[i] = x[i] * (1.0 / size); 
}

static void
fftutl_bit_inv (Complex * x, int l)
{
  int n = ipow2 (l);
  int j, k, i;

  j = 0;
  for (i = 0; i < n - 1; i++)
    {
      if (i <= j)
	{
	  Complex tmp =x[i];
	  x[i] = x[j];
	  x[j] = tmp;
	}
      k = n / 2;
      while (k <= j)
	{
	  j = j - k;
	  k /= 2;
	}
      j = j + k;
    }
}

static void
fftb_loop (Complex * x_top, int depth, Complex base_w)
{
  int i;
  for (i = depth; i > 0; i--) /* 深さ */
    {
      const int nest_enrollee = ipow2 (depth - i);
      const int n = ipow2 (i);
      const int n2 = n / 2;

      Complex w1 = 1.0;
      int k;
      for (k = 0; k < n2; k++)
	{
	  int j;
	  for (j = 0; j < nest_enrollee; j++)
	    {
	      /* アドレスのバイアス */
	      Complex *x = x_top + j * n;

	      const int k0 = k + n2 * 0;
	      const int k1 = k + n2 * 1;

	      const Complex x0 = x[k0];
	      const Complex x1 = x[k1];

	      x[k0] = x0 + x1;
	      x[k1] = (x0 - x1) * w1;
	    }
	  w1 *= base_w;
	}
      base_w = base_w * base_w;
    }
}

static void
fftb_nest (Complex * x, int depth, Complex base_w)
{
  if(depth <= CACHE_SIZE_L)
    {
      /* 扱うデータが，一次データキャッシュより小さくなった． */
      fftb_loop (x, depth, base_w);
    }
  else
    {
      const int n = ipow2 (depth);
      const int n2 = n / 2;

      Complex x0, x1, w1 = 1.0;
      int i;
      for (i = 0; i < n2; i++)
	{
	  const int i0 = i + n2 * 0;
	  const int i1 = i + n2 * 1;

	  x0 = x[i0];
	  x1 = x[i1];

	  x[i0] = x0 + x1;
	  x[i1] = (x0 - x1)* w1;

	  w1 *= base_w;
	}
      
      { /* 二進ツリーのようにデータを分割 */ 
	Complex w2 = base_w * base_w;
	fftb_nest (x + n2 * 0, depth - 1, w2);
	fftb_nest (x + n2 * 1, depth - 1, w2);
      }
    }
}

void
fftb (Complex * x, int depth, int f)
{
  Complex base_w = c_arg((Real) f * 2.0 * M_PI / ipow2 (depth));

  /* 入力正順，出力ビット逆順 */
  fftb_nest (x, depth, base_w);

  if (f < 0)
    fftutl_div_n (x, ipow2 (depth));

  fftutl_bit_inv (x, depth);

  return;
}


/*** 
  fftb() の三角関数の計算部を改善
***/

static void
fftb2__loop (Complex * x_top, int depth, Complex base_w)
{
  int i, j ,k;
  for (i = depth; i > 1; i--) /* 深さ */
    {
      const int nest_enrollee = ipow2 (depth - i);
      const int n = ipow2 (i);
      const int n2 = n / 2;

      Complex x0, x1;
      Complex w1_prev = 1.0;
      Complex w1;
      Real ss;

      /* equiv w1 = basw_w * base_w */
      RE(w1) = 1 - 2 * IM(base_w) * IM(base_w);
      IM(w1) = 2 * IM(base_w) * RE(base_w);

      ss = 2 * IM(w1);
      base_w = w1;

      { /* k = 0,1; */
	for (j = 0; j < nest_enrollee; j++)
	  {
	    /* アドレスのバイアス */
	    Complex *x = x_top + j * n;	  

	    const int k0 = n2 * 0;
	    const int k1 = n2 * 1;

	    x0 = x[k0];
	    x1 = x[k1];

	    x[k0] = x0 + x1;
	    x[k1] = x0 - x1;

	    /* 1 */
	    x+=1;

	    x0 = x[k0];
	    x1 = x[k1];

	    x[k0] = x0 + x1;
	    x[k1] = w1 * (x0 -x1);
	  }
	}
      for (k = 2; k < n2; k+=2)
	{
	  RE(w1_prev) -= ss * IM(w1);
	  IM(w1_prev) += ss * RE(w1);
	  RE(w1) -= ss * IM(w1_prev);
	  IM(w1) += ss * RE(w1_prev);

	  for (j = 0; j < nest_enrollee; j++)
	    {
	      /* アドレスのバイアス */
	      Complex *x = x_top + j * n;

	      const int k0 = k + n2 * 0;
	      const int k1 = k + n2 * 1;

	      x0 = x[k0];
	      x1 = x[k1];

	      x[k0] = x0 + x1;
	      x[k1] = w1_prev * (x0 - x1);

	      /* 1 */
	      x+=1;

	      x0 = x[k0];
	      x1 = x[k1];

	      x[k0] = x0 + x1;
	      x[k1] = w1 * (x0 - x1);
	    }
	}
    }

  for (; i > 0; i--) /* 深さ */
    {
      const int nest_enrollee = ipow2 (depth - i);
      const int n = ipow2 (i);
      const int n2 = n / 2;

      for (j = 0; j < nest_enrollee; j++)
	{
	  /* アドレスのバイアス */
	  Complex *x = x_top + j * n;
	  
	  const int k0 = n2 * 0;
	  const int k1 = n2 * 1;

	  const Complex x0 = x[k0];
	  const Complex x1 = x[k1];

	  x[k0] = x0 + x1;
	  x[k1] = x0 - x1;
	}
    }
}

static void
fftb2__nest (Complex * x, int depth, Complex base_w)
{
  if(depth <= CACHE_SIZE_L)
    fftb2__loop (x, depth, base_w);
  else
    {
      const int n = ipow2 (depth);
      const int n2 = n / 2;
      int i;

      Complex x0, x1;
      Complex w1_prev = 1.0;
      Complex w1; Real ss;

      /* equiv w1 = basw_w * base_w */
      RE(w1) = 1 - 2 * IM(base_w) * IM(base_w);
      IM(w1) = 2 * IM(base_w) * RE(base_w);

      ss = 2 * IM(w1);
      base_w = w1;

      {
	x0 = x[0];
	x1 = x[n2];

	x[0] = x0 + x1;
	x[n2] = x0 - x1;

	x0 = x[1];
	x1 = x[n2+1];

	x[1] = x0 + x1;
	x[n2+1] = w1 * (x0 - x1);
      }

      for (i = 2; i < n2; i+=2)
	{
	  /* 0 */
	  int i0 = i + n2 * 0;
	  int i1 = i + n2 * 1;

	  x0 = x[i0];
	  x1 = x[i1];

	  RE(w1_prev) -= ss * IM(w1);
	  IM(w1_prev) += ss * RE(w1);

	  x[i0] = x0 + x1;
	  x[i1] = w1_prev * (x0 - x1);

	  /* 1 */
	  i0+=1;
	  i1+=1;

	  x0 = x[i0];
	  x1 = x[i1];

	  RE(w1) -= ss * IM(w1_prev);
	  IM(w1) += ss * RE(w1_prev);

	  x[i0] = x0 +x1;
	  x[i1] = w1 * (x0 - x1);
	}
      
      fftb2__nest (x + n2 * 0, depth - 1, base_w);
      fftb2__nest (x + n2 * 1, depth - 1, base_w);
    }
}

void
fftb2 (Complex * x, int depth, int f)
{
  Complex base_w = c_arg((Real) f * M_PI / ipow2 (depth));

  /* 入力正順，出力ビット逆順 */
  fftb2__nest (x, depth, base_w);

  if (f < 0)
    fftutl_div_n (x, ipow2 (depth));

  fftutl_bit_inv (x, depth);

  return;
}

/* EOF */