diff --git a/singular_values.c b/singular_values.c index bc06c61..3907c0f 100644 --- a/singular_values.c +++ b/singular_values.c @@ -9,7 +9,33 @@ #define MAX_ELEMENTS 14000 //#define DRAW_PICTURE 1 -void mpq_quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e) +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); @@ -46,12 +72,12 @@ void initialize_triangle_generators(mat_workspace *ws, mat *gen, mpq_t s, mpq_t mpq_set_ui(*mat_ref(gen[1], 0, 0), 1, 1); mpq_set_ui(*mat_ref(gen[1], 1, 0), 0, 1); mpq_set_ui(*mat_ref(gen[1], 2, 0), 0, 1); - mpq_quartic(*mat_ref(gen[1], 0, 1), t, 0, 0, 1, -1, 2); - mpq_quartic(*mat_ref(gen[1], 1, 1), t, 0, 0, -1, 2, -2); - mpq_quartic(*mat_ref(gen[1], 2, 1), t, 0, 0, 1, -3, 3); - mpq_quartic(*mat_ref(gen[1], 0, 2), t, 0, 0, 1, 0, 3); - mpq_quartic(*mat_ref(gen[1], 1, 2), t, 0, 0, -1, 1, -1); - mpq_quartic(*mat_ref(gen[1], 2, 2), t, 0, 0, 1, -2, 1); + quartic(*mat_ref(gen[1], 0, 1), t, 0, 0, 1, -1, 2); + quartic(*mat_ref(gen[1], 1, 1), t, 0, 0, -1, 2, -2); + quartic(*mat_ref(gen[1], 2, 1), t, 0, 0, 1, -3, 3); + quartic(*mat_ref(gen[1], 0, 2), t, 0, 0, 1, 0, 3); + quartic(*mat_ref(gen[1], 1, 2), t, 0, 0, -1, 1, -1); + quartic(*mat_ref(gen[1], 2, 2), t, 0, 0, 1, -2, 1); mat_pseudoinverse(ws, gen[3], gen[0]); // p^{-1} mat_pseudoinverse(ws, gen[4], gen[1]); // q^{-1} @@ -79,95 +105,152 @@ char *print_word(groupelement_t *g, char *str) return str; } +void enumerate(groupelement_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 < MAX_ELEMENTS; i++) { + if(group[i].length % 2 != 0) + continue; + if(!group[i].inverse) + continue; + + int parent = group[i].parent->id; + int grandparent = group[i].parent->parent->id; + int letter; + + if(group[parent].letter == 1 && group[i].letter == 2) + letter = 0; // p = bc + else if(group[parent].letter == 2 && group[i].letter == 0) + letter = 1; // q = ca + else if(group[parent].letter == 0 && group[i].letter == 1) + letter = 2; // r = ab + if(group[parent].letter == 2 && group[i].letter == 1) + letter = 3; // p^{-1} = cb + else if(group[parent].letter == 0 && group[i].letter == 2) + letter = 4; // q^{-1} = ac + else if(group[parent].letter == 1 && group[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(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t) +{ + 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]); + + 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))); + } + + 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[]) { - groupelement_t *group; - mat_workspace *ws; - mat *matrices; - mat tmp; - mat gen[6]; - char buf[100], buf2[100], buf3[100]; - mpq_t s,t; - mpq_t det, tr, trinv; + mpq_t s, t, tmp; + mpz_t accuracy; + double t_; + mat *matrices; + groupelement_t *group; + int index; - // allocate stuff - group = malloc(MAX_ELEMENTS*sizeof(groupelement_t)); - ws = mat_workspace_init(3); - matrices = malloc(MAX_ELEMENTS*sizeof(mat)); - for(int i = 0; i < MAX_ELEMENTS; i++) - mat_init(matrices[i], 3); - for(int i = 0; i < 6; i++) - mat_init(gen[i], 3); - mat_init(tmp, 3); + mpq_inits(s, t, tmp, NULL); + mpz_init(accuracy); + group = malloc(MAX_ELEMENTS*sizeof(groupelement_t)); + matrices = malloc(MAX_ELEMENTS*sizeof(mat)); + for(int i = 0; i < MAX_ELEMENTS; i++) + mat_init(matrices[i], 3); - mpq_inits(s, t, det, tr, trinv, NULL); + // mpq_set_str(t, argv[1], 10); + mpz_set_ui(accuracy, 100); + for(int i = 0; ; i++) { + mpq_set(t, tmp); + continued_fraction_approximation(tmp, atof(argv[1]), i); + if(mpz_cmp(mpq_numref(tmp),accuracy) > 0 && mpz_cmp(mpq_denref(tmp),accuracy) > 0) + break; + } + mpq_canonicalize(t); + gmp_fprintf(stdout, "\"t = %Qd = %.3f\"\n", mpq_get_d(t), t); - mpq_set_ui(s,1,1); - double t_ = atof(argv[1]); - mpq_set_ui(t,(int)(t_*100),100); - mpq_canonicalize(s); - mpq_canonicalize(t); + if(argc > 2 && strcmp(argv[2],"p") == 0) { + gmp_fprintf(stdout, "%Qd\n", t); + return 0; + } - // the real action - generate_triangle_group(group, MAX_ELEMENTS, 3, 3, 4); - initialize_triangle_generators(ws, gen, s, t); + generate_triangle_group(group, MAX_ELEMENTS, 3, 3, 4); - mat_identity(matrices[0]); - for(int i = 1; i < MAX_ELEMENTS; i++) { - if(group[i].length % 2 != 0) - continue; - if(!group[i].inverse) - continue; +// for(int i = 0; i < 10; i++) { +// mpq_set_ui(t,100+i,100); +// mpq_canonicalize(t); - int parent = group[i].parent->id; - int grandparent = group[i].parent->parent->id; - int letter; + enumerate(group, matrices, s, t); + //printf("%f %f\n", mpq_get_d(t), max_slope(group, matrices, s, t, &index)); + output_invariants(group, matrices, s, t); +// } - if(group[parent].letter == 1 && group[i].letter == 2) - letter = 0; // p = bc - else if(group[parent].letter == 2 && group[i].letter == 0) - letter = 1; // q = ca - else if(group[parent].letter == 0 && group[i].letter == 1) - letter = 2; // r = ab - if(group[parent].letter == 2 && group[i].letter == 1) - letter = 3; // p^{-1} = cb - else if(group[parent].letter == 0 && group[i].letter == 2) - letter = 4; // q^{-1} = ac - else if(group[parent].letter == 1 && group[i].letter == 0) - letter = 5; // r^{-1} = ba - - mat_multiply(ws, matrices[i], matrices[grandparent], gen[letter]); - } - - for(int i = 0; i < MAX_ELEMENTS; i++) { - if(group[i].length % 2 != 0) - continue; - if(!group[i].inverse) - continue; - - mat_trace(tr, matrices[i]); - mat_trace(trinv, matrices[group[i].inverse->id]); - - double lambda1, lambda2, lambda3; -// int realevs = gsl_poly_solve_cubic(-mpq_get_d(tr), mpq_get_d(trinv), -1, &lambda3, &lambda2, &lambda1); -// if(realevs != 3) -// continue; -// if(lambda1 < 0 || lambda2 < 0 || lambda3 < 0) -// continue; - -// gmp_printf("%d %d %s %Qd %Qd\n", i, group[i].length, print_word(&group[i], buf), tr, trinv); - printf("%d %d %s %f %f %f\n", i, group[i].length, print_word(&group[i], buf), log(mpq_get_d(tr)), log(mpq_get_d(trinv))); - - } - - // free stuff - for(int i = 0; i < MAX_ELEMENTS; i++) - mat_clear(matrices[i]); - for(int i = 0; i < 6; i++) - mat_clear(gen[i]); - mat_clear(tmp); - mpq_clears(s, t, det, tr, trinv, NULL); - mat_workspace_clear(ws); - - return 0; + for(int i = 0; i < MAX_ELEMENTS; i++) + mat_clear(matrices[i]); + free(matrices); + free(group); + mpq_clears(s, t, tmp, NULL); + mpz_clear(accuracy); } diff --git a/singular_values.plt b/singular_values.plt index ce251a7..cd4894e 100644 --- a/singular_values.plt +++ b/singular_values.plt @@ -3,7 +3,7 @@ if(!exists("logs")) logs = log(1.0) file = sprintf("< ./singular_values %f", exp(logt)) #title = sprintf("s = %f, t = %f", exp(logs), exp(logt)) -title = sprintf("t = %.3f", floor(exp(logt)*100)/100.0) +title = sprintf("t = %.3f", exp(logt)) # print title set zeroaxis @@ -23,7 +23,7 @@ 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 4:5 w p pt 7 ps 0.5 lc 1 t title, \ +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 ""