#ifndef KLEIN_H
#warning Please include klein/klein.h before klein/parm.h
#endif

#ifdef KLEIN_PARM_H
#error file included twice
#endif
#define KLEIN_PARM_H


typedef void (*function_t)(float *, int *, unsigned char *);

struct parm
{
	unsigned char *grid;
	unsigned char m, n;
	function_t f;
};


void klein_parametrize( struct klein * klein, struct parm parm );

#ifdef KLEIN_IMPLEMENT
#ifdef TEST
#include <assert.h>
#endif

static inline uint64_t __factorial(uint64_t n)
{
	if (n == 1)
		return 1;

	return n * __factorial(n - 1);
}

static inline uint64_t __face(int n)
{
	if (n == 2)
		return 1;

	return (1 << (n - 3)) * __factorial(n) / __factorial(n - 2);
}

void klein_parametrize( struct klein * klein, struct parm parm)
{
	unsigned long i, j,  o, p, n;
	uint64_t k, size, q = 0;
	int *face;

#ifdef TEST
	assert(__face(2) == 1);
	assert(__face(3) == 6);
	assert(__face(4) == 24);
#endif

	klein->dim = parm.n;
	klein->vertex_size = 0;
	{
		uint64_t test = 0;

		for (o = 0; o < parm.m; o++)
			for (p = 0; p < o; p++)
				test += 1;


		for (o = 0; o < parm.m; o++)
			for (p = 0; p < o; p++)
				klein->vertex_size += (uint64_t)parm.grid[p] * parm.grid[o] * 6 * __face(parm.m)/test;

	}

	size = klein->vertex_size*klein->dim;

	klein->vertex = malloc(size * sizeof(float));
	face = malloc(parm.m * sizeof(int));


	for (o = 0; o < parm.m; o++)
	{
		for (p = 0; p < o; p++)
		{
			for (k = 0; k < ((uint64_t)1 << (parm.m - 2)); k++)
			{
				unsigned char skip = 0;
				for (n = 0; n < parm.m; n++)
				{
					if (n == o || n == p)
						skip++;

					face[n] = (k & ((uint64_t)1 << (n - skip))) ? parm.grid[n] : 0;
				}

				for (i = 0; i < parm.grid[p]; i++)
				{
					for (j = 0; j < parm.grid[o]; j++)
					{
						face[p] = i;
						face[o] = j;
						parm.f(&klein->vertex[q], face, parm.grid);
						q += parm.n;

						face[p] = i + 1;
						face[o] = j;
						parm.f(&klein->vertex[q], face, parm.grid);
						q += parm.n;

						face[p] = i + 1;
						face[o] = j + 1;
						parm.f(&klein->vertex[q], face, parm.grid);
						q += parm.n;

						face[p] = i;
						face[o] = j;
						parm.f(&klein->vertex[q], face, parm.grid);
						q += parm.n;

						face[p] = i;
						face[o] = j + 1;
						parm.f(&klein->vertex[q], face, parm.grid);
						q += parm.n;

						face[p] = i + 1;
						face[o] = j + 1;
						parm.f(&klein->vertex[q], face, parm.grid);
						q += parm.n;
					}
				}
			}
		}
	}

#ifdef TEST
	assert(q == size);
#endif
}
#endif