From 0763056ccb1a5bfa09117938eab87cba4b1a815a Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Tue, 14 Jun 2022 15:41:43 +0200 Subject: [PATCH] compute complex max slope --- Makefile | 4 ++-- complex_anosov.c | 48 +++++++++++++++++++++++++++++++++++++++++------ enumerate.c | 49 ++++++++++++++++++++++++++++++++++++++++++++++++ enumerate.h | 4 ++++ 4 files changed, 97 insertions(+), 8 deletions(-) diff --git a/Makefile b/Makefile index daf8750..3abece2 100644 --- a/Makefile +++ b/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= diff --git a/complex_anosov.c b/complex_anosov.c index 7ecc5cc..1c7e379 100644 --- a/complex_anosov.c +++ b/complex_anosov.c @@ -9,6 +9,7 @@ #include #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", - 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)); + solve_characteristic_polynomial_d(solver, poly, + gaussian_sqrt5_real(traces[i].tr), + gaussian_sqrt5_imag(traces[i].tr), + gaussian_sqrt5_real(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); diff --git a/enumerate.c b/enumerate.c index f779988..bd847e9 100644 --- a/enumerate.c +++ b/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; +} diff --git a/enumerate.h b/enumerate.h index cca03e9..0e3b171 100644 --- a/enumerate.h +++ b/enumerate.h @@ -4,6 +4,8 @@ #include "mat.h" #include "coxeter.h" +#include + 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