//N*Nの配列 (N=2^n)　g_array[N][N] を2次元フーリエ変換したいとき（g_array[][]は用意していません）
//行と列の数が異なる場合はもう少し面倒くさいです。

/*
#include <stdio.h>
#include <stdlib.h>
#include "fft2dc.h"
のように読み込んでやる。
*/

void inputdata2d(int n1, int n2, double **a)
{
    int j1, j2;

    for (j1 = 0; j1 < n1; j1++) {
        for (j2 = 0; j2 < n2; j2++) {
            a[j1][2*j2] = g_array[j1][j2];//実部
			a[j1][2*j2+1] = 0;//虚部
        }
    }
}

void outputdata2d(int n1, int n2, double **a)
{
	    int j1, j2;

    for (j1 = 0; j1 < n1; j1++) {
        for (j2 = 0; j2 < n2; j2++) {
			g_array[j1][j2] = a[j1][2*j2];
		}
	}
}

int main(void)
{
	int *ip,  m, i;
	int j1, j2 ;
	double **a, *w;

	a = alloc_2d_double(N, 2*N);

    ip = alloc_1d_int(2 + (int) sqrt(N + 0.5));
    m = N/2 + 2*N;
    w = alloc_1d_double(m);
    ip[0] = 0;
    

	inputdata2d(N, N, a);//g_array[][] を複素数の配列 a[][] に拡張　サイズN*N→N*2N

	cdft2d(N, 2*N, 1, a, ip, w);//a[j1][j2]をFourier変換 → a[k1][k2]

	
    cdft2d(N, 2*N, -1, a, ip, w);//逆フーリエ変換　a[k1][k2] → a[j1][j2]　開始

	for (j1 = 0; j1 < N; j1++) {
		for (j2 = 0; j2 < 2 * N; j2++) {
			a[j1][j2] *= 1.0 / (N * N);
		}
	}//逆フーリエ変換終了
	
	outputdata2d(N, N, a );

    free_1d_double(w);
    free_1d_int(ip);
    free_2d_double(a);//作った配列は消します
}