E.2 プログラムリスト : fft.c

#include <math.h>
#include "complex.h"

/* 高速フーリエ変換 : fft.c */
void fft(complex *f, complex *ww, int *jun, int m, int ift)
{
  complex *F, *WW, c, w, ctmp;
  int *JUN, i, j, k, nnn, nh, jj, ia, ib, ic, n, nn;
  int mm, kban, kzou, ja;
  double ang, a, b;

  F = f;     --F;
  WW = ww;   --WW;
  JUN = jun; --JUN;
  n = pow(2,m);
  nn = n / 2;
  ang = 2. * 3.14159265 / (double) n;
  if(ift == 0) w = cmplx(cos(ang), -sin(ang));
  else         w = cmplx(cos(ang),  sin(ang));
  nnn = n / 4;
  WW[1] = cmplx(1.0 , 0.0);
  for( i = 1 ; i <=nnn ; ++i ){
    WW[i+1] =  cmul( &w, &WW[i] );
    a = real( &WW[i+1] );
    b = aimag( &WW[i+1] );
    WW[nn+1-i] = cmplx(-a, b);
  }
  for( j = 1 ; j <= m ; j++ ){
    nh = n / ( 1 << j );              /* ２のj乗 */
    jj = 1 << ( j-1 );                /* ２のj-1乗 */
    ia = 0;
    ib = nh;
    for( k = 1 ; k <= jj ; ++k ){
      for( i = 1 ; i <= nh ; ++i ){
        c = cadd( &F[i+ia], &F[i+ib] );
        ic = jj * ( i-1 ) + 1;
        ctmp = csub( &F[i+ia], &F[i+ib]);
        F[i+ib] = cmul( &ctmp, &WW[ic] );
        F[i+ia] = c;
      }
      ia = ia + 2 * nh;
      ib = ib + 2 * nh;
    }
  }
  if( ift != 0 ){
    ctmp = cmplx( (double) n, 0.0 );
    for( i = 1 ; i <= n; ++i )
      F[i] = cdiv( &F[i], &ctmp );
  }
  JUN[1] = 1;
  mm = m - 1;
  for( i = 0 ; i <= mm ; i++ ){
    kban = 1 << i;                    /* ２のi乗 */
    kzou = 1 << ( mm - i );           /* ２のmm-1乗 */
    for( j = 1; j <= kban; ++j ){
      ja = kban + j;
      JUN[ja] = JUN[j] + kzou;
    }
  }
  for( i = 2 ; i <= n-1 ; ++i ){
    if( i >= JUN[i] ) continue;
    w = F[ JUN[i] ];
    F[ JUN[i] ] = F[i];
    F[i] = w;
  }
  return;
}