From 2f14076337de5decd3fb786ebb1534a050dbd625 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 2 Mar 2022 10:48:04 -0600 Subject: [PATCH] redo hyperbolic pictures --- hyperbolic.c | 205 +++++++++++++++++++++------------------------------ linalg.c | 10 ++- linalg.h | 2 +- 3 files changed, 89 insertions(+), 128 deletions(-) diff --git a/hyperbolic.c b/hyperbolic.c index 4185fc3..4bd525c 100644 --- a/hyperbolic.c +++ b/hyperbolic.c @@ -198,6 +198,16 @@ void draw_triangle(point *p, gsl_matrix *frame, const char *arguments) #endif } +void draw_dot(point p, gsl_matrix *frame, const char *arguments) +{ + double x = coord(p, 0, frame); + double y = coord(p, 1, frame); + +#ifdef POINCARE + printf("\n", CONV(x), CONV(y), arguments); +#endif +} + void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments) { char buffer[100]; @@ -217,11 +227,16 @@ void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments) #endif } -void compute_word(workspace_t *ws, gsl_matrix *result, gsl_matrix **gen, const char *word, int modifier) +void compute_word(workspace_t *ws, gsl_matrix *result, gsl_matrix **gen, const char *word, int modifier, int inverse) { gsl_matrix_set_identity(result); - for(int i = 0; word[i] != 0; i++) - multiply_right(result, gen[(word[i]-'a'+modifier)%3], ws); + if(inverse) { + for(int i = 0; word[i] != 0; i++) + multiply_left(gen[(word[i]-'a'+modifier)%3], result, ws); + } else { + for(int i = 0; word[i] != 0; i++) + multiply_right(result, gen[(word[i]-'a'+modifier)%3], ws); + } } int main(int argc, const char *argv[]) @@ -230,20 +245,22 @@ int main(int argc, const char *argv[]) gsl_matrix **matrices; gsl_matrix *cartan; gsl_matrix *gen[3]; - gsl_matrix *coxeter[3]; - gsl_matrix *coxeter2[3]; - gsl_matrix *coxeter_eigenvectors[3]; - gsl_matrix *coxeter_eigenvectors2[3]; + gsl_matrix **special; + gsl_matrix **special_eigenvectors; + point *special_attracting, *special_repelling; gsl_matrix *frame; workspace_t *ws; + int p,q,r; + int elements, nspecial; if(argc < 5) { - fprintf(stderr, "Usage: %s

\n", argv[0]); + fprintf(stderr, "Usage: %s

...\n", argv[0]); exit(1); } - int elements = atoi(argv[4]); - int p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]); + nspecial = argc - 5; + elements = atoi(argv[4]); + p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]); group = malloc(elements*sizeof(groupelement_t)); matrices = malloc(elements*sizeof(gsl_matrix*)); @@ -252,79 +269,44 @@ int main(int argc, const char *argv[]) cartan = gsl_matrix_alloc(3, 3); frame = gsl_matrix_alloc(3, 3); LOOP(i) gen[i] = gsl_matrix_alloc(3, 3); - LOOP(i) coxeter[i] = gsl_matrix_alloc(3, 3); - LOOP(i) coxeter_eigenvectors[i] = gsl_matrix_alloc(3, 3); - LOOP(i) coxeter2[i] = gsl_matrix_alloc(3, 3); - LOOP(i) coxeter_eigenvectors2[i] = gsl_matrix_alloc(3, 3); ws = workspace_alloc(3); + special = malloc(3*nspecial*sizeof(gsl_matrix*)); + special_eigenvectors = malloc(3*nspecial*sizeof(gsl_matrix*)); + special_attracting = malloc(3*nspecial*sizeof(point)); + special_repelling = malloc(3*nspecial*sizeof(point)); + for(int i = 0; i < 3*nspecial; i++) { + special[i] = gsl_matrix_alloc(3, 3); + special_eigenvectors[i] = gsl_matrix_alloc(3, 3); + } generate_triangle_group(group, elements, p, q, r); cartan_matrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0); initialize_triangle_generators(gen, cartan); - int pos = diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form + diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form gsl_matrix_set_identity(matrices[0]); for(int i = 1; i < elements; i++) multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]); -// LOOP(i) multiply_many(ws, coxeter[i], 3, gen[i%3], gen[(i+1)%3], gen[(i+2)%3]); // coxeter -// LOOP(i) multiply_many(ws, coxeter[i], 4, gen[i%3], gen[(i+1)%3], gen[i%3], gen[(i+2)%3]); // abcb - LOOP(i) compute_word(ws, coxeter[i], gen, "abcb", i); -/* LOOP(i) multiply_many(ws, coxeter[i], 10, - gen[i%3], gen[(i+1)%3], gen[(i+2)%3], - gen[i%3], gen[(i+1)%3], gen[(i+2)%3], - gen[i%3], gen[(i+1)%3], gen[(i+2)%3], - gen[(i+1)%3]); // (abc)^3 b */ - LOOP(i) eigenvectors(coxeter[i], coxeter_eigenvectors[i], ws); - LOOP(i) eigenvectors(coxeter2[i], coxeter_eigenvectors2[i], ws); + for(int i = 0; i < nspecial; i++) + LOOP(j) { + compute_word(ws, special[3*i+j], gen, argv[i+5], j, 0); + eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws); + special_attracting[3*i+j] = column(special_eigenvectors[3*i+j], 0); + + // repelling = attracting of inverse + compute_word(ws, special[3*i+j], gen, argv[i+5], j, 1); + eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws); + special_repelling[3*i+j] = column(special_eigenvectors[3*i+j], 0); + } /* - for(int i = 0; i < elements; i++) { - printf("%4d: ", i); - for(groupelement_t *cur = &group[i]; cur->parent; cur = cur->parent) - fputc(cur->letter+'a', stdout); - fputc('\n', stdout); - } - - return 0;*/ - - point coxeter_attracting[3]; - point coxeter_repelling[3]; - point coxeter_axes[3]; - point coxeter2_attracting[3]; - point coxeter2_repelling[3]; - point coxeter2_axes[3]; - point edge_midpoints[3]; - point reflection_lines[3]; - point triangle_points[3]; - point transformed[3]; - point transformed2[3]; - point transformed3[3]; - point transformed4[3]; - point transformed5[3]; - point transformed6[3]; - point center; - - LOOP(i) coxeter_attracting[i] = column(coxeter_eigenvectors[i], 0); - LOOP(i) coxeter_repelling[i] = column(coxeter_eigenvectors[i], 2); - LOOP(i) coxeter_axes[i] = incidence(coxeter_attracting[i], coxeter_repelling[i]); - LOOP(i) edge_midpoints[i] = incidence(coxeter_axes[(i+1)%3], coxeter_axes[(i+2)%3]); - LOOP(i) reflection_lines[i] = row(cartan, i); - LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]); - LOOP(i) coxeter2_attracting[i] = column(coxeter_eigenvectors2[i], 0); - LOOP(i) coxeter2_repelling[i] = column(coxeter_eigenvectors2[i], 2); - LOOP(i) coxeter2_axes[i] = incidence(coxeter2_attracting[i], coxeter2_repelling[i]); - - print_svg_header(); - // let's correct the frame of reference by using hyperbolic transformations center = apply(frame, triangle_points[2]); double angle = atan2(center.x[1], center.x[0]); double boost = atanh(-sqrt(center.x[0]*center.x[0]+center.x[1]*center.x[1])/center.x[2]); gsl_matrix *frame_correction = gsl_matrix_alloc(3, 3); gsl_matrix_set_identity(frame_correction); - - /* gsl_matrix_set(frame_correction, 0, 0, cos(angle-M_PI/2)); gsl_matrix_set(frame_correction, 0, 1, sin(angle-M_PI/2)); gsl_matrix_set(frame_correction, 1, 0, -sin(angle-M_PI/2)); @@ -338,80 +320,57 @@ int main(int argc, const char *argv[]) gsl_matrix_set(frame_correction, 2, 0, sinh(-boost)); gsl_matrix_set(frame_correction, 2, 2, cosh(-boost)); multiply_left(frame_correction, frame, ws); - */ gsl_matrix_free(frame_correction); + */ -// int indices[10] = {0, 1, 6, 10, 30, 46, 124, 185, 484, 717}; -// int indices[10] = {0, 1, 4, 10, 22}; + // the actual drawing + point transformed[3]; + point reflection_lines[3]; + point triangle_points[3]; + + print_svg_header(); for(int k = 0; k < elements; k++) { if(group[k].length % 2) continue; + LOOP(i) reflection_lines[i] = row(cartan, i); + LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]); LOOP(i) transformed[i] = apply(matrices[k], triangle_points[i]); -// draw_triangle(transformed, frame, "black,fill=black!10,line width=0pt"); draw_triangle(transformed, frame, "fill:#cfcfcf;"); -// draw_triangle(transformed, frame, "fill:#000000;"); } + // draw special elements for(int k = 0; k < elements; k++) { // if(group[k].length % 2) // continue; - LOOP(i) transformed[i] = apply(matrices[k], edge_midpoints[(i+2)%3]); - LOOP(i) transformed2[i] = apply(matrices[k], coxeter_repelling[i]); - LOOP(i) transformed3[i] = apply(matrices[k], coxeter_repelling[(i+1)%3]); - LOOP(i) transformed4[i] = apply(matrices[k], coxeter_attracting[i%3]); - LOOP(i) transformed5[i] = apply(matrices[k], coxeter2_repelling[i]); - LOOP(i) transformed6[i] = apply(matrices[k], coxeter2_attracting[i%3]); - //LOOP(i) draw_line(transformed2[i], transformed4[i], frame, "red"); - draw_line(transformed2[0], transformed4[0], frame, "fill:none;stroke:red;stroke-width:1;"); -// draw_line(transformed2[1], transformed4[1], frame, "fill:none;stroke:blue;stroke-width:1;"); -// draw_line(transformed2[2], transformed4[2], frame, "fill:none;stroke:darkgreen;stroke-width:1;"); - -// draw_line(transformed5[0], transformed6[0], frame, "fill:none;stroke:blue;stroke-width:1;"); -// draw_line(transformed5[1], transformed6[1], frame, "fill:none;stroke:blue;stroke-width:1;"); -// draw_line(transformed5[2], transformed6[2], frame, "fill:none;stroke:blue;stroke-width:1;"); -// draw_line(transformed2[1], transformed4[1], frame, "fill:none;stroke:darkgreen;stroke-width:1;"); -// draw_line(transformed2[2], transformed4[2], frame, "fill:none;stroke:blue;stroke-width:1;"); -// LOOP(i) draw_line(transformed[i], transformed3[i], frame, "red"); - -// LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]); -// draw_line(transformed[1], transformed[2], frame, "red"); + for(int i = 0; i < 3*nspecial; i++) { +// draw_dot(apply(matrices[k], special_repelling[i]), frame, "fill:red;stroke-width:0;"); +// draw_dot(apply(matrices[k], special_attracting[i]), frame, "fill:blue;stroke-width:0;"); + draw_line(apply(matrices[k], special_repelling[i]), + apply(matrices[k], special_attracting[i]), + frame, "fill:none;stroke:red;stroke-width:1;"); + } } - /* - draw_line(apply(matrices[0], coxeter_repelling[0]), - apply(matrices[0], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[1], coxeter_repelling[0]), - apply(matrices[1], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[2], coxeter_repelling[0]), - apply(matrices[2], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[3], coxeter_repelling[0]), - apply(matrices[3], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[4], coxeter_repelling[0]), - apply(matrices[4], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[5], coxeter_repelling[0]), - apply(matrices[5], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[14], coxeter_repelling[0]), - apply(matrices[14], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[15], coxeter_repelling[0]), - apply(matrices[15], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[22], coxeter_repelling[0]), - apply(matrices[22], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - draw_line(apply(matrices[23], coxeter_repelling[0]), - apply(matrices[23], coxeter_attracting[0]), - frame, "fill:none;stroke:red;stroke-width:1;"); - */ - print_svg_footer(); + + // clean up + free(group); + for(int i = 0; i < elements; i++) + gsl_matrix_free(matrices[i]); + free(matrices); + gsl_matrix_free(cartan); + gsl_matrix_free(frame); + LOOP(i) gsl_matrix_free(gen[i]); + workspace_free(ws); + for(int i = 0; i < 3*nspecial; i++) { + gsl_matrix_free(special[i]); + gsl_matrix_free(special_eigenvectors[i]); + } + free(special); + free(special_eigenvectors); + free(special_attracting); + free(special_repelling); } diff --git a/linalg.c b/linalg.c index ba0a30e..7d06f43 100644 --- a/linalg.c +++ b/linalg.c @@ -145,11 +145,11 @@ double determinant(gsl_matrix *g, workspace_t *ws) return gsl_linalg_LU_det(ws->tmp, s); } -void jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws) +int jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws) { gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv); gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); - gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_VAL_DESC); + gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); int real = 1; for(int i = 0; i < ws->n; i++) @@ -165,6 +165,8 @@ void jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws) evs[2*i] = GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i)); evs[2*i+1] = GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)); } + + return real; } void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) @@ -182,8 +184,8 @@ void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) real = 0; if(!real) { - printf("We have non-real eigenvalues!\n"); - exit(1); + fprintf(stderr,"We have non-real eigenvalues!\n"); +// exit(1); } for(int i = 0; i < ws->n; i++) diff --git a/linalg.h b/linalg.h index 11cbe80..1077090 100644 --- a/linalg.h +++ b/linalg.h @@ -38,7 +38,7 @@ void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...); void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws); void initialize(gsl_matrix *g, double *data, int x, int y); void rotation_matrix(gsl_matrix *g, double *vector); -void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws); +int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws); double trace(gsl_matrix *g); double determinant(gsl_matrix *g, workspace_t *ws); void eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);