From 7f0cb0878725441eb1debd4ce1fde554c83167b3 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sat, 23 Oct 2021 16:04:12 -0500 Subject: [PATCH] combine singular_values and special_elements and test Charlie's prediction --- Makefile | 15 +- billiard_words.hs | 10 +- coxeter.c | 6 +- singular_values.c | 359 ++-------------------------------------------- special_element.c | 236 +++++++----------------------- 5 files changed, 80 insertions(+), 546 deletions(-) diff --git a/Makefile b/Makefile index a88cd22..0269d46 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -HEADERS=linalg.h mat.h coxeter.h +HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h #SPECIAL_OPTIONS=-O0 -g -D_DEBUG SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline @@ -16,14 +16,14 @@ convert: convert.hs billiard_words: billiard_words.hs ghc --make -dynamic billiard_words.hs -singular_values: singular_values.o coxeter.o mat.o - gcc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o -lm -lgmp -lmps +singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o + gcc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o enumerate_triangle_group.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 -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 -lgmp -lmps -lgsl -lcblas +special_element: special_element.o coxeter.o linalg.o mat.o enumerate_triangle_group.o + gcc $(OPTIONS) -o special_element coxeter.o linalg.o special_element.o mat.o enumerate_triangle_group.o -lm -lgmp -lmps -lgsl -lcblas singular_values.o: singular_values.c $(HEADERS) gcc $(OPTIONS) -c singular_values.c @@ -34,6 +34,9 @@ singular_values_mpi.o: singular_values_mpi.c $(HEADERS) special_element.o: special_element.c $(HEADERS) gcc $(OPTIONS) -c special_element.c +enumerate_triangle_group.o: enumerate_triangle_group.c $(HEADERS) + gcc $(OPTIONS) -c enumerate_triangle_group.c + linalg.o: linalg.c $(HEADERS) gcc $(OPTIONS) -c linalg.c @@ -44,4 +47,4 @@ mat.o: mat.c $(HEADERS) gcc $(OPTIONS) -c mat.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 + 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 diff --git a/billiard_words.hs b/billiard_words.hs index 4d959d3..9e88f51 100644 --- a/billiard_words.hs +++ b/billiard_words.hs @@ -1,18 +1,22 @@ import Data.List import Data.Ord import Text.Printf +import System.Environment -main = listWordsUpToLength 200 +main = do + argv <- getArgs + listWordsUpToLength $ read $ argv !! 0 listWordsUpToLength :: Int -> IO () listWordsUpToLength n = do - 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 "%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] 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..200], q <- [0..200], 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], p /= 0 || q /= 0] where sl ((p,q),_) = fromIntegral p / fromIntegral q diff --git a/coxeter.c b/coxeter.c index b8ca58f..128d9cd 100644 --- a/coxeter.c +++ b/coxeter.c @@ -9,9 +9,9 @@ group_t *coxeter_init_triangle(int p, int q, int r, int nmax) { int coxeter_matrix[9] = { - 1, p, r, - p, 1, q, - r, q, 1}; + 1, r, q, + r, 1, p, + q, p, 1}; return coxeter_init(3, coxeter_matrix, nmax); } diff --git a/singular_values.c b/singular_values.c index fe0cff0..e613e69 100644 --- a/singular_values.c +++ b/singular_values.c @@ -1,17 +1,12 @@ #include "coxeter.h" #include "linalg.h" #include "mat.h" - -#include -#include +#include "enumerate_triangle_group.h" #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 OUTPUT_POINTS -//#define OUTPUT_POINTS - struct result { int id; int count; @@ -67,326 +62,6 @@ static int compare_result_by_slope(const void *a_, const void *b_) return slopea > slopeb ? -1 : slopea < slopeb ? 1 : 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; - - 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); -} - -void initialize_triangle_generators(mat_workspace *ws, mat *gen, 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); - - quartic(rho1, s, 0, 0, 1, -1, 1); - quartic(rho2, s, 0, 0, 1, -1, 1); - quartic(rho3, s, 0, 0, 1, 0, 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, 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, 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); -} - -void output_invariants(group_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); - - for(int i = 0; i < group->size; i++) { - if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse) - continue; - - mat_trace(tr, matrices[i]); - mat_trace(trinv, matrices[group->elements[i].inverse->id]); - - 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->elements[i].length, print_word(&group->elements[i], buf), tr, trinv, log(evs[0]), -log(evs[2])); - } - - mpq_clears(tr, trinv, NULL); -} - -/* -double max_slope(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t, int *index) -{ - double max = 0; - double slope; - - mpq_t tr, trinv; - char buf[100]; - - mpq_inits(tr, trinv, NULL); - - for(int i = 0; i < MAX_ELEMENTS; i++) { - if(group[i].length % 2 != 0 || !group[i].inverse) - continue; - - mat_trace(tr, matrices[i]); - mat_trace(trinv, matrices[group[i].inverse->id]); - - slope = log(mpq_get_d(trinv))/log(mpq_get_d(tr)); - if(slope > max) - { - *index = i; - max = slope; - } - } - - mpq_clears(tr, trinv, NULL); - - return max; -} -*/ - int main(int argc, char *argv[]) { mpq_t s, q, t, tmp; @@ -407,7 +82,7 @@ int main(int argc, char *argv[]) struct result **distinct_invariants; if(argc < 4) { - fprintf(stderr, "Usage: %s \n", argv[0]); + fprintf(stderr, "Usage: %s \n", argv[0]); exit(1); } @@ -435,9 +110,10 @@ int main(int argc, char *argv[]) // get approximate s and q values sapprox = atof(argv[2]); - tapprox = atof(argv[3]); - tqfactor = pow((sapprox*sapprox-sapprox+1)*(sapprox*sapprox-sapprox+1)*(sapprox*sapprox+1), 1/6.0); - qapprox = tapprox/tqfactor; + qapprox = atof(argv[3]); +// tapprox = atof(argv[3]); +// tqfactor = pow((sapprox*sapprox-sapprox+1)*(sapprox*sapprox-sapprox+1)*(sapprox*sapprox+1), 1/6.0); +// qapprox = tapprox/tqfactor; for(int i = 0; ; i++) { continued_fraction_approximation(tmp, sapprox, i); @@ -459,16 +135,14 @@ int main(int argc, char *argv[]) tqfactor = pow((mpq_get_d(s)*mpq_get_d(s)-mpq_get_d(s)+1)*(mpq_get_d(s)*mpq_get_d(s)-mpq_get_d(s)+1)*(mpq_get_d(s)*mpq_get_d(s)+1), 1/6.0); -#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 + // order of the triangle reflection generators: a, b, c + // order of the rotation orders: bc, ac, ab DEBUG("Generate group\n"); group = coxeter_init_triangle(4, 3, 3, nmax); DEBUG("Compute matrices\n"); - enumerate(group, matrices, s, q); + enumerate(group, matrices, 4, 3, 3, s, q); n = 0; DEBUG("Compute traces\n"); @@ -538,6 +212,8 @@ int main(int argc, char *argv[]) qsort(distinct_invariants, nuniq, sizeof(struct result*), compare_result_by_slope); + gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor); + // printf("- 0 0 - - - - 0.5\n"); int cumulative = 0; double slope; @@ -572,19 +248,6 @@ int main(int argc, char *argv[]) } // printf("- 0 %d - - - - 2.0\n", cumulative); -#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)); -#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++) { mpq_clear(invariants[i].tr); diff --git a/special_element.c b/special_element.c index 84c382c..98c60b0 100644 --- a/special_element.c +++ b/special_element.c @@ -1,10 +1,7 @@ - #include "coxeter.h" #include "linalg.h" #include "mat.h" - -#include -#include +#include "enumerate_triangle_group.h" #define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0); #define DEBUG(msg, ...) @@ -14,180 +11,30 @@ #define HEIGHT 300 #define IDX(i,j) (((i)-1)*HEIGHT + ((j)-1)) -int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues) +double mpq_log(mpq_t m_op) { - 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; + static double logB = log(ULONG_MAX); - mpq_inits(coeff1, coeff2, zero, NULL); - mpq_set(coeff1, trinv); - mpq_sub(coeff2, zero, tr); + // Undefined logs (should probably return NAN in second case?) + if (mpz_get_ui(mpq_numref(m_op)) == 0 || mpz_sgn(mpq_numref(m_op)) < 0) + return -INFINITY; - 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); + // Log of numerator + double lognum = log(mpq_numref(m_op)->_mp_d[abs(mpq_numref(m_op)->_mp_size) - 1]); + lognum += (abs(mpq_numref(m_op)->_mp_size)-1) * logB; - 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; - } + // Subtract log of denominator, if it exists + if (abs(mpq_denref(m_op)->_mp_size) > 0) + { + lognum -= log(mpq_denref(m_op)->_mp_d[abs(mpq_denref(m_op)->_mp_size)-1]); + lognum -= (abs(mpq_denref(m_op)->_mp_size)-1) * logB; } - - cplx_vfree(roots); - mpq_clears(coeff1, coeff2, zero, NULL); - - return result; -} - -// 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) -{ - mpq_t s,sinv,q,x,y; - mpq_t zero, one, two; - mpq_t tmp; - - 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); - - // 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); - - // 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); - - // 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); - - // 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); - } - - LOOP(i,3) mat_pseudoinverse(ws, gen[i+3], gen[i]); - - mpq_inits(s,sinv,q,x,y,zero,one,two,tmp,NULL); -} - -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, mpq_t m, 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, m, 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); + return lognum; } int main(int argc, char *argv[]) { - mpq_t m, t, tmp; - double s; + mpq_t m, t, s, q, tmp, tmp2; mat_workspace *ws; mat gen[6]; mps_context *solver; @@ -203,7 +50,7 @@ int main(int argc, char *argv[]) DEBUG("Allocate\n"); - mpq_inits(m, t, tmp, tr, trinv, NULL); + mpq_inits(m, t, s, q, tmp, tmp2, tr, trinv, NULL); ws = mat_workspace_init(3); for(int i = 0; i < 6; i++) mat_init(gen[i], 3); @@ -218,15 +65,18 @@ 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); - for(int i = STARTX; i <= WIDTH; i++) { - for(int j = 1; j <= HEIGHT; j++) { - for(int w = 1; 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)); + 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)); DEBUG("Compute matrix\n"); - initialize_triangle_generators(ws, gen, m, t); + initialize_triangle_generators(ws, gen, 6, 4, 3, s, q); mat_identity(element); mat_identity(inverse); @@ -272,25 +122,39 @@ int main(int argc, char *argv[]) x = log(fabs(evs[0])); y = -log(fabs(evs[2])); - slope = y/x > 1 ? y/x : x/y; - if(slope > max_slope[IDX(i,j)]) { - max_slope[IDX(i,j)] = slope; - max_slope_index[IDX(i,j)] = w; + 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 > max_slope[IDX(i,j)]) { +// max_slope[IDX(i,j)] = slope; +// max_slope_index[IDX(i,j)] = w; +// } + // gmp_printf("%Qd %Qd %f %f %f\n", tr, trinv, x, y, y/x); -// gmp_printf("%.5f %.5f %.7f %.9f\n", mpq_get_d(t), mpq_get_d(m), s, slope); +// 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); } - printf("%.5f %.5f %d %.9f\n", (double)i/DENOMINATOR, (double)j/DENOMINATOR, max_slope_index[IDX(i,j)], max_slope[IDX(i,j)]); +// 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, tmp, tr, trinv, NULL); + mpq_clears(m, t, s, q, tmp, tmp2, tr, trinv, NULL); mat_workspace_clear(ws); for(int i = 0; i < 6; i++) mat_clear(gen[i]);