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



fftw_complex *init ( int n )
{
  int i;
  fftw_complex *in;
  in = fftw_malloc ( sizeof ( fftw_complex ) * n );
  
   for ( i = 0; i < n; i++ )
   { 
      if (i >= (n/4-1) && i <= (3*n/4-1)) 
        {
          in[i][0] = 1.;
          in[i][1] = 0.;
        } 
      else
        {
          in[i][0] = 0.;
          in[i][1] = 0.;
        }       
   } 
  return in;
}


void print_cplx (fftw_complex x[], int n)
{
  int i;
  for ( i = 0; i < n; i++ ) 
  {
      printf( "i = %d ;  x[i] = (%f + %fi); \n", i, x[i][0], x[i][1]);
  }
}

void print_results (fftw_complex x[], fftw_complex y[], fftw_complex z[], int n)
{
  int i;
  for ( i = 0; i < n; i++ ) 
  {
      printf( "i = %d ;  in[i] = (%f + %fi); out[i] = (%f + %fi); newout[i] = (%f + %fi);\n", i, x[i][0], x[i][1], y[i][0], y[i][1], z[i][0], z[i][1]);
  }
}


void norm_cplx (fftw_complex x[], int n)
{
  int i;
  for ( i = 0; i < n; i++ ) 
  {
     x[i][0] = x[i][0]/n;
     x[i][1] = x[i][1]/n;
  }
}


void check_results(fftw_complex x[], fftw_complex y[], int n)
{
   int i;
   double diffr, diffi, modulex;
   double prec = 1.0e-7;
   double cdiff;
   double maxdiff = 0.;
   double maxmod = FLT_MIN;
   
   for( i = 0; i < n; i++ ) 
   { 
      modulex = sqrt(pow(x[i][0],2.)+pow(x[i][1],2.));
      if (maxmod < modulex) 
      {
	maxmod = modulex;
      }   
   }  
   
   for( i = 0; i < n; i++ ) 
   {  
      cdiff = 0.;
      diffr = fabs(x[i][0] - y[i][0])/maxmod;
      diffi = fabs(x[i][1] - y[i][1])/maxmod;
      if ( cdiff < diffr) 
      {
         cdiff = diffr;
      }
      if( cdiff < diffi) 
      {
         cdiff = diffi;
      }  
      
      if (maxdiff < cdiff )
      {
	maxdiff = cdiff;
      }
   }
   printf("******************************************\n");
   printf("max diff = %E\n",cdiff);
   if(maxdiff > prec)
   {
      printf("Results are incorrect! Error Max = %E \n", maxdiff);
   }
   else 
   {  
      printf("Results are correct! Error Max= %E \n", maxdiff);
   }   
   printf("******************************************\n");
}

int main ( void )

{
/* Declare variabiles */  
  const ptrdiff_t n = 64;
  fftw_complex *in;
  fftw_complex *out;
  fftw_complex *newout;
  fftw_plan plan_backward;
  fftw_plan plan_forward;
/* Create arrays. */
  in = fftw_malloc ( sizeof ( fftw_complex ) * n );
  out = fftw_malloc ( sizeof ( fftw_complex ) * n );
  newout = fftw_malloc ( sizeof ( fftw_complex ) * n );  

/* Initialize data */

  in = init(n);

  
/* Create plans. */
  plan_forward = fftw_plan_dft_1d ( n, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
  plan_backward = fftw_plan_dft_1d ( n, out, newout, FFTW_BACKWARD, FFTW_ESTIMATE );

/* Print initial solution */
  print_cplx(in, n);

/* Compute transform (as many times as desired) */
  fftw_execute ( plan_forward );
  
/* Compute anti-transform */
  fftw_execute ( plan_backward );

/* Normalization */
  norm_cplx(newout, n);
  
  
/* Print results */
  printf("******************************************\n");
  print_results(in, out, newout, n);

/* Check results */  
  printf("******************************************\n");
  check_results(in, newout, n);

/* deallocate and destroy plans */
  fftw_destroy_plan ( plan_forward );
  fftw_destroy_plan ( plan_backward );
  fftw_free ( in );
  fftw_free ( newout );
  fftw_free ( out );


 
}