#include <math.h>
#include <malloc.h>
#include <string.h>

#include "coxeter.h"

#define LOOP(i,n) for(int i = 0; i < (n); i++)
#define MIN(x,y) ((x) > (y) ? (y) : (x))
#define MAX(x,y) ((x) > (y) ? (x) : (y))

group_t *coxeter_init_triangle(int p, int q, int r, int nmax)
{
	int coxeter_matrix[9] = {
		1, r, q,
		r, 1, p,
		q, p, 1};

	return coxeter_init(3, coxeter_matrix, nmax);
}

group_t *coxeter_init(int rank, int *coxeter_matrix, int nmax)
{
	int n;
	group_t *group;
	double *schlaefli_matrix;
	double *vectors;
	double *cur_vec;
	groupelement_t *cur, *cur_inv;
	int word[100];  // reasonable estimate?

	// allocate
	group = malloc(sizeof(group_t));
	group->coxeter_matrix = malloc(rank*rank*sizeof(int));
	schlaefli_matrix = malloc(rank*rank*sizeof(double));
	vectors = malloc(rank*nmax*sizeof(double));
	cur_vec = malloc(rank*sizeof(double));
	group->elements = malloc(nmax*sizeof(groupelement_t));
	group->lists = malloc(2*nmax*rank*sizeof(groupelement_t*));
	memset(group->lists, 0, 2*nmax*rank*sizeof(groupelement_t*));
	LOOP(i, nmax) group->elements[i].left = group->lists + i*rank;
	LOOP(i, nmax) group->elements[i].right = group->lists + (nmax+i)*rank;
	LOOP(i, nmax) group->elements[i].letter = -1;
	LOOP(i, nmax) group->elements[i].id = i;

	// copy coxeter matrix
	group->rank = rank;
	memcpy(group->coxeter_matrix, coxeter_matrix, rank*rank*sizeof(int));

	// generate Schläfli matrix
	LOOP(i,rank) LOOP(j,rank) {
		if(group->coxeter_matrix[i*rank+j] == -1)
			schlaefli_matrix[i*rank+j] = -2;
		else
			schlaefli_matrix[i*rank+j] = -2*cos(M_PI/group->coxeter_matrix[i*rank+j]);
	}

	// identity element
	group->elements[0].length = 0;
	group->elements[0].letter = -1;
	LOOP(i, rank) vectors[i] = 1.0;
	n = 1;

	// elements and left multiplication
	for(int i = 0; n < nmax && i < n; i++) {
		LOOP(j, rank) {
			if(n >= nmax)
				break;
			if(vectors[i*rank+j] < 0) // this generator decreases length
				continue;
			if(group->elements[i].left[j] != 0) // we already added this element
				continue;

			LOOP(k, rank) {
				if(k == j)
					vectors[n*rank+k] = -vectors[i*rank+k];
				else
					vectors[n*rank+k] = vectors[i*rank+k] - vectors[i*rank+j]*schlaefli_matrix[k*rank+j];
			}

			group->elements[n].length = group->elements[i].length + 1;

			// if s_k * w is shorter than w (the new element), find out what it is and update "left"
			LOOP(k, rank) {
				if(vectors[n*rank+k] > 0)
					continue;

				if(group->elements[n].letter == -1)
					group->elements[n].letter = k;

				// get w
				LOOP(l, rank) cur_vec[l] = vectors[n*rank+l];

				// apply s_k
				LOOP(l, rank)
					if(l != k)
						cur_vec[l] -= cur_vec[k]*schlaefli_matrix[l*rank+k];
				cur_vec[k] = -cur_vec[k];

				// find a reduced word for s_k w
				LOOP(m, group->elements[i].length) { // s_k w should have same length as element i
					int p;

					// find a generator s_p decreasing the word length
					for(p = 0; p < rank; p++)
						if(cur_vec[p] < 0)
							break;

					if(p == rank) // this should not happen
						fprintf(stderr, "Uh oh!\n");

					word[m] = p;

					// apply s_p
					LOOP(l, rank)
						if(l != p)
							cur_vec[l] -= cur_vec[p]*schlaefli_matrix[l*rank+p];
					cur_vec[p] = -cur_vec[p];
				}

				// find the element corresponding to the word
				groupelement_t *cur = &group->elements[0];
				LOOP(m, group->elements[i].length) {
					cur = cur->left[word[group->elements[i].length - 1 - m]];
				}

				cur->left[k] = &group->elements[n];
				group->elements[n].left[k] = cur;
			}

			n++;
		}
	}

	group->size = n;

	// parent
	LOOP(i, n) {
		if(i == 0)
			group->elements[i].parent = NULL;
		else
			group->elements[i].parent = group->elements[i].left[group->elements[i].letter];
	}

	// right multiplication
	LOOP(i, n) {
		LOOP(j, rank) {
			if(group->elements[i].right[j])
				continue;

			int k = 0;
			for(groupelement_t *cur = &group->elements[i]; cur->parent; cur = cur->parent)
				word[k++] = cur->letter;
			cur = group->elements[0].left[j];
			for(int k = group->elements[i].length - 1; k >= 0; k--)
				if(cur)
					cur = cur->left[word[k]];

			if(cur) {
				group->elements[i].right[j] = cur;
				cur->right[j] = &group->elements[i];
			}
		}
	}

	// inverse
	LOOP(i, n) {
		cur_inv = &group->elements[0];
		for(groupelement_t *cur = &group->elements[i]; cur->parent; cur = cur->parent) {
			if(cur_inv == NULL)
				break;
			cur_inv = cur_inv->left[cur->letter];
		}
		group->elements[i].inverse = cur_inv;
	}

	// free
	free(schlaefli_matrix);
	free(vectors);

	return group;
}

void coxeter_clear(group_t *g)
{
	free(g->coxeter_matrix);
	free(g->elements);
	free(g->lists);
	free(g);
}

int coxeter_snprint(char *str, int buflen, groupelement_t *g)
{
	int n = MIN(g->length, buflen-1); // number of characters without null byte

	str[n] = 0;
	int i = n-1;
	while(g->parent) {
		str[i--] = 'a' + g->letter;
		g = g->parent;
	}

	return n;
}