From 244784794d1a355143e53bf588ffc0f2414ab41d Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Tue, 14 Jun 2022 14:22:22 +0200 Subject: [PATCH] compute complex traces --- .gitignore | 18 ++---- complex_anosov.c | 141 ++++++++++++++++++++++++++++++++++++++++++++++- coxeter.c | 16 ++++++ coxeter.h | 2 + enumerate.c | 98 ++++++++++++++++++++++++++++++++ enumerate.h | 9 +++ generators.c | 103 +++++++++++++++++++++++----------- generators.h | 6 +- mat.h | 4 +- qext.c | 102 +++++++++++++++++++++++++++++++--- qext.h | 6 +- 11 files changed, 446 insertions(+), 59 deletions(-) diff --git a/.gitignore b/.gitignore index 0680d5f..c6c8671 100644 --- a/.gitignore +++ b/.gitignore @@ -1,17 +1,7 @@ *.o -triangle_group/singular_values -.#* -singular_values -singular_values_mpi -singular_values_barbot -output/ -special_element -max_slope_picture/generate -convert -billiard_words -*.pnm -*.png -*.hi +.\#* +\#*\# gmon.out -restart core +output/ +complex_anosov diff --git a/complex_anosov.c b/complex_anosov.c index 62c2301..7ecc5cc 100644 --- a/complex_anosov.c +++ b/complex_anosov.c @@ -5,9 +5,148 @@ #include "qext.h" #include +#include +#include + +#define LOOP(i,n) for(int i = 0; i < (n); i++) + +/* + Elements up to length 0: 1 + Elements up to length 1: 4 + Elements up to length 2: 10 + Elements up to length 3: 22 + Elements up to length 4: 46 + Elements up to length 5: 91 + Elements up to length 6: 175 + Elements up to length 7: 334 + Elements up to length 8: 634 + Elements up to length 9: 1198 + Elements up to length 10: 2260 + Elements up to length 11: 4261 + Elements up to length 12: 8029 + Elements up to length 13: 15124 + Elements up to length 14: 28486 + Elements up to length 15: 53650 + Elements up to length 16: 101038 + Elements up to length 17: 190279 + Elements up to length 18: 358339 + Elements up to length 19: 674830 + Elements up to length 20: 1270846 + Elements up to length 21: 2393266 + Elements up to length 22: 4507012 + Elements up to length 23: 8487625 +*/ + +static double gaussian_sqrt5_real(NUMBER x) +{ + double result = 0.0; + + mpq_t tmp; + mpq_init(tmp); + + // a_0 + sqrt(5)a_1 + 4a_2 + 2sqrt(5)a_3 + mpq_set_si(tmp, 4, 1); + mpq_mul(tmp, tmp, x->a[2]); + mpq_add(tmp, tmp, x->a[0]); + result = mpq_get_d(tmp); + mpq_set_si(tmp, 2, 1); + mpq_mul(tmp, tmp, x->a[3]); + mpq_add(tmp, tmp, x->a[1]); + result += mpq_get_d(tmp)*sqrt(5); + + mpq_clear(tmp); + + return result; +} + +static double gaussian_sqrt5_imag(NUMBER x) +{ + double result = 0.0; + + mpq_t tmp; + mpq_init(tmp); + + // a_1 + 2sqrt(5)a_2 + 14a_3 + mpq_set_si(tmp, 14, 1); + mpq_mul(tmp, tmp, x->a[3]); + mpq_add(tmp, tmp, x->a[1]); + result = mpq_get_d(tmp); + mpq_set_si(tmp, 2, 1); + mpq_mul(tmp, tmp, x->a[2]); + result += mpq_get_d(tmp)*sqrt(5); + + mpq_clear(tmp); + + return result; +} int main(int argc, char *argv[]) { - printf("initial main function.\n"); + char buf[100]; + int n = atoi(argv[1]); + + mpq_t qreal, qimag; + mpq_inits(qreal, qimag, NULL); + mpq_set_si(qreal, 50, 10); + mpq_set_si(qimag, 1, 10); + mpq_canonicalize(qreal); + mpq_canonicalize(qimag); + + /* + int length = 0; + LOOP(i, n) { + if(group->elements[i].length > length) { + printf("Elements up to length %d: %d\n", length, i); + length = group->elements[i].length; + } + } + return 0; + */ + + /* + LOOP(i, n) { + groupelement_t *cur = &group->elements[i]; + groupelement_t *other; + + cur->conjugacy_class = cur; // start with itself and reduce if possible + + LOOP(j, 3) { + if(cur->left[j] && cur->left[j]->right[j]) { + other = cur->left[j]->right[j]; + if(other->id < cur->id) + cur->conjugacy_class = other->conjugacy_class; + } + if(cur->right[j] && cur->right[j]->left[j]) { + other = cur->right[j]->left[j]; + if(other->id < cur->id) + cur->conjugacy_class = other->conjugacy_class; + } + } + } + */ + + group_t *group = coxeter_init_triangle(5, 5, 5, n); + + mat gen[3]; + LOOP(i, 3) mat_init(gen[i], 3, QT_GAUSS_SQRT5); + generators_triangle_reflection_group_555_complex(gen, 2, 2, 2, qreal, qimag); + + struct tracedata *traces; + int nuniq = enumerate_coxeter_group_traces(group, gen, &traces); + + LOOP(i, nuniq) { + printf("%d %f %f %f %f\n", + traces[i].id, + gaussian_sqrt5_real(traces[i].tr), + gaussian_sqrt5_imag(traces[i].tr), + gaussian_sqrt5_real(traces[i].trinv), + gaussian_sqrt5_imag(traces[i].trinv)); + } + + enumerate_tracedata_clear(traces, nuniq); + LOOP(i, 3) mat_clear(gen[i]); + coxeter_clear(group); + mpq_clears(qreal, qimag, NULL); + return 0; } diff --git a/coxeter.c b/coxeter.c index 128d9cd..38bc13e 100644 --- a/coxeter.c +++ b/coxeter.c @@ -5,6 +5,8 @@ #include "coxeter.h" #define LOOP(i,n) for(int i = 0; i < (n); i++) +#define MIN(x,y) ((x) > (y) ? (y) : (x)) +#define MAX(x,y) ((x) > (y) ? (x) : (y)) group_t *coxeter_init_triangle(int p, int q, int r, int nmax) { @@ -185,3 +187,17 @@ void coxeter_clear(group_t *g) free(g->lists); free(g); } + +int coxeter_snprint(char *str, int buflen, groupelement_t *g) +{ + int n = MIN(g->length, buflen-1); // number of characters without null byte + + str[n] = 0; + int i = n-1; + while(g->parent) { + str[i--] = 'a' + g->letter; + g = g->parent; + } + + return n; +} diff --git a/coxeter.h b/coxeter.h index ead28b5..89de856 100644 --- a/coxeter.h +++ b/coxeter.h @@ -13,6 +13,7 @@ typedef struct _groupelement { struct _groupelement *inverse; struct _groupelement **left; struct _groupelement **right; + struct _groupelement *conjugacy_class; } groupelement_t; typedef struct _group { @@ -26,5 +27,6 @@ typedef struct _group { group_t *coxeter_init(int rank, int *coxeter_matrix, int nmax); group_t *coxeter_init_triangle(int p, int q, int r, int nmax); void coxeter_clear(group_t *g); +int coxeter_snprint(char *str, int buflen, groupelement_t *g); #endif diff --git a/enumerate.c b/enumerate.c index 90ea826..f779988 100644 --- a/enumerate.c +++ b/enumerate.c @@ -1,5 +1,7 @@ #include "enumerate.h" +#include + #define LOOP(i,n) for(int i = 0; i < (n); i++) void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices) @@ -22,3 +24,99 @@ void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices) mat_workspace_clear(ws); } + +static int compare_tracedata(const void *a_, const void *b_) +{ + int d = 0; + + struct tracedata **a = (struct tracedata **)a_; + struct tracedata **b = (struct tracedata **)b_; + + d = CMP((*a)->tr,(*b)->tr); + if(d == 0) { + d = CMP((*a)->trinv, (*b)->trinv); + } + + return d; +} + +static int compare_tracedata_id(const void *a, const void *b) +{ + int ida = (*(struct tracedata **)a)->id; + int idb = (*(struct tracedata **)b)->id; + + return ida > idb ? 1 : ida < idb ? -1 : 0; +} + +int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out) +{ + TYPE t = GETTYPE(gen[0]->x[0]); + int n = group->size; + + mat *matrices = malloc(n*sizeof(mat)); + struct tracedata *traces = malloc(n*sizeof(struct tracedata)); + struct tracedata **distinct_traces = malloc(n*sizeof(struct tracedata*)); + + LOOP(i, n) mat_init(matrices[i], 3, t); + enumerate_coxeter_group(group, gen, matrices); + + LOOP(i, n) { + INIT(traces[i].tr, t); + INIT(traces[i].trinv, t); + } + LOOP(i, n) { + if(!group->elements[i].inverse) + continue; + mat_trace(traces[i].tr, matrices[i]); + mat_trace(traces[i].trinv, matrices[group->elements[i].inverse->id]); + traces[i].id = i; + distinct_traces[i] = &traces[i]; + } + + qsort(distinct_traces, n, sizeof(struct tracedata*), compare_tracedata); + + int nuniq = 0; + LOOP(i, n) { + if(i == 0 || compare_tracedata(&distinct_traces[i], &distinct_traces[nuniq-1]) != 0) { + distinct_traces[nuniq] = distinct_traces[i]; + nuniq++; + } else { + int oldlength = group->elements[distinct_traces[nuniq-1]->id].length; + int newlength = group->elements[distinct_traces[i]->id].length; + if(newlength < oldlength) + distinct_traces[nuniq-1]->id = distinct_traces[i]->id; + } + } + + qsort(distinct_traces, nuniq, sizeof(struct tracedata*), compare_tracedata_id); + + struct tracedata *unique_traces = malloc(nuniq*sizeof(struct tracedata)); + LOOP(i, nuniq) { + INIT(unique_traces[i].tr, t); + INIT(unique_traces[i].trinv, t); + SET(unique_traces[i].tr, distinct_traces[i]->tr); + SET(unique_traces[i].trinv, distinct_traces[i]->trinv); + unique_traces[i].id = distinct_traces[i]->id; + } + + LOOP(i, n) mat_clear(matrices[i]); + LOOP(i, n) { + CLEAR(traces[i].tr); + CLEAR(traces[i].trinv); + } + free(matrices); + free(traces); + free(distinct_traces); + + *traces_out = unique_traces; + return nuniq; +} + +void enumerate_tracedata_clear(struct tracedata *traces, int n) +{ + LOOP(i, n) { + CLEAR(traces[i].tr); + CLEAR(traces[i].trinv); + } + free(traces); +} diff --git a/enumerate.h b/enumerate.h index 6e1a83c..cca03e9 100644 --- a/enumerate.h +++ b/enumerate.h @@ -4,6 +4,15 @@ #include "mat.h" #include "coxeter.h" +struct tracedata { + int id; + NUMBER tr; + NUMBER trinv; +}; + void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices); +int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out); + +void enumerate_tracedata_clear(struct tracedata *traces, int n); #endif diff --git a/generators.c b/generators.c index 5639e2d..da958af 100644 --- a/generators.c +++ b/generators.c @@ -1,33 +1,17 @@ #include "generators.h" +#include "mat.h" #define LOOP(i,n) for(int i = 0; i < (n); i++) -void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, mpq_t q) +void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER q) { - NUMBER b1,c1,a2,c2,a3,b3; - mpq_t qinv; + NUMBER tmp; + NUMBER qinv; TYPE t = GETTYPE(rho1); - INIT(b1,t);INIT(c1,t);INIT(a2,t);INIT(c2,t);INIT(a3,t);INIT(b3,t); - - // qinv = q^{-1} - mpq_init(qinv); - mpq_inv(qinv, q); - - // c1 = rho2 q, a2 = rho3 q, b3 = rho1 q, b1 = c2 = a3 = 1/q - SET_ZERO(c1); - 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); - SET_Q(b1, qinv); - SET_Q(c2, qinv); - SET_Q(a3, qinv); + INIT(tmp, t); + INIT(qinv, t); + INV(qinv, q); LOOP(i, 3) { mat_init(gen[i], 3, t); @@ -37,18 +21,21 @@ void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, } } - NEG(*mat_ref(gen[0], 1, 0), b1); - NEG(*mat_ref(gen[0], 2, 0), c1); - NEG(*mat_ref(gen[1], 0, 1), a2); - NEG(*mat_ref(gen[1], 2, 1), c2); - NEG(*mat_ref(gen[2], 0, 2), a3); - NEG(*mat_ref(gen[2], 1, 2), b3); + MUL(tmp, q, rho1); + NEG(*mat_ref(gen[2], 1, 2), tmp); + MUL(tmp, q, rho2); + NEG(*mat_ref(gen[0], 2, 0), tmp); + MUL(tmp, q, rho3); + NEG(*mat_ref(gen[1], 0, 1), tmp); - CLEAR(b1);CLEAR(c1);CLEAR(a2);CLEAR(c2);CLEAR(a3);CLEAR(b3); - mpq_clear(qinv); + NEG(*mat_ref(gen[1], 2, 1), qinv); + NEG(*mat_ref(gen[2], 0, 2), qinv); + NEG(*mat_ref(gen[0], 1, 0), qinv); + + CLEAR(tmp); + CLEAR(qinv); } - int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t t) { int p[3] = {p1, p2, p3}; @@ -103,9 +90,59 @@ int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q } } - generators_triangle_reflection_generic(gen, rho[0], rho[1], rho[2], t); + NUMBER param; + INIT(param, type); + SET_Q(param, t); + + generators_triangle_reflection_generic(gen, rho[0], rho[1], rho[2], param); LOOP(i, 3) CLEAR(rho[i]); + CLEAR(param); return 1; } + +int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, int q3, mpq_t treal, mpq_t timag) +{ + + int q[3] = {q1, q2, q3}; + + NUMBER rho[3]; + LOOP(i, 3) INIT(rho[i], QT_GAUSS_SQRT5); + + // compute rho1, rho2, rho3 + LOOP(i, 3) { + if(q[i] == 1) { + // 4cos(pi/5)^2 = 3/2 + 1/2*sqrt(5) + mpq_set_si(rho[i]->a[0], 3, 2); + mpq_set_si(rho[i]->a[1], 7, 12); + mpq_set_si(rho[i]->a[2], 0, 1); + mpq_set_si(rho[i]->a[3], -1, 24); + } else if(q[i] == 2) { + // 4cos(pi/5)^2 = 3/2 - 1/2*sqrt(5) + mpq_set_si(rho[i]->a[0], 3, 2); + mpq_set_si(rho[i]->a[1], -7, 12); + mpq_set_si(rho[i]->a[2], 0, 1); + mpq_set_si(rho[i]->a[3], 1, 24); + } else { + return 0; + } + } + + NUMBER param; + INIT(param, QT_GAUSS_SQRT5); + mpq_set(param->a[0], treal); + mpq_set_si(param->a[1], -1, 6); + mpq_mul(param->a[1], param->a[1], timag); + mpq_set_si(param->a[2], 0, 1); + mpq_set_si(param->a[3], 1, 12); + mpq_mul(param->a[3], param->a[3], timag); + + generators_triangle_reflection_generic(gen, rho[0], rho[1], rho[2], param); + + LOOP(i, 3) CLEAR(rho[i]); + CLEAR(param); + + return 1; + +} diff --git a/generators.h b/generators.h index ca97d52..853eb3f 100644 --- a/generators.h +++ b/generators.h @@ -3,7 +3,11 @@ #include "mat.h" -void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, mpq_t q); +void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER q); int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t q); +#ifdef QEXT +int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, int q3, mpq_t real, mpq_t imag); +#endif + #endif diff --git a/mat.h b/mat.h index fcc2b0b..18db837 100644 --- a/mat.h +++ b/mat.h @@ -24,7 +24,8 @@ #define SUB qext_sub #define NEG qext_neg #define MUL qext_mul -// #define DIV qext_div +#define DIV qext_div +#define INV qext_inv #define CMP qext_cmp #define PRINT qext_print #define SNPRINT qext_snprint @@ -46,6 +47,7 @@ #define NEG mpq_neg #define MUL mpq_mul #define DIV mpq_div +#define INV mpq_inv #define PRINT(x) gmp_printf("%Qd", x) #define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x) #define TYPE int diff --git a/qext.c b/qext.c index ebccf0c..0864aa2 100644 --- a/qext.c +++ b/qext.c @@ -1,5 +1,7 @@ #include "qext.h" +#include + int qext_rank; mpq_t *qext_coefficient; @@ -9,13 +11,20 @@ mpq_t *qext_coefficient; #define LOOP(i,n) for(int i = 0; i < (n); i++) #define RANK(x) ((x)->type->rank) -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; +struct qext_type *QT_TRIVIAL = &(struct qext_type) { + .rank = 1, + .integer_coeffs = (const int[]) {-1, 1} +}; -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; +struct qext_type *QT_SQRT5 = &(struct qext_type) { + .rank = 2, + .integer_coeffs = (const int[]) {-5, 0, 1} +}; + +struct qext_type *QT_GAUSS_SQRT5 = &(struct qext_type) { + .rank = 4, + .integer_coeffs = (const int[]) {36, 0, -8, 0, 1} +}; static void qext_init_type(struct qext_type *type) { @@ -161,7 +170,86 @@ void qext_mul(qext_number out, qext_number x, qext_number y) mpq_clear(tmp); } +void qext_inv(qext_number out, qext_number in) +{ + qext_number one; + qext_init(one, in->type); + qext_set_int(one, 1); + qext_div(out, one, in); + qext_clear(one); +} + void qext_div(qext_number out, qext_number x, qext_number y) { - // todo: implement this by solving a linear system + int rk = RANK(x); + struct qext_type *type = x->type; + + mpq_t tmp; + mpq_t result[rk]; + mpq_t matrix[rk*rk]; + + mpq_init(tmp); + LOOP(i, rk) mpq_init(result[i]); + LOOP(i, rk*rk) mpq_init(matrix[i]); + + // initialize a matrix which represents multiplication by y + // that means the i-th column is y*a^i + LOOP(i, rk) + mpq_set(matrix[i*rk], y->a[i]); + LOOP(j, rk-1) { // columns + mpq_mul(matrix[j+1], matrix[(rk-1)*rk+j], type->coeffs[0]); + LOOP(i, rk-1) { // rows + mpq_mul(matrix[(i+1)*rk+(j+1)], matrix[(rk-1)*rk+j], type->coeffs[i+1]); + mpq_add(matrix[(i+1)*rk+(j+1)], matrix[(i+1)*rk+(j+1)], matrix[i*rk+j]); + } + } + + // initialize result as x + LOOP(i, rk) mpq_set(result[i], x->a[i]); + + // use Gaussian elimination to solve the system matrix * out = result + LOOP(j, rk) + { + // find nonzero entry + int row = j; + while(row < rk && mpq_sgn(matrix[row*rk+j]) == 0) row++; + + // if the entire column is zero, the matrix was not invertible + // this could happen if the polynomial is not irreducible / we're not in a field + assert(row != rk); + + // permute rows + if(row != j) { + mpq_swap(result[row], result[j]); + LOOP(k, rk) { + mpq_swap(matrix[row*rk+k], matrix[j*rk+k]); + } + } + + // normalize + LOOP(k, rk) + if(k != j) + mpq_div(matrix[j*rk+k], matrix[j*rk+k], matrix[j*rk+j]); + mpq_div(result[j], result[j], matrix[j*rk+j]); + mpq_set_ui(matrix[j*rk+j], 1, 1); + + // subtract + LOOP(i, rk) { + if(i == j) + continue; + mpq_mul(tmp, matrix[i*rk+j], result[j]); + mpq_sub(result[i], result[i], tmp); + for(int k = j+1; k < rk; k++) { + mpq_mul(tmp, matrix[i*rk+j], matrix[j*rk+k]); + mpq_sub(matrix[i*rk+k], matrix[i*rk+k], tmp); + } + mpq_set_ui(matrix[i*rk+j], 0, 1); // it isn't strictly necessary to do this as we won't use this entry again, but helps debugging + } + } + + LOOP(i, rk) mpq_set(out->a[i], result[i]); + + mpq_clear(tmp); + LOOP(i, rk) mpq_clear(result[i]); + LOOP(i, rk*rk) mpq_clear(matrix[i]); } diff --git a/qext.h b/qext.h index 78a5622..7181353 100644 --- a/qext.h +++ b/qext.h @@ -5,8 +5,8 @@ struct qext_type { int rank; - const int *integer_coeffs; - mpq_t *coeffs; + const int *integer_coeffs; // actually the rank+1 coefficients of the polynomial + mpq_t *coeffs; // components of a^rank }; struct qext_number_internal { @@ -18,6 +18,7 @@ typedef struct qext_number_internal qext_number[1]; extern struct qext_type *QT_TRIVIAL; extern struct qext_type *QT_SQRT5; +extern struct qext_type *QT_GAUSS_SQRT5; struct qext_type *qext_newtype(int rank, const int *coeffs); void qext_init(qext_number x, struct qext_type *type); @@ -29,6 +30,7 @@ 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_inv(qext_number out, qext_number in); 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);