/*------------------------------------------------------------------------------

    File    :	FastFT_float.c

    Author  :	Henryk Richter
                (based on Sources by Stéphane TAVENARD)

    $VER:   FFT_float.c  0.1  (17/09/1997)

    FFT support routines for double precision floating point values

------------------------------------------------------------------------------*/

#define FFT_Buggs

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

#ifdef FFT_Buggs
#include "FastFT_float_asm.h"
#endif

#include "FastFT_float.h"

static WORD INDEX[ FFT_RANGE_MAX ];
static FLOAT SINCOS [ FFT_RANGE_MAX  ]; // FFT_RANGE_MAX / 2 * 2
static FLOAT SINCOS2 [ FFT_RANGE_MAX  ];

static int range = 0;
static int power = 0;

static void build_index( int power )
{
   WORD n, i, j, l, m, x;

   x = power - 1;
   n = 1<<(x);
   for( i=0; i<n; i++ ) {
      INDEX[ i ] = 0;
      l = 1;
      m = n;
      for( j=0; j<x; j++ ) {
	 m = m>>1;
	 if( i & l ) INDEX[ i ] += m;
	 l = l<<1;
      }
   }
}

static void build_sin_cos( int power )
{
   WORD i, n, j;
   double p;

   n = 1<<(power-1);
   p = PI;
   for( i=0; i<n; i++ ) {
      j = INDEX [i];
      SINCOS[ 2*i ]     = (float)sin( 2 * p * (double)j / (double)n);
      SINCOS[ 2*i + 1 ] = (float)cos( 2 * p * (double)j / (double)n);

      SINCOS2[ 2*i ]     = (float)sin( p * (double)i / (double)n);
      SINCOS2[ 2*i + 1 ] = (float)cos( p * (double)i / (double)n);
   }
}

void FastFT_FLOAT_forward( FLOAT *x_real, FLOAT *x_imag,
			FLOAT *energy, FLOAT *phi, int n )
{
   if( n > FFT_RANGE_MAX ) {
      fprintf( stderr, "FFT_DOUBLE_forward: n is out of range (%ld>%ld)\n",
		       n, FFT_RANGE_MAX );
      return;
   }
   if( n != range ) {
      range = 1;
      power = 0;
      while( range < n ) {
	 range = range << 1;
	 power++;
      }
      if( range != n ) {
	 range = 0;
	 power = 0;
	 fprintf( stderr, "FFT_DOUBLE_forward: n is not a power of 2 (%ld)\n", n );
	 return;
      }
      build_index( power );
      build_sin_cos( power );
   }

   ASM_FastFT_main_loop( x_real, x_imag, SINCOS, SINCOS2, power );
   ASM_FastFT_energy_phi( x_real, x_imag, energy, phi, n );
}