#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fft.h"

#define NOERR 0

static siftFlg = 1 ;	/* DC成分の位置  0 シフトしない  1 シフトする */

static double *sinTable ;
static double *cosTable ;
static double *subBuf ;
static tableFlg = 1 ;	/* 0のときはテーブルを作成しない */

int fft_2( char *work, double real[], double imag[], int ex, int inv )
{
	int i, length ;

	length = 1 << ex ;

	makeTable( work, length, inv ) ;
	tableFlg = 0 ;

	for( i=0 ; i<length ; i++ )
	{
		if( inv == -1 )
			sift( &real[i*length], &imag[i*length], length ) ;
		fft_1( work, &real[i*length], &imag[i*length], ex, inv ) ;
		if( inv == 1 )
			sift( &real[i*length], &imag[i*length], length ) ;
	}
	matrixReverce( real, length ) ;
	matrixReverce( imag, length ) ;
	for( i=0 ; i<length ; i++ )
	{
		if( inv == -1 )
			sift( &real[i*length], &imag[i*length], length ) ;
		fft_1( work, &real[i*length], &imag[i*length], ex, inv ) ;
		if( inv == 1 )
			sift( &real[i*length], &imag[i*length], length ) ;
	}
	matrixReverce( real, length ) ;
	matrixReverce( imag, length ) ;

	tableFlg = 1 ;
	return NOERR ;
}

static int sift( double real[], double imag[], int length )
{
	int i ;
	double temp ;

	if( siftFlg == 0 )return NOERR ;

	for( i=0 ; i<(length >> 1) ; i++ )
	{
		temp = real[i] ;
		real[i] = real[ i + (length >> 1) ] ;
		real[ i + (length >> 1) ] = temp ;

		temp = imag[i] ;
		imag[i] = imag[ i + (length >> 1) ] ;
		imag[ i + (length >> 1) ] = temp ;
	}

	return NOERR ;
}

static int matrixReverce( double a[], int length )
{
	int i, j ;
	double temp ;

	for( i=0 ; i<length ; i++ )
	{
		for( j=0 ; j<length ; j++ )
		{
			if( i < j )
			{
				temp = a[ j + i*length ] ;
				a[ j + i*length ] = a[ i + j*length ] ;
				a[ i + j*length ] = temp ;
			}
		}
	}

	return NOERR ;
}

int fft_1( char *work, double real[], double imag[], int ex, int inv )
{
	int i, j, k, w, j1, j2 ;
	int length, length2, count, temp ;
	double ar, ai, br, bi, nrml ;

	length = 1 << ex ;

	if( tableFlg )makeTable( work, length, inv ) ;

	count = 1 ;
	length2 = length ;
	for( i=0 ; i<ex ; i++ )
	{
		length2 >>= 1 ;
		temp = 0 ;
		for( j=0 ; j<count ; j++ )
		{
			w = 0 ;
			for( k=0 ; k<length2 ; k++ )
			{
				j1 = temp + k ;
				j2 = j1 + length2 ;
				ar = real[j1] ;
				ai = imag[j1] ;
				br = real[j2] ;
				bi = imag[j2] ;
				real[j1] = ar + br ;
				imag[j1] = ai + bi ;
				ar = ar - br ;
				ai = ai - bi ;
				real[j2] = ar*cosTable[w] - ai*sinTable[w] ;
				imag[j2] = ar*sinTable[w] + ai*cosTable[w] ;
				w += count ;
			}
			temp += (2*length2) ;
		}
		count <<= 1 ;
	}
	bitReverce( subBuf, real, ex ) ;
	bitReverce( subBuf, imag, ex ) ;
	nrml = 1.0 / sqrt( (double)length ) ;
	for( i=0 ; i<length ; i++ )
	{
		real[i] *= nrml ;
		imag[i] *= nrml ;
	}

	return NOERR ;
}

static int makeTable( char *work, int length, int inv )
{
	int i ;

	sinTable = (double *)work ;
	cosTable = sinTable + length ;
	subBuf   = cosTable + length ;

	if( inv == 1 )
	{
		for( i=0 ; i<length ; i++ )
		{
			sinTable[i] = sin( -_PI * 2 * i / length ) ;
			cosTable[i] = cos( -_PI * 2 * i / length ) ;
		}
	}
	else
	{
		for( i=0 ; i<length ; i++ )
		{
			sinTable[i] = sin( _PI * 2 * i / length ) ;
			cosTable[i] = cos( _PI * 2 * i / length ) ;
		}
	}

	return NOERR ;
}

static bitReverce( double *buf, double a[], int ex )
{
	int i, j, bit, length ;

	length = 1 << ex ;
	for( i=0 ; i<length ; i++ )
	{
		bit = 0 ;
		for( j=0 ; j<ex ; j++ )
		{
			if( i & (1 << j) )
			{
				bit |= (1 << (ex - j - 1)) ;
			}
		}
		buf[i] = a[bit] ;
	}
	for( i=0 ; i<length ; i++ )a[i] = buf[i] ;

	return NOERR ;
}