/*-------------------------------------------------------------*/
/* Cut this and save to fftlib.c			       */
/* Program by Gregory F. Shay  1988			*/
#include <math.h>
#include <stdio.h>

#define MAXN 2048

float wr[MAXN],wi[MAXN],xinr[MAXN],xini[MAXN],xoutr[MAXN],xouti[MAXN];
int bits[MAXN];
int LOGN,N;

fft_main(forward)
	int forward;
{
	int i,j,k,flag;
	char dtype;

flag = 0;
/* input size of array */
scanf("%c\n",&dtype);
if ((dtype != 'r') && (dtype != 'c'))
	{
	fprintf(stderr,"Incompatible input data file, expect real or complex\n");
	flag = 1;
	}
scanf("%d \n",&LOGN);
if ((LOGN<2) || (LOGN>14) )
	{
	flag = 1;
	fprintf(stderr,"Error, size of array, %d, out of range.\n",LOGN);
	}

N = 1 << LOGN;
if (N>MAXN) 
	{
	fprintf(stderr,"Error, maximum size of array is set to %d\n",MAXN);
	flag = 1;
	}

if (flag == 0)
  {
/* Initialize omega array and bit reversal array */

initw();
bitrev();

/* Read in data */
if (dtype == 'r')
	for(i=0;i<N;i++)
	{
	scanf("%f ",&xinr[i]);
	xini[i] = 0.;
	}
if (dtype == 'c')
	for(i=0;i<N;i++)
	{
	scanf("%f ",&xinr[i]);
	scanf("%f ",&xini[i]);
	}

/* Perform fft */

fft(xinr,xini,xoutr,xouti,forward);

/* Output data */

printf("c\n");
printf("%d \n",LOGN);
for(i=0;i<N;i++)
	printf("%g %g\n",xoutr[i],xouti[i]);

  } /* end of if flag == 0 clause */
} /* End of main */

initw()
{
	float pi,theta;
	int i,j,k;

pi = 3.14159265;
theta = 2.*pi/((float)N);
for(k=0;k<= N-1;k++)
	{
	wr[k]=cos(theta*k);
	wi[k]=sin(theta*k);
	}
}

bitrev()
{
	int i,j,m;
bits[0]=0;
m=1;
for(i=0;i<=LOGN-1;i++)
	{
	for(j=0;j<=m-1;j++)
		{
		bits[j] <<= 1;
		bits[j+m]=bits[j]+1;
		}
	m <<= 1;
	}
} /* End of bitrev */

fft(xinr,xini,xoutr,xouti,forward)
	float xinr[],xini[],xoutr[],xouti[];
	int forward;
{
	int i,j,k,istart,c,d;
	float ar,ai,br,bi,fn;
	int joff,jpnt,ioff,jwpnt,jbk;

for(i=0;i<=N-1;i++)
	{
	xoutr[i]=xinr[bits[i]];
	xouti[i]=xini[bits[i]];
	}
joff=N/2;
jpnt=N/2;
jbk=2;
ioff=1;
for(i=1;i<=LOGN;i++)
	{
	for(istart = 0;istart<=N-1;istart+=jbk)
		{
		jwpnt=0;
		for(k=istart;k<=istart+ioff-1;k++)
			{
			c=k+ioff;
			d=jwpnt+joff;
			if (forward == 1)
			   {
			   ar=xoutr[k]+xoutr[c]*wr[jwpnt]-xouti[c]*wi[jwpnt];
			   ai=xouti[k]+xoutr[c]*wi[jwpnt]+xouti[c]*wr[jwpnt];
			   br=xoutr[k]+xoutr[c]*wr[d]-xouti[c]*wi[d];
			   bi=xouti[k]+xoutr[c]*wi[d]+xouti[c]*wr[d];
			   }
			else
			   {
			   ar=xoutr[k]+xoutr[c]*wr[jwpnt]+xouti[c]*wi[jwpnt];
			   ai=xouti[k]+xouti[c]*wr[jwpnt]-xoutr[c]*wi[jwpnt];
			   br=xoutr[k]+xoutr[c]*wr[d]+xouti[c]*wi[d];
			   bi=xouti[k]+xouti[c]*wr[d]-xoutr[c]*wi[d];
			   }
			xoutr[k]=ar;
			xouti[k]=ai;
			xoutr[c]=br;
			xouti[c]=bi;
			jwpnt += jpnt;
			} /* End of k loop */
		} /* End of istart loop */
	jpnt >>= 1;
	ioff = jbk;
	jbk <<= 1;
	} /* End of i loop */

if(forward == 1)
	fn = (float)N/2.;
else
	fn = 2.;   

for (i=0;i<=N-1;i++)
	{
	xoutr[i] = xoutr[i]/fn;
	xouti[i] = xouti[i]/fn;
	}
	

} /* End of fft subroutine */