From ab2a83b0bcee754519393500fcab37e585405d28 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sat, 26 May 2018 09:19:03 +0200 Subject: [PATCH 1/4] first try of hyperbolic triangle group --- Makefile | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Makefile b/Makefile index d50572b..ecd48d8 100644 --- a/Makefile +++ b/Makefile @@ -7,7 +7,7 @@ SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS) -all: singular_values element_path limit_set +all: singular_values element_path limit_set hyperbolic singular_values: singular_values.o triangle.o linalg.o gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas @@ -18,6 +18,9 @@ element_path: element_path.o linalg.o limit_set: limit_set.o linalg.o triangle.o gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o -lm -lgsl -lcblas +hyperbolic: hyperbolic.o triangle.o linalg.o + gcc $(OPTIONS) -o hyperbolic hyperbolic.o triangle.o linalg.o -lm -lgsl -lcblas + singular_values.o: singular_values.c $(HEADERS) gcc $(OPTIONS) -c singular_values.c @@ -33,5 +36,8 @@ triangle.o: triangle.c $(HEADERS) limit_set.o: limit_set.c $(HEADERS) gcc $(OPTIONS) -c limit_set.c +hyperbolic.o: hyperbolic.c $(HEADERS) + gcc $(OPTIONS) -c hyperbolic.c + clean: - rm -f singular_values element_path limit_set triangle.o linalg.o singular_values.o element_path.o limit_set.o + rm -f singular_values element_path limit_set hyperbolic triangle.o linalg.o singular_values.o element_path.o limit_set.o hyperbolic.o From aea8496d6364d05b00801241e388eb0dc863735d Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Thu, 31 May 2018 00:13:53 +0200 Subject: [PATCH 2/4] these boxes could work! --- Makefile | 4 +-- linalg.c | 75 +++++++++++++++++++++++++++++++++++++++++++++++++++--- linalg.h | 27 ++++++++++++-------- triangle.c | 1 + 4 files changed, 91 insertions(+), 16 deletions(-) diff --git a/Makefile b/Makefile index ecd48d8..f651f9b 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,8 @@ HEADERS=triangle.h linalg.h queue.h -#SPECIAL_OPTIONS=-O0 -g -D_DEBUG +SPECIAL_OPTIONS=-O0 -g -D_DEBUG #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline -SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline +#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline #SPECIAL_OPTIONS= OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS) diff --git a/linalg.c b/linalg.c index 66ef82d..69b857f 100644 --- a/linalg.c +++ b/linalg.c @@ -18,6 +18,7 @@ workspace_t *workspace_alloc(int n) workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t)); result->n = n; result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n); + result->work_symmv = gsl_eigen_symmv_alloc(n); result->work_sv = gsl_vector_alloc(n); result->eval_complex = gsl_vector_complex_alloc(n); result->evec_complex = gsl_matrix_complex_alloc(n, n); @@ -34,6 +35,7 @@ workspace_t *workspace_alloc(int n) void workspace_free(workspace_t *workspace) { gsl_eigen_nonsymmv_free(workspace->work_nonsymmv); + gsl_eigen_symmv_free(workspace->work_symmv); gsl_vector_free(workspace->work_sv); gsl_vector_complex_free(workspace->eval_complex); gsl_matrix_complex_free(workspace->evec_complex); @@ -67,6 +69,33 @@ void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out) gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, out); } +void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws) +{ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, ws->stack[0]); + gsl_matrix_memcpy(a, ws->stack[0]); +} + +void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws) +{ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, ws->stack[0]); + gsl_matrix_memcpy(b, ws->stack[0]); +} + +void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...) +{ + va_list args; + va_start(args, n); + + gsl_matrix_set_identity(out); + + for(int i = 0; i < n; i++) { + gsl_matrix *cur = va_arg(args, gsl_matrix *); + multiply_right(out, cur, ws); + } + + va_end(args); +} + void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws) { gsl_matrix_memcpy(ws->tmp, g); @@ -139,14 +168,15 @@ void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws) void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) { + gsl_matrix_memcpy(ws->stack[0], g); gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); - int r = gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); + int r = gsl_eigen_nonsymmv(ws->stack[0], ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); ERROR(r, "gsl_eigen_nonsymmv failed!\n"); gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); int real = 1; - for(int i = 0; i < 3; i++) + for(int i = 0; i < ws->n; i++) if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) real = 0; @@ -155,8 +185,45 @@ void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) exit(1); } - for(int i = 0; i < 3; i++) - for(int j = 0; j < 3; j++) + for(int i = 0; i < ws->n; i++) + for(int j = 0; j < ws->n; j++) gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j))); } + +void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws) +{ + gsl_matrix_memcpy(ws->stack[0], g); + int r = gsl_eigen_symmv (ws->stack[0], eval, evec, ws->work_symmv); + ERROR(r, "gsl_eigen_symmv failed!\n"); + + gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC); +} + +// returns number of positive directions and matrix transforming TO diagonal basis +int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws) +{ + gsl_matrix_memcpy(ws->stack[0], A); + int r = gsl_eigen_symmv (ws->stack[0], ws->eval_real, cob, ws->work_symmv); + ERROR(r, "gsl_eigen_symmv failed!\n"); + + gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC); + + gsl_matrix_transpose(cob); + + int positive = 0; + + for(int i = 0; i < ws->n; i++) { + if(gsl_vector_get(ws->eval_real, i) > 0) + positive++; + + for(int j = 0; j < ws->n; j++) + *gsl_matrix_ptr(cob, i, j) *= sqrt(fabs(gsl_vector_get(ws->eval_real, i))); + } + + return positive; + +// printf("Eigenvalues: %.10f, %.10f, %.10f\n", gsl_vector_get(ws->eval_real, 0), gsl_vector_get(ws->eval_real, 1), gsl_vector_get(ws->eval_real, 2)); + +// return 0; +} diff --git a/linalg.h b/linalg.h index b6cdb1f..11cbe80 100644 --- a/linalg.h +++ b/linalg.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -13,16 +14,17 @@ #define FCMP(x, y) gsl_fcmp(x, y, 1e-10) typedef struct _workspace { - int n; - gsl_eigen_nonsymmv_workspace *work_nonsymmv; - gsl_vector *work_sv; - gsl_vector_complex *eval_complex; - gsl_matrix_complex *evec_complex; - gsl_vector *eval_real; - gsl_matrix *evec_real; - gsl_matrix *tmp; - gsl_permutation *permutation; - gsl_matrix *stack[20]; + int n; + gsl_eigen_nonsymmv_workspace *work_nonsymmv; + gsl_eigen_symmv_workspace *work_symmv; + gsl_vector *work_sv; + gsl_vector_complex *eval_complex; + gsl_matrix_complex *evec_complex; + gsl_vector *eval_real; + gsl_matrix *evec_real; + gsl_matrix *tmp; + gsl_permutation *permutation; + gsl_matrix *stack[20]; } workspace_t; workspace_t *workspace_alloc(int n); @@ -30,6 +32,9 @@ void workspace_free(workspace_t *workspace); void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws); void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws); void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out); +void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws); +void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws); +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); @@ -37,5 +42,7 @@ void 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); +void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws); +int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws); #endif diff --git a/triangle.c b/triangle.c index 5b23d10..770c7fe 100644 --- a/triangle.c +++ b/triangle.c @@ -28,6 +28,7 @@ int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int int k[3] = {k1, k2, k3}; memset(group, 0, size*sizeof(groupelement_t)); + group[0].letter = -1; n = 1; queue_init(&q); queue_put(&q, 0); From 442d14ce0c7fb818cc981e21e23ca9276ef6c2d8 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Thu, 31 May 2018 00:17:32 +0200 Subject: [PATCH 3/4] hyperbolic tiling --- hyperbolic.c | 265 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 265 insertions(+) create mode 100644 hyperbolic.c diff --git a/hyperbolic.c b/hyperbolic.c new file mode 100644 index 0000000..53c2450 --- /dev/null +++ b/hyperbolic.c @@ -0,0 +1,265 @@ +#include +#include +#include + +#include "triangle.h" +#include "linalg.h" + +#define POINCARE 1 +#define LOOP(i) for(int i = 0; i < 3; i++) + +void cartan_matrix(gsl_matrix *cartan, double a1, double a2, double a3, double s) +{ + gsl_matrix_set(cartan, 0, 0, -2); + gsl_matrix_set(cartan, 0, 1, 2*s*cos(a3)); + gsl_matrix_set(cartan, 0, 2, 2/s*cos(a2)); + + gsl_matrix_set(cartan, 1, 0, 2/s*cos(a3)); + gsl_matrix_set(cartan, 1, 1, -2); + gsl_matrix_set(cartan, 1, 2, 2*s*cos(a1)); + + gsl_matrix_set(cartan, 2, 0, 2*s*cos(a2)); + gsl_matrix_set(cartan, 2, 1, 2/s*cos(a1)); + gsl_matrix_set(cartan, 2, 2, -2); +} + +void initialize_triangle_generators(gsl_matrix **gen, gsl_matrix *cartan) +{ + LOOP(i) { + gsl_matrix_set_identity(gen[i]); + LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j); + } +} + +typedef struct { double x[3]; } point; + +point apply(gsl_matrix *g, point x) +{ + point out; + LOOP(i) out.x[i] = 0.0; + LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(g, i, j) * x.x[j]; + return out; +} + +point incidence(point a, point b) +{ + point c; + LOOP(i) c.x[i] = a.x[(i+1)%3]*b.x[(i+2)%3] - a.x[(i+2)%3]*b.x[(i+1)%3]; + return c; +} + +double applyForm(point x, point y, gsl_matrix *form) +{ + double out = 0.0; + LOOP(i) LOOP(j) out += x.x[i] * gsl_matrix_get(form, i, j) * y.x[j]; + return out; +} + +void boundary(point x, point y, point *a, point *b, gsl_matrix *form) +{ + double formX = applyForm(x, x, form); + double formY = applyForm(y, y, form); + double formXY = applyForm(x, y, form); + double t1 = (- formXY + sqrt(formXY * formXY - formX * formY)) / formX; + double t2 = (- formXY - sqrt(formXY * formXY - formX * formY)) / formX; + LOOP(i) a->x[i] = t1 * x.x[i] + y.x[i]; + LOOP(i) b->x[i] = t2 * x.x[i] + y.x[i]; +} + +point column(gsl_matrix *m, int j) +{ + point out; + LOOP(i) out.x[i] = gsl_matrix_get(m, i, j); + return out; +} + +point row(gsl_matrix *m, int i) +{ + point out; + LOOP(j) out.x[j] = gsl_matrix_get(m, i, j); + return out; +} + +double coord(point p, int k, gsl_matrix *frame) +{ + double tmp[3] = {0, 0, 0}; + double factor = 1.0; + + LOOP(i) LOOP(j) tmp[i] += gsl_matrix_get(frame, i, j)*p.x[j]; + +#ifdef POINCARE + double norm2 = (tmp[0]*tmp[0] + tmp[1]*tmp[1]) / (tmp[2]*tmp[2]); + if(fabs(norm2) < 0.5) // just to avoid dividing by 0 + factor = 1.0 / (1.0 + sqrt(1.0 - norm2)); + else if(fabs(norm2) < 1.0) + factor = (1.0 - sqrt(1.0 - norm2)) / norm2; + else + factor = 1.0; +#endif + + return factor*tmp[k]/tmp[2]; +} + +void print_tex_header() { + char *header = + "\\documentclass{standalone}\n" + "\\usepackage[utf8]{inputenc}\n" + "\\usepackage{tikz}\n" + "\\begin{document}\n" + "\\begin{tikzpicture}[scale=5]\n"; + + printf("%s", header); +} + +void print_tex_footer() { + char *footer = + "\\draw[thick] (0,0) circle (1.0);\n" + "\\end{tikzpicture}\n" + "\\end{document}\n"; + + printf("%s", footer); +} + +char *poincare_arc(double x1, double x2, double y1, double y2, char *buffer) +{ + double t = (1 - x1*y1 - x2*y2)/(x1*y2 - x2*y1)/2; + double m1 = x1/2 + y1/2 + t*(y2 - x2); // center of circle + double m2 = x2/2 + y2/2 + t*(x1 - y1); + + double r = sqrt((x1-m1)*(x1-m1) + (x2-m2)*(x2-m2)); + double phix = atan2(x2-m2,x1-m1); + double phiy = atan2(y2-m2,y1-m1); + + if(phix - phiy > M_PI) + phiy += 2*M_PI; + else if(phiy - phix > M_PI) + phix += 2*M_PI; + + sprintf(buffer, "%f:%f:%f", phix/M_PI*180, phiy/M_PI*180, r); + return buffer; +} + +void draw_triangle(point *p, gsl_matrix *frame, const char *arguments) +{ +#ifdef POINCARE + char buffer1[100], buffer2[100], buffer3[100]; + + double x1 = coord(p[0], 0, frame); + double y1 = coord(p[0], 1, frame); + double x2 = coord(p[1], 0, frame); + double y2 = coord(p[1], 1, frame); + double x3 = coord(p[2], 0, frame); + double y3 = coord(p[2], 1, frame); + + printf("\\draw[%s] (%f, %f) arc (%s) arc (%s) arc (%s);\n", arguments, x1, y1, + poincare_arc(x1, y1, x2, y2, buffer1), + poincare_arc(x2, y2, x3, y3, buffer2), + poincare_arc(x3, y3, x1, y1, buffer3)); + +#else + printf("\\draw[%s] (%f, %f) -- (%f, %f) -- (%f, %f) -- cycle;\n", arguments, + coord(p[0], 0, frame), coord(p[0], 1, frame), + coord(p[1], 0, frame), coord(p[1], 1, frame), + coord(p[2], 0, frame), coord(p[2], 1, frame)); +#endif +} + +void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments) +{ + char buffer[100]; + + double x1 = coord(p1, 0, frame); + double y1 = coord(p1, 1, frame); + double x2 = coord(p2, 0, frame); + double y2 = coord(p2, 1, frame); + +#ifdef POINCARE + printf("\\draw[%s] (%f, %f) arc (%s);\n", arguments, x1, y1, + poincare_arc(x1, y1, x2, y2, buffer)); +#else + printf("\\draw[%s] (%f, %f) -- (%f, %f);\n", arguments, x1, y1, x2, y2); +#endif +} + +int main() +{ + groupelement_t *group; + gsl_matrix **matrices; + gsl_matrix *cartan; + gsl_matrix *gen[3]; + gsl_matrix *coxeter[3]; + gsl_matrix *coxeter_eigenvectors[3]; + gsl_matrix *frame; + workspace_t *ws; + int elements = 5000; + int p = 5, q = 5, r = 5; + + group = malloc(elements*sizeof(groupelement_t)); + matrices = malloc(elements*sizeof(gsl_matrix*)); + for(int i = 0; i < elements; i++) + matrices[i] = gsl_matrix_alloc(3, 3); + 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); + ws = workspace_alloc(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 + + 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]); + LOOP(i) eigenvectors(coxeter[i], coxeter_eigenvectors[i], ws); + + /* + 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 reflection_lines[3]; + point triangle_points[3]; + point transformed[3]; + point transformed2[3]; + + LOOP(i) coxeter_attracting[i] = column(coxeter_eigenvectors[i], 0); + LOOP(i) coxeter_repelling[i] = column(coxeter_eigenvectors[i], 2); + LOOP(i) reflection_lines[i] = row(cartan, i); + LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]); + + print_tex_header(); + +// int indices[10] = {0, 1, 6, 10, 30, 46, 124, 185, 484, 717}; +// int indices[10] = {0, 1, 4, 10, 22}; + + for(int k = 0; k < elements; k++) { + if(group[k].length % 2) + continue; + + LOOP(i) transformed[i] = apply(matrices[k], triangle_points[i]); + draw_triangle(transformed, frame, "black,fill=black!10,line width=0pt"); + } + + for(int k = 0; k < elements; k++) { + if(group[k].length % 2) + continue; + + LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]); + LOOP(i) transformed2[i] = apply(matrices[k], coxeter_repelling[i]); + LOOP(i) draw_line(transformed[i], transformed2[i], frame, "red"); + } + + print_tex_footer(); +} From 9019e4865cec4945076f6b49a161ffcdca16eb4e Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Fri, 15 Jun 2018 19:33:59 +0200 Subject: [PATCH 4/4] played around with boxes --- hyperbolic.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hyperbolic.c b/hyperbolic.c index 53c2450..b7430a3 100644 --- a/hyperbolic.c +++ b/hyperbolic.c @@ -192,7 +192,7 @@ int main() gsl_matrix *frame; workspace_t *ws; int elements = 5000; - int p = 5, q = 5, r = 5; + int p = 3, q = 5, r = 7; group = malloc(elements*sizeof(groupelement_t)); matrices = malloc(elements*sizeof(gsl_matrix*));