diff --git a/singular_values.c b/singular_values.c index e63a15f..58013ad 100644 --- a/singular_values.c +++ b/singular_values.c @@ -9,6 +9,9 @@ //#define DEBUG(msg, ...) fprintf(stderr, msg, ##__VA_ARGS__) #define DEBUG(msg, ...) +#define OUTPUT_POINTS +//#define OUTPUT_POINTS + struct result { mpq_t tr; mpq_t trinv; @@ -23,7 +26,7 @@ static int compare_result(const void *a_, const void *b_) d = mpq_cmp((*a)->tr,(*b)->tr); if(d == 0) - mpq_cmp((*a)->trinv, (*b)->trinv); + d = mpq_cmp((*a)->trinv, (*b)->trinv); return d; } @@ -356,11 +359,12 @@ int main(int argc, char *argv[]) int index; mps_context *solver; int acc = 100; - int n; + int n, nuniq, nmax; int retval; double evs[3]; double max_slope; - int nmax; + char buf[100]; + char buf2[100]; struct result *invariants; struct result **distinct_invariants; @@ -413,15 +417,18 @@ 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); -// gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor); +#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(3, 3, 4, nmax); + group = coxeter_init_triangle(4, 3, 3, nmax); DEBUG("Compute matrices\n"); enumerate(group, matrices, s, q); + n = 0; DEBUG("Compute traces\n"); for(int i = 0; i < nmax; i++) { if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse) @@ -429,21 +436,25 @@ int main(int argc, char *argv[]) mat_trace(invariants[i].tr, matrices[i]); mat_trace(invariants[i].trinv, matrices[group->elements[i].inverse->id]); + + distinct_invariants[n++] = &invariants[i]; } DEBUG("Get unique traces\n"); - qsort(distinct_invariants, nmax, sizeof(struct result*), compare_result); - n = 0; - for(int i = 0; i < nmax; i++) { - if(i == 0 || compare_result(&distinct_invariants[i], &distinct_invariants[n-1]) != 0) - distinct_invariants[n++] = distinct_invariants[i]; + qsort(distinct_invariants, n, sizeof(struct result*), compare_result); + + 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 < n; i++) { + 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"); @@ -459,16 +470,22 @@ int main(int argc, char *argv[]) if(fabs(evs[0]) < fabs(evs[1])) SWAP(double, evs[0], evs[1]); - double x = log(evs[0]); - double y = -log(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)) + if(y/x > max_slope && (x > 0.1 || y > 0.1)) { + max_slope_index = distinct_invariants[i] - invariants; max_slope = y/x; + } - gmp_printf("%Qd %Qd %f %f %f\n", distinct_invariants[i]->tr, distinct_invariants[i]->trinv, x, y, 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); +#endif } -// printf("%.3f %.3f %f\n", mpq_get_d(s), mpq_get_d(q)*tqfactor, max_slope); +#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); diff --git a/singular_values.plt b/singular_values.plt index af86668..f85674c 100644 --- a/singular_values.plt +++ b/singular_values.plt @@ -1,7 +1,8 @@ if(!exists("logt")) logt = log(1.78) if(!exists("logs")) logs = log(0.5) -file = sprintf("< ./singular_values 1000 %f %f", exp(logs), exp(logt)) +#file = sprintf("< ./singular_values 713698 %f %f", exp(logs), exp(logt)) +file = sprintf("< ./singular_values 1621 %f %f", exp(logs), exp(logt)) set zeroaxis set samples 1000