From 1a85a5ef144f85e01faadfe04eb3cefa2070951d Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sat, 11 Jun 2022 17:36:32 +0200 Subject: [PATCH] update QEXT to new dynamic version --- enumerate_triangle_group.c | 23 ++++++----- mat.c | 12 +++--- mat.h | 22 +++++----- qext.c | 82 +++++++++++++++++++++++++------------- qext.h | 16 ++++++-- singular_values_barbot.c | 15 +++---- 6 files changed, 106 insertions(+), 64 deletions(-) diff --git a/enumerate_triangle_group.c b/enumerate_triangle_group.c index d4b08b3..2c03e7a 100644 --- a/enumerate_triangle_group.c +++ b/enumerate_triangle_group.c @@ -73,7 +73,7 @@ void continued_fraction_approximation(mpq_t out, double in, int level) void quartic(NUMBER out, NUMBER in, int a, int b, int c, int d, int e) { NUMBER tmp; - INIT(tmp); + INIT(tmp, GETTYPE(in)); SET_INT(out, a); @@ -100,8 +100,9 @@ void quartic(NUMBER out, NUMBER in, int a, int b, int c, int d, int e) // that is the function x^2 y^2 - 4 x^3 - 4 y^3 - 27 + 18 xy void discriminant(NUMBER out, NUMBER x, NUMBER y) { + TYPE type = GETTYPE(out); NUMBER x2, x3, y2, y3, tmp; - INIT(x2);INIT(x3);INIT(y2);INIT(y3);INIT(tmp); + INIT(x2, type);INIT(x3, type);INIT(y2, type);INIT(y3, type);INIT(tmp, type); MUL(x2, x, x); MUL(x3, x2, x); @@ -134,15 +135,16 @@ void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NU mat_workspace *ws; mat r1,r2,r3; NUMBER b1,c1,a2,c2,a3,b3; + TYPE type = GETTYPE(rho1); mpq_t sinv, qinv; - ws = mat_workspace_init(3); - INIT(b1);INIT(c1);INIT(a2);INIT(c2);INIT(a3);INIT(b3); + ws = mat_workspace_init(3, type); + INIT(b1, type);INIT(c1, type);INIT(a2, type);INIT(c2, type);INIT(a3, type);INIT(b3, type); mpq_init(sinv); mpq_init(qinv); - mat_init(r1, 3); - mat_init(r2, 3); - mat_init(r3, 3); + mat_init(r1, 3, type); + mat_init(r2, 3, type); + mat_init(r3, 3, type); // sinv = s^{-1} mpq_inv(sinv, s); @@ -250,7 +252,7 @@ void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3 void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s_, mpq_t q_) { NUMBER rho, c, one, s; - INIT(rho);INIT(c);INIT(one);INIT(s); + INIT(rho, QT_SQRT5);INIT(c, QT_SQRT5);INIT(one, QT_SQRT5);INIT(s, QT_SQRT5); // c = - (1+sqrt(5))/2 mpq_set_si(c->a[0], -1, 2); @@ -292,8 +294,9 @@ void enumerate_triangle_rotation_subgroup(group_t *group, mat *gen, mat *matrice char buf[100], buf2[100], buf3[100]; // allocate stuff - ws = mat_workspace_init(3); - mat_init(tmp, 3); + TYPE type = GETTYPE(gen[0]->x[0]); + ws = mat_workspace_init(3, type); + mat_init(tmp, 3, type); // initialize_triangle_generators(ws, gen, p1, p2, p3, s, q); diff --git a/mat.c b/mat.c index 391ea22..0e7aa37 100644 --- a/mat.c +++ b/mat.c @@ -1,11 +1,11 @@ #include "mat.h" -mat_workspace *mat_workspace_init(int n) +mat_workspace *mat_workspace_init(int n, TYPE t) { mat_workspace *ws = (mat_workspace*)malloc(sizeof(mat_workspace)); - mat_init(ws->tmp_mat, n); - INIT(ws->tmp_num); - INIT(ws->tmp_num2); + mat_init(ws->tmp_mat, n, t); + INIT(ws->tmp_num, t); + INIT(ws->tmp_num2, t); return ws; } @@ -17,11 +17,11 @@ void mat_workspace_clear(mat_workspace *ws) free(ws); } -void mat_init(mat m, int n) +void mat_init(mat m, int n, TYPE t) { m->n = n; m->x = malloc(n*n*sizeof(NUMBER)); - LOOP(i,n) LOOP(j,n) INIT(m->x[i+j*n]); + LOOP(i,n) LOOP(j,n) INIT(m->x[i+j*n], t); } void mat_get(NUMBER out, mat m, int i, int j) diff --git a/mat.h b/mat.h index 5040fbe..fcc2b0b 100644 --- a/mat.h +++ b/mat.h @@ -22,28 +22,34 @@ #define SET_ONE(x) qext_set_int((x),1) #define ADD qext_add #define SUB qext_sub +#define NEG qext_neg #define MUL qext_mul // #define DIV qext_div #define CMP qext_cmp #define PRINT qext_print #define SNPRINT qext_snprint +#define TYPE struct qext_type* +#define GETTYPE(x) ((x)->type) -#endif - -#ifdef QEXT_TRIVIAL +#else #define NUMBER mpq_t -#define INIT mpq_init +#define INIT(x,t) mpq_init(x) #define CLEAR mpq_clear #define SET mpq_set #define SET_INT(x,y) mpq_set_si(x,y,1) +#define SET_Q SET #define SET_ZERO(x) SET_INT(x,0) #define SET_ONE(x) SET_INT(x,1) #define ADD mpq_add #define SUB mpq_sub +#define NEG mpq_neg #define MUL mpq_mul #define DIV mpq_div #define PRINT(x) gmp_printf("%Qd", x) +#define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x) +#define TYPE int +#define GETTYPE(x) 0 #endif @@ -52,7 +58,7 @@ struct _mat{ int n; NUMBER *x; -} ; +}; typedef struct _mat mat[1]; @@ -62,9 +68,9 @@ typedef struct _mat_workspace { NUMBER tmp_num2; } mat_workspace; -mat_workspace *mat_workspace_init(int n); +mat_workspace *mat_workspace_init(int n, TYPE t); void mat_workspace_clear(mat_workspace *ws); -void mat_init(mat m, int n); +void mat_init(mat m, int n, TYPE t); void mat_get(NUMBER out, mat m, int i, int j); void mat_set(mat m, int i, int j, NUMBER x); NUMBER *mat_ref(mat m, int i, int j); @@ -73,10 +79,8 @@ void mat_identity(mat m); void mat_copy(mat to, mat from); void mat_clear(mat m); int mat_same(mat m1, mat m2); -static void mat_multiply_outofplace(mat_workspace *ws, mat out, mat in1, mat in2); void mat_multiply(mat_workspace *ws, mat out, mat in1, mat in2); void mat_det(mat_workspace *ws, NUMBER out, mat in); -static void mat_pseudoinverse_outofplace(mat_workspace *ws, mat out, mat in); void mat_pseudoinverse(mat_workspace *ws, mat out, mat in); void mat_trace(NUMBER out, mat in); void mat_print(mat in); diff --git a/qext.c b/qext.c index 66a0d7f..ebccf0c 100644 --- a/qext.c +++ b/qext.c @@ -7,62 +7,90 @@ mpq_t *qext_coefficient; #include #define LOOP(i,n) for(int i = 0; i < (n); i++) +#define RANK(x) ((x)->type->rank) -void qext_setup(int rank, const int *coeffs) +const int qext_type_trivial_coeffs[] = {-1, 1}; +struct qext_type qext_type_trivial_v = {1, qext_type_trivial_coeffs, NULL}; +struct qext_type *QT_TRIVIAL = &qext_type_trivial_v; + +const int qext_type_sqrt5_coeffs[] = {-5, 0, 1}; +struct qext_type qext_type_sqrt5_v = {2, qext_type_sqrt5_coeffs, NULL}; +struct qext_type *QT_SQRT5 = &qext_type_sqrt5_v; + +static void qext_init_type(struct qext_type *type) { - 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]); + type->coeffs = malloc(type->rank*sizeof(mpq_t)); + LOOP(i, type->rank) { + mpq_init(type->coeffs[i]); + mpq_set_si(type->coeffs[i], + -type->integer_coeffs[i], + type->integer_coeffs[type->rank]); } } -void qext_init(qext_number x) +struct qext_type *qext_newtype(int rank, const int *coeffs) { - x->rk = qext_rank; - x->a = malloc(x->rk*sizeof(mpq_t)); - LOOP(i, x->rk) mpq_init(x->a[i]); + struct qext_type *type = malloc(sizeof(struct qext_type)); + type->rank = rank; + type->integer_coeffs = coeffs; + qext_init_type(type); + + return type; +} + +void qext_init(qext_number x, struct qext_type *type) +{ + if(type->coeffs == NULL) // uninitialized default type + qext_init_type(type); + + x->type = type; + x->a = malloc(RANK(x)*sizeof(mpq_t)); + LOOP(i, RANK(x)) mpq_init(x->a[i]); } void qext_clear(qext_number x) { - LOOP(i, x->rk) mpq_clear(x->a[i]); + LOOP(i, RANK(x)) mpq_clear(x->a[i]); free(x->a); } void qext_set(qext_number x, qext_number y) { - LOOP(i, x->rk) mpq_set(x->a[i], y->a[i]); + LOOP(i, RANK(x)) mpq_set(x->a[i], y->a[i]); } void qext_set_int(qext_number x, int y) { mpq_set_si(x->a[0], y, 1); - for(int i = 1; i < x->rk; i++) + for(int i = 1; i < RANK(x); i++) 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++) + for(int i = 1; i < RANK(x); 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]); + LOOP(i, RANK(x)) mpq_add(result->a[i], x->a[i], y->a[i]); } void qext_sub(qext_number result, qext_number x, qext_number y) { - LOOP(i, x->rk) mpq_sub(result->a[i], x->a[i], y->a[i]); + LOOP(i, RANK(x)) mpq_sub(result->a[i], x->a[i], y->a[i]); +} + +void qext_neg(qext_number result, qext_number x) +{ + LOOP(i, RANK(x)) mpq_neg(result->a[i], x->a[i]); } void qext_print(qext_number x) { - LOOP(i, x->rk) { + LOOP(i, RANK(x)) { if(i == 0) gmp_printf("%Qd", x->a[i]); else if(i == 1) @@ -76,7 +104,7 @@ void qext_snprint(char *str, size_t size, qext_number x) { int len = 0; - LOOP(i, x->rk) { + LOOP(i, RANK(x)) { if(i == 0) len += gmp_snprintf(str+len, size-len, "%Qd", x->a[i]); else if(i == 1) @@ -90,7 +118,7 @@ int qext_cmp(qext_number x, qext_number y) { int cmp; - LOOP(i, x->rk) { + LOOP(i, RANK(x)) { cmp = mpq_cmp(x->a[i], y->a[i]); if(cmp) return cmp; @@ -101,11 +129,11 @@ int qext_cmp(qext_number x, qext_number y) void qext_mul(qext_number out, qext_number x, qext_number y) { - mpq_t result[2*x->rk-1]; + mpq_t result[2*RANK(x)-1]; mpq_t tmp; mpq_init(tmp); - LOOP(i, 2*x->rk-1) { + LOOP(i, 2*RANK(x)-1) { mpq_init(result[i]); mpq_set_ui(result[i], 0, 1); } @@ -115,21 +143,21 @@ 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, x->rk) LOOP(j, x->rk) { + LOOP(i, RANK(x)) LOOP(j, RANK(x)) { mpq_mul(tmp, x->a[i], y->a[j]); mpq_add(result[i+j], result[i+j], tmp); } - for(int i = x->rk-2; i >= 0; i--) { - LOOP(j, x->rk) { - mpq_mul(tmp, qext_coefficient[j], result[x->rk+i]); + for(int i = RANK(x)-2; i >= 0; i--) { + LOOP(j, RANK(x)) { + mpq_mul(tmp, x->type->coeffs[j], result[RANK(x)+i]); mpq_add(result[i+j], result[i+j], tmp); } } - LOOP(i, x->rk) mpq_set(out->a[i], result[i]); + LOOP(i, RANK(x)) mpq_set(out->a[i], result[i]); - LOOP(i, 2*x->rk-1) mpq_clear(result[i]); + LOOP(i, 2*RANK(x)-1) mpq_clear(result[i]); mpq_clear(tmp); } diff --git a/qext.h b/qext.h index fed4cf9..78a5622 100644 --- a/qext.h +++ b/qext.h @@ -3,21 +3,31 @@ #include +struct qext_type { + int rank; + const int *integer_coeffs; + mpq_t *coeffs; +}; + struct qext_number_internal { - int rk; + struct qext_type *type; mpq_t *a; }; typedef struct qext_number_internal qext_number[1]; -void qext_setup(int rank, const int *coeffs); -void qext_init(qext_number x); +extern struct qext_type *QT_TRIVIAL; +extern struct qext_type *QT_SQRT5; + +struct qext_type *qext_newtype(int rank, const int *coeffs); +void qext_init(qext_number x, struct qext_type *type); 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_neg(qext_number result, qext_number x); 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); diff --git a/singular_values_barbot.c b/singular_values_barbot.c index 850eb06..eca49b8 100644 --- a/singular_values_barbot.c +++ b/singular_values_barbot.c @@ -146,8 +146,8 @@ 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); + INIT(disc, QT_SQRT5); + INIT(zero, QT_SQRT5); SET_ZERO(zero); for(int i = 0; i < nuniq; i++) { discriminant(disc, invariants[i]->tr, invariants[i]->trinv); @@ -289,12 +289,12 @@ int init_node(const void *_g, void *_n) g->id_list = (int*)(g+1); // pointers get scrambled by transmission, reconstruct n->matrices = malloc(g->nmax*sizeof(mat)); for(int i = 0; i < g->nmax; i++) - mat_init(n->matrices[i], 3); + mat_init(n->matrices[i], 3, QT_SQRT5); n->invariants = malloc(g->nmax*sizeof(struct result)); n->distinct_invariants = malloc(g->nmax*sizeof(struct result)); // we won't need that many, but just in case for(int i = 0; i < g->nmax; i++) { - INIT(n->invariants[i].tr); - INIT(n->invariants[i].trinv); + INIT(n->invariants[i].tr, QT_SQRT5); + INIT(n->invariants[i].trinv, QT_SQRT5); n->invariants[i].id = i; } @@ -370,7 +370,7 @@ int do_computation(const void *_g, void *_n, const void *_in, void *_out) DEBUG("Compute matrices\n"); mat gen[6]; for(int i = 0; i < 6; i++) - mat_init(gen[i], 3); + mat_init(gen[i], 3, QT_SQRT5); generators_triangle_rotation_555_barbot(gen, s, q); enumerate_triangle_rotation_subgroup(n->group, gen, n->matrices); @@ -420,9 +420,6 @@ 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,