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


/* Killer Radix-2 DIF FFT, IFFT and Postscrambler by David J. Savinkoff	*/
/* 2002   */
/* To compile with UNIX  cc multiply.c -lm                              */
/* To compile with Linux  gcc -W -Wall -lm -o multiply multiply.c       */
#define SIZE 256	/* 256 point FFT. */
#define L 8		/* SIZE = 2^L	  */
#define V (2*M_PI/SIZE)

/*	  Using double instead of float may yield little improvement.	*/

float		real_[SIZE],	imag_[SIZE];

/* Note that this is actually an unscaled IFFT because I prefer results */
/* consistent with the Fourier Series, not The Laplace Transform.       */
void fft(void)
{
int	aye, bee, scramble, mask,
	i1, i2, i3, i4,
	i, k, m, s;

float	temp, b1, b2, csin_, sine_, x;

	/* FFT in-place algorithm					*/

	i1 = SIZE>>1;	i2 = 1;

	for( i=L; i--; )
	  {   i3 = 0;	i4 = i1;	scramble = 0;
	      for( k=i2; k--; )
		{   x  =  (float)scramble*V;
		    csin_ =  cos(x);	sine_ = sin(x);
		    /* Ripple Counter with MSB swapped for LSB  etc..	*/
		    mask = SIZE>>1;
		    do
		     {    bee = scramble;
			  scramble ^=  mask >>= 1;
		     }while( scramble < bee);

		    for( m=i3; m<i4; m++ )
		      {   s  = m+i1;
			  b1 = csin_*real_[s] - sine_*imag_[s];
			  b2 = csin_*imag_[s] + sine_*real_[s];
			  real_[s] =  real_[m]-b1; imag_[s] =  imag_[m]-b2;
			  real_[m] += b1;	   imag_[m] += b2;
		      }
		    i3 += i1<<1;	i4 += i1<<1;
		}
	      i1 >>= 1;		i2 <<= 1;
	  }
	/* Complex in-place postscrambler				*/
	for(aye=scramble=0; aye<SIZE; aye++)
	  {  if(scramble>aye)
	      {  temp=imag_[aye]; imag_[aye]=imag_[scramble];
			imag_[scramble]=temp;
		 temp=real_[aye]; real_[aye]=real_[scramble];
			real_[scramble]=temp;
	      }
	      /* Ripple Counter with MSB swapped for LSB  etc..		*/
	      mask = SIZE;
	      do
	      {	 bee = scramble;
		 scramble ^=  mask >>= 1;
	      }while( scramble < bee);
	  }
}

/* IFFT differs in scaling and -sin(x)					*/
/* Note that this is actually a scaled FFT because I prefer results     */
/* consistent with the Fourier Series, not The Laplace Transform.       */
void ifft(void)
{
int	aye, bee, scramble, mask,
	i1, i2, i3, i4,
	i, k, m, s;

float	temp, b1, b2, csin_, sine_, x;

	/* IFFT in-place algorithm					*/

	i1 = SIZE>>1;	i2 = 1;

	for( i=L; i--; )
	  {   i3 = 0;	i4 = i1;	scramble = 0;
	      for( k=i2; k--; )
		{   x  =  (float)scramble*V;
		    csin_ =  cos(x);	sine_ = sin(x);
		    /* Ripple Counter with MSB swapped for LSB  etc..	*/
		    mask = SIZE>>1;
		    do
		     {    bee = scramble;
			  scramble ^=  mask >>= 1;
		     }while( scramble < bee);

		    for( m=i3; m<i4; m++ )
		      {   s  = m+i1;
			  b1 = csin_*real_[s] + sine_*imag_[s];
			  b2 = csin_*imag_[s] - sine_*real_[s];
			  real_[s] =  real_[m]-b1; imag_[s] =  imag_[m]-b2;
			  real_[m] += b1;	   imag_[m] += b2;
		      }
		    i3 += i1<<1;	i4 += i1<<1;
		}
	      i1 >>= 1;		i2 <<= 1;
	  }
	/* Complex in-place postscrambler				*/
	for(aye=scramble=0; aye<SIZE; aye++)
	  {  if(scramble>aye)
	      {  temp=real_[aye]; real_[aye]=real_[scramble];
			real_[scramble]=temp;
		 temp=imag_[aye]; imag_[aye]=imag_[scramble];
			imag_[scramble]=temp;
	      }
	      /* Ripple Counter with MSB swapped for LSB  etc..		*/
	      mask = SIZE;
	      do
	      {	 bee = scramble;
		 scramble ^=  mask >>= 1;
	      }while( scramble < bee);
	  }
	/* Scale Complex output time function				*/
	for( aye=SIZE; aye--; )
	  {   real_[aye] /= SIZE;
	      imag_[aye] /= SIZE;
	  }
}


char	buff[SIZE+2];

void enter(void)
{
char  *	ascii_pointer;
float *	float_pointer;

	printf("\nEnter number A: ");
	fgets(buff, SIZE+2, stdin);
        for(ascii_pointer=buff; *ascii_pointer; ascii_pointer++)
		;
	*--ascii_pointer = 0;
        for(float_pointer = real_; ascii_pointer > buff; )
	{	*float_pointer++ = atof(--ascii_pointer);
		*ascii_pointer = 0;
	}

	printf("\rEnter number B: ");
	fgets(buff, SIZE+2, stdin);
        for(ascii_pointer=buff; *ascii_pointer; ascii_pointer++)
		;
	*--ascii_pointer = 0;
        for(float_pointer = imag_; ascii_pointer > buff; )
	{	*float_pointer++ = atof(--ascii_pointer);
		*ascii_pointer = 0;
	}
}

void dual_convolute(void)
{
int     k, retro_k;
float   temp1, temp2;

  for( k = retro_k = SIZE/2; k; )
  {
    /* Calculate Real stuff */
    temp1 = (imag_[k]+imag_[retro_k])*(real_[k]+real_[retro_k]) +
            (real_[k]-real_[retro_k])*(imag_[k]-imag_[retro_k]);

    /* Calculate Imaginary stuff */
    temp2 = (imag_[k]+imag_[retro_k])*(imag_[k]-imag_[retro_k]) -
            (real_[k]+real_[retro_k])*(real_[k]-real_[retro_k]);

    imag_[k] = temp2 * 0.25;
    imag_[retro_k] = temp2 * -0.25;

    /* Assign Real stuff */
    real_[k--] = real_[retro_k++] = temp1 * 0.25;
  }
  /* Calculate more Real stuff */
  real_[0] *= imag_[0];
  /* Calculate more Imaginary stuff */
  imag_[0] = 0.0;
}
void complex_display(void)
{
int     n;

  for(n=0; n<20; n++)
   printf("( %g + j%g ),", real_[n], imag_[n]);

}
void convolution_display(void)
{
int     n;

  for(n=0; n<SIZE; n++)
   printf("%d,", (int)(real_[n]+0.5));

}
void multiplication_display(void)
{
int     n;
float * summation;
float * digit_out;
float   carry_cell;

summation = real_;                  digit_out = imag_;

carry_cell = *summation;

  for(n=SIZE; n--; )
  {
    *digit_out++ = (float)((long)(carry_cell + 0.5) % 10);
    carry_cell   = (float)((long)(carry_cell + 0.5) / 10) + *++summation;
  }
  for(n=SIZE; n--; )
   printf("%d", (int)(imag_[n]) );
}
int main(void)
{

        enter();
/*	printf("\nEnter result\n");  */
/*	complex_display();  */

	fft();
/*	printf("\nFFT result\n");  */
/*	complex_display();  */

	dual_convolute();

	ifft();

	printf("\nMultiplication result:\n");
	multiplication_display();
	printf("\nConvolution result:\n");
        convolution_display();

	getchar();			/*	wait for user		*/
	return 1;
}