reorganize singular_values, do one precomputation and then only compute unique traces

Makefile

@@ -1,7 +1,7 @@
 HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h
 
 #SPECIAL_OPTIONS=-O0 -g -D_DEBUG
-SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
+SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline
 #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL
 #SPECIAL_OPTIONS=

billiard_words.hs

@@ -9,17 +9,21 @@ main = do
 
 listWordsUpToLength :: Int -> IO ()
 listWordsUpToLength n = do
-  putStrLn $ unlines [printf "%s %d/%d %f %s"
+  putStr $ unlines [printf "%s %d/%d %f"
                              w
                              (p `div` gcd p q)
                              (q `div` gcd p q)
-                             (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double)
-                             (slopeWord "bca" (orthogonalSlope (p,q))) |
+                             (atan (sqrt 3 / (2*q_/p_ + 1))) |
                       ((p,q),w) <- wordlist (n `div` 2, n `div` 2),
+                      let p_ = fromIntegral p :: Double,
+                      let q_ = fromIntegral q :: Double,
                       length w <= n,
                       let x = 2*q + p,
                       let y = 2*p + q]
 
+--                             (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double) |
+--                             (slopeWord "bca" (orthogonalSlope (p,q))) |
+
 wordlist :: (Int,Int) -> [((Int,Int),String)]
 wordlist (pmax,qmax) = nub $
                      sortBy (comparing sl)

coxeter.h

@@ -7,6 +7,7 @@ typedef struct _groupelement {
 	int id;
 	int letter;
 	struct _groupelement *parent;
+	int need_to_compute; // optimization flag; if 0, will skip matrix computation
 
 	int length;
 	struct _groupelement *inverse;

enumerate_triangle_group.c

@@ -1,13 +1,12 @@
 #include "enumerate_triangle_group.h"
 #include "linalg.h"
 
-int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues)
+int solve_characteristic_polynomial(mps_context *solv, mps_monomial_poly *poly, mpq_t tr, mpq_t trinv, double *eigenvalues)
 {
 	mpq_t coeff1, coeff2, zero;
 	cplx_t *roots;
 	double radii[3];
 	double *radii_p[3];
-	mps_monomial_poly *poly;
 	mps_boolean errors;
 	int result = 0;
 
@@ -15,7 +14,6 @@ int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, do
 	mpq_set(coeff1, trinv);
 	mpq_sub(coeff2, zero, tr);
 
-	poly = mps_monomial_poly_new(solv, 3);
 	mps_monomial_poly_set_coefficient_int(solv, poly, 0, -1, 0);
 	mps_monomial_poly_set_coefficient_q(solv, poly, 1, coeff1, zero);
 	mps_monomial_poly_set_coefficient_q(solv, poly, 2, coeff2, zero);
@@ -215,7 +213,7 @@ char *print_word(groupelement_t *g, char *str)
   return str;
 }
 
-void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, mpq_t t)
+void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, mpq_t q)
 {
 	mat_workspace *ws;
 	mat tmp;
@@ -228,7 +226,7 @@ void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, m
 		mat_init(gen[i], 3);
 	mat_init(tmp, 3);
 
-	initialize_triangle_generators(ws, gen, p1, p2, p3, s, t);
+	initialize_triangle_generators(ws, gen, p1, p2, p3, s, q);
 
 	mat_identity(matrices[0]);
 	for(int i = 1; i < group->size; i++) {
@@ -236,6 +234,8 @@ void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, m
 			continue;
 		if(!group->elements[i].inverse)
 			continue;
+		if(!group->elements[i].need_to_compute)
+			continue;
 
 		int parent = group->elements[i].parent->id;
 		int grandparent = group->elements[i].parent->parent->id;

enumerate_triangle_group.h

@@ -6,7 +6,7 @@
 
 #include <mps/mps.h>
 
-int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues);
+int solve_characteristic_polynomial(mps_context *solv, mps_monomial_poly *poly, mpq_t tr, mpq_t trinv, double *eigenvalues);
 void continued_fraction_approximation(mpq_t out, double in, int level);
 void quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e);
 // p1,p2,p3 are only allowed to be 2,3,4,6

singular_values.c

@@ -3,9 +3,12 @@
 #include "mat.h"
 #include "enumerate_triangle_group.h"
 
+#include <time.h>
+
 #define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
-//#define DEBUG(msg, ...) fprintf(stderr, msg, ##__VA_ARGS__)
-#define DEBUG(msg, ...)
+
+#define DEBUG(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__);
+//#define DEBUG(msg, ...)
 
 struct result {
 	int id;
@@ -31,7 +34,17 @@ static int compare_result(const void *a_, const void *b_)
 	return d;
 }
 
-static int compare_result_with_id(const void *a_, const void *b_)
+static int compare_result_by_id(const void *a_, const void *b_)
+{
+	int d = 0;
+
+	struct result **a = (struct result **)a_;
+	struct result **b = (struct result **)b_;
+
+	return (*a)->id - (*b)->id;
+}
+
+static int compare_result_by_tr_trinv_id(const void *a_, const void *b_)
 {
 	int d = 0;
 
@@ -62,128 +75,73 @@ static int compare_result_by_slope(const void *a_, const void *b_)
 	return slopea > slopeb ? -1 : slopea < slopeb ? 1 : 0;
 }
 
-int main(int argc, char *argv[])
+struct timespec starttime;
+static void start_timer()
 {
-	mpq_t s, q, t, tmp;
-	double sapprox, tapprox, qapprox, tqfactor;
-	mat *matrices;
-	group_t *group;
-	int index;
+	clock_gettime(CLOCK_MONOTONIC, &starttime);
+}
+
+static double runtime()
+{
+	struct timespec curtime;
+	double diff;
+	clock_gettime(CLOCK_MONOTONIC, &curtime);
+	return (curtime.tv_sec - starttime.tv_sec) + (curtime.tv_nsec - starttime.tv_nsec) / 1e9;
+}
+
+int compute_invariants(group_t *group, mat *matrices, struct result **invariants, int *n, int unique)
+{
+	mpq_t tmp;
 	mps_context *solver;
-	int acc = 100;
-	int n, nuniq, nmax;
+	mps_monomial_poly *poly;
+	int index;
+	int ntraces = *n, nuniq;
 	int retval;
 	double evs[3];
+	int max_slope_id;
 	double max_slope;
 	char buf[100];
-	char buf2[100];
 
-	struct result *invariants;
-	struct result **distinct_invariants;
-
-	if(argc < 4) {
-		fprintf(stderr, "Usage: %s <N> <s> <q>\n", argv[0]);
-		exit(1);
+//	DEBUG("Compute traces\n");
+	for(int i = 0; i < ntraces; i++) {
+		int id = invariants[i]->id;
+		int invid = group->elements[id].inverse->id;
+		mat_trace(invariants[i]->tr, matrices[id]);
+		mat_trace(invariants[i]->trinv, matrices[invid]);
 	}
 
-	nmax = atoi(argv[1]);
+	if(!unique)
+		nuniq = ntraces;
+	else {
+//		DEBUG("Get unique traces\n");
+		qsort(invariants, ntraces, sizeof(struct result*), compare_result);
 
-	DEBUG("Allocate\n");
-
-	mpq_inits(s, q, t, tmp, NULL);
-	matrices = malloc(nmax*sizeof(mat));
-	for(int i = 0; i < nmax; i++)
-		mat_init(matrices[i], 3);
-	invariants = malloc(nmax*sizeof(struct result));
-	distinct_invariants = malloc(nmax*sizeof(struct result));
-	for(int i = 0; i < nmax; i++) {
-		mpq_init(invariants[i].tr);
-		mpq_init(invariants[i].trinv);
-		distinct_invariants[i] = &invariants[i];
+		nuniq = 0;
+		for(int i = 0; i < ntraces; i++) {
+			if(i == 0 || compare_result(&invariants[i], &invariants[nuniq-1]) != 0) {
+				invariants[nuniq] = invariants[i];
+				invariants[nuniq]->count = 1;
+				nuniq++;
+			} else {
+				invariants[nuniq-1]->count++;
+				int oldlength = group->elements[invariants[nuniq-1]->id].length;
+				int newlength = group->elements[invariants[i]->id].length;
+				if(newlength < oldlength)
+					invariants[nuniq-1]->id = invariants[i]->id;
+			}
+		}
 	}
 
+	DEBUG("Solve characteristic polynomials\n");
 	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);
 
-	DEBUG("Approximate parameters\n");
-
-	// get approximate s and q values
-	sapprox = atof(argv[2]);
-	qapprox = atof(argv[3]);
-//	tapprox = atof(argv[3]);
-//	tqfactor = pow((sapprox*sapprox-sapprox+1)*(sapprox*sapprox-sapprox+1)*(sapprox*sapprox+1), 1/6.0);
-//	qapprox = tapprox/tqfactor;
-
-	for(int i = 0; ; i++) {
-		continued_fraction_approximation(tmp, sapprox, i);
-		if(fabs(mpq_get_d(t)-sapprox) < 1e-10
-		   || (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0))
-			break;
-		mpq_set(s, tmp);
-	}
-	mpq_canonicalize(s);
-
-	for(int i = 0; ; i++) {
-		continued_fraction_approximation(tmp, qapprox, i);
-		if(fabs(mpq_get_d(t)-qapprox) < 1e-10
-		   || (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0))
-			break;
-		mpq_set(q, tmp);
-	}
-	mpq_canonicalize(q);
-
-	tqfactor = pow((mpq_get_d(s)*mpq_get_d(s)-mpq_get_d(s)+1)*(mpq_get_d(s)*mpq_get_d(s)-mpq_get_d(s)+1)*(mpq_get_d(s)*mpq_get_d(s)+1), 1/6.0);
-
-	// group
-	// order of the triangle reflection generators: a, b, c
-	// order of the rotation orders: bc, ac, ab
-	DEBUG("Generate group\n");
-	group = coxeter_init_triangle(4, 3, 3, nmax);
-
-	DEBUG("Compute matrices\n");
-	enumerate(group, matrices, 4, 3, 3, s, q);
-
-	n = 0;
-	DEBUG("Compute traces\n");
-	for(int i = 0; i < nmax; i++) {
-		if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse)
-			continue;
-
-		invariants[i].id = i;
-		mat_trace(invariants[i].tr, matrices[i]);
-		mat_trace(invariants[i].trinv, matrices[group->elements[i].inverse->id]);
-
-		distinct_invariants[n++] = &invariants[i];
-	}
-
-	DEBUG("Get unique traces\n");
-
-	qsort(distinct_invariants, n, sizeof(struct result*), compare_result);
-
-	nuniq = 0;
-	for(int i = 0; i < n; i++) {
-		if(i == 0 || compare_result(&distinct_invariants[i], &distinct_invariants[nuniq-1]) != 0) {
-			distinct_invariants[nuniq] = distinct_invariants[i];
-			distinct_invariants[nuniq]->count = 1;
-			nuniq++;
-		} else {
-			distinct_invariants[nuniq-1]->count++;
-			int oldlength = group->elements[distinct_invariants[nuniq-1]->id].length;
-			int newlength = group->elements[distinct_invariants[i]->id].length;
-			if(newlength < oldlength)
-				distinct_invariants[nuniq-1]->id = distinct_invariants[i]->id;
-		}
-
-//		gmp_printf("%d %d %s\n", i, nuniq-1, print_word(&group->elements[i], buf));
-	}
-
 	max_slope = 0;
-	int max_slope_index;
-
-	DEBUG("Solve characteristic polynomials\n");
 	for(int i = 0; i < nuniq; i++) {
-		retval = solve_characteristic_polynomial(solver, distinct_invariants[i]->tr, distinct_invariants[i]->trinv, evs);
+		retval = solve_characteristic_polynomial(solver, poly, invariants[i]->tr, invariants[i]->trinv, evs);
+		retval = 0;evs[0] = 2;evs[1] = 1;evs[2] = 0.5; // fake solving the polynomial for memory leak test
 		if(retval == 1) {
 			fprintf(stderr, "Error! Could not solve polynomial.\n");
 			continue;
@@ -201,52 +159,204 @@ int main(int argc, char *argv[])
 		double x = log(fabs(evs[0]));
 		double y = -log(fabs(evs[2]));
 
-		distinct_invariants[i]->x = x;
-		distinct_invariants[i]->y = y;
+		invariants[i]->x = x;
+		invariants[i]->y = y;
 
-		if(y/x > max_slope && (x > 0.1 || y > 0.1)) {
-			max_slope_index = distinct_invariants[i] - invariants;
+		if(y/x > max_slope + 1e-12 && (x > 0.1 || y > 0.1)) {
+			max_slope_id = invariants[i]->id;
 			max_slope = y/x;
+		} else if(y/x > max_slope - 1e-12 && (x > 0.1 || y > 0.1)) {
+//			DEBUG("%s didn't quite make it\n",
+//			      print_word(&group->elements[invariants[i]->id], buf));
+		}
+	}
+	mps_context_free(solver);
+
+	qsort(invariants, nuniq, sizeof(struct result*), compare_result_by_id);
+
+	*n = nuniq;
+	return max_slope_id;
+}
+
+long check_memory_usage(mat *matrices, int n)
+{
+	mpq_t x;
+	long total;
+
+	for(int i = 0; i < n; i++)
+	{
+		LOOP(j,3) LOOP(k,3) {
+			total += mpq_numref(M(matrices[i], j, k))->_mp_size;
+			total += mpq_denref(M(matrices[i], j, k))->_mp_size;
 		}
 	}
 
-	qsort(distinct_invariants, nuniq, sizeof(struct result*), compare_result_by_slope);
+	return total;
+}
 
-	gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor);
+int main(int argc, char *argv[])
+{
+	mpq_t s, q, t, tmp;
+	int p1, p2, p3;
+	int sstart, send, sdenom, qstart, qend, qdenom;
+	mat *matrices;
+	group_t *group;
+	int nmax, n;
+	int max_slope_id;
+	char buf[100];
+	char buf2[100];
+	struct result *invariants;
+	struct result **distinct_invariants;
 
-//	printf("- 0 0 - - - - 0.5\n");
-	int cumulative = 0;
-	double slope;
-	for(int i = 0; i < nuniq; i++) {
-		slope = distinct_invariants[i]->y/distinct_invariants[i]->x;
+	start_timer();
+
+	mpq_inits(s, q, t, tmp, NULL);
+	if(argc < 11) {
+		fprintf(stderr, "Usage: %s <N> <p1> <p2> <p3> <s start> <s end> <s denom> <q start> <q end> <q denom>\n", argv[0]);
+		exit(1);
+	}
+	nmax = atoi(argv[1]);
+	p1 = atoi(argv[2]);
+	p2 = atoi(argv[3]);
+	p3 = atoi(argv[4]);
+	sstart = atoi(argv[5]);
+	send = atoi(argv[6]);
+	sdenom = atoi(argv[7]);
+	qstart = atoi(argv[8]);
+	qend = atoi(argv[9]);
+	qdenom = atoi(argv[10]);
+
+	DEBUG("Allocate\n");
+	matrices = malloc(nmax*sizeof(mat));
+	for(int i = 0; i < nmax; i++)
+		mat_init(matrices[i], 3);
+	invariants = malloc(nmax*sizeof(struct result));
+	distinct_invariants = malloc(nmax*sizeof(struct result));
+	for(int i = 0; i < nmax; i++) {
+		mpq_init(invariants[i].tr);
+		mpq_init(invariants[i].trinv);
+	}
+
+	// order of the triangle reflection generators: a, b, c
+	// order of the rotation orders: bc, ac, ab
+	DEBUG("Generate group\n");
+	group = coxeter_init_triangle(p1, p2, p3, nmax);
+
+	// first run; compute all matrices
+	for(int i = 0; i < group->size; i++)
+		group->elements[i].need_to_compute = 1;
+
+	// use very generic values for the pilot run unless sstart=send and qstart=qend
+	if(sstart == send && qstart == qend) {
+		mpq_set_ui(s, sstart, sdenom);
+		mpq_set_ui(q, qstart, qdenom);
+		DEBUG("Single run for s = %d/%d, q = %d/%d\n", sstart, sdenom, qstart, qdenom);
+	} else {
+		mpq_set_ui(s, 4, 100);
+		mpq_set_ui(q, 7, 100);
+		DEBUG("Initial run for s = %d/%d, q = %d/%d\n", 4, 100, 7, 100);
+	}
+
+	DEBUG("Compute matrices\n");
+	enumerate(group, matrices, p1, p2, p3, s, q);
+
+	// prepare array of ids
+	n = 0;
+	for(int i = 0; i < group->size; i++)
+	{
+		if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse)
+			continue;
+		invariants[i].id = i;
+		distinct_invariants[n++] = &invariants[i];
+	}
+
+	DEBUG("Compute invariants\n");
+	max_slope_id = compute_invariants(group, matrices, distinct_invariants, &n, 1);
+
+	// prepare for next time; don't need to change ids in distinct_invariants!
+	for(int i = 0; i < group->size; i++)
+		group->elements[i].need_to_compute = 0;
+	group->elements[0].need_to_compute = 1;
+	int multiplication_count = 1;
+	for(int i = 0; i < n; i++) {
+		groupelement_t *cur = &group->elements[distinct_invariants[i]->id];
+		while(cur->need_to_compute == 0) {
+			cur->need_to_compute = 1;
+			multiplication_count++;
+			cur = cur->parent->parent;  // also need to compute its even-length ancestors
+		}
+		cur = group->elements[distinct_invariants[i]->id].inverse;
+		while(cur->need_to_compute == 0) {
+			cur->need_to_compute = 1;
+			multiplication_count++;
+			cur = cur->parent->parent;
+		}
+
+	}
+
+	DEBUG("Would have needed %d matrix multiplications for %d unique traces up to reflection length %d\n", multiplication_count, n, group->elements[group->size-1].length);
+
+	if(sstart != send || qstart != qend) {
+		for(int sloop = sstart; sloop <= send; sloop++) {
+			for(int qloop = qstart; qloop <= qend; qloop++) {
+				DEBUG("Loop for s = %d/%d, q = %d/%d\n", sloop, sdenom, qloop, qdenom);
+				mpq_set_ui(s, sloop, sdenom);
+				mpq_set_ui(q, qloop, qdenom);
+				DEBUG("Compute matrices\n");
+				enumerate(group, matrices, p1, p2, p3, s, q);
+				DEBUG("Compute invariants\n");
+				max_slope_id = compute_invariants(group, matrices, distinct_invariants, &n, 0);
+				// output
+				gmp_printf("%Qd %Qd %s\n", s, q,
+				           print_word(&group->elements[max_slope_id], buf));
+				fflush(stdout);
+			}
+		}
+	} else {
+		// output
+		for(int i = 0; i < n; i++) {
+			double slope = distinct_invariants[i]->y/distinct_invariants[i]->x;
+
+			// exclude tr = trinv = 2/1/0/-1/3
+			mpq_set_si(tmp, 2, 1);
+			if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 &&
+			   mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0)
+				continue;
+			mpq_set_si(tmp, 1, 1);
+			if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 &&
+			   mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0)
+				continue;
+			mpq_set_si(tmp, 0, 1);
+			if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 &&
+			   mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0)
+				continue;
+			mpq_set_si(tmp, -1, 1);
+			if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 &&
+			   mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0)
+				continue;
+			mpq_set_si(tmp, 3, 1);
+			if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 &&
+			   mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0)
+				continue;
+
+			gmp_printf("%d %d %s %f\n",
+			           distinct_invariants[i]->id, distinct_invariants[i]->count,
+			           print_word(&group->elements[distinct_invariants[i]->id], buf),
+			           slope
+				);
+
+			/*
+			gmp_printf("%d %d %d %Qd %Qd %f %f %f %f %f %s\n",
+			           distinct_invariants[i]->id, distinct_invariants[i]->count, cumulative,
+			           distinct_invariants[i]->tr, distinct_invariants[i]->trinv,
+			           log(fabs(mpq_get_d(distinct_invariants[i]->tr))), log(fabs(mpq_get_d(distinct_invariants[i]->trinv))),
+			           distinct_invariants[i]->x, distinct_invariants[i]->y, slope,
+			           print_word(&group->elements[distinct_invariants[i]->id], buf)
+				);
+			*/
+		}
 
-		mpq_set_si(tmp, 1, 1);
-		if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) {
-			continue;
-		}
-		mpq_set_si(tmp, 0, 1);
-		if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) {
-			continue;
-		}
-		mpq_set_si(tmp, -1, 1);
-		if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) {
-			continue;
-		}
-		mpq_set_si(tmp, 3, 1);
-		if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) {
-			continue;
-		}
-
-		cumulative += distinct_invariants[i]->count;
-		gmp_printf("%d %d %d %Qd %Qd %f %f %f %f %f %s\n",
-		           distinct_invariants[i]->id, distinct_invariants[i]->count, cumulative,
-		           distinct_invariants[i]->tr, distinct_invariants[i]->trinv,
-		           log(fabs(mpq_get_d(distinct_invariants[i]->tr))), log(fabs(mpq_get_d(distinct_invariants[i]->trinv))),
-		           distinct_invariants[i]->x, distinct_invariants[i]->y, slope,
-		           print_word(&group->elements[distinct_invariants[i]->id], buf)
-			);
 	}
-//	printf("- 0 %d - - - - 2.0\n", cumulative);
 
 	DEBUG("Clean up\n");
 	for(int i = 0; i < nmax; i++) {
@@ -260,5 +370,4 @@ int main(int argc, char *argv[])
 	free(matrices);
 	coxeter_clear(group);
 	mpq_clears(s, q, t, tmp, NULL);
-	mps_context_free(solver);
 }

special_element.c

@@ -33,6 +33,7 @@ int main(int argc, char *argv[])
 	mat_workspace *ws;
 	mat gen[6];
 	mps_context *solver;
+	mps_monomial_poly *poly;
 	mat element, inverse;
 	int letter1, letter2, letter;
 	mpq_t tr, trinv;
@@ -68,6 +69,7 @@ int main(int argc, char *argv[])
 	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);
 
@@ -112,7 +114,7 @@ int main(int argc, char *argv[])
 		mat_trace(tr, element);
 		mat_trace(trinv, inverse);
 
-		retval = solve_characteristic_polynomial(solver, tr, trinv, evs);
+		retval = solve_characteristic_polynomial(solver, poly, tr, trinv, evs);
 		if(retval == 1) {
 			fprintf(stderr, "Error! Could not solve polynomial.\n");
 			return 1;
@@ -150,7 +152,7 @@ int main(int argc, char *argv[])
 		}
 
 		if(mode != 2)
-			gmp_printf("%s %.9f\n", argv[w+3], slope);
+			gmp_printf("%s %.9f %Qd %Qd\n", argv[w+3], slope, tr, trinv);
 	}
 
 	if(mode != 0)