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

//==============================================================================
//	1次元FFT
//		引数:複素数, サイズ, 指数, フラグ[0:FFT / 1:逆FFT]
//==============================================================================

void fft1(complex c[], int n, int iter, int flag)
{
	int i, it, j, j1, j2, k, p, q, xp, xp2;
	double arg, sign, w;
	complex t, ww;

	// 指定サイズチェック
	if (n < 2) {
		fprintf(stderr, "Error : n < 2  in fft1()\n");
		return;
	}
	if (iter <= 0) {
		iter = 0;
		i = n;
		while((i /= 2) != 0) iter++;
	}
	j = 1;
	for (i = 0; i < iter; i++) j *= 2;
	if (n != j) {
		fprintf(stderr, "Error : n != 2 ^ k  in fft1()\n");
		return;
	}

	sign = (flag == 1)? 1.: -1.;
	xp2 = n;
	for (it = 0; it < iter; it++) {
		xp = xp2;
		xp2 = xp / 2;
		w = M_PI / xp2;
		for (k = 0; k < xp2; k++) {
			arg = k * w;
			ww = cos(arg) + I * sign * sin(arg);
			i = k - xp;
			for (j = xp, p = i+xp, q = i+xp+xp2;
			     j <= n; j += xp, p += xp, q += xp) {
				t = c[p] - c[q];
				c[p] = c[p] + c[q];
				c[q] = t * ww;
			}
		}
	}
	j1 = n / 2;
	j2 = n - 1;
	j = 1;
	for (i = 1, p = 0; i <= j2; i++, p++) {
		if (i < j) {
			t = c[p];
			c[p] = c[j-1];
			c[j-1] = t;
		}
		k = j1;
		while (k < j) {
			j -= k;
			k /= 2;
		}
		j += k;
	}
	if (flag == 0) return;
	w = n;
	for (i = 0, p = 0; i < n; i++, p++) c[p] = c[p] / w;
	return;
}


//==============================================================================
//	2次元FFT
//		引数:複素数, サイズ(横), サイズ(縦), フラグ[0:FFT / 1:逆FFT]
//==============================================================================

void fft2(complex c[], int m, int n, int flag)
{
	int i, iter, j, k, p;
	complex *w;

	// サイズチェック
	if (m < 2) {
		fprintf(stderr, "Error : Illegal parameter in fft2()\n");
		return;
	}
	iter = 0;
	i = m;
	while ((i /= 2) != 0) iter++;
	j = 1;
	for (i = 1; i <= iter; i++) j *= 2;
	if (m != j) {
		fprintf(stderr, "Error : m != 2 ^ k  in fft2()\n");
		return;
	}

	w = (complex *)malloc(m * sizeof(complex));
	if (!w) {
		fprintf(stderr, "Error : Out of memory  in fft2()\n");
		return;
	}

	for (j = 0; j < n; j++, k += n) {
		for (i = 0, p = j; i < m; i++, p += n) w[i] = c[p];
		fft1(w, n, iter, flag);
		for (i = 0, p = j; i < m; i++, p += n) c[p] = w[i];
	}
	for (i = k = 0; i < m; i++, k += n) {
		for (j = 0, p = k; j < m; j++) w[j] = c[p++];
		fft1(w, n, iter, flag);
		for (j = 0, p = k; j < m; j++) c[p++] = w[j];
	}

	free(w);
}