#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, ...)

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;
	mps_monomial_poly *poly;
	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 = 0;
	int max_slope_index = 0;
	double min_slope = INFINITY;
	int min_slope_index = 0;
	char *env;
	int mode;

	if(argc < 2) {
		fprintf(stderr,
		        "Usage: %s <s> <q> <word1> <word2> ...\n"
		        "Computes jordan slopes of a list of group elements for a fixed representation.\n"
		        "s,q: representation in the Hitchin component, given as rational numbers, e.g. 2/7\n"
		        "word1, word2, ...: elements in the triangle rotation group, given as reflection group words\n"
		        "output: word - jordan slope pairs\n"
		        "+max slope index, max slope value, max slope word, min slope index, min slope value, min slope word\n"
		        "controlled by environment variable MAXIMUM=no/yes/only, default yes\n",
		        argv[0]);
		exit(0);
	}

	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);

	solver = mps_context_new();
	poly = mps_monomial_poly_new(solver, 3);
	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);

	env = getenv("MAXIMUM");
	if(!env || strcmp(env, "yes") == 0) {
		mode = 1; // yes
	} else if(strcmp(env, "no") == 0) {
		mode = 0; // no
	} else if(strcmp(env, "only") == 0) {
		mode = 2; // only
	}

	for(int w = 0; w < argc - 3; w++) {
		initialize_triangle_generators(ws, gen, 6, 4, 3, s, q);

		mat_identity(element);
		mat_identity(inverse);
		for(int k = 0; k < strlen(argv[w+3]); k+=2) {
			letter1 = argv[w+3][k] - 'a';
			letter2 = argv[w+3][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);
		}

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

		retval = solve_characteristic_polynomial(solver, poly, 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) {
			max_slope = slope;
			max_slope_index = w;
		}

		if(slope < min_slope) {
			min_slope = slope;
			min_slope_index = w;
		}

		if(mode != 2)
			gmp_printf("%s %.9f %Qd %Qd\n", argv[w+3], slope, tr, trinv);
	}

	if(mode != 0)
		printf("%d %.9f %s %d %.9f %s\n",
		       max_slope_index, max_slope, argv[max_slope_index+3],
		       min_slope_index, min_slope, argv[min_slope_index+3]);
	fflush(stdout);

	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);
}