2D FFT (IFFT) source code in C99

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

#ifndef M_PI
#define M_PI      3.14159265358979323846
#endif
#define FFT2K     (2048)
#define Rotor(a)  (double _Complex)( cos(a) + I*sin(a) )

enum { FFT=DFT, IFFT=IDFT };

void FFT2k_TwoDimensionalC99( double _Complex **d, int flag ) { // ∵d[2048][2048]
#define elements  FFT2K
int     x, y, e;    // flag   IFFT(2*M_PI/e..), FFT(-2*M_PI/e..)
double  unit = (flag==IFFT) ? 2*M_PI/elements : -2*M_PI/elements;
double _Complex rotor[elements]; for(e=0;e<elements;e++) rotor[e]=Rotor(unit*e);

    for( y=0; y<elements; y++ ) {  // x軸方向の1次元FFT(IFFT)をy行分行う
        int     i, j, k, m, mh, ir;// ir:irev
        double _Complex  w, t;     // t: temp
        for( i=0, j=1; j<elements-1; j++ ) { // scrambler, for initial placement
            for( k=elements>>1; k>(i^=k); k>>=1 );
            if( j<i ){ t=d[y][j]; d[y][j]=d[y][i]; d[y][i]=t; } //exchange
        }
        for( mh=1; (m=mh<<1)<=elements; mh=m ) {
            for( ir=0, i=0; i<elements; i+=m ) {
                for( w=rotor[ir], k=elements>>2; k>(ir^=k); k>>=1 );//for next ir
                for( j=i; j<mh+i; j++ ) {
                         k  = j + mh;
                         t  = d[y][j] - d[y][k];
                    d[y][j] = d[y][j] + d[y][k];    // a ->- (a + b) ->- a
                    d[y][k] = t * w;                // b ->- (a-b)*w ->- b
        }   }   }
        if(flag==IFFT) for( i=0; i<elements; i++ ) d[y][i] /= elements;
    } // X_END
    for( x=0; x<elements; x++ ) {  // y軸方向の1次元FFT(IFFT)をx列分行う
        int     i, j, k, m, mh, ir;// ir:irev
        double _Complex  w, t;     // t: temp
        for( i=0, j=1; j<elements-1; j++ ) { // scrambler, for initial placement
            for( k=elements>>1; k>(i^=k); k>>=1 );
            if( j<i ){ t=d[j][x]; d[j][x]=d[i][x]; d[i][x]=t; } //exchange
        }
        for( mh=1; (m=mh<<1)<=elements; mh=m ) {
            for( ir=0, i=0; i<elements; i+=m ) {
                for( w=rotor[ir], k=elements>>2; k>(ir^=k); k>>=1 );//for next ir
                for( j=i; j<mh+i; j++ ) {
                         k  = j + mh;
                         t  = d[j][x] - d[k][x];
                    d[j][x] = d[j][x] + d[k][x];    // a ->- (a + b) ->- a
                    d[k][x] = t * w;                // b ->- (a-b)*w ->- b
        }   }   }
        if(flag==IFFT) for( i=0; i<elements; i++ ) d[i][x] /= elements;
    } // Y_END
}


////////////////////////////////////////////////////////////////////////
// How to use FFT2k_TwoDimensionalC99(...)  2019/04/01
//
void dft_TwoDimensional( COMPLEX **s, int lx, int ly, int flag ) {
int     w = lx>ly ? lx : ly;
COMPLEX **spreadsheet = NULL;

    double _Complex **ss = (double _Complex **)s;
    if( lx==2048 && ly==2048 ) return FFT2k_TwoDimensionalC99(ss,flag);
    if( (spreadsheet=(COMPLEX**)alloc2Darray(sizeof(COMPLEX),w,w)) ) {
       dft_TwoDimensionalWithTable( spreadsheet, w, s, lx, ly, flag );
       free(spreadsheet);
    } else dft_TwoDimensionalStandard( s, lx, ly, flag );
}