/*  FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X
    (A COMPLEX ARRAY OF DIMENSION N). SOURCE: Ferziger; Numerical
        methods for engineering applications.
 
    X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM.
    N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<4294967296).
    K = 1 FOR FORWARD TRANSFORM.
    K = -1 FOR INVERSE TRANSFORM. */

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

static int log2n(int n){
    int m=1,y=0;
    while(m<n) m+=m,y++;
    if(m!=n) {
        fprintf(stderr,"N=%d is not a power of 2.\n",n);
        exit(1);
    }
    return y;
}
static int l2n;
static double gain;
static complex double u[32];
fft(complex double x[],int n,int k){
    l2n=log2n(n);
    gain=1.0/n;
    u[0]=cexp(-k*2.0*M_PI*I*gain);
    int i;
    for(i=1;i<l2n;i++) u[i]=u[i-1]*u[i-1];
    int sby2=n;
    int stage;
    for(stage=0;stage<l2n;stage++){
        complex double w=1.0;
        complex double v=u[stage];
        int s=sby2;
        sby2/=2;
        int l;
        for(l=0;l<sby2;l++){
            int i;
            for(i=0;i<n;i+=s){
                int p=i+l;
                int q=p+sby2;
                complex double t=x[p]+x[q];
                x[q]=(x[p]-x[q])*w;
                x[p]=t;
            }
            w*=v;
        }
    }
    for(i=0;i<n;i++){
        int index=i;
        int jndex=0;
        int j;
        for(j=0;j<l2n;j++){
            jndex+=jndex;
            int itemp=index/2;
            if(itemp+itemp != index) jndex++;
            index=itemp;
        }
        j=jndex;
        if(j>=i){
            complex double xtemp=x[j];
            x[j]=x[i];
            x[i]=xtemp;
        }
    }
    if(k<0){
        int i;
        for(i=0;i<n;i++) x[i]*=gain;
    }
}