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

double PI;
typedef double complex cplx;

void _fft(cplx buf[], cplx out[], int n, int step)
{
	if (step < n) {
		_fft(out, buf, n, step * 2);
		_fft(out + step, buf + step, n, step * 2);

		for (int i = 0; i < n; i += 2 * step) {
			cplx t = cexp(-I * PI * i / n) * out[i + step];
			buf[i / 2]     = out[i] + t;
			buf[(i + n)/2] = out[i] - t;
		}
	}
}

void fft(cplx buf[], int n)
{
	cplx out[n];
	for (int i = 0; i < n; i++) out[i] = buf[i];
	_fft(buf, out, n, 1);
}

/* pad array length to power of two */
cplx *pad_two(double g[], int len, int *ns)
{
	int n = 1;
	if (*ns) n = *ns;
	else while (n < len) n *= 2;

	cplx *buf = calloc(sizeof(cplx), n);
	for (int i = 0; i < len; i++) buf[i] = g[i];
	*ns = n;
	return buf;
}

void deconv(double g[], int lg, double f[], int lf, double out[], int row_len) {
	int ns = 0;
	cplx *g2 = pad_two(g, lg, &ns);
	cplx *f2 = pad_two(f, lf, &ns);

	fft(g2, ns);
	fft(f2, ns);

	cplx h[ns];
	for (int i = 0; i < ns; i++) h[i] = g2[i] / f2[i];
	fft(h, ns);

	for (int i = 0; i < ns; i++) {
		if (cabs(creal(h[i])) < 1e-10)
			h[i] = 0;
	}

	for (int i = 0; i > lf - lg - row_len; i--)
		out[-i] = h[(i + ns) % ns]/32;
	free(g2);
	free(f2);
}

double* unpack2(void *m, int rows, int len, int to_len)
{
	double *buf = calloc(sizeof(double), rows * to_len);
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < len; j++)
			buf[i * to_len + j] = ((double(*)[len])m)[i][j];
	return buf;
}

void pack2(double * buf, int rows, int from_len, int to_len, void *out)
{
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < to_len; j++)
			((double(*)[to_len])out)[i][j] = buf[i * from_len + j] / 4;
}

void deconv2(void *g, int row_g, int col_g, void *f, int row_f, int col_f, void *out) {
	double *g2 = unpack2(g, row_g, col_g, col_g);
	double *f2 = unpack2(f, row_f, col_f, col_g);

	double ff[(row_g - row_f + 1) * col_g];
	deconv(g2, row_g * col_g, f2, row_f * col_g, ff, col_g);
	pack2(ff, row_g - row_f + 1, col_g, col_g - col_f + 1, out);

	free(g2);
	free(f2);
}

double* unpack3(void *m, int x, int y, int z, int to_y, int to_z)
{
	double *buf = calloc(sizeof(double), x * to_y * to_z);
	for (int i = 0; i < x; i++)
		for (int j = 0; j < y; j++) {
			for (int k = 0; k < z; k++)
				buf[(i * to_y + j) * to_z + k] =
					((double(*)[y][z])m)[i][j][k];
		}
	return buf;
}

void pack3(double * buf, int x, int y, int z, int to_y, int to_z, void *out)
{
	for (int i = 0; i < x; i++)
		for (int j = 0; j < to_y; j++)
			for (int k = 0; k < to_z; k++)
				((double(*)[to_y][to_z])out)[i][j][k] =
					buf[(i * y + j) * z + k] / 4;
}

void deconv3(void *g, int gx, int gy, int gz, void *f, int fx, int fy, int fz, void *out) {
	double *g2 = unpack3(g, gx, gy, gz, gy, gz);
	double *f2 = unpack3(f, fx, fy, fz, gy, gz);

	double ff[(gx - fx + 1) * gy * gz];
	deconv(g2, gx * gy * gz, f2, fx * gy * gz, ff, gy * gz);
	pack3(ff, gx - fx + 1, gy, gz, gy - fy + 1, gz - fz + 1, out);

	free(g2);
	free(f2);
}

int main()
{
	PI = atan2(1,1) * 4;
	double h[2][3][4] = {
		{{-6, -8, -5,  9}, {-7, 9, -6, -8}, { 2, -7,  9,  8}},
		{{ 7,  4,  4, -6}, { 9, 9,  4, -4}, {-3,  7, -2, -3}}
	};
	int hx = 2, hy = 3, hz = 4;
	double f[3][2][3] = {	{{-9,  5, -8}, { 3,  5,  1}},
				{{-1, -7,  2}, {-5, -6,  6}},
				{{ 8,  5,  8}, {-2, -6, -4}} };
	int fx = 3, fy = 2, fz = 3;
	double g[4][4][6] = {
		{	{ 54,  42,  53, -42,  85, -72}, { 45,-170,  94, -36,  48,  73},
			{-39,  65,-112, -16, -78, -72}, {  6, -11,  -6,  62,  49,   8} },
		{ 	{-57,  49, -23,   52, -135,  66},{-23, 127, -58,   -5, -118,  64},
			{ 87, -16,  121,  23,  -41, -12},{-19,  29,   35,-148,  -11,  45} },
		{	{-55, -147, -146, -31,  55,  60},{-88,  -45,  -28,  46, -26,-144},
			{-12, -107,  -34, 150, 249,  66},{ 11,  -15,  -34,  27, -78, -50} },
		{	{ 56,  67, 108,   4,  2,-48},{ 58,  67,  89,  32, 32, -8},
			{-42, -31,-103, -30,-23, -8},{  6,   4, -26, -10, 26, 12}
		}
	};
	int gx = 4, gy = 4, gz = 6;

	double h2[gx - fx + 1][gy - fy + 1][gz - fz + 1];
	deconv3(g, gx, gy, gz, f, fx, fy, fz, h2);
	printf("deconv3(g, f):\n");
	for (int i = 0; i < gx - fx + 1; i++) {
		for (int j = 0; j < gy - fy + 1; j++) {
			for (int k = 0; k < gz - fz + 1; k++)
				printf("%g ", h2[i][j][k]);
			printf("\n");
		}
		if (i < gx - fx) printf("\n");
	}

	double f2[gx - hx + 1][gy - hy + 1][gz - hz + 1];
	deconv3(g, gx, gy, gz, h, hx, hy, hz, f2);
	printf("\ndeconv3(g, h):\n");
	for (int i = 0; i < gx - hx + 1; i++) {
		for (int j = 0; j < gy - hy + 1; j++) {
			for (int k = 0; k < gz - hz + 1; k++)
				printf("%g ", f2[i][j][k]);
			printf("\n");
		}
		if (i < gx - hx) printf("\n");
	}
}

/* two-D case; since task doesn't require showing it, it's commented out */
/*
int main()
{
	PI = atan2(1,1) * 4;
	double h[][6] = { 	{-8, 1, -7, -2, -9, 4},
				{4, 5, -5, 2, 7, -1},
				{-6, -3, -3, -6, 9, 5} };
	int hr = 3, hc = 6;

	double f[][5] = {	{-5, 2, -2, -6, -7},
				{9, 7, -6, 5, -7},
				{1, -1, 9, 2, -7},
				{5, 9, -9, 2, -5},
				{-8, 5, -2, 8, 5} };
	int fr = 5, fc = 5;
	double g[][10] = {
			{40, -21, 53, 42, 105, 1, 87, 60, 39, -28},
			{-92, -64, 19, -167, -71, -47, 128, -109, 40, -21},
			{58, 85, -93, 37, 101, -14, 5, 37, -76, -56},
			{-90, -135, 60, -125, 68, 53, 223, 4, -36, -48},
			{78, 16, 7, -199, 156, -162, 29, 28, -103, -10},
			{-62, -89, 69, -61, 66, 193, -61, 71, -8, -30},
			{48, -6, 21, -9, -150, -22, -56, 32, 85, 25}	};
	int gr = 7, gc = 10;

	double h2[gr - fr + 1][gc - fc + 1];
	deconv2(g, gr, gc, f, fr, fc, h2);
	for (int i = 0; i < gr - fr + 1; i++) {
		for (int j = 0; j < gc - fc + 1; j++)
			printf(" %g", h2[i][j]);
		printf("\n");
	}

	double f2[gr - hr + 1][gc - hc + 1];
	deconv2(g, gr, gc, h, hr, hc, f2);
	for (int i = 0; i < gr - hr + 1; i++) {
		for (int j = 0; j < gc - hc + 1; j++)
			printf(" %g", f2[i][j]);
		printf("\n");
	}
}*/