#include "coxeter.h"
#include "linalg.h"
#include "mat.h"
#include "enumerate_triangle_group.h"

#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
#define DEBUG(msg, ...)
#define DENOMINATOR 300
#define WIDTH 135
#define STARTX 121
#define HEIGHT 300
#define IDX(i,j) (((i)-1)*HEIGHT + ((j)-1))

double mpq_log(mpq_t m_op)
{
	static double logB = log(ULONG_MAX);

        // Undefined logs (should probably return NAN in second case?)
	if (mpz_get_ui(mpq_numref(m_op)) == 0 || mpz_sgn(mpq_numref(m_op)) < 0)
		return -INFINITY;

	// Log of numerator
	double lognum = log(mpq_numref(m_op)->_mp_d[abs(mpq_numref(m_op)->_mp_size) - 1]);
	lognum += (abs(mpq_numref(m_op)->_mp_size)-1) * logB;

	// Subtract log of denominator, if it exists
	if (abs(mpq_denref(m_op)->_mp_size) > 0)
	{
		lognum -= log(mpq_denref(m_op)->_mp_d[abs(mpq_denref(m_op)->_mp_size)-1]);
		lognum -= (abs(mpq_denref(m_op)->_mp_size)-1) * logB;
	}
	return lognum;
}

int main(int argc, char *argv[])
{
	mpq_t m, t, s, q, tmp, tmp2;
	mat_workspace *ws;
	mat gen[6];
	mps_context *solver;
	mat element, inverse;
	int letter1, letter2, letter;
	mpq_t tr, trinv;
	double x, y, slope;
	int retval;
	double evs[3];
	char buf[100];
	double *max_slope;
	int *max_slope_index;

	DEBUG("Allocate\n");

	mpq_inits(m, t, s, q, tmp, tmp2, tr, trinv, NULL);
	ws = mat_workspace_init(3);
	for(int i = 0; i < 6; i++)
		mat_init(gen[i], 3);
	mat_init(element, 3);
	mat_init(inverse, 3);
	max_slope = malloc(sizeof(double)*WIDTH*HEIGHT);
	max_slope_index = malloc(sizeof(int)*WIDTH*HEIGHT);
	memset(max_slope_index, 0, sizeof(int)*WIDTH*HEIGHT);
	memset(max_slope, 0, sizeof(int)*WIDTH*HEIGHT);

	solver = mps_context_new();
	mps_context_set_output_prec(solver, 20); // relative precision
	mps_context_set_output_goal(solver, MPS_OUTPUT_GOAL_APPROXIMATE);

	mpq_set_str(s, argv[1], 10);
	mpq_set_str(q, argv[2], 10);

//	for(int i = STARTX; i <= WIDTH; i++) {
//		for(int j = 1; j <= HEIGHT; j++) {
			for(int w = 3; w < argc; w++) {
//				mpq_set_ui(t, j, DENOMINATOR);
//				mpq_set_ui(m, i, DENOMINATOR); // 414/1000 ~ sqrt(2)-1 <-> s=1
//				s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m));

				DEBUG("Compute matrix\n");
				initialize_triangle_generators(ws, gen, 6, 4, 3, s, q);

				mat_identity(element);
				mat_identity(inverse);
				for(int k = 0; k < strlen(argv[w]); k+=2) {
					letter1 = argv[w][k] - 'a';
					letter2 = argv[w][k+1] - 'a';

					if(letter1 == 1 && letter2 == 2)
						letter = 0; // p = bc
					else if(letter1 == 2 && letter2 == 0)
						letter = 1; // q = ca
					else if(letter1 == 0 && letter2 == 1)
						letter = 2; // r = ab
					else if(letter1 == 2 && letter2 == 1)
						letter = 3; // p^{-1} = cb
					else if(letter1 == 0 && letter2 == 2)
						letter = 4; // q^{-1} = ac
					else if(letter1 == 1 && letter2 == 0)
						letter = 5; // r^{-1} = ba

					mat_multiply(ws, element, element, gen[letter]);
					mat_multiply(ws, inverse, gen[(letter+3)%6], inverse);
				}

				DEBUG("Compute traces\n");

				mat_trace(tr, element);
				mat_trace(trinv, inverse);

				DEBUG("Solve characteristic polynomials\n");
				retval = solve_characteristic_polynomial(solver, tr, trinv, evs);
				if(retval == 1) {
					fprintf(stderr, "Error! Could not solve polynomial.\n");
					return 1;
				}

				if(fabs(evs[0]) < fabs(evs[1]))
					SWAP(double, evs[0], evs[1]);
				if(fabs(evs[1]) < fabs(evs[2]))
					SWAP(double, evs[1], evs[2]);
				if(fabs(evs[0]) < fabs(evs[1]))
					SWAP(double, evs[0], evs[1]);

				x = log(fabs(evs[0]));
				y = -log(fabs(evs[2]));

				if(x > DBL_MAX || y > DBL_MAX) {
					mpq_abs(tmp, tr);
					mpq_abs(tmp2, trinv);
					slope = mpq_log(tmp)/mpq_log(tmp2);
				} else {
					slope = y/x;
				}

				if(slope < 1)
					slope = 1/slope;


//				if(slope > max_slope[IDX(i,j)]) {
//					max_slope[IDX(i,j)] = slope;
//					max_slope_index[IDX(i,j)] = w;
//				}

//	gmp_printf("%Qd %Qd %f %f %f\n", tr, trinv, x, y, y/x);
//				gmp_printf("%s %.5f %.5f %Qd %Qd %.9f %.9f %.9f\n", argv[w], mpq_get_d(s), mpq_get_d(q),
//				           tr, trinv,
//				           x, y, slope);
				gmp_printf("%s %.9f\n", argv[w], slope);

			}

//			printf("%.5f %.5f %d %.9f\n", (double)i/DENOMINATOR, (double)j/DENOMINATOR, max_slope_index[IDX(i,j)], max_slope[IDX(i,j)]);
			fflush(stdout);
//		}
//	}

	DEBUG("Clean up\n");
	mpq_clears(m, t, s, q, tmp, tmp2, tr, trinv, NULL);
	mat_workspace_clear(ws);
	for(int i = 0; i < 6; i++)
		mat_clear(gen[i]);
	mat_clear(element);
	mat_clear(inverse);
	mps_context_free(solver);
	free(max_slope);
	free(max_slope_index);
}