From 273528130032b66ddbb204bc8e865e25b51a24c2 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sat, 5 Feb 2022 17:51:45 -0600 Subject: [PATCH] reorganize singular_values, do one precomputation and then only compute unique traces --- Makefile | 2 +- billiard_words.hs | 10 +- coxeter.h | 1 + enumerate_triangle_group.c | 10 +- enumerate_triangle_group.h | 2 +- singular_values.c | 399 +++++++++++++++++++++++-------------- special_element.c | 6 +- 7 files changed, 273 insertions(+), 157 deletions(-) diff --git a/Makefile b/Makefile index 0269d46..0c7d7ba 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ 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 +SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL #SPECIAL_OPTIONS= diff --git a/billiard_words.hs b/billiard_words.hs index 9a60079..bedd44e 100644 --- a/billiard_words.hs +++ b/billiard_words.hs @@ -9,17 +9,21 @@ main = do listWordsUpToLength :: Int -> IO () listWordsUpToLength n = do - putStrLn $ unlines [printf "%s %d/%d %f %s" + putStr $ unlines [printf "%s %d/%d %f" w (p `div` gcd p q) (q `div` gcd p q) - (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double) - (slopeWord "bca" (orthogonalSlope (p,q))) | + (atan (sqrt 3 / (2*q_/p_ + 1))) | ((p,q),w) <- wordlist (n `div` 2, n `div` 2), + let p_ = fromIntegral p :: Double, + let q_ = fromIntegral q :: Double, length w <= n, let x = 2*q + p, let y = 2*p + q] +-- (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double) | +-- (slopeWord "bca" (orthogonalSlope (p,q))) | + wordlist :: (Int,Int) -> [((Int,Int),String)] wordlist (pmax,qmax) = nub $ sortBy (comparing sl) diff --git a/coxeter.h b/coxeter.h index 33ce018..ead28b5 100644 --- a/coxeter.h +++ b/coxeter.h @@ -7,6 +7,7 @@ typedef struct _groupelement { int id; int letter; struct _groupelement *parent; + int need_to_compute; // optimization flag; if 0, will skip matrix computation int length; struct _groupelement *inverse; diff --git a/enumerate_triangle_group.c b/enumerate_triangle_group.c index fc820d3..66e82f9 100644 --- a/enumerate_triangle_group.c +++ b/enumerate_triangle_group.c @@ -1,13 +1,12 @@ #include "enumerate_triangle_group.h" #include "linalg.h" -int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues) +int solve_characteristic_polynomial(mps_context *solv, mps_monomial_poly *poly, 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; @@ -15,7 +14,6 @@ int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, do 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); @@ -215,7 +213,7 @@ 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 t) +void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, mpq_t q) { mat_workspace *ws; mat tmp; @@ -228,7 +226,7 @@ void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, m mat_init(gen[i], 3); mat_init(tmp, 3); - initialize_triangle_generators(ws, gen, p1, p2, p3, s, t); + initialize_triangle_generators(ws, gen, p1, p2, p3, s, q); mat_identity(matrices[0]); for(int i = 1; i < group->size; i++) { @@ -236,6 +234,8 @@ void enumerate(group_t *group, mat *matrices, int p1, int p2, int p3, mpq_t s, m continue; if(!group->elements[i].inverse) continue; + if(!group->elements[i].need_to_compute) + continue; int parent = group->elements[i].parent->id; int grandparent = group->elements[i].parent->parent->id; diff --git a/enumerate_triangle_group.h b/enumerate_triangle_group.h index 5d8a85e..7488861 100644 --- a/enumerate_triangle_group.h +++ b/enumerate_triangle_group.h @@ -6,7 +6,7 @@ #include -int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues); +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 diff --git a/singular_values.c b/singular_values.c index e613e69..75e267a 100644 --- a/singular_values.c +++ b/singular_values.c @@ -3,9 +3,12 @@ #include "mat.h" #include "enumerate_triangle_group.h" +#include + #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 DEBUG(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__); +//#define DEBUG(msg, ...) struct result { int id; @@ -31,7 +34,17 @@ static int compare_result(const void *a_, const void *b_) return d; } -static int compare_result_with_id(const void *a_, const void *b_) +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; @@ -62,128 +75,73 @@ static int compare_result_by_slope(const void *a_, const void *b_) return slopea > slopeb ? -1 : slopea < slopeb ? 1 : 0; } -int main(int argc, char *argv[]) +struct timespec starttime; +static void start_timer() { - mpq_t s, q, t, tmp; - double sapprox, tapprox, qapprox, tqfactor; - mat *matrices; - group_t *group; - int index; + clock_gettime(CLOCK_MONOTONIC, &starttime); +} + +static double runtime() +{ + struct timespec curtime; + double diff; + clock_gettime(CLOCK_MONOTONIC, &curtime); + return (curtime.tv_sec - starttime.tv_sec) + (curtime.tv_nsec - starttime.tv_nsec) / 1e9; +} + +int compute_invariants(group_t *group, mat *matrices, struct result **invariants, int *n, int unique) +{ + mpq_t tmp; mps_context *solver; - int acc = 100; - int n, nuniq, nmax; + mps_monomial_poly *poly; + int index; + int ntraces = *n, nuniq; int retval; double evs[3]; + int max_slope_id; double max_slope; char buf[100]; - char buf2[100]; - struct result *invariants; - struct result **distinct_invariants; - - if(argc < 4) { - fprintf(stderr, "Usage: %s \n", argv[0]); - exit(1); +// 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]); } - nmax = atoi(argv[1]); + if(!unique) + nuniq = ntraces; + else { +// DEBUG("Get unique traces\n"); + qsort(invariants, ntraces, sizeof(struct result*), compare_result); - DEBUG("Allocate\n"); - - mpq_inits(s, q, t, tmp, NULL); - matrices = malloc(nmax*sizeof(mat)); - for(int i = 0; i < nmax; i++) - mat_init(matrices[i], 3); - invariants = malloc(nmax*sizeof(struct result)); - distinct_invariants = malloc(nmax*sizeof(struct result)); - for(int i = 0; i < nmax; i++) { - mpq_init(invariants[i].tr); - mpq_init(invariants[i].trinv); - distinct_invariants[i] = &invariants[i]; + nuniq = 0; + for(int i = 0; i < ntraces; i++) { + if(i == 0 || compare_result(&invariants[i], &invariants[nuniq-1]) != 0) { + invariants[nuniq] = invariants[i]; + invariants[nuniq]->count = 1; + nuniq++; + } else { + invariants[nuniq-1]->count++; + 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); - DEBUG("Approximate parameters\n"); - - // get approximate s and q values - sapprox = atof(argv[2]); - 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); - if(fabs(mpq_get_d(t)-sapprox) < 1e-10 - || (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0)) - break; - mpq_set(s, tmp); - } - mpq_canonicalize(s); - - for(int i = 0; ; i++) { - continued_fraction_approximation(tmp, qapprox, i); - if(fabs(mpq_get_d(t)-qapprox) < 1e-10 - || (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0)) - break; - mpq_set(q, tmp); - } - mpq_canonicalize(q); - - 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); - - // 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, 4, 3, 3, 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) - continue; - - invariants[i].id = i; - 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, 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]; - distinct_invariants[nuniq]->count = 1; - nuniq++; - } else { - distinct_invariants[nuniq-1]->count++; - int oldlength = group->elements[distinct_invariants[nuniq-1]->id].length; - int newlength = group->elements[distinct_invariants[i]->id].length; - if(newlength < oldlength) - distinct_invariants[nuniq-1]->id = distinct_invariants[i]->id; - } - -// gmp_printf("%d %d %s\n", i, nuniq-1, print_word(&group->elements[i], buf)); - } - max_slope = 0; - int max_slope_index; - - DEBUG("Solve characteristic polynomials\n"); for(int i = 0; i < nuniq; i++) { - retval = solve_characteristic_polynomial(solver, distinct_invariants[i]->tr, distinct_invariants[i]->trinv, evs); + retval = solve_characteristic_polynomial(solver, poly, invariants[i]->tr, invariants[i]->trinv, evs); + retval = 0;evs[0] = 2;evs[1] = 1;evs[2] = 0.5; // fake solving the polynomial for memory leak test if(retval == 1) { fprintf(stderr, "Error! Could not solve polynomial.\n"); continue; @@ -201,52 +159,204 @@ int main(int argc, char *argv[]) double x = log(fabs(evs[0])); double y = -log(fabs(evs[2])); - distinct_invariants[i]->x = x; - distinct_invariants[i]->y = y; + invariants[i]->x = x; + invariants[i]->y = y; - if(y/x > max_slope && (x > 0.1 || y > 0.1)) { - max_slope_index = distinct_invariants[i] - invariants; + if(y/x > max_slope + 1e-12 && (x > 0.1 || y > 0.1)) { + max_slope_id = invariants[i]->id; max_slope = y/x; + } else if(y/x > max_slope - 1e-12 && (x > 0.1 || y > 0.1)) { +// DEBUG("%s didn't quite make it\n", +// print_word(&group->elements[invariants[i]->id], buf)); + } + } + mps_context_free(solver); + + qsort(invariants, nuniq, sizeof(struct result*), compare_result_by_id); + + *n = nuniq; + return max_slope_id; +} + +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; } } - qsort(distinct_invariants, nuniq, sizeof(struct result*), compare_result_by_slope); + return total; +} - gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor); +int main(int argc, char *argv[]) +{ + mpq_t s, q, t, tmp; + int p1, p2, p3; + int sstart, send, sdenom, qstart, qend, qdenom; + mat *matrices; + group_t *group; + int nmax, n; + int max_slope_id; + char buf[100]; + char buf2[100]; + struct result *invariants; + struct result **distinct_invariants; -// printf("- 0 0 - - - - 0.5\n"); - int cumulative = 0; - double slope; - for(int i = 0; i < nuniq; i++) { - slope = distinct_invariants[i]->y/distinct_invariants[i]->x; + start_timer(); + + mpq_inits(s, q, t, tmp, NULL); + if(argc < 11) { + fprintf(stderr, "Usage: %s \n", argv[0]); + exit(1); + } + nmax = atoi(argv[1]); + p1 = atoi(argv[2]); + p2 = atoi(argv[3]); + p3 = atoi(argv[4]); + sstart = atoi(argv[5]); + send = atoi(argv[6]); + sdenom = atoi(argv[7]); + qstart = atoi(argv[8]); + qend = atoi(argv[9]); + qdenom = atoi(argv[10]); + + DEBUG("Allocate\n"); + matrices = malloc(nmax*sizeof(mat)); + for(int i = 0; i < nmax; i++) + mat_init(matrices[i], 3); + invariants = malloc(nmax*sizeof(struct result)); + distinct_invariants = malloc(nmax*sizeof(struct result)); + for(int i = 0; i < nmax; i++) { + mpq_init(invariants[i].tr); + mpq_init(invariants[i].trinv); + } + + // 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(p1, p2, p3, nmax); + + // first run; compute all matrices + for(int i = 0; i < group->size; i++) + group->elements[i].need_to_compute = 1; + + // use very generic values for the pilot run unless sstart=send and qstart=qend + if(sstart == send && qstart == qend) { + mpq_set_ui(s, sstart, sdenom); + mpq_set_ui(q, qstart, qdenom); + DEBUG("Single run for s = %d/%d, q = %d/%d\n", sstart, sdenom, qstart, qdenom); + } else { + mpq_set_ui(s, 4, 100); + mpq_set_ui(q, 7, 100); + DEBUG("Initial run for s = %d/%d, q = %d/%d\n", 4, 100, 7, 100); + } + + DEBUG("Compute matrices\n"); + enumerate(group, matrices, p1, p2, p3, s, q); + + // prepare array of ids + n = 0; + for(int i = 0; i < group->size; i++) + { + if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse) + continue; + invariants[i].id = i; + distinct_invariants[n++] = &invariants[i]; + } + + DEBUG("Compute invariants\n"); + max_slope_id = compute_invariants(group, matrices, distinct_invariants, &n, 1); + + // prepare for next time; don't need to change ids in distinct_invariants! + for(int i = 0; i < group->size; i++) + group->elements[i].need_to_compute = 0; + group->elements[0].need_to_compute = 1; + int multiplication_count = 1; + for(int i = 0; i < n; i++) { + groupelement_t *cur = &group->elements[distinct_invariants[i]->id]; + while(cur->need_to_compute == 0) { + cur->need_to_compute = 1; + multiplication_count++; + cur = cur->parent->parent; // also need to compute its even-length ancestors + } + cur = group->elements[distinct_invariants[i]->id].inverse; + while(cur->need_to_compute == 0) { + cur->need_to_compute = 1; + multiplication_count++; + cur = cur->parent->parent; + } + + } + + DEBUG("Would have needed %d matrix multiplications for %d unique traces up to reflection length %d\n", multiplication_count, n, group->elements[group->size-1].length); + + if(sstart != send || qstart != qend) { + for(int sloop = sstart; sloop <= send; sloop++) { + for(int qloop = qstart; qloop <= qend; qloop++) { + DEBUG("Loop for s = %d/%d, q = %d/%d\n", sloop, sdenom, qloop, qdenom); + mpq_set_ui(s, sloop, sdenom); + mpq_set_ui(q, qloop, qdenom); + DEBUG("Compute matrices\n"); + enumerate(group, matrices, p1, p2, p3, s, q); + DEBUG("Compute invariants\n"); + max_slope_id = compute_invariants(group, matrices, distinct_invariants, &n, 0); + // output + gmp_printf("%Qd %Qd %s\n", s, q, + print_word(&group->elements[max_slope_id], buf)); + fflush(stdout); + } + } + } else { + // output + for(int i = 0; i < n; i++) { + double slope = distinct_invariants[i]->y/distinct_invariants[i]->x; + + // exclude tr = trinv = 2/1/0/-1/3 + mpq_set_si(tmp, 2, 1); + if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, 1, 1); + if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, 0, 1); + if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, -1, 1); + if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) + continue; + mpq_set_si(tmp, 3, 1); + if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && + mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) + continue; + + gmp_printf("%d %d %s %f\n", + distinct_invariants[i]->id, distinct_invariants[i]->count, + print_word(&group->elements[distinct_invariants[i]->id], buf), + slope + ); + + /* + gmp_printf("%d %d %d %Qd %Qd %f %f %f %f %f %s\n", + distinct_invariants[i]->id, distinct_invariants[i]->count, cumulative, + distinct_invariants[i]->tr, distinct_invariants[i]->trinv, + log(fabs(mpq_get_d(distinct_invariants[i]->tr))), log(fabs(mpq_get_d(distinct_invariants[i]->trinv))), + distinct_invariants[i]->x, distinct_invariants[i]->y, slope, + print_word(&group->elements[distinct_invariants[i]->id], buf) + ); + */ + } - mpq_set_si(tmp, 1, 1); - if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) { - continue; - } - mpq_set_si(tmp, 0, 1); - if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) { - continue; - } - mpq_set_si(tmp, -1, 1); - if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) { - continue; - } - mpq_set_si(tmp, 3, 1); - if(mpq_cmp(distinct_invariants[i]->tr, tmp) == 0 && mpq_cmp(distinct_invariants[i]->trinv, tmp) == 0) { - continue; - } - - cumulative += distinct_invariants[i]->count; - gmp_printf("%d %d %d %Qd %Qd %f %f %f %f %f %s\n", - distinct_invariants[i]->id, distinct_invariants[i]->count, cumulative, - distinct_invariants[i]->tr, distinct_invariants[i]->trinv, - log(fabs(mpq_get_d(distinct_invariants[i]->tr))), log(fabs(mpq_get_d(distinct_invariants[i]->trinv))), - distinct_invariants[i]->x, distinct_invariants[i]->y, slope, - print_word(&group->elements[distinct_invariants[i]->id], buf) - ); } -// printf("- 0 %d - - - - 2.0\n", cumulative); DEBUG("Clean up\n"); for(int i = 0; i < nmax; i++) { @@ -260,5 +370,4 @@ int main(int argc, char *argv[]) free(matrices); coxeter_clear(group); mpq_clears(s, q, t, tmp, NULL); - mps_context_free(solver); } diff --git a/special_element.c b/special_element.c index 3a80b27..de84489 100644 --- a/special_element.c +++ b/special_element.c @@ -33,6 +33,7 @@ int main(int argc, char *argv[]) mat_workspace *ws; mat gen[6]; mps_context *solver; + mps_monomial_poly *poly; mat element, inverse; int letter1, letter2, letter; mpq_t tr, trinv; @@ -68,6 +69,7 @@ int main(int argc, char *argv[]) mat_init(inverse, 3); 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); @@ -112,7 +114,7 @@ int main(int argc, char *argv[]) mat_trace(tr, element); mat_trace(trinv, inverse); - retval = solve_characteristic_polynomial(solver, tr, trinv, evs); + retval = solve_characteristic_polynomial(solver, poly, tr, trinv, evs); if(retval == 1) { fprintf(stderr, "Error! Could not solve polynomial.\n"); return 1; @@ -150,7 +152,7 @@ int main(int argc, char *argv[]) } if(mode != 2) - gmp_printf("%s %.9f\n", argv[w+3], slope); + gmp_printf("%s %.9f %Qd %Qd\n", argv[w+3], slope, tr, trinv); } if(mode != 0)