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

#include "triangle.h"
#include "linalg.h"

#define POINCARE 1
#define LOOP(i) for(int i = 0; i < 3; i++)

void cartan_matrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
{
	gsl_matrix_set(cartan, 0, 0, -2);
	gsl_matrix_set(cartan, 0, 1, 2*s*cos(a3));
	gsl_matrix_set(cartan, 0, 2, 2/s*cos(a2));

	gsl_matrix_set(cartan, 1, 0, 2/s*cos(a3));
	gsl_matrix_set(cartan, 1, 1, -2);
	gsl_matrix_set(cartan, 1, 2, 2*s*cos(a1));

	gsl_matrix_set(cartan, 2, 0, 2*s*cos(a2));
	gsl_matrix_set(cartan, 2, 1, 2/s*cos(a1));
	gsl_matrix_set(cartan, 2, 2, -2);
}

void initialize_triangle_generators(gsl_matrix **gen, gsl_matrix *cartan)
{
	LOOP(i) {
		gsl_matrix_set_identity(gen[i]);
		LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j);
	}
}

typedef struct { double x[3]; } point;

point apply(gsl_matrix *g, point x)
{
	point out;
	LOOP(i) out.x[i] = 0.0;
	LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(g, i, j) * x.x[j];
	return out;
}

point incidence(point a, point b)
{
	point c;
	LOOP(i) c.x[i] = a.x[(i+1)%3]*b.x[(i+2)%3] - a.x[(i+2)%3]*b.x[(i+1)%3];
	return c;
}

double applyForm(point x, point y, gsl_matrix *form)
{
	double out = 0.0;
	LOOP(i) LOOP(j) out += x.x[i] * gsl_matrix_get(form, i, j) * y.x[j];
	return out;
}

void boundary(point x, point y, point *a, point *b, gsl_matrix *form)
{
	double formX = applyForm(x, x, form);
	double formY = applyForm(y, y, form);
	double formXY = applyForm(x, y, form);
	double t1 = (- formXY + sqrt(formXY * formXY - formX * formY)) / formX;
	double t2 = (- formXY - sqrt(formXY * formXY - formX * formY)) / formX;
	LOOP(i) a->x[i] = t1 * x.x[i] + y.x[i];
	LOOP(i) b->x[i] = t2 * x.x[i] + y.x[i];
}

point column(gsl_matrix *m, int j)
{
	point out;
	LOOP(i) out.x[i] = gsl_matrix_get(m, i, j);
	return out;
}

point row(gsl_matrix *m, int i)
{
	point out;
	LOOP(j) out.x[j] = gsl_matrix_get(m, i, j);
	return out;
}

double coord(point p, int k, gsl_matrix *frame)
{
	double tmp[3] = {0, 0, 0};
	double factor = 1.0;

	LOOP(i) LOOP(j) tmp[i] += gsl_matrix_get(frame, i, j)*p.x[j];

#ifdef POINCARE
	double norm2 = (tmp[0]*tmp[0] + tmp[1]*tmp[1]) / (tmp[2]*tmp[2]);
	if(fabs(norm2) < 0.5) // just to avoid dividing by 0
		factor = 1.0 / (1.0 + sqrt(1.0 - norm2));
	else if(fabs(norm2) < 1.0)
		factor = (1.0 - sqrt(1.0 - norm2)) / norm2;
	else
		factor = 1.0;
#endif

	return factor*tmp[k]/tmp[2];
}

void print_tex_header() {
	char *header =
		"\\documentclass{standalone}\n"
		"\\usepackage[utf8]{inputenc}\n"
		"\\usepackage{tikz}\n"
		"\\begin{document}\n"
		"\\begin{tikzpicture}[scale=5]\n";

	printf("%s", header);
}

void print_tex_footer() {
	char *footer =
		"\\draw[thick] (0,0) circle (1.0);\n"
		"\\end{tikzpicture}\n"
		"\\end{document}\n";

	printf("%s", footer);
}

char *poincare_arc(double x1, double x2, double y1, double y2, char *buffer)
{
	double t = (1 - x1*y1 - x2*y2)/(x1*y2 - x2*y1)/2;
	double m1 = x1/2 + y1/2 + t*(y2 - x2);  // center of circle
	double m2 = x2/2 + y2/2 + t*(x1 - y1);

	double r = sqrt((x1-m1)*(x1-m1) + (x2-m2)*(x2-m2));
	double phix = atan2(x2-m2,x1-m1);
	double phiy = atan2(y2-m2,y1-m1);

	if(phix - phiy > M_PI)
		phiy += 2*M_PI;
	else if(phiy - phix > M_PI)
		phix += 2*M_PI;

	sprintf(buffer, "%f:%f:%f", phix/M_PI*180, phiy/M_PI*180, r);
	return buffer;
}

void draw_triangle(point *p, gsl_matrix *frame, const char *arguments)
{
#ifdef POINCARE
	char buffer1[100], buffer2[100], buffer3[100];

	double x1 = coord(p[0], 0, frame);
	double y1 = coord(p[0], 1, frame);
	double x2 = coord(p[1], 0, frame);
	double y2 = coord(p[1], 1, frame);
	double x3 = coord(p[2], 0, frame);
	double y3 = coord(p[2], 1, frame);

	printf("\\draw[%s] (%f, %f) arc (%s) arc (%s) arc (%s);\n", arguments, x1, y1,
		   poincare_arc(x1, y1, x2, y2, buffer1),
		   poincare_arc(x2, y2, x3, y3, buffer2),
		   poincare_arc(x3, y3, x1, y1, buffer3));

#else
	printf("\\draw[%s] (%f, %f) -- (%f, %f) -- (%f, %f) -- cycle;\n", arguments,
		   coord(p[0], 0, frame), coord(p[0], 1, frame),
		   coord(p[1], 0, frame), coord(p[1], 1, frame),
		   coord(p[2], 0, frame), coord(p[2], 1, frame));
#endif
}

void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments)
{
	char buffer[100];

	double x1 = coord(p1, 0, frame);
	double y1 = coord(p1, 1, frame);
	double x2 = coord(p2, 0, frame);
	double y2 = coord(p2, 1, frame);

#ifdef POINCARE
	printf("\\draw[%s] (%f, %f) arc (%s);\n", arguments, x1, y1,
		   poincare_arc(x1, y1, x2, y2, buffer));
#else
	printf("\\draw[%s] (%f, %f) -- (%f, %f);\n", arguments, x1, y1, x2, y2);
#endif
}

int main()
{
	groupelement_t *group;
	gsl_matrix **matrices;
	gsl_matrix *cartan;
	gsl_matrix *gen[3];
	gsl_matrix *coxeter[3];
	gsl_matrix *coxeter_eigenvectors[3];
	gsl_matrix *frame;
	workspace_t *ws;
	int elements = 5000;
	int p = 5, q = 5, r = 5;

	group = malloc(elements*sizeof(groupelement_t));
	matrices = malloc(elements*sizeof(gsl_matrix*));
	for(int i = 0; i < elements; i++)
		matrices[i] = gsl_matrix_alloc(3, 3);
	cartan = gsl_matrix_alloc(3, 3);
	frame = gsl_matrix_alloc(3, 3);
	LOOP(i) gen[i] = gsl_matrix_alloc(3, 3);
	LOOP(i) coxeter[i] = gsl_matrix_alloc(3, 3);
	LOOP(i) coxeter_eigenvectors[i] = gsl_matrix_alloc(3, 3);
	ws = workspace_alloc(3);

	generate_triangle_group(group, elements, p, q, r);
	cartan_matrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0);
	initialize_triangle_generators(gen, cartan);
	int pos = diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form

	gsl_matrix_set_identity(matrices[0]);
	for(int i = 1; i < elements; i++)
		multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);

	LOOP(i) multiply_many(ws, coxeter[i], 3, gen[i%3], gen[(i+1)%3], gen[(i+2)%3]);
	LOOP(i) eigenvectors(coxeter[i], coxeter_eigenvectors[i], ws);

	/*
	for(int i = 0; i < elements; i++) {
		printf("%4d: ", i);
		for(groupelement_t *cur = &group[i]; cur->parent; cur = cur->parent)
			fputc(cur->letter+'a', stdout);
		fputc('\n', stdout);
	}

	return 0;*/

	point coxeter_attracting[3];
	point coxeter_repelling[3];
	point reflection_lines[3];
	point triangle_points[3];
	point transformed[3];
	point transformed2[3];

	LOOP(i) coxeter_attracting[i] = column(coxeter_eigenvectors[i], 0);
	LOOP(i) coxeter_repelling[i] = column(coxeter_eigenvectors[i], 2);
	LOOP(i) reflection_lines[i] = row(cartan, i);
	LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]);

	print_tex_header();

//	int indices[10] = {0, 1, 6, 10, 30, 46, 124, 185, 484, 717};
//	int indices[10] = {0, 1, 4, 10, 22};

	for(int k = 0; k < elements; k++) {
		if(group[k].length % 2)
			continue;

		LOOP(i) transformed[i] = apply(matrices[k], triangle_points[i]);
		draw_triangle(transformed, frame, "black,fill=black!10,line width=0pt");
	}

	for(int k = 0; k < elements; k++) {
		if(group[k].length % 2)
			continue;

		LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]);
		LOOP(i) transformed2[i] = apply(matrices[k], coxeter_repelling[i]);
		LOOP(i) draw_line(transformed[i], transformed2[i], frame, "red");
	}

	print_tex_footer();
}