diff --git a/Makefile b/Makefile index 4faffa1..5766967 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,8 @@ -HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h +HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h qext.h -#SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT_SQRT5 +SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT #SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline -SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT_SQRT5 +#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL #SPECIAL_OPTIONS= @@ -19,8 +19,8 @@ billiard_words: billiard_words.hs singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o parallel.o mpicc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps -singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o - mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps +singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o qext.o + mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o qext.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 @@ -55,5 +55,8 @@ mat.o: mat.c $(HEADERS) parallel.o: parallel.c $(HEADERS) gcc $(OPTIONS) -c parallel.c +qext.o: qext.c $(HEADERS) + gcc $(OPTIONS) -c qext.c + clean: - rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.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 parallel.o + rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.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 parallel.o qext.o diff --git a/enumerate_triangle_group.c b/enumerate_triangle_group.c index 673424f..d4b08b3 100644 --- a/enumerate_triangle_group.c +++ b/enumerate_triangle_group.c @@ -129,33 +129,38 @@ void discriminant(NUMBER out, NUMBER x, NUMBER y) CLEAR(x2);CLEAR(x3);CLEAR(y2);CLEAR(y3);CLEAR(tmp); } -void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER s, NUMBER q) +void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, mpq_t s, mpq_t q) { mat_workspace *ws; mat r1,r2,r3; NUMBER b1,c1,a2,c2,a3,b3; - NUMBER sinv; + mpq_t sinv, qinv; ws = mat_workspace_init(3); - INIT(sinv);INIT(b1);INIT(c1);INIT(a2);INIT(c2);INIT(a3);INIT(b3); + INIT(b1);INIT(c1);INIT(a2);INIT(c2);INIT(a3);INIT(b3); + mpq_init(sinv); + mpq_init(qinv); mat_init(r1, 3); mat_init(r2, 3); mat_init(r3, 3); // sinv = s^{-1} - SET_INT(sinv, 1); - DIV(sinv, sinv, s); + mpq_inv(sinv, s); + mpq_inv(qinv, q); // c1 = rho2 q, a2 = rho3 q, b3 = rho1 q, b1 = c2 = a3 = 1/q - MUL(c1, rho2, q); - MUL(a2, rho3, q); - MUL(b3, rho1, q); + SET_Q(c1, q); + SET_Q(a2, q); + SET_Q(b3, q); + MUL(c1, c1, rho2); + MUL(a2, a2, rho3); + MUL(b3, b3, rho1); SET_INT(b1, 1); SET_INT(c2, 1); SET_INT(a3, 1); - DIV(b1, b1, q); - DIV(c2, c2, q); - DIV(a3, a3, q); + SET_Q(b1, qinv); + SET_Q(c2, qinv); + SET_Q(a3, qinv); mat_zero(r1); mat_zero(r2); @@ -183,17 +188,17 @@ void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NU // gen[0] = diag(1,s^{-1},s) SET_INT(*mat_ref(gen[0], 0, 0), 1); - mat_set(gen[0], 1, 1, sinv); - mat_set(gen[0], 2, 2, s); + SET_Q (*mat_ref(gen[0], 1, 1), sinv); + SET_Q (*mat_ref(gen[0], 2, 2), s); // gen[1] = diag(s,1,s^{-1}) - mat_set(gen[1], 0, 0, s); + SET_Q (*mat_ref(gen[1], 0, 0), s); SET_INT(*mat_ref(gen[1], 1, 1), 1); - mat_set(gen[1], 2, 2, sinv); + SET_Q (*mat_ref(gen[1], 2, 2), sinv); - // gen[3] = diag(s^{-1},s,1) - mat_set(gen[2], 0, 0, sinv); - mat_set(gen[2], 1, 1, s); + // gen[2] = diag(s^{-1},s,1) + SET_Q (*mat_ref(gen[2], 0, 0), sinv); + SET_Q (*mat_ref(gen[2], 1, 1), s); SET_INT(*mat_ref(gen[2], 2, 2), 1); // gen[0] = r2 * gen[0] * r3 @@ -214,7 +219,9 @@ void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NU mat_pseudoinverse(ws, gen[5], gen[2]); mat_workspace_clear(ws); - CLEAR(sinv);CLEAR(b1);CLEAR(c1);CLEAR(a2);CLEAR(c2);CLEAR(a3);CLEAR(b3); + CLEAR(b1);CLEAR(c1);CLEAR(a2);CLEAR(c2);CLEAR(a3);CLEAR(b3); + mpq_clear(sinv); + mpq_clear(qinv); mat_clear(r1); mat_clear(r2); mat_clear(r3); @@ -239,31 +246,29 @@ void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3 } #endif -#ifdef QEXT_SQRT5 +#ifdef QEXT void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s_, mpq_t q_) { - NUMBER rho, c, one, s, q; - INIT(rho);INIT(c);INIT(one);INIT(s);INIT(q); + NUMBER rho, c, one, s; + INIT(rho);INIT(c);INIT(one);INIT(s); // c = - (1+sqrt(5))/2 - mpq_set_si(c[0], -1, 2); - mpq_set_si(c[1], -1, 2); + mpq_set_si(c->a[0], -1, 2); + mpq_set_si(c->a[1], -1, 2); SET_ONE(one); - mpq_set(s[0], s_); - mpq_set_ui(s[1], 0, 1); - mpq_set(q[0], q_); - mpq_set_ui(q[1], 0, 1); + SET_Q(s, s_); + // rho = s^2 + cs + 1 SET(rho, one); MUL(rho, rho, s); ADD(rho, rho, c); MUL(rho, rho, s); ADD(rho, rho, one); - generators_triangle_rotation_generic(gen, rho, rho, rho, s, q); + generators_triangle_rotation_generic(gen, rho, rho, rho, s_, q_); - CLEAR(rho);CLEAR(c);CLEAR(one);CLEAR(s);CLEAR(q); + CLEAR(rho);CLEAR(c);CLEAR(one);CLEAR(s); } #endif diff --git a/enumerate_triangle_group.h b/enumerate_triangle_group.h index 47b46fc..bb78d73 100644 --- a/enumerate_triangle_group.h +++ b/enumerate_triangle_group.h @@ -10,15 +10,15 @@ 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); -// geneartors -void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER s, NUMBER q); +// generators +void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, mpq_t s, mpq_t q); #ifdef QEXT_TRIVIAL void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3, mpq_t s, mpq_t q); #endif -#ifdef QEXT_SQRT5 -void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s, mpq_t q); +#ifdef QEXT +void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s_, mpq_t q_); #endif // general functions diff --git a/mat.h b/mat.h index 85ac648..5040fbe 100644 --- a/mat.h +++ b/mat.h @@ -8,119 +8,25 @@ // library for matrix computations in variable rings (based on GMP types) -#ifdef QEXT_SQRT5 +#ifdef QEXT -typedef mpq_t NUMBER [2]; -#define INIT(x) do { mpq_init((x)[0]); mpq_init((x)[1]); } while(0) -#define CLEAR(x) do { mpq_clear((x)[0]); mpq_clear((x)[1]); } while(0) -#define SET(x,y) do { mpq_set((x)[0],(y)[0]); mpq_set((x)[1],(y)[1]); } while(0) -#define SET_INT(x,y) do { mpq_set_si((x)[0],(y),1); mpq_set_si((x)[1],0,1); } while(0) -#define SET_ZERO(x) SET_INT(x,0) -#define SET_ONE(x) SET_INT(x,1) -#define ADD(x,y,z) do { mpq_add((x)[0],(y)[0],(z)[0]); mpq_add((x)[1],(y)[1],(z)[1]); } while(0) -#define SUB(x,y,z) do { mpq_sub((x)[0],(y)[0],(z)[0]); mpq_sub((x)[1],(y)[1],(z)[1]); } while(0) -#define MUL multiply_sqrt5 -#define DIV divide_sqrt5 -#define CMP compare_sqrt5 -#define PRINT(x) gmp_printf("%Qd + %Qd*sqrt(5)", (x)[0], (x)[1]) -#define SPRINT(buf, x) gmp_sprintf(buf, "%Qd+%Qd*sqrt(5)", (x)[0], (x)[1]) +#include "qext.h" -static void multiply_sqrt5(NUMBER out, NUMBER a, NUMBER b) -{ - // could be in place!!! - - mpq_t tmp, result[2]; - mpq_inits(tmp, result[0], result[1], 0); - - mpq_mul(result[0], a[1], b[1]); - mpq_set_si(tmp, 5, 1); - mpq_mul(result[0], result[0], tmp); - mpq_mul(tmp, a[0], b[0]); - mpq_add(result[0], result[0], tmp); - - mpq_mul(result[1], a[1], b[0]); - mpq_mul(tmp, a[0], b[1]); - mpq_add(result[1], result[1], tmp); - - mpq_set(out[0], result[0]); - mpq_set(out[1], result[1]); - - mpq_clears(tmp, result[0], result[1], 0); -} - -static void divide_sqrt5(NUMBER out, NUMBER a, NUMBER b) -{ - mpq_t denom, num, tmp, neg5, result[2]; - mpq_inits(denom, num, tmp, neg5, result[0], result[1], 0); - mpq_set_si(neg5, -5, 1); - - mpq_mul(denom, b[1], b[1]); - mpq_mul(denom, denom, neg5); - mpq_mul(tmp, b[0], b[0]); - mpq_add(denom, denom, tmp); - - mpq_mul(num, a[1], b[1]); - mpq_mul(num, num, neg5); - mpq_mul(tmp, a[0], b[0]); - mpq_add(num, num, tmp); - - mpq_div(result[0], num, denom); - - mpq_mul(num, a[1], b[0]); - mpq_mul(tmp, a[0], b[1]); - mpq_sub(num, num, tmp); - - mpq_div(result[1], num, denom); - - mpq_set(out[0], result[0]); - mpq_set(out[1], result[1]); - - mpq_clears(denom, num, tmp, neg5, result[0], result[1], 0); -} - -static int compare_sqrt5(NUMBER x, NUMBER y) -{ - int result; - mpq_t p, q, p2, q2, c5; - mpq_inits(p, q, p2, q2, c5, 0); - - mpq_sub(p, x[0], y[0]); - mpq_sub(q, x[1], y[1]); - - /* - want to know if p + sqrt(5) q > 0 - if p>0 and q>0: always true - if p>0 and q<0: equivalent to |p|^2 > 5 |q|^2 - if p<0 and q>0: equivalent to |p|^2 < 5 |q|^2 - if p<0 and q<0: always false - */ - - if(mpq_sgn(p) > 0 && mpq_sgn(q) > 0) { - result = 1; - goto done; - } - - if(mpq_sgn(p) < 0 && mpq_sgn(q) < 0) { - result = -1; - goto done; - } - - mpq_mul(p2, p, p); - mpq_mul(q2, q, q); - mpq_set_si(c5, 5, 1); - mpq_mul(q2, q2, c5); - - if(mpq_sgn(p) > 0) - result = mpq_cmp(p2, q2); // this can't be zero, or else |p/q| = sqrt(5), but it's irrational - else if(mpq_sgn(p) < 0) - result = -mpq_cmp(p2, q2); // this can be zero if p = 0 and q = 0 - else // p = 0 - return mpq_sgn(q); - -done: - mpq_clears(p, q, p2, q2, c5, 0); - return result; -} +#define NUMBER qext_number +#define INIT qext_init +#define CLEAR qext_clear +#define SET qext_set +#define SET_INT qext_set_int +#define SET_Q qext_set_q +#define SET_ZERO(x) qext_set_int((x),0) +#define SET_ONE(x) qext_set_int((x),1) +#define ADD qext_add +#define SUB qext_sub +#define MUL qext_mul +// #define DIV qext_div +#define CMP qext_cmp +#define PRINT qext_print +#define SNPRINT qext_snprint #endif diff --git a/qext.c b/qext.c index fe6933e..66a0d7f 100644 --- a/qext.c +++ b/qext.c @@ -1,15 +1,26 @@ #include "qext.h" -const int qext_coefficients[QEXT_RANK+1] = {36, 0, -8, 0, 1}; +int qext_rank; +mpq_t *qext_coefficient; #include #include #define LOOP(i,n) for(int i = 0; i < (n); i++) +void qext_setup(int rank, const int *coeffs) +{ + qext_rank = rank; + qext_coefficient = malloc(rank*sizeof(mpq_t)); + LOOP(i, rank) { + mpq_init(qext_coefficient[i]); + mpq_set_si(qext_coefficient[i], -coeffs[i], coeffs[rank]); + } +} + void qext_init(qext_number x) { - x->rk = QEXT_RANK; + x->rk = qext_rank; x->a = malloc(x->rk*sizeof(mpq_t)); LOOP(i, x->rk) mpq_init(x->a[i]); } @@ -32,6 +43,13 @@ void qext_set_int(qext_number x, int y) mpq_set_ui(x->a[i], 0, 1); } +void qext_set_q(qext_number x, mpq_t y) +{ + mpq_set(x->a[0], y); + for(int i = 1; i < x->rk; i++) + mpq_set_ui(x->a[i], 0, 1); +} + void qext_add(qext_number result, qext_number x, qext_number y) { LOOP(i, x->rk) mpq_add(result->a[i], x->a[i], y->a[i]); @@ -83,21 +101,13 @@ int qext_cmp(qext_number x, qext_number y) void qext_mul(qext_number out, qext_number x, qext_number y) { - static initialized = 0; - static mpq_t result[2*x->rk-1]; - static mpq_t coefficient[x->rk]; - static mpq_t tmp; + mpq_t result[2*x->rk-1]; + mpq_t tmp; - if(!initialized) { - run_before = 1; - mpq_init(tmp); - LOOP(i, 2*x->rk-1) { - mpq_init(result[i]); - } - LOOP(i, x->rk) { - mpq_init(coefficient[i]); - mpq_set_si(coefficient[i], -qext_coefficients[i], qext_coefficients[x->rk]); - } + mpq_init(tmp); + LOOP(i, 2*x->rk-1) { + mpq_init(result[i]); + mpq_set_ui(result[i], 0, 1); } // rk*rk multiplications + (rk-1)*rk multiplications with fixed coefficients @@ -105,8 +115,6 @@ void qext_mul(qext_number out, qext_number x, qext_number y) // degree 2: 4+2 // degree 4: 16+12, including 6 times 0* - LOOP(i, 2*x->rk-1) mpq_set_ui(result[i], 0, 1); - LOOP(i, x->rk) LOOP(j, x->rk) { mpq_mul(tmp, x->a[i], y->a[j]); mpq_add(result[i+j], result[i+j], tmp); @@ -114,18 +122,18 @@ void qext_mul(qext_number out, qext_number x, qext_number y) for(int i = x->rk-2; i >= 0; i--) { LOOP(j, x->rk) { - mpq_mul(tmp, coefficient[j], result[x->rk+i]); + mpq_mul(tmp, qext_coefficient[j], result[x->rk+i]); mpq_add(result[i+j], result[i+j], tmp); } } LOOP(i, x->rk) mpq_set(out->a[i], result[i]); - /* LOOP(i, 2*x->rk-1) mpq_clear(result[i]); - LOOP(i, x->rk) mpq_clear(coefficient[i]); mpq_clear(tmp); - */ } -void qext_div(qext_number out, qext_number x, qext_number y); +void qext_div(qext_number out, qext_number x, qext_number y) +{ + // todo: implement this by solving a linear system +} diff --git a/qext.h b/qext.h index 2e81d32..fed4cf9 100644 --- a/qext.h +++ b/qext.h @@ -3,8 +3,6 @@ #include -#define QEXT_RANK 4 - struct qext_number_internal { int rk; mpq_t *a; @@ -12,26 +10,18 @@ struct qext_number_internal { typedef struct qext_number_internal qext_number[1]; +void qext_setup(int rank, const int *coeffs); void qext_init(qext_number x); void qext_clear(qext_number x); void qext_set(qext_number x, qext_number y); void qext_set_int(qext_number x, int y); +void qext_set_q(qext_number x, mpq_t y); void qext_add(qext_number result, qext_number x, qext_number y); void qext_sub(qext_number result, qext_number x, qext_number y); void qext_mul(qext_number out, qext_number x, qext_number y); void qext_div(qext_number out, qext_number x, qext_number y); int qext_cmp(qext_number x, qext_number y); void qext_print(qext_number x); -void qext_sprint(qext_number x); - -static void qext_set_zero(qext_number x) -{ - qext_set_int(x, 0); -} - -static void qext_set_one(qext_number x) -{ - qext_set_int(x, 1); -} +void qext_snprint(char *str, size_t len, qext_number x); #endif diff --git a/singular_values_barbot.c b/singular_values_barbot.c index 2c668e8..850eb06 100644 --- a/singular_values_barbot.c +++ b/singular_values_barbot.c @@ -145,12 +145,15 @@ int invariants_trace_loxodromic(group_t *group, mat *matrices, struct result **i // check if loxodromic NUMBER disc, zero; + double disc_real; INIT(disc); INIT(zero); SET_ZERO(zero); for(int i = 0; i < nuniq; i++) { discriminant(disc, invariants[i]->tr, invariants[i]->trinv); - invariants[i]->disc_sign = CMP(disc, zero); + disc_real = mpq_get_d(disc->a[0]) + sqrt(5)*mpq_get_d(disc->a[1]); + gmp_printf("%Qd %Qd %f\n", disc->a[0], disc->a[1], disc_real); + invariants[i]->disc_sign = disc_real > 0 ? 1 : disc_real < 0 ? -1 : 0; } CLEAR(disc); CLEAR(zero); @@ -417,6 +420,9 @@ int main(int argc, char *argv[]) g->qend = atoi(argv[9]); g->qdenom = atoi(argv[10]); + const int coefficients[] = {-5, 0, 1}; + qext_setup(2, coefficients); + // initialize parallel_context *ctx = parallel_init(); parallel_set_datasize_and_callbacks(ctx, init_node, do_computation, destroy_node, @@ -537,13 +543,13 @@ int main(int argc, char *argv[]) mpq_clear(tmp); */ - SPRINT(buf, n.distinct_invariants[i]->tr); - SPRINT(buf2, n.distinct_invariants[i]->trinv); - printf("%d %s %f %f %d\n", + SNPRINT(buf, sizeof(buf), n.distinct_invariants[i]->tr); + SNPRINT(buf2, sizeof(buf2), n.distinct_invariants[i]->trinv); + printf("%d %s %d\n", n.distinct_invariants[i]->id, print_word(&n.group->elements[n.distinct_invariants[i]->id], buf3), - mpq_get_d(n.distinct_invariants[i]->tr[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->tr[1]), - mpq_get_d(n.distinct_invariants[i]->trinv[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->trinv[1]), +// mpq_get_d(n.distinct_invariants[i]->tr[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->tr[1]), +// mpq_get_d(n.distinct_invariants[i]->trinv[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->trinv[1]), n.distinct_invariants[i]->disc_sign); }