2D DFT(IDFT) source code in C99

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

#ifndef M_PI
#define M_PI      3.14159265358979323846
#endif
#define Rotor(a)  (double _Complex)( cos(a) + I*sin(a) )
#define cZERO     (double _Complex)( 0. + I*0. )

enum { DFT, IDFT };

void dft_TwoDimensional( double _Complex **d, int lx, int ly, int flag ) {
int     x, y, u, v, elements;
double  unit, theta, basis = (flag==DFT) ? -2*M_PI : 2*M_PI;
double _Complex temp[(lx>ly?lx:ly)]; // 動的作業用バッファー領域 毎回初期化する事

    ////////////////////////////////////////////////////////////////////////////////////
    ///   .........  START, Two Dimensional DFT,IDFT Operation .........
    ///  処理手順： 作業エリア確保 => 畳み込み演算 => 作業エリアから結果の書き戻し
    ///
    for( elements=lx, unit=basis/elements, y=0; y<ly; y++ ) {
        for( u=0; u<elements; u++ ) { theta=unit*u; // Euler's formula: eｉθ = cosθ +ｉsinθ
          //for( temp[u]=cZERO, v=0; v<elements; v++ ) temp[u] += d[y][v] *  cexp(theta*v);
            for( temp[u]=cZERO, v=0; v<elements; v++ ) temp[u] += d[y][v] * Rotor(theta*v);
        }
        if(flag==DFT) for( u=0; u<elements; u++ ) d[y][u] = temp[u];
        else          for( u=0; u<elements; u++ ) d[y][u] = temp[u] / elements;
    }// ............. End of Y direction, Σ dft_oneDimensional( X elements ) //

    for( elements=ly, unit=basis/elements, x=0; x<lx; x++ ) {
        for( u=0; u<elements; u++ ) { theta=unit*u; // Euler's formula: eｉθ = cosθ +ｉsinθ
          //for( temp[u]=cZERO, v=0; v<elements; v++ ) temp[u] += d[v][x] *  cexp(theta*v);
            for( temp[u]=cZERO, v=0; v<elements; v++ ) temp[u] += d[v][x] * Rotor(theta*v);
        }
        if(flag==DFT) for( u=0; u<elements; u++ ) d[u][x] = temp[u];
        else          for( u=0; u<elements; u++ ) d[u][x] = temp[u] / elements;
    }// ............. End of X direction, Σ dft_oneDimensional( Y elements ) //
}