compute complex max slope

Makefile

@@ -1,8 +1,8 @@
 HEADERS=mat.h coxeter.h enumerate.h generators.h qext.h
 
-SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
+#SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
 #SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline
-#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
+SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
 #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL
 #SPECIAL_OPTIONS=
 

complex_anosov.c

@@ -9,6 +9,7 @@
 #include <math.h>
 
 #define LOOP(i,n) for(int i = 0; i < (n); i++)
+#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
 
 /*
   Elements up to length 0: 1
@@ -134,15 +135,50 @@ int main(int argc, char *argv[])
 	struct tracedata *traces;
 	int nuniq = enumerate_coxeter_group_traces(group, gen, &traces);
 
+	mps_context *solver = mps_context_new();
+	mps_monomial_poly *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);
+
+	double ev_real[3], ev_imag[3], ev_abs2[3];
+	double max_slope = 0;
+	int max_slope_id = 0;
+
 	LOOP(i, nuniq) {
-		printf("%d %f %f %f %f\n",
+		solve_characteristic_polynomial_d(solver, poly,
-		       traces[i].id,
 		                                  gaussian_sqrt5_real(traces[i].tr),
 		                                  gaussian_sqrt5_imag(traces[i].tr),
 		                                  gaussian_sqrt5_real(traces[i].trinv),
-		       gaussian_sqrt5_imag(traces[i].trinv));
+		                                  gaussian_sqrt5_imag(traces[i].trinv),
+		                                  ev_real, ev_imag);
+
+		LOOP(j, 3) ev_abs2[j] = ev_real[j]*ev_real[j] + ev_imag[j]*ev_imag[j];
+
+		if(fabs(ev_abs2[0]) < fabs(ev_abs2[1]))
+			SWAP(double, ev_abs2[0], ev_abs2[1]);
+		if(fabs(ev_abs2[1]) < fabs(ev_abs2[2]))
+			SWAP(double, ev_abs2[1], ev_abs2[2]);
+		if(fabs(ev_abs2[0]) < fabs(ev_abs2[1]))
+			SWAP(double, ev_abs2[0], ev_abs2[1]);
+
+		if(log(ev_abs2[0]) < 1e-3) // we regard this as a finite order element
+			continue;
+
+		double slope = - log(ev_abs2[0]) / log(ev_abs2[2]);
+		if(slope > max_slope) {
+			max_slope = slope;
+			max_slope_id = traces[i].id;
 		}
 
+		printf("%d %f %f %f\n",
+		       traces[i].id, log(ev_abs2[0]), log(ev_abs2[1]), log(ev_abs2[2]));
+	}
+
+	coxeter_snprint(buf, sizeof(buf), &group->elements[max_slope_id]);
+	gmp_printf("q = %Qd + i*%Qd\tElements: %d\tTraces: %d\tMaximal slope: %f at %s\n", qreal, qimag, n, nuniq, max_slope, buf);
+
+	mps_monomial_poly_free(solver, MPS_POLYNOMIAL(poly));
+	mps_context_free(solver);
 	enumerate_tracedata_clear(traces, nuniq);
 	LOOP(i, 3) mat_clear(gen[i]);
 	coxeter_clear(group);

enumerate.c

@@ -120,3 +120,52 @@ void enumerate_tracedata_clear(struct tracedata *traces, int n)
 	}
 	free(traces);
 }
+
+// returns squares of absolute values
+int solve_characteristic_polynomial_d(mps_context *solv, mps_monomial_poly *poly, double tr_real, double tr_imag, double trinv_real, double trinv_imag, double *eigenvalues_real, double *eigenvalues_imag)
+{
+//        mpq_t coeff1, coeff2, zero;
+        cplx_t *roots;
+        double radii[3];
+        double *radii_p[3];
+        mps_boolean errors;
+        int result = 0;
+
+        /*
+        mpq_inits(coeff1, coeff2, zero, NULL);
+        mpq_set(coeff1, trinv);
+        mpq_sub(coeff2, zero, tr);
+        */
+
+        mps_monomial_poly_set_coefficient_d(solv, poly, 0, -1, 0);
+        mps_monomial_poly_set_coefficient_d(solv, poly, 1, trinv_real, trinv_imag);
+        mps_monomial_poly_set_coefficient_d(solv, poly, 2, -tr_real, -tr_imag);
+//        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_d(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_real[i] = cplx_Re(roots[i]);
+                        eigenvalues_imag[i] = cplx_Im(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;
+}

enumerate.h

@@ -4,6 +4,8 @@
 #include "mat.h"
 #include "coxeter.h"
 
+#include <mps/mps.h>
+
 struct tracedata {
 	int id;
 	NUMBER tr;
@@ -15,4 +17,6 @@ int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **
 
 void enumerate_tracedata_clear(struct tracedata *traces, int n);
 
+int solve_characteristic_polynomial_d(mps_context *solv, mps_monomial_poly *poly, double tr_real, double tr_imag, double trinv_real, double trinv_imag, double *eigenvalues_real, double *eigenvalues_imag);
+
 #endif