From 2b384447ca40e33a4e1e8dd889c4a0ec92e0b196 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sun, 30 Jan 2022 19:09:56 -0600 Subject: [PATCH] missing files from previous commit + clean up of special_element program + billiard picture script --- billiard_picture.sh | 20 +++ billiard_words.hs | 21 ++- enumerate_triangle_group.c | 254 +++++++++++++++++++++++++++++++++++++ enumerate_triangle_group.h | 17 +++ special_element.c | 173 ++++++++++++------------- 5 files changed, 394 insertions(+), 91 deletions(-) create mode 100755 billiard_picture.sh create mode 100644 enumerate_triangle_group.c create mode 100644 enumerate_triangle_group.h diff --git a/billiard_picture.sh b/billiard_picture.sh new file mode 100755 index 0000000..be0a731 --- /dev/null +++ b/billiard_picture.sh @@ -0,0 +1,20 @@ +#!/bin/bash + +trap 'exit 130' INT + +wordlength=30 +sdenom=100 +sstart=1 +send=100 +qdenom=100 +qstart=1 +qend=100 # 1/sqrt(2) = 0.7071... + +words="$(./billiard_words $wordlength | awk '{print $1}')" + +for s in $(seq $sstart $send); do + for q in $(seq $qstart $qend); do + echo -n "$s/$sdenom $q/$qdenom " + MAXIMUM=only ./special_element $s/$sdenom $q/$qdenom $words + done +done diff --git a/billiard_words.hs b/billiard_words.hs index 9e88f51..2ae728c 100644 --- a/billiard_words.hs +++ b/billiard_words.hs @@ -9,14 +9,23 @@ main = do listWordsUpToLength :: Int -> IO () listWordsUpToLength n = do - putStrLn $ unlines [printf "%s %d/%d %f" w (x `div` gcd x y) (y `div` gcd x y) (fromIntegral x / fromIntegral y :: Double) | ((p,q),w) <- wordlist (n`div`2) (n`div`2), length w <= n, let x = 2*q + p, let y = 2*p + q] --- putStrLn $ unlines [printf "%d/%d\t%d/%d\t%.7f\t%d\t%s" p q (x `div` gcd x y) (y `div` gcd x y) (sqrt 3 / (1 + 2*fromIntegral q / fromIntegral p) :: Double) (length w) w | ((p,q),w) <- wordlist (n`div`2) (n`div`2), length w <= n, let x = 2*q + p, let y = 2*p + q] - - - -- putStrLn $ unlines [printf "%d/%d\t%.5f\t%.5f\t%d\t%s" p q (fromIntegral p / fromIntegral q :: Double) (sqrt 3 / (1 + 2*fromIntegral q / fromIntegral p) :: Double) (length w) w | ((p,q),w) <- wordlist (n`div`2) (n`div`2), length w <= n] + putStrLn $ unlines [printf "%s %d/%d %f" + w + (x `div` gcd x y) + (y `div` gcd x y) + (fromIntegral x / fromIntegral y :: Double) | + ((p,q),w) <- wordlist (n `div` 2) (n `div` 2), + length w <= n, + let x = 2*q + p, + let y = 2*p + q] wordlist :: Int -> Int -> [((Int,Int),String)] -wordlist pmax qmax = nub $ sortBy (comparing sl) [((p `div` gcd p q, q `div` gcd p q), slopeWord "bca" p q) | p <- [0..pmax], q <- [0..qmax], p /= 0 || q /= 0] +wordlist pmax qmax = nub $ + sortBy (comparing sl) + [((p `div` gcd p q, q `div` gcd p q), slopeWord "bca" p q) | + p <- [0..pmax], + q <- [0..qmax], + q /= 0] -- use p /= 0 || q /= 0 for more symmetric output where sl ((p,q),_) = fromIntegral p / fromIntegral q diff --git a/enumerate_triangle_group.c b/enumerate_triangle_group.c new file mode 100644 index 0000000..6262f4b --- /dev/null +++ b/enumerate_triangle_group.c @@ -0,0 +1,254 @@ +#include "enumerate_triangle_group.h" +#include "linalg.h" + +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; + + if(in < 0) { + mpq_init(tmp); + mpq_set_ui(tmp, 0, 1); + continued_fraction_approximation(out, -in, level); + mpq_sub(out, tmp, out); + mpq_clear(tmp); + return; + } + + if(level == 0) { + mpq_set_si(out, (signed long int)round(in), 1); // floor(in) + } else { + continued_fraction_approximation(out, 1/(in - floor(in)), level - 1); + mpq_init(tmp); + mpq_set_ui(tmp, 1, 1); + mpq_div(out, tmp, out); // out -> 1/out + mpq_set_si(tmp, (signed long int)in, 1); // floor(in) + mpq_add(out, out, tmp); + mpq_clear(tmp); + } +} + +void quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e) +{ + mpq_t tmp; + mpq_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); + + mpq_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) +{ + mat r1,r2,r3; + mpq_t rho1, rho2, rho3; + mpq_t b1,c1,a2,c2,a3,b3; + mpq_t sinv; + + 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_set_ui(sinv, 1, 1); + mpq_div(sinv, sinv, s); + + // 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); + + 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); + + // 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); + + 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]); + + /* + 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_clear(r1); + mat_clear(r2); + mat_clear(r3); +} + +char *print_word(groupelement_t *g, char *str) +{ + int i = g->length - 1; + + str[g->length] = 0; + while(g->parent) { + str[i--] = 'a' + g->letter; + g = g->parent; + } + + return str; +} + +void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, mpq_t t) +{ + 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, t); + + mat_identity(matrices[0]); + for(int i = 1; i < group->size; i++) { + if(group->elements[i].length % 2 != 0) + continue; + if(!group->elements[i].inverse) + continue; + + int parent = group->elements[i].parent->id; + int grandparent = group->elements[i].parent->parent->id; + int letter; + + if(group->elements[parent].letter == 1 && group->elements[i].letter == 2) + letter = 0; // p = bc + else if(group->elements[parent].letter == 2 && group->elements[i].letter == 0) + letter = 1; // q = ca + else if(group->elements[parent].letter == 0 && group->elements[i].letter == 1) + letter = 2; // r = ab + if(group->elements[parent].letter == 2 && group->elements[i].letter == 1) + letter = 3; // p^{-1} = cb + else if(group->elements[parent].letter == 0 && group->elements[i].letter == 2) + letter = 4; // q^{-1} = ac + else if(group->elements[parent].letter == 1 && group->elements[i].letter == 0) + letter = 5; // r^{-1} = ba + + mat_multiply(ws, matrices[i], matrices[grandparent], gen[letter]); + } + + // 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 new file mode 100644 index 0000000..5d8a85e --- /dev/null +++ b/enumerate_triangle_group.h @@ -0,0 +1,17 @@ +#ifndef ENUMERATE_TRIANGLE_GROUP_H +#define ENUMERATE_TRIANGLE_GROUP_H + +#include "mat.h" +#include "coxeter.h" + +#include + +int solve_characteristic_polynomial(mps_context *solv, 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); +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); + +#endif diff --git a/special_element.c b/special_element.c index 98c60b0..87ea78f 100644 --- a/special_element.c +++ b/special_element.c @@ -5,11 +5,6 @@ #define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0); #define DEBUG(msg, ...) -#define DENOMINATOR 300 -#define WIDTH 135 -#define STARTX 121 -#define HEIGHT 300 -#define IDX(i,j) (((i)-1)*HEIGHT + ((j)-1)) double mpq_log(mpq_t m_op) { @@ -45,10 +40,25 @@ int main(int argc, char *argv[]) int retval; double evs[3]; char buf[100]; - double *max_slope; - int *max_slope_index; + double max_slope = 0; + int max_slope_index = 0; + double min_slope = INFINITY; + int min_slope_index = 0; + char *env; + int mode; - DEBUG("Allocate\n"); + if(argc < 2) { + fprintf(stderr, + "Usage: %s ...\n" + "Computes jordan slopes of a list of group elements for a fixed representation.\n" + "s,q: representation in the Hitchin component, given as rational numbers, e.g. 2/7\n" + "word1, word2, ...: elements in the triangle rotation group, given as reflection group words\n" + "output: word - jordan slope pairs\n" + "+max slope index, max slope value, max slope word, min slope index, min slope value, min slope word\n" + "controlled by environment variable MAXIMUM=no/yes/only, default yes\n", + argv[0]); + exit(0); + } mpq_inits(m, t, s, q, tmp, tmp2, tr, trinv, NULL); ws = mat_workspace_init(3); @@ -56,10 +66,6 @@ int main(int argc, char *argv[]) mat_init(gen[i], 3); mat_init(element, 3); mat_init(inverse, 3); - max_slope = malloc(sizeof(double)*WIDTH*HEIGHT); - max_slope_index = malloc(sizeof(int)*WIDTH*HEIGHT); - memset(max_slope_index, 0, sizeof(int)*WIDTH*HEIGHT); - memset(max_slope, 0, sizeof(int)*WIDTH*HEIGHT); solver = mps_context_new(); mps_context_set_output_prec(solver, 20); // relative precision @@ -68,92 +74,91 @@ int main(int argc, char *argv[]) mpq_set_str(s, argv[1], 10); mpq_set_str(q, argv[2], 10); -// for(int i = STARTX; i <= WIDTH; i++) { -// for(int j = 1; j <= HEIGHT; j++) { - for(int w = 3; w < argc; w++) { -// mpq_set_ui(t, j, DENOMINATOR); -// mpq_set_ui(m, i, DENOMINATOR); // 414/1000 ~ sqrt(2)-1 <-> s=1 -// s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m)); + env = getenv("MAXIMUM"); + if(!env || strcmp(env, "yes") == 0) { + mode = 1; // yes + } else if(strcmp(env, "no") == 0) { + mode = 0; // no + } else if(strcmp(env, "only") == 0) { + mode = 2; // only + } - DEBUG("Compute matrix\n"); - initialize_triangle_generators(ws, gen, 6, 4, 3, s, q); + for(int w = 0; w < argc - 3; w++) { + initialize_triangle_generators(ws, gen, 4, 4, 4, s, q); - mat_identity(element); - mat_identity(inverse); - for(int k = 0; k < strlen(argv[w]); k+=2) { - letter1 = argv[w][k] - 'a'; - letter2 = argv[w][k+1] - 'a'; + mat_identity(element); + mat_identity(inverse); + for(int k = 0; k < strlen(argv[w+3]); k+=2) { + letter1 = argv[w+3][k] - 'a'; + letter2 = argv[w+3][k+1] - 'a'; - if(letter1 == 1 && letter2 == 2) - letter = 0; // p = bc - else if(letter1 == 2 && letter2 == 0) - letter = 1; // q = ca - else if(letter1 == 0 && letter2 == 1) - letter = 2; // r = ab - else if(letter1 == 2 && letter2 == 1) - letter = 3; // p^{-1} = cb - else if(letter1 == 0 && letter2 == 2) - letter = 4; // q^{-1} = ac - else if(letter1 == 1 && letter2 == 0) - letter = 5; // r^{-1} = ba + if(letter1 == 1 && letter2 == 2) + letter = 0; // p = bc + else if(letter1 == 2 && letter2 == 0) + letter = 1; // q = ca + else if(letter1 == 0 && letter2 == 1) + letter = 2; // r = ab + else if(letter1 == 2 && letter2 == 1) + letter = 3; // p^{-1} = cb + else if(letter1 == 0 && letter2 == 2) + letter = 4; // q^{-1} = ac + else if(letter1 == 1 && letter2 == 0) + letter = 5; // r^{-1} = ba - mat_multiply(ws, element, element, gen[letter]); - mat_multiply(ws, inverse, gen[(letter+3)%6], inverse); - } + mat_multiply(ws, element, element, gen[letter]); + mat_multiply(ws, inverse, gen[(letter+3)%6], inverse); + } - DEBUG("Compute traces\n"); + mat_trace(tr, element); + mat_trace(trinv, inverse); - mat_trace(tr, element); - mat_trace(trinv, inverse); + retval = solve_characteristic_polynomial(solver, tr, trinv, evs); + if(retval == 1) { + fprintf(stderr, "Error! Could not solve polynomial.\n"); + return 1; + } - DEBUG("Solve characteristic polynomials\n"); - retval = solve_characteristic_polynomial(solver, tr, trinv, evs); - if(retval == 1) { - fprintf(stderr, "Error! Could not solve polynomial.\n"); - return 1; - } + 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]); - 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]); + x = log(fabs(evs[0])); + y = -log(fabs(evs[2])); - x = log(fabs(evs[0])); - y = -log(fabs(evs[2])); + if(x > DBL_MAX || y > DBL_MAX) { + mpq_abs(tmp, tr); + mpq_abs(tmp2, trinv); + slope = mpq_log(tmp)/mpq_log(tmp2); + } else { + slope = y/x; + } - if(x > DBL_MAX || y > DBL_MAX) { - mpq_abs(tmp, tr); - mpq_abs(tmp2, trinv); - slope = mpq_log(tmp)/mpq_log(tmp2); - } else { - slope = y/x; - } + if(slope < 1) + slope = 1/slope; - if(slope < 1) - slope = 1/slope; + if(slope > max_slope) { + max_slope = slope; + max_slope_index = w; + } + if(slope < min_slope) { + min_slope = slope; + min_slope_index = w; + } -// if(slope > max_slope[IDX(i,j)]) { -// max_slope[IDX(i,j)] = slope; -// max_slope_index[IDX(i,j)] = w; -// } + if(mode != 2) + gmp_printf("%s %.9f\n", argv[w+3], slope); + } -// gmp_printf("%Qd %Qd %f %f %f\n", tr, trinv, x, y, y/x); -// gmp_printf("%s %.5f %.5f %Qd %Qd %.9f %.9f %.9f\n", argv[w], mpq_get_d(s), mpq_get_d(q), -// tr, trinv, -// x, y, slope); - gmp_printf("%s %.9f\n", argv[w], slope); + if(mode != 0) + printf("%d %.9f %s %d %.9f %s\n", + max_slope_index, max_slope, argv[max_slope_index+3], + min_slope_index, min_slope, argv[min_slope_index+3]); + fflush(stdout); - } - -// printf("%.5f %.5f %d %.9f\n", (double)i/DENOMINATOR, (double)j/DENOMINATOR, max_slope_index[IDX(i,j)], max_slope[IDX(i,j)]); - fflush(stdout); -// } -// } - - DEBUG("Clean up\n"); mpq_clears(m, t, s, q, tmp, tmp2, tr, trinv, NULL); mat_workspace_clear(ws); for(int i = 0; i < 6; i++) @@ -161,6 +166,4 @@ int main(int argc, char *argv[]) mat_clear(element); mat_clear(inverse); mps_context_free(solver); - free(max_slope); - free(max_slope_index); }