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



void FFT(double complex *l,double complex *r,int p)
 { int N,i1; double complex *w,r1inter,r2inter;
     N=pow(2,p);
     w=(double complex *)malloc(N*sizeof(double complex));
     w[0]=1+0*I;
     /*mise en place de la liste des puissances de w*/
     double complex c;
     c=cos(2.*M_PI/N)-sin(2.*M_PI/N)*I;
     for(i1=0;i1<=(N-2);i1++)
     {w[i1+1]=w[i1]*c;
      
      }
      
     
     
    int k,w1,d,k1,z,s;
      
      
   /*inversions des bits*/   
   for(k=0;k<=(N-1) ;k++)
     {w1=0; z=k; s=0; d=N/2;
                       
      for(k1=0; k1<=(p-1);k1++)
         {
                 
           if (z==0){ break; }
                         
            else if(z!=0)
             {w1=z%2; s=s+w1*d;
              d=d/2;
             }
              z=(z-w1)/2;
     
          }
      
      r[k]=l[s];
      
      
      }

     
     /*la partie calcul proprement dite*/
   {    int d1,d2,i,j,b,f; 
        d1=N/2; d2=1;
        for (j=1; j<=p; j++)
            {
             for(i=0; i<=(d1-1); i++)
               {
                for(b=0; b<=(d2-1); b++)
  
                 {
                         
                  r1inter=r[i*2*d2+b];
                  r2inter=r[i*2*d2+b+d2];       
                  r[i*2*d2+b]=(r1inter+r2inter*w[b])*0.5; 
                  r[i*2*d2+b+d2]=(r1inter-r2inter*w[b])*0.5; 
                  
                  }
                  
                 
                }
                
	    d1=d1/2; d2=d2*2; 
        
        
             }

                     
    }            
free(w);           
}           


int main(void)
{  
int lambda;          
double complex *testprime;
double complex *result1prime;


testprime=(double complex *)malloc(8*sizeof(double complex));
result1prime=(double complex *)malloc(8*sizeof(double complex));


testprime[0]=40.;
testprime[1]=23.;
testprime[2]=1.;
testprime[3]=54.;
testprime[4]=1.;
testprime[5]=8.;   
testprime[6]=1.;
testprime[7]=4.;

       
FFT(testprime,result1prime,3);

for(lambda=0; lambda<=7; lambda++)
{printf("%lg + %lg*i \n",creal(result1prime[lambda]),cimag(result1prime[lambda]));}



free(testprime);
free(result1prime);
free(result2prime);                    
return(0);            
}