From b5bd9d398fb52b73ddfe523025092a0030569a8f Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 31 Mar 2021 22:45:25 -0500 Subject: [PATCH] draw character variety picture --- Makefile | 10 +- singular_values.c | 281 +++++++++++++++++++++++----------------------- 2 files changed, 148 insertions(+), 143 deletions(-) diff --git a/Makefile b/Makefile index 5710721..813e839 100644 --- a/Makefile +++ b/Makefile @@ -7,14 +7,20 @@ SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS) -all: singular_values +all: singular_values special_element singular_values: singular_values.o coxeter.o linalg.o mat.o gcc $(OPTIONS) -o singular_values coxeter.o linalg.o singular_values.o mat.o -lm -lgsl -lcblas -lgmp -lmps -lpthread +special_element: special_element.o coxeter.o linalg.o mat.o + gcc $(OPTIONS) -o special_element coxeter.o linalg.o special_element.o mat.o -lm -lgsl -lcblas -lgmp -lmps -lpthread + singular_values.o: singular_values.c $(HEADERS) gcc $(OPTIONS) -c singular_values.c +special_element.o: special_element.c $(HEADERS) + gcc $(OPTIONS) -c special_element.c + linalg.o: linalg.c $(HEADERS) gcc $(OPTIONS) -c linalg.c @@ -25,4 +31,4 @@ mat.o: mat.c $(HEADERS) gcc $(OPTIONS) -c mat.c clean: - rm -f singular_values coxeter.o linalg.o singular_values.o mat.o + rm -f singular_values special_element coxeter.o linalg.o singular_values.o mat.o special_element.o diff --git a/singular_values.c b/singular_values.c index 58013ad..25600d7 100644 --- a/singular_values.c +++ b/singular_values.c @@ -6,11 +6,26 @@ #include #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, ...) do { print_time(); fprintf(stderr, msg, ##__VA_ARGS__); } while (0); +//#define DEBUG(msg, ...) -#define OUTPUT_POINTS //#define OUTPUT_POINTS +#define OUTPUT_SUMMARY + +struct timespec starttime; + +void print_time() +{ + double diff; + struct timespec current; + + clock_gettime(CLOCK_REALTIME, ¤t); + + diff = (current.tv_sec - starttime.tv_sec) + (current.tv_nsec - starttime.tv_nsec)*1e-9; + + fprintf(stderr, "[%.3f] ", diff); +} + struct result { mpq_t tr; @@ -124,99 +139,75 @@ void quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e) mpq_clear(tmp); } -void initialize_triangle_generators(mat_workspace *ws, mat *gen, mpq_t s, mpq_t q) +// this version is only for the (4,4,4) group +void initialize_triangle_generators(mat_workspace *ws, mat *gen, mpq_t m, mpq_t t) { - mat r1,r2,r3; - mpq_t rho1, rho2, rho3; - mpq_t b1,c1,a2,c2,a3,b3; - mpq_t sinv; + mpq_t s,sinv,q,x,y; + mpq_t zero, one, two; + mpq_t tmp; - mpq_inits(sinv,rho1,rho2,rho3,b1,c1,a2,c2,a3,b3,NULL); - mat_init(r1, 3); - mat_init(r2, 3); - mat_init(r3, 3); + mpq_inits(s,sinv,q,x,y,zero,one,two,tmp,NULL); + mpq_set_ui(zero, 0, 1); + mpq_set_ui(one, 1, 1); + mpq_set_ui(two, 2, 1); - mpq_set_ui(sinv, 1, 1); - mpq_div(sinv, sinv, s); + // s = (1-m^2)/2m + mpq_mul(s, m, m); + mpq_sub(s, one, s); + mpq_div(s, s, m); + mpq_div(s, s, two); + mpq_div(sinv, one, s); - quartic(rho1, s, 0, 0, 1, -1, 1); - quartic(rho2, s, 0, 0, 1, -1, 1); - quartic(rho3, s, 0, 0, 1, 0, 1); + // q = (1+m^2)/(1-m^2) = 2/(1-m^2) - 1 + mpq_mul(q, m, m); + mpq_sub(q, one, q); + mpq_div(q, two, q); + mpq_sub(q, q, one); - mpq_mul(c1, rho2, q); - mpq_mul(a2, rho3, q); - mpq_mul(b3, rho1, q); - mpq_set_ui(b1, 1, 1); - mpq_set_ui(c2, 1, 1); - mpq_set_ui(a3, 1, 1); - mpq_div(b1, b1, q); - mpq_div(c2, c2, q); - mpq_div(a3, a3, q); + // x = -tq, y = -q/t + mpq_mul(x, q, t); + mpq_sub(x, zero, x); + mpq_div(y, q, t); + mpq_sub(y, zero, y); - // actually, we want minus everything - mat_zero(r1); - mat_zero(r2); - mat_zero(r3); - mpq_set_si(*mat_ref(r1, 0, 0), -1, 1); - mpq_set_si(*mat_ref(r1, 1, 1), 1, 1); - mpq_set_si(*mat_ref(r1, 2, 2), 1, 1); - mpq_set_si(*mat_ref(r2, 0, 0), 1, 1); - mpq_set_si(*mat_ref(r2, 1, 1), -1, 1); - mpq_set_si(*mat_ref(r2, 2, 2), 1, 1); - mpq_set_si(*mat_ref(r3, 0, 0), 1, 1); - mpq_set_si(*mat_ref(r3, 1, 1), 1, 1); - mpq_set_si(*mat_ref(r3, 2, 2), -1, 1); + // q^2 = xy = 1 + 1/s^2 + // [ -s s*y 0] + // [ -s*x s*x*y - 1/s 0] + // [ -s*y s*y^2 - x 1] + LOOP(i,3) { + mat_zero(gen[i]); + mpq_sub(tmp, zero, s); + mat_set(gen[i%3], i%3, i%3, tmp); + mpq_mul(tmp, s, y); + mat_set(gen[i%3], i%3, (i+1)%3, tmp); + mpq_mul(tmp, s, x); + mpq_sub(tmp, zero, tmp); + mat_set(gen[i%3], (i+1)%3, i%3, tmp); + mpq_mul(tmp, s, x); + mpq_mul(tmp, tmp, y); + mpq_sub(tmp, tmp, sinv); + mat_set(gen[i%3], (i+1)%3, (i+1)%3, tmp); + mpq_mul(tmp, s, y); + mpq_sub(tmp, zero, tmp); + mat_set(gen[i%3], (i+2)%3, i%3, tmp); + mpq_mul(tmp, s, y); + mpq_mul(tmp, tmp, y); + mpq_sub(tmp, tmp, x); + mat_set(gen[i%3], (i+2)%3, (i+1)%3, tmp); + mat_set(gen[i%3], (i+2)%3, (i+2)%3, one); + } - mpq_set(*mat_ref(r1, 1, 0), b1); - mpq_set(*mat_ref(r1, 2, 0), c1); - mpq_set(*mat_ref(r2, 0, 1), a2); - mpq_set(*mat_ref(r2, 2, 1), c2); - mpq_set(*mat_ref(r3, 0, 2), a3); - mpq_set(*mat_ref(r3, 1, 2), b3); - - mat_zero(gen[0]); - mat_zero(gen[1]); - mat_zero(gen[2]); - - mpq_set_ui(*mat_ref(gen[0], 0, 0), 1, 1); - mat_set(gen[0], 1, 1, sinv); - mat_set(gen[0], 2, 2, s); - - mat_set(gen[1], 0, 0, s); - mpq_set_ui(*mat_ref(gen[1], 1, 1), 1, 1); - mat_set(gen[1], 2, 2, sinv); - - mat_set(gen[2], 0, 0, sinv); - mat_set(gen[2], 1, 1, s); - mpq_set_ui(*mat_ref(gen[2], 2, 2), 1, 1); - - mat_multiply(ws, gen[0], r2, gen[0]); - mat_multiply(ws, gen[0], gen[0], r3); - mat_multiply(ws, gen[1], r3, gen[1]); - mat_multiply(ws, gen[1], gen[1], r1); - mat_multiply(ws, gen[2], r1, gen[2]); - mat_multiply(ws, gen[2], gen[2], r2); - - mat_pseudoinverse(ws, gen[3], gen[0]); - mat_pseudoinverse(ws, gen[4], gen[1]); - mat_pseudoinverse(ws, gen[5], gen[2]); + LOOP(i,3) mat_pseudoinverse(ws, gen[i+3], gen[i]); + // debug output /* - mat_print(r1); - mat_print(r2); - mat_print(r3); + gmp_printf("m = %Qd, s = %Qd, t = %Qd, q = %Qd, x = %Qd, y = %Qd\n", m, s, t, q, x, y); mat_print(gen[0]); mat_print(gen[1]); mat_print(gen[2]); - mat_print(gen[3]); - mat_print(gen[4]); - mat_print(gen[5]); */ - mpq_clears(sinv,rho1,rho2,rho3,b1,c1,a2,c2,a3,b3,NULL); - mat_clear(r1); - mat_clear(r2); - mat_clear(r3); + mpq_inits(s,sinv,q,x,y,zero,one,two,tmp,NULL); } char *print_word(groupelement_t *g, char *str) @@ -232,7 +223,7 @@ char *print_word(groupelement_t *g, char *str) return str; } -void enumerate(group_t *group, mat *matrices, mpq_t s, mpq_t t) +void enumerate(group_t *group, mat *matrices, mpq_t m, mpq_t t) { mat_workspace *ws; mat tmp; @@ -245,7 +236,7 @@ void enumerate(group_t *group, mat *matrices, mpq_t s, mpq_t t) mat_init(gen[i], 3); mat_init(tmp, 3); - initialize_triangle_generators(ws, gen, s, t); + initialize_triangle_generators(ws, gen, m, t); mat_identity(matrices[0]); for(int i = 1; i < group->size; i++) { @@ -352,8 +343,8 @@ double max_slope(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t, int *in int main(int argc, char *argv[]) { - mpq_t s, q, t, tmp; - double sapprox, tapprox, qapprox, tqfactor; + mpq_t m, t, tmp; + double s; mat *matrices; group_t *group; int index; @@ -369,11 +360,13 @@ int main(int argc, char *argv[]) struct result *invariants; struct result **distinct_invariants; + clock_gettime(CLOCK_REALTIME, &starttime); + nmax = atoi(argv[1]); DEBUG("Allocate\n"); - mpq_inits(s, q, t, tmp, NULL); + mpq_inits(m, t, tmp, NULL); matrices = malloc(nmax*sizeof(mat)); for(int i = 0; i < nmax; i++) mat_init(matrices[i], 3); @@ -389,6 +382,7 @@ int main(int argc, char *argv[]) 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 @@ -420,81 +414,86 @@ int main(int argc, char *argv[]) #ifdef OUTPUT_POINTS gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor); #endif + */ // group DEBUG("Generate group\n"); - group = coxeter_init_triangle(4, 3, 3, nmax); + group = coxeter_init_triangle(4, 4, 4, nmax); - DEBUG("Compute matrices\n"); - enumerate(group, matrices, s, q); + fprintf(stderr, "max word length = %d\n", group->elements[nmax-1].length); - 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; + for(int ttick = 45; ttick <= 65; ttick++) { + for(int mtick = 45; mtick < 65; mtick++) { + mpq_set_ui(t, ttick, 100); + mpq_set_ui(m, mtick, 100); // 414/1000 ~ sqrt(2)-1 <-> s=1 + s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m)); - mat_trace(invariants[i].tr, matrices[i]); - mat_trace(invariants[i].trinv, matrices[group->elements[i].inverse->id]); + DEBUG("Compute matrices\n"); + enumerate(group, matrices, m, t); - distinct_invariants[n++] = &invariants[i]; - } + 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; - DEBUG("Get unique traces\n"); + mat_trace(invariants[i].tr, matrices[i]); + mat_trace(invariants[i].trinv, matrices[group->elements[i].inverse->id]); - qsort(distinct_invariants, n, sizeof(struct result*), compare_result); + distinct_invariants[n++] = &invariants[i]; - 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]; - } +// gmp_printf("%Qd %Qd %d %s\n", invariants[i].tr, invariants[i].trinv, i, print_word(&group->elements[i], buf)); + } - max_slope = 0; - int max_slope_index; + DEBUG("Get unique traces\n"); - 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); - if(retval == 1) { - fprintf(stderr, "Error! Could not solve polynomial.\n"); - continue; - } else if(retval == 2) { - continue; - } + qsort(distinct_invariants, n, sizeof(struct result*), compare_result); - 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])) + 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]; + } + + 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); + 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]); - double x = log(fabs(evs[0])); - double y = -log(fabs(evs[2])); + double x = log(fabs(evs[0])); + double y = -log(fabs(evs[2])); - if(y/x > max_slope && (x > 0.1 || y > 0.1)) { - max_slope_index = distinct_invariants[i] - invariants; - max_slope = y/x; - } + if(y/x > max_slope && (x > 0.1 || y > 0.1)) { + max_slope_index = distinct_invariants[i] - invariants; + max_slope = y/x; + } #ifdef OUTPUT_POINTS - gmp_printf("%Qd %Qd %f %f %f %d\n", distinct_invariants[i]->tr, distinct_invariants[i]->trinv, x, y, y/x); + gmp_printf("%Qd %Qd %f %f %f\n", distinct_invariants[i]->tr, distinct_invariants[i]->trinv, x, y, y/x); #endif - } + } #ifdef OUTPUT_SUMMARY - fprintf(stdout, "%.3f %.3f %f %s\n", mpq_get_d(s), mpq_get_d(q)*tqfactor, max_slope, print_word(&group->elements[max_slope_index], buf)); +// fprintf(stdout, "%.5f %.5f %.5f %f %s\n", mpq_get_d(t), mpq_get_d(m), s, max_slope, print_word(&group->elements[max_slope_index], buf)); + fprintf(stdout, "%.5f %.5f %.5f %f %s\n", mpq_get_d(t), mpq_get_d(m), s, max_slope, print_word(&group->elements[max_slope_index], buf)); #endif - -// output_invariants(group, matrices, s, q, solver); - -// for(int i = 0; i < 10; i++) { -// mpq_set_ui(t,100+i,100); -// mpq_canonicalize(t); - - //printf("%f %f\n", mpq_get_d(t), max_slope(group, matrices, s, t, &index)); -// } + } + } DEBUG("Clean up\n"); for(int i = 0; i < nmax; i++) { @@ -507,6 +506,6 @@ int main(int argc, char *argv[]) mat_clear(matrices[i]); free(matrices); coxeter_clear(group); - mpq_clears(s, q, t, tmp, NULL); + mpq_clears(m, t, tmp, NULL); mps_context_free(solver); }