use multiprecision poly solver

Makefile

@@ -1,8 +1,8 @@
 HEADERS=triangle.h linalg.h queue.h mat.h
 
 #SPECIAL_OPTIONS=-O0 -g -D_DEBUG
-#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
+SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
-SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
+#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 #SPECIAL_OPTIONS=
 
 OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
@@ -10,7 +10,7 @@ OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTI
 all: singular_values
 
 singular_values: singular_values.o triangle.o linalg.o mat.o
-	gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o mat.o -lm -lgsl -lcblas -lgmp
+	gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o mat.o -lm -lgsl -lcblas -lgmp -lmps -lpthread
 
 singular_values.o: singular_values.c $(HEADERS)
 	gcc $(OPTIONS) -c singular_values.c

singular_values.c

@@ -3,12 +3,60 @@
 #include "mat.h"
 
 #include <gsl/gsl_poly.h>
+#include <mps/mps.h>
 
 //#define MAX_ELEMENTS 2800000
 //#define MAX_ELEMENTS 720000
-#define MAX_ELEMENTS 3000
+#define MAX_ELEMENTS 1000
 //#define DRAW_PICTURE 1
 
+#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
+
+int solve_characteristic_polynomial(mps_context *solv, 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;
+
+	mpq_inits(coeff1, coeff2, zero, NULL);
+	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);
+	mps_monomial_poly_set_coefficient_int(solv, poly, 3, 1, 0);
+
+	mps_context_set_input_poly(solv, (mps_polynomial*)poly);
+	mps_mpsolve(solv);
+
+	roots = cplx_valloc(3);
+	for(int i = 0; i < 3; i++)
+		radii_p[i] = &(radii[i]);
+	mps_context_get_roots_d(solv, &roots, radii_p);
+	errors = mps_context_has_errors(solv);
+
+	if(errors) {
+		result = 1;
+	} else {
+		for(int i = 0; i < 3; i++) {
+			eigenvalues[i] = cplx_Re(roots[i]);
+			if(fabs(cplx_Im(roots[i])) > radii[i]) // non-real root
+				result = 2;
+		}
+	}
+
+	cplx_vfree(roots);
+	mpq_clears(coeff1, coeff2, zero, NULL);
+
+	return result;
+}
+
 void continued_fraction_approximation(mpq_t out, double in, int level)
 {
 	mpq_t tmp;
@@ -214,10 +262,12 @@ void enumerate(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t)
 	mat_workspace_clear(ws);
 }
 
-void output_invariants(groupelement_t *group, mat *matrices, mpq_t s, mpq_t q)
+void output_invariants(groupelement_t *group, mat *matrices, mpq_t s, mpq_t q, mps_context *solver)
 {
 	mpq_t tr, trinv;
 	char buf[100];
+	double evs[3];
+	int retval;
 
 	mpq_inits(tr, trinv, NULL);
 
@@ -228,7 +278,22 @@ void output_invariants(groupelement_t *group, mat *matrices, mpq_t s, mpq_t q)
 		mat_trace(tr, matrices[i]);
 		mat_trace(trinv, matrices[group[i].inverse->id]);
 
-		gmp_printf("%d %d %s %Qd %Qd %f %f\n", i, group[i].length, print_word(&group[i], buf), tr, trinv, log(mpq_get_d(tr)), log(mpq_get_d(trinv)));
+		retval = solve_characteristic_polynomial(solver, tr, trinv, evs);
+		if(retval == 1) {
+			fprintf(stderr, "Error! Could not solve polynomial.\n");
+			continue;
+		} else if(retval == 2) {
+			continue;
+		}
+
+		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]);
+
+		gmp_printf("%d %d %s %Qd %Qd %f %f\n", i, group[i].length, print_word(&group[i], buf), tr, trinv, log(evs[0]), -log(evs[2]));
 	}
 
 	mpq_clears(tr, trinv, NULL);
@@ -271,6 +336,7 @@ int main(int argc, char *argv[])
 	mat *matrices;
 	groupelement_t *group;
 	int index;
+	mps_context *solver;
 
 	mpq_inits(s, q, t, tmp, NULL);
 	group = malloc(MAX_ELEMENTS*sizeof(groupelement_t));
@@ -278,6 +344,10 @@ int main(int argc, char *argv[])
 	for(int i = 0; i < MAX_ELEMENTS; i++)
 		mat_init(matrices[i], 3);
 
+	solver = mps_context_new();
+	mps_context_set_output_prec(solver, 20); // relative precision
+	mps_context_set_output_goal(solver, MPS_OUTPUT_GOAL_APPROXIMATE);
+
 	int acc = 100;
 
 	sapprox = atof(argv[1]);
@@ -315,7 +385,7 @@ int main(int argc, char *argv[])
 
 		enumerate(group, matrices, s, q);
 		//printf("%f %f\n", mpq_get_d(t), max_slope(group, matrices, s, t, &index));
-		output_invariants(group, matrices, s, q);
+		output_invariants(group, matrices, s, q, solver);
 //	}
 
 	for(int i = 0; i < MAX_ELEMENTS; i++)
@@ -323,4 +393,5 @@ int main(int argc, char *argv[])
 	free(matrices);
 	free(group);
 	mpq_clears(s, q, t, tmp, NULL);
+	mps_context_free(solver);
 }

singular_values.plt

@@ -6,9 +6,9 @@ file = sprintf("< ./singular_values %f %f", exp(logs), exp(logt))
 set zeroaxis
 set samples 1000
 set size square
-set xrange [0:20]
+set xrange [0:30]
-set yrange [0:20]
+set yrange [0:30]
-set trange [0:20]
+set trange [0:30]
 set grid
 set parametric
 
@@ -20,9 +20,14 @@ set parametric
 tr(a,b) = exp(a) + exp(b-a) + exp(-b)
 trinv(a,b) = exp(-a) + exp(a-b) + exp(b)
 
+#plot file using 6:7 w p pt 7 ps 0.5 lc 1 t columnheader,
+#	 log(tr(t,t*2)),log(trinv(t,2*t)) w l lw 2 t "", \
+#	 log(tr(t,t/2)),log(trinv(t,t/2)) w l lw 2 t ""
+
 plot file using 6:7 w p pt 7 ps 0.5 lc 1 t columnheader, \
-	 log(tr(t,t*2)),log(trinv(t,2*t)) w l lw 2 t "", \
+	 t,2*t w l lw 2 t "", \
-	 log(tr(t,t/2)),log(trinv(t,t/2)) w l lw 2 t ""
+	 t,t/2 w l lw 2 t ""
+
 
 #plot for[i=-10:10] log(tr(t,t*exp(log(2)*i/10.0))),log(trinv(t,t*exp(log(2)*i/10.0))) w l lw 2 t ""