From 721b139307d0da5bd9be6729965ac63b4d298f11 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 13 Apr 2022 19:23:38 -0500 Subject: [PATCH] compute whether a representation in the barbot component of the (5,5,5) group is (almost) all loxodromic --- .gitignore | 1 + Makefile | 12 +- enumerate_triangle_group.c | 220 +++++++++------ enumerate_triangle_group.h | 21 +- mat.c | 14 +- mat.h | 134 ++++++++- singular_values_barbot.c | 545 +++++++++++++++++++++++++++++++++++++ 7 files changed, 843 insertions(+), 104 deletions(-) create mode 100644 singular_values_barbot.c diff --git a/.gitignore b/.gitignore index ebfcc8f..0680d5f 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ triangle_group/singular_values .#* singular_values singular_values_mpi +singular_values_barbot output/ special_element max_slope_picture/generate diff --git a/Makefile b/Makefile index eec1e19..96200e3 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,8 @@ HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h -#SPECIAL_OPTIONS=-O0 -g -D_DEBUG +SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT_SQRT5 #SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline -SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline +#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT_SQRT5 #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL #SPECIAL_OPTIONS= @@ -19,6 +19,9 @@ billiard_words: billiard_words.hs singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o parallel.o mpicc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps +singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o + mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps + #singular_values_mpi: singular_values_mpi.o coxeter.o mat.o # mpicc $(OPTIONS) -o singular_values_mpi coxeter.o singular_values_mpi.o mat.o -lm -lgmp -lmps @@ -28,6 +31,9 @@ special_element: special_element.o coxeter.o linalg.o mat.o enumerate_triangle_g singular_values.o: singular_values.c $(HEADERS) gcc $(OPTIONS) -c singular_values.c +singular_values_barbot.o: singular_values_barbot.c $(HEADERS) + gcc $(OPTIONS) -c singular_values_barbot.c + #singular_values_mpi.o: singular_values_mpi.c $(HEADERS) # mpicc $(OPTIONS) -c singular_values_mpi.c @@ -50,4 +56,4 @@ parallel.o: parallel.c $(HEADERS) gcc $(OPTIONS) -c parallel.c clean: - rm -f singular_values special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_mpi.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words enumerate_triangle_group.o parallel.o + rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.o singular_values_mpi.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words enumerate_triangle_group.o parallel.o diff --git a/enumerate_triangle_group.c b/enumerate_triangle_group.c index 66e82f9..673424f 100644 --- a/enumerate_triangle_group.c +++ b/enumerate_triangle_group.c @@ -70,100 +70,131 @@ void continued_fraction_approximation(mpq_t out, double in, int level) } } -void quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e) +void quartic(NUMBER out, NUMBER in, int a, int b, int c, int d, int e) { - mpq_t tmp; - mpq_init(tmp); + NUMBER tmp; + INIT(tmp); - mpq_set_si(out, a, 1); - mpq_mul(out, out, in); - mpq_set_si(tmp, b, 1); - mpq_add(out, out, tmp); - mpq_mul(out, out, in); - mpq_set_si(tmp, c, 1); - mpq_add(out, out, tmp); - mpq_mul(out, out, in); - mpq_set_si(tmp, d, 1); - mpq_add(out, out, tmp); - mpq_mul(out, out, in); - mpq_set_si(tmp, e, 1); - mpq_add(out, out, tmp); + SET_INT(out, a); - mpq_clear(tmp); + MUL(out, out, in); + SET_INT(tmp, b); + ADD(out, out, tmp); + + MUL(out, out, in); + SET_INT(tmp, c); + ADD(out, out, tmp); + + MUL(out, out, in); + SET_INT(tmp, d); + ADD(out, out, tmp); + + MUL(out, out, in); + SET_INT(tmp, e); + ADD(out, out, tmp); + + CLEAR(tmp); } -// p1,p2,p3 are only allowed to be 2,3,4,6 -void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2, int p3, mpq_t s, mpq_t q) +// discriminant of the polynomial z^3 - x z^2 + y z - 1 +// that is the function x^2 y^2 - 4 x^3 - 4 y^3 - 27 + 18 xy +void discriminant(NUMBER out, NUMBER x, NUMBER y) { - mat r1,r2,r3; - mpq_t rho1, rho2, rho3; - mpq_t b1,c1,a2,c2,a3,b3; - mpq_t sinv; + NUMBER x2, x3, y2, y3, tmp; + INIT(x2);INIT(x3);INIT(y2);INIT(y3);INIT(tmp); - mpq_inits(sinv,rho1,rho2,rho3,b1,c1,a2,c2,a3,b3,NULL); + MUL(x2, x, x); + MUL(x3, x2, x); + MUL(y2, y, y); + MUL(y3, y2, y); + + MUL(out, x2, y2); + + SET_INT(tmp, -4); + MUL(tmp, tmp, x3); + ADD(out, out, tmp); + + SET_INT(tmp, -4); + MUL(tmp, tmp, y3); + ADD(out, out, tmp); + + SET_INT(tmp, -27); + ADD(out, out, tmp); + + SET_INT(tmp, 18); + MUL(tmp, tmp, x); + MUL(tmp, tmp, y); + ADD(out, out, tmp); + + CLEAR(x2);CLEAR(x3);CLEAR(y2);CLEAR(y3);CLEAR(tmp); +} + +void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER s, NUMBER q) +{ + mat_workspace *ws; + mat r1,r2,r3; + NUMBER b1,c1,a2,c2,a3,b3; + NUMBER sinv; + + ws = mat_workspace_init(3); + INIT(sinv);INIT(b1);INIT(c1);INIT(a2);INIT(c2);INIT(a3);INIT(b3); mat_init(r1, 3); mat_init(r2, 3); mat_init(r3, 3); // sinv = s^{-1} - mpq_set_ui(sinv, 1, 1); - mpq_div(sinv, sinv, s); - - // rho_i = s^2 + 2s cos(2 pi / p_i) + 1 - // coefficient 2 is the value for p=infinity, not sure if that would even work - quartic(rho1, s, 0, 0, 1, p1 == 2 ? -2 : p1 == 3 ? -1 : p1 == 4 ? 0 : p1 == 6 ? 1 : 2, 1); - quartic(rho2, s, 0, 0, 1, p2 == 2 ? -2 : p2 == 3 ? -1 : p2 == 4 ? 0 : p2 == 6 ? 1 : 2, 1); - quartic(rho3, s, 0, 0, 1, p3 == 2 ? -2 : p3 == 3 ? -1 : p3 == 4 ? 0 : p3 == 6 ? 1 : 2, 1); + SET_INT(sinv, 1); + DIV(sinv, sinv, s); // c1 = rho2 q, a2 = rho3 q, b3 = rho1 q, b1 = c2 = a3 = 1/q - 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); + MUL(c1, rho2, q); + MUL(a2, rho3, q); + MUL(b3, rho1, q); + SET_INT(b1, 1); + SET_INT(c2, 1); + SET_INT(a3, 1); + DIV(b1, b1, q); + DIV(c2, c2, q); + DIV(a3, a3, q); 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); + SET_INT(*mat_ref(r1, 0, 0), -1); + SET_INT(*mat_ref(r1, 1, 1), 1); + SET_INT(*mat_ref(r1, 2, 2), 1); + SET_INT(*mat_ref(r2, 0, 0), 1); + SET_INT(*mat_ref(r2, 1, 1), -1); + SET_INT(*mat_ref(r2, 2, 2), 1); + SET_INT(*mat_ref(r3, 0, 0), 1); + SET_INT(*mat_ref(r3, 1, 1), 1); + SET_INT(*mat_ref(r3, 2, 2), -1); - 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); + SET(*mat_ref(r1, 1, 0), b1); + SET(*mat_ref(r1, 2, 0), c1); + SET(*mat_ref(r2, 0, 1), a2); + SET(*mat_ref(r2, 2, 1), c2); + SET(*mat_ref(r3, 0, 2), a3); + SET(*mat_ref(r3, 1, 2), b3); mat_zero(gen[0]); mat_zero(gen[1]); mat_zero(gen[2]); // gen[0] = diag(1,s^{-1},s) - mpq_set_ui(*mat_ref(gen[0], 0, 0), 1, 1); + SET_INT(*mat_ref(gen[0], 0, 0), 1); mat_set(gen[0], 1, 1, sinv); mat_set(gen[0], 2, 2, s); // gen[1] = diag(s,1,s^{-1}) mat_set(gen[1], 0, 0, s); - mpq_set_ui(*mat_ref(gen[1], 1, 1), 1, 1); + SET_INT(*mat_ref(gen[1], 1, 1), 1); mat_set(gen[1], 2, 2, sinv); // gen[3] = diag(s^{-1},s,1) 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); + SET_INT(*mat_ref(gen[2], 2, 2), 1); // gen[0] = r2 * gen[0] * r3 // gen[1] = r3 * gen[1] * r1 @@ -182,24 +213,60 @@ void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2, mat_pseudoinverse(ws, gen[4], gen[1]); mat_pseudoinverse(ws, gen[5], gen[2]); - /* - mat_print(r1); - mat_print(r2); - mat_print(r3); - 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_workspace_clear(ws); + CLEAR(sinv);CLEAR(b1);CLEAR(c1);CLEAR(a2);CLEAR(c2);CLEAR(a3);CLEAR(b3); mat_clear(r1); mat_clear(r2); mat_clear(r3); } +#ifdef QEXT_TRIVIAL +// p1,p2,p3 are only allowed to be 2,3,4,6 +void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3, mpq_t s, mpq_t q) +{ + mpq_t rho1, rho2, rho3; + mpq_inits(rho1,rho2,rho3,NULL); + + // rho_i = s^2 + 2s cos(2 pi / p_i) + 1 + // coefficient 2 is the value for p=infinity, not sure if that would even work + quartic(rho1, s, 0, 0, 1, p1 == 2 ? -2 : p1 == 3 ? -1 : p1 == 4 ? 0 : p1 == 6 ? 1 : 2, 1); + quartic(rho2, s, 0, 0, 1, p2 == 2 ? -2 : p2 == 3 ? -1 : p2 == 4 ? 0 : p2 == 6 ? 1 : 2, 1); + quartic(rho3, s, 0, 0, 1, p3 == 2 ? -2 : p3 == 3 ? -1 : p3 == 4 ? 0 : p3 == 6 ? 1 : 2, 1); + + generators_triangle_rotation_generic(gen, rho1, rho2, rho3, s, q); + + mpq_clears(rho1,rho2,rho3,NULL); +} +#endif + +#ifdef QEXT_SQRT5 +void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s_, mpq_t q_) +{ + NUMBER rho, c, one, s, q; + INIT(rho);INIT(c);INIT(one);INIT(s);INIT(q); + + // c = - (1+sqrt(5))/2 + mpq_set_si(c[0], -1, 2); + mpq_set_si(c[1], -1, 2); + SET_ONE(one); + + mpq_set(s[0], s_); + mpq_set_ui(s[1], 0, 1); + mpq_set(q[0], q_); + mpq_set_ui(q[1], 0, 1); + + SET(rho, one); + MUL(rho, rho, s); + ADD(rho, rho, c); + MUL(rho, rho, s); + ADD(rho, rho, one); + + generators_triangle_rotation_generic(gen, rho, rho, rho, s, q); + + CLEAR(rho);CLEAR(c);CLEAR(one);CLEAR(s);CLEAR(q); +} +#endif + char *print_word(groupelement_t *g, char *str) { int i = g->length - 1; @@ -213,20 +280,17 @@ char *print_word(groupelement_t *g, char *str) return str; } -void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, mpq_t q) +void enumerate_triangle_rotation_subgroup(group_t *group, mat *gen, mat *matrices) { mat_workspace *ws; mat tmp; - mat gen[6]; char buf[100], buf2[100], buf3[100]; // allocate stuff ws = mat_workspace_init(3); - for(int i = 0; i < 6; i++) - mat_init(gen[i], 3); mat_init(tmp, 3); - initialize_triangle_generators(ws, gen, p1, p2, p3, s, q); +// initialize_triangle_generators(ws, gen, p1, p2, p3, s, q); mat_identity(matrices[0]); for(int i = 1; i < group->size; i++) { @@ -258,8 +322,6 @@ void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, m } // free stuff - for(int i = 0; i < 6; i++) - mat_clear(gen[i]); mat_clear(tmp); mat_workspace_clear(ws); } diff --git a/enumerate_triangle_group.h b/enumerate_triangle_group.h index 7488861..47b46fc 100644 --- a/enumerate_triangle_group.h +++ b/enumerate_triangle_group.h @@ -6,12 +6,25 @@ #include +// rational only functions int solve_characteristic_polynomial(mps_context *solv, mps_monomial_poly *poly, mpq_t tr, mpq_t trinv, double *eigenvalues); void continued_fraction_approximation(mpq_t out, double in, int level); -void quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e); -// p1,p2,p3 are only allowed to be 2,3,4,6 -void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2, int p3, mpq_t s, mpq_t q); + +// geneartors +void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER s, NUMBER q); + +#ifdef QEXT_TRIVIAL +void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3, mpq_t s, mpq_t q); +#endif + +#ifdef QEXT_SQRT5 +void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s, mpq_t q); +#endif + +// general functions +void quartic(NUMBER out, NUMBER in, int a, int b, int c, int d, int e); +void discriminant(NUMBER out, NUMBER x, NUMBER y); char *print_word(groupelement_t *g, char *str); -void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, mpq_t t); +void enumerate_triangle_rotation_subgroup(group_t *group, mat *gen, mat *matrices); #endif diff --git a/mat.c b/mat.c index 7ed0b5a..391ea22 100644 --- a/mat.c +++ b/mat.c @@ -78,7 +78,7 @@ static void mat_multiply_outofplace(mat_workspace *ws, mat out, mat in1, mat in2 LOOP(i,n) LOOP(j,n) { SET_ZERO(M(out,i,j)); LOOP(k,n) { - MULTIPLY(*tmp, M(in1,i,k), M(in2,k,j)); + MUL(*tmp, M(in1,i,k), M(in2,k,j)); ADD(M(out,i,j), M(out,i,j), *tmp); } } @@ -103,12 +103,12 @@ void mat_det(mat_workspace *ws, NUMBER out, mat in) SET_ZERO(out); LOOP(i,n) { - MULTIPLY(*tmp, M(in,0,i), M(in,1,(i+1)%3)); - MULTIPLY(*tmp, *tmp, M(in,2,(i+2)%3)); + MUL(*tmp, M(in,0,i), M(in,1,(i+1)%3)); + MUL(*tmp, *tmp, M(in,2,(i+2)%3)); ADD(out, out, *tmp); - MULTIPLY(*tmp, M(in,0,(i+2)%3), M(in,1,(i+1)%3)); - MULTIPLY(*tmp, *tmp, M(in,2,i)); + MUL(*tmp, M(in,0,(i+2)%3), M(in,1,(i+1)%3)); + MUL(*tmp, *tmp, M(in,2,i)); SUB(out, out, *tmp); } } @@ -121,8 +121,8 @@ static void mat_pseudoinverse_outofplace(mat_workspace *ws, mat out, mat in) int n = 3; LOOP(i,n) LOOP(j,n) { - MULTIPLY(*tmp, M(in,(i+1)%3,(j+1)%3), M(in,(i+2)%3,(j+2)%3)); - MULTIPLY(*tmp2, M(in,(i+1)%3,(j+2)%3), M(in,(i+2)%3,(j+1)%3)); + MUL(*tmp, M(in,(i+1)%3,(j+1)%3), M(in,(i+2)%3,(j+2)%3)); + MUL(*tmp2, M(in,(i+1)%3,(j+2)%3), M(in,(i+2)%3,(j+1)%3)); SUB(M(out,j,i), *tmp, *tmp2); } } diff --git a/mat.h b/mat.h index 5c548e3..85ac648 100644 --- a/mat.h +++ b/mat.h @@ -8,27 +8,139 @@ // library for matrix computations in variable rings (based on GMP types) -/* - needed features: - x multiply matrices - - inverse - x pseudoinverse - x set - - eigenvalues - */ +#ifdef QEXT_SQRT5 + +typedef mpq_t NUMBER [2]; +#define INIT(x) do { mpq_init((x)[0]); mpq_init((x)[1]); } while(0) +#define CLEAR(x) do { mpq_clear((x)[0]); mpq_clear((x)[1]); } while(0) +#define SET(x,y) do { mpq_set((x)[0],(y)[0]); mpq_set((x)[1],(y)[1]); } while(0) +#define SET_INT(x,y) do { mpq_set_si((x)[0],(y),1); mpq_set_si((x)[1],0,1); } while(0) +#define SET_ZERO(x) SET_INT(x,0) +#define SET_ONE(x) SET_INT(x,1) +#define ADD(x,y,z) do { mpq_add((x)[0],(y)[0],(z)[0]); mpq_add((x)[1],(y)[1],(z)[1]); } while(0) +#define SUB(x,y,z) do { mpq_sub((x)[0],(y)[0],(z)[0]); mpq_sub((x)[1],(y)[1],(z)[1]); } while(0) +#define MUL multiply_sqrt5 +#define DIV divide_sqrt5 +#define CMP compare_sqrt5 +#define PRINT(x) gmp_printf("%Qd + %Qd*sqrt(5)", (x)[0], (x)[1]) +#define SPRINT(buf, x) gmp_sprintf(buf, "%Qd+%Qd*sqrt(5)", (x)[0], (x)[1]) + +static void multiply_sqrt5(NUMBER out, NUMBER a, NUMBER b) +{ + // could be in place!!! + + mpq_t tmp, result[2]; + mpq_inits(tmp, result[0], result[1], 0); + + mpq_mul(result[0], a[1], b[1]); + mpq_set_si(tmp, 5, 1); + mpq_mul(result[0], result[0], tmp); + mpq_mul(tmp, a[0], b[0]); + mpq_add(result[0], result[0], tmp); + + mpq_mul(result[1], a[1], b[0]); + mpq_mul(tmp, a[0], b[1]); + mpq_add(result[1], result[1], tmp); + + mpq_set(out[0], result[0]); + mpq_set(out[1], result[1]); + + mpq_clears(tmp, result[0], result[1], 0); +} + +static void divide_sqrt5(NUMBER out, NUMBER a, NUMBER b) +{ + mpq_t denom, num, tmp, neg5, result[2]; + mpq_inits(denom, num, tmp, neg5, result[0], result[1], 0); + mpq_set_si(neg5, -5, 1); + + mpq_mul(denom, b[1], b[1]); + mpq_mul(denom, denom, neg5); + mpq_mul(tmp, b[0], b[0]); + mpq_add(denom, denom, tmp); + + mpq_mul(num, a[1], b[1]); + mpq_mul(num, num, neg5); + mpq_mul(tmp, a[0], b[0]); + mpq_add(num, num, tmp); + + mpq_div(result[0], num, denom); + + mpq_mul(num, a[1], b[0]); + mpq_mul(tmp, a[0], b[1]); + mpq_sub(num, num, tmp); + + mpq_div(result[1], num, denom); + + mpq_set(out[0], result[0]); + mpq_set(out[1], result[1]); + + mpq_clears(denom, num, tmp, neg5, result[0], result[1], 0); +} + +static int compare_sqrt5(NUMBER x, NUMBER y) +{ + int result; + mpq_t p, q, p2, q2, c5; + mpq_inits(p, q, p2, q2, c5, 0); + + mpq_sub(p, x[0], y[0]); + mpq_sub(q, x[1], y[1]); + + /* + want to know if p + sqrt(5) q > 0 + if p>0 and q>0: always true + if p>0 and q<0: equivalent to |p|^2 > 5 |q|^2 + if p<0 and q>0: equivalent to |p|^2 < 5 |q|^2 + if p<0 and q<0: always false + */ + + if(mpq_sgn(p) > 0 && mpq_sgn(q) > 0) { + result = 1; + goto done; + } + + if(mpq_sgn(p) < 0 && mpq_sgn(q) < 0) { + result = -1; + goto done; + } + + mpq_mul(p2, p, p); + mpq_mul(q2, q, q); + mpq_set_si(c5, 5, 1); + mpq_mul(q2, q2, c5); + + if(mpq_sgn(p) > 0) + result = mpq_cmp(p2, q2); // this can't be zero, or else |p/q| = sqrt(5), but it's irrational + else if(mpq_sgn(p) < 0) + result = -mpq_cmp(p2, q2); // this can be zero if p = 0 and q = 0 + else // p = 0 + return mpq_sgn(q); + +done: + mpq_clears(p, q, p2, q2, c5, 0); + return result; +} + +#endif + +#ifdef QEXT_TRIVIAL #define NUMBER mpq_t #define INIT mpq_init #define CLEAR mpq_clear #define SET mpq_set -#define SET_ZERO(x) mpq_set_ui(x,0,1) -#define SET_ONE(x) mpq_set_ui(x,1,1) +#define SET_INT(x,y) mpq_set_si(x,y,1) +#define SET_ZERO(x) SET_INT(x,0) +#define SET_ONE(x) SET_INT(x,1) #define ADD mpq_add #define SUB mpq_sub -#define MULTIPLY mpq_mul +#define MUL mpq_mul #define DIV mpq_div #define PRINT(x) gmp_printf("%Qd", x) +#endif + #define M(m,i,j) ((m)->x[(i)+(m)->n*(j)]) struct _mat{ diff --git a/singular_values_barbot.c b/singular_values_barbot.c new file mode 100644 index 0000000..af49dff --- /dev/null +++ b/singular_values_barbot.c @@ -0,0 +1,545 @@ +#include "coxeter.h" +#include "linalg.h" +#include "mat.h" +#include "enumerate_triangle_group.h" +#include "parallel.h" + +#include + +#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0); + +//#define DEBUG(msg, ...) fprintf(stderr, "[%003d%10.3f] " msg, mpi_rank(0), runtime(), ##__VA_ARGS__) +#define DEBUG(msg, ...) +#define INFO(msg, ...) fprintf(stderr, "[%003d%10.3f] " msg, mpi_rank(0), runtime(), ##__VA_ARGS__) + +struct result { + int id; + NUMBER tr; + NUMBER trinv; + int disc_sign; +}; + +// we want as much as possible to be node data, except if it is only known to the main node +// (command line arguments) or should only be computed once (id list) + +struct global_data { + // command line arguments + unsigned int nmax; + unsigned int p1, p2, p3; + unsigned int sstart, send, sdenom; + unsigned int qstart, qend, qdenom; + + unsigned int *id_list; + unsigned int id_list_length; +}; + +struct node_data { + group_t *group; + mat* matrices; + struct result *invariants; + struct result **distinct_invariants; + int distinct_invariants_length; + mps_context *solver; +}; + +struct input_data { + unsigned int snum, sden; + unsigned int qnum, qden; +}; + +struct output_data { + int has_non_loxodromic; +}; + +static int compare_result(const void *a_, const void *b_) +{ + int d = 0; + + struct result **a = (struct result **)a_; + struct result **b = (struct result **)b_; + + d = CMP((*a)->tr,(*b)->tr); + if(d == 0) { + d = CMP((*a)->trinv, (*b)->trinv); + } + + return d; +} + +static int compare_result_by_id(const void *a_, const void *b_) +{ + int d = 0; + + struct result **a = (struct result **)a_; + struct result **b = (struct result **)b_; + + return (*a)->id - (*b)->id; +} + +static int compare_result_by_tr_trinv_id(const void *a_, const void *b_) +{ + int d = 0; + + struct result **a = (struct result **)a_; + struct result **b = (struct result **)b_; + + d = CMP((*a)->tr,(*b)->tr); + if(d == 0) { + d = CMP((*a)->trinv, (*b)->trinv); + if(d == 0) { + d = (*b)->id - (*a)->id; + } + } + + return d; +} + +/* +static int compare_result_by_slope(const void *a_, const void *b_) +{ + int d = 0; + + struct result **a = (struct result **)a_; + struct result **b = (struct result **)b_; + + double slopea = (*a)->x / (*a)->y; + double slopeb = (*b)->x / (*b)->y; + + return slopea > slopeb ? -1 : slopea < slopeb ? 1 : 0; +} +*/ + +int invariants_trace_loxodromic(group_t *group, mat *matrices, struct result **invariants, int *n, int unique) +{ + int ntraces = *n, nuniq; // ntraces is the number of traces we are asked to compute, nuniq is the number of unique ones after we eliminate duplicates + + // compute the traces + for(int i = 0; i < ntraces; i++) { + int id = invariants[i]->id; + int invid = group->elements[id].inverse->id; + mat_trace(invariants[i]->tr, matrices[id]); + mat_trace(invariants[i]->trinv, matrices[invid]); + } + + // throw out duplicates if unique == 1 + if(!unique) + nuniq = ntraces; + else { + qsort(invariants, ntraces, sizeof(struct result*), compare_result); + + nuniq = 0; + for(int i = 0; i < ntraces; i++) { + if(i == 0 || compare_result(&invariants[i], &invariants[nuniq-1]) != 0) { + invariants[nuniq] = invariants[i]; + nuniq++; + } else { + int oldlength = group->elements[invariants[nuniq-1]->id].length; + int newlength = group->elements[invariants[i]->id].length; + if(newlength < oldlength) + invariants[nuniq-1]->id = invariants[i]->id; + } + } + } + + // check if loxodromic + NUMBER disc, zero; + INIT(disc); + INIT(zero); + SET_ZERO(zero); + for(int i = 0; i < nuniq; i++) { + discriminant(disc, invariants[i]->tr, invariants[i]->trinv); + invariants[i]->disc_sign = CMP(disc, zero); + } + CLEAR(disc); + CLEAR(zero); + + // sort by ID again + qsort(invariants, nuniq, sizeof(struct result*), compare_result_by_id); + + *n = nuniq; + return 0; +} + +/* +int invariants_trace_slope(group_t *group, mat *matrices, struct result **invariants, int *n, int unique) +{ + mpq_t tmp; + mps_context *solver; + mps_monomial_poly *poly; + int index; + int ntraces = *n, nuniq; + int retval; + double evs[3]; + char buf[100]; + +// DEBUG("Compute traces\n"); + for(int i = 0; i < ntraces; i++) { + int id = invariants[i]->id; + int invid = group->elements[id].inverse->id; + mat_trace(invariants[i]->tr, matrices[id]); + mat_trace(invariants[i]->trinv, matrices[invid]); + } + + if(!unique) + nuniq = ntraces; + else { +// DEBUG("Get unique traces\n"); + qsort(invariants, ntraces, sizeof(struct result*), compare_result); + + nuniq = 0; + for(int i = 0; i < ntraces; i++) { + if(i == 0 || compare_result(&invariants[i], &invariants[nuniq-1]) != 0) { + invariants[nuniq] = invariants[i]; + nuniq++; + } else { + int oldlength = group->elements[invariants[nuniq-1]->id].length; + int newlength = group->elements[invariants[i]->id].length; + if(newlength < oldlength) + invariants[nuniq-1]->id = invariants[i]->id; + } + } + } + + DEBUG("Solve characteristic polynomials\n"); + solver = mps_context_new(); + 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); + + for(int i = 0; i < nuniq; i++) { + retval = solve_characteristic_polynomial(solver, poly, invariants[i]->tr, 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])); + + invariants[i]->x = x; + invariants[i]->y = y; + invariants[i]->slope = y/x; + } + mps_context_free(solver); + + qsort(invariants, nuniq, sizeof(struct result*), compare_result_by_id); + + *n = nuniq; + return 0; +} +*/ + +/* +long check_memory_usage(mat *matrices, int n) +{ + mpq_t x; + long total; + + for(int i = 0; i < n; i++) + { + LOOP(j,3) LOOP(k,3) { + total += mpq_numref(M(matrices[i], j, k))->_mp_size; + total += mpq_denref(M(matrices[i], j, k))->_mp_size; + } + } + + return total; +} +*/ + + +void destroy_node(const void *_g, void *_n) +{ + struct global_data *g = (struct global_data *)_g; + struct node_data *n = (struct node_data *)_n; + + for(int i = 0; i < g->nmax; i++) { + CLEAR(n->invariants[i].tr); + CLEAR(n->invariants[i].trinv); + } + free(n->invariants); + free(n->distinct_invariants); + for(int i = 0; i < g->nmax; i++) + mat_clear(n->matrices[i]); + free(n->matrices); + coxeter_clear(n->group); +} + +int init_node(const void *_g, void *_n) +{ + struct global_data *g = (struct global_data *)_g; + struct node_data *n = (struct node_data *)_n; + + DEBUG("Allocate\n"); + g->id_list = (int*)(g+1); // pointers get scrambled by transmission, reconstruct + n->matrices = malloc(g->nmax*sizeof(mat)); + for(int i = 0; i < g->nmax; i++) + mat_init(n->matrices[i], 3); + n->invariants = malloc(g->nmax*sizeof(struct result)); + n->distinct_invariants = malloc(g->nmax*sizeof(struct result)); // we won't need that many, but just in case + for(int i = 0; i < g->nmax; i++) { + INIT(n->invariants[i].tr); + INIT(n->invariants[i].trinv); + n->invariants[i].id = i; + } + + // order of the triangle reflection generators: a, b, c + // order of the rotation orders: bc, ac, ab + DEBUG("Generate group\n"); + n->group = coxeter_init_triangle(g->p1, g->p2, g->p3, g->nmax); + + return 0; +} + +int process_output(group_t *group, mat *matrices, struct result **invariants, int invariants_length, struct output_data *out) +{ + out->has_non_loxodromic = 0; + for(int i = 0; i < invariants_length; i++) { + if(invariants[i]->disc_sign <= 0 && invariants[i]->id != 0 && invariants[i]->id != 4 && invariants[i]->id != 22) { + out->has_non_loxodromic = 1; + } + } +} + +int do_computation(const void *_g, void *_n, const void *_in, void *_out) +{ + struct global_data *g = (struct global_data *)_g; + struct node_data *n = (struct node_data *)_n; + struct input_data *in = (struct input_data *)_in; + struct output_data *out = (struct output_data *)_out; + + mpq_t s, q; + + mpq_inits(s, q, NULL); + mpq_set_ui(s, in->snum, in->sden); + mpq_set_ui(q, in->qnum, in->qden); + + INFO("Computing represention with s = %d/%d and q = %d/%d.\n", + in->snum, in->sden, + in->qnum, in->qden); + + // we need to compute all the elements pointed to in id_list, and all their suffixes or prefixes + // I can imagine a smarter way of doing this which checks if there is a shorter route to the element + for(int i = 0; i < n->group->size; i++) + n->group->elements[i].need_to_compute = 0; + n->group->elements[0].need_to_compute = 1; + + int needed_elements = 1; + for(int i = 0; i < g->id_list_length; i++) + { + int id = g->id_list[i]; + n->distinct_invariants[i] = &n->invariants[id]; + groupelement_t *cur = &n->group->elements[id]; + while(cur->need_to_compute == 0) { + cur->need_to_compute = 1; + needed_elements++; + cur = cur->parent->parent; // also need to compute its even-length ancestors + } + cur = n->group->elements[id].inverse; + while(cur->need_to_compute == 0) { + cur->need_to_compute = 1; + needed_elements++; + cur = cur->parent->parent; + } + } + n->distinct_invariants_length = g->id_list_length; + + DEBUG("Need to compute %d elements to get %d traces up to reflection length %d\n", + needed_elements, g->id_list_length, n->group->elements[n->group->size-1].length); + + DEBUG("Compute matrices\n"); + mat gen[6]; + for(int i = 0; i < 6; i++) + mat_init(gen[i], 3); + + generators_triangle_rotation_555_barbot(gen, s, q); + enumerate_triangle_rotation_subgroup(n->group, gen, n->matrices); + + for(int i = 0; i < 6; i++) + mat_clear(gen[i]); + + DEBUG("Compute invariants\n"); + invariants_trace_loxodromic( + n->group, n->matrices, + n->distinct_invariants, &n->distinct_invariants_length, 1); + +// DEBUG("Find max slopes\n"); + process_output(n->group, n->matrices, n->distinct_invariants, n->distinct_invariants_length, out); + + mpq_clears(s, q, NULL); + + return 0; +} + +int main(int argc, char *argv[]) +{ + char buf[1000]; + char buf2[1000]; + char buf3[1000]; + struct global_data *g; + struct node_data n; + + start_timer(); + + // parse command line arguments + if(argc < 11) { + fprintf(stderr, "Usage: %s [restart]\n", argv[0]); + exit(1); + } + int nmax = atoi(argv[1]); + g = (struct global_data*)malloc(sizeof(struct global_data) + nmax*sizeof(int)); + g->id_list = (int*)(g+1); + g->nmax = nmax; + g->p1 = atoi(argv[2]); + g->p2 = atoi(argv[3]); + g->p3 = atoi(argv[4]); + g->sstart = atoi(argv[5]); + g->send = atoi(argv[6]); + g->sdenom = atoi(argv[7]); + g->qstart = atoi(argv[8]); + g->qend = atoi(argv[9]); + g->qdenom = atoi(argv[10]); + + // initialize + parallel_context *ctx = parallel_init(); + parallel_set_datasize_and_callbacks(ctx, init_node, do_computation, destroy_node, + sizeof(struct global_data) + g->nmax*sizeof(int), + sizeof(struct node_data), + sizeof(struct input_data), + sizeof(struct output_data)); + if(ctx->mpi_mode == 1 && ctx->rank != 0) { + // worker mode + parallel_work(ctx); + parallel_destroy(ctx); + exit(0); + } + + init_node(g, &n); + + // use very generic values for the pilot run unless sstart=send and qstart=qend + struct input_data pilot_input; + struct output_data pilot_output; + if(g->sstart == g->send && g->qstart == g->qend) { + pilot_input.snum = g->sstart; + pilot_input.sden = g->sdenom; + pilot_input.qnum = g->qstart; + pilot_input.qden = g->qdenom; + DEBUG("Single run for s = %d/%d, q = %d/%d\n", g->sstart, g->sdenom, g->qstart, g->qdenom); + } else { + pilot_input.snum = 4; + pilot_input.sden = 100; + pilot_input.qnum = 7; + pilot_input.qden = 100; + DEBUG("Initial run for s = %d/%d, q = %d/%d\n", 4, 100, 7, 100); + } + + g->id_list_length = 0; + for(int i = 0; i < n.group->size; i++) + if(n.group->elements[i].length % 2 == 0 && n.group->elements[i].inverse) + g->id_list[g->id_list_length++] = i; + + do_computation(g, &n, &pilot_input, &pilot_output); + + for(int i = 0; i < n.distinct_invariants_length; i++) + g->id_list[i] = n.distinct_invariants[i]->id; + g->id_list_length = n.distinct_invariants_length; + + if(g->sstart != g->send || g->qstart != g->qend) { + + struct input_data *inputs = malloc((g->send - g->sstart + 1)*(g->qend - g->qstart + 1)*sizeof(struct input_data)); + struct output_data *outputs = malloc((g->send - g->sstart + 1)*(g->qend - g->qstart + 1)*sizeof(struct input_data)); + + int njobs = 0; + for(int sloop = g->sstart; sloop <= g->send; sloop++) { + for(int qloop = g->qstart; qloop <= g->qend; qloop++) { + inputs[njobs].sden = g->sdenom; + inputs[njobs].qden = g->qdenom; + inputs[njobs].snum = sloop; + inputs[njobs].qnum = qloop; + njobs++; + } + } + + if(argc >= 12) + parallel_run(ctx, g, inputs, outputs, njobs, argv[11]); + else + parallel_run(ctx, g, inputs, outputs, njobs, NULL); + +// DEBUG("Loop for s = %d/%d, q = %d/%d\n", sloop, g->sdenom, qloop, g->qdenom); + + for(int i = 0; i < njobs; i++) + { + /* + gmp_printf("%d/%d %d/%d %d %s %f\n", + inputs[i].snum, inputs[i].sden, inputs[i].qnum, inputs[i].qden, + outputs[i].max_slope_id, + print_word(&n.group->elements[outputs[i].max_slope_id], buf), + outputs[i].max_slope); + */ + + + printf("%d/%d %d/%d %d\n", + inputs[i].snum, inputs[i].sden, inputs[i].qnum, inputs[i].qden, + outputs[i].has_non_loxodromic); + } + + free(inputs); + free(outputs); + } else { + // output + for(int i = 0; i < n.distinct_invariants_length; i++) { + // exclude tr = trinv = 2/1/0/-1/3 + + /* + mpq_t tmp; + mpq_init(tmp); + mpq_set_si(tmp, 2, 1); + if(mpq_cmp(n.distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(n.distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, 1, 1); + if(mpq_cmp(n.distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(n.distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, 0, 1); + if(mpq_cmp(n.distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(n.distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, -1, 1); + if(mpq_cmp(n.distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(n.distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, 3, 1); + if(mpq_cmp(n.distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(n.distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_clear(tmp); + */ + + SPRINT(buf, n.distinct_invariants[i]->tr); + SPRINT(buf2, n.distinct_invariants[i]->trinv); + printf("%d %s %f %f %d\n", + n.distinct_invariants[i]->id, + print_word(&n.group->elements[n.distinct_invariants[i]->id], buf3), + mpq_get_d(n.distinct_invariants[i]->tr[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->tr[1]), + mpq_get_d(n.distinct_invariants[i]->trinv[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->trinv[1]), + n.distinct_invariants[i]->disc_sign); + + } + + } + + destroy_node(g, &n); + free(g); + parallel_destroy(ctx); +}