From 2568b4ef7134f3202f04193d794b2a269d3b8c03 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 29 Jul 2020 20:23:15 -0400 Subject: [PATCH] use multiprecision poly solver --- Makefile | 6 ++-- singular_values.c | 79 ++++++++++++++++++++++++++++++++++++++++++--- singular_values.plt | 15 ++++++--- 3 files changed, 88 insertions(+), 12 deletions(-) diff --git a/Makefile b/Makefile index 544235f..e9218e5 100644 --- a/Makefile +++ b/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 -flto -funroll-loops -Winline +SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline +#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 diff --git a/singular_values.c b/singular_values.c index 04d3068..d655792 100644 --- a/singular_values.c +++ b/singular_values.c @@ -3,12 +3,60 @@ #include "mat.h" #include +#include //#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); } diff --git a/singular_values.plt b/singular_values.plt index 03613ce..567ff08 100644 --- a/singular_values.plt +++ b/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 yrange [0:20] -set trange [0:20] +set xrange [0:30] +set yrange [0:30] +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 "", \ - log(tr(t,t/2)),log(trinv(t,t/2)) w l lw 2 t "" + t,2*t 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 ""