florian/hyperbolic_tilings
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compare commits

base: florian:c62044d6374c39107748b7bbc63acf28eb4ab41d
Branches Tags
florian:main
..
compare: florian:bea92930cac0baf4ab94c039ecff93d5d160f9cb
Branches Tags
florian:main

No commits in common. "c62044d6374c39107748b7bbc63acf28eb4ab41d" and "bea92930cac0baf4ab94c039ecff93d5d160f9cb" have entirely different histories.
c62044d637 ... bea92930ca

10 changed files with 182 additions and 137 deletions
Show Stats Expand all files Collapse all files

3
.gitignore vendored
View File

@ -1,3 +0,0 @@
*.o
hyperbolic
output/

46
Makefile
View File

@ -7,19 +7,61 @@ SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
all: hyperbolic
all: discreteness singular_values nondiscrete traces element_path limit_set hyperbolic ellipticity
singular_values: singular_values.o triangle.o linalg.o
gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas
nondiscrete: nondiscrete.o triangle.o linalg.o
gcc $(OPTIONS) -o nondiscrete triangle.o linalg.o nondiscrete.o -lm -lgsl -lcblas
traces: traces.o triangle.o linalg.o
gcc $(OPTIONS) -o traces triangle.o linalg.o traces.o -lm -lgsl -lcblas
discreteness: discreteness.o triangle.o linalg.o
gcc $(OPTIONS) -o discreteness triangle.o linalg.o discreteness.o -lm -lgsl -lcblas
element_path: element_path.o linalg.o
gcc $(OPTIONS) -o element_path element_path.o linalg.o -lm -lgsl -lcblas
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
ellipticity: ellipticity.o triangle.o linalg.o
gcc $(OPTIONS) -o ellipticity ellipticity.o triangle.o linalg.o -lm -lgsl -lcblas
singular_values.o: singular_values.c $(HEADERS)
gcc $(OPTIONS) -c singular_values.c
nondiscrete.o: nondiscrete.c $(HEADERS)
gcc $(OPTIONS) -c nondiscrete.c
traces.o: traces.c $(HEADERS)
gcc $(OPTIONS) -c traces.c
discreteness.o: discreteness.c $(HEADERS)
gcc $(OPTIONS) -c discreteness.c
element_path.o: element_path.c $(HEADERS)
gcc $(OPTIONS) -c element_path.c
linalg.o: linalg.c $(HEADERS)
gcc $(OPTIONS) -c linalg.c
triangle.o: triangle.c $(HEADERS)
gcc $(OPTIONS) -c triangle.c
limit_set.o: limit_set.c $(HEADERS)
gcc $(OPTIONS) -c limit_set.c
hyperbolic.o: hyperbolic.c $(HEADERS)
gcc $(OPTIONS) -c hyperbolic.c
ellipticity.o: ellipticity.c $(HEADERS)
gcc $(OPTIONS) -c ellipticity.c
clean:
rm -f hyperbolic triangle.o linalg.o hyperbolic.o
rm -f singular_values nondiscrete traces discreteness element_path limit_set hyperbolic ellipticity triangle.o linalg.o singular_values.o nondiscrete.o traces.o discreteness.o element_path.o limit_set.o hyperbolic.o ellipticity.o

226
hyperbolic.c
View File

@ -198,16 +198,6 @@ 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("<circle cx=\"%f\" cy=\"%f\" r=\"6\" style=\"%s\"/>\n", CONV(x), CONV(y), arguments);
#endif
}
void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments)
{
char buffer[100];
@ -227,16 +217,11 @@ 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, int inverse)
void compute_word(workspace_t *ws, gsl_matrix *result, gsl_matrix **gen, const char *word, int modifier)
{
gsl_matrix_set_identity(result);
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);
}
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[])
@ -245,22 +230,20 @@ int main(int argc, const char *argv[])
gsl_matrix **matrices;
gsl_matrix *cartan;
gsl_matrix *gen[3];
gsl_matrix **special;
gsl_matrix **special_eigenvectors;
point *special_attracting, *special_repelling, *special_rotation;
gsl_matrix *coxeter[3];
gsl_matrix *coxeter2[3];
gsl_matrix *coxeter_eigenvectors[3];
gsl_matrix *coxeter_eigenvectors2[3];
gsl_matrix *frame;
workspace_t *ws;
int p,q,r;
int elements, nspecial, nspecial_hyp, nspecial_rot;
if(argc < 5) {
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements> <word1> <word2> ...\n", argv[0]);
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements>\n", argv[0]);
exit(1);
}
nspecial = argc - 5;
elements = atoi(argv[4]);
p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]);
int elements = atoi(argv[4]);
int p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]);
group = malloc(elements*sizeof(groupelement_t));
matrices = malloc(elements*sizeof(gsl_matrix*));
@ -269,57 +252,79 @@ 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));
special_rotation = 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);
diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form
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]);
nspecial_hyp = nspecial_rot = 0;
for(int i = 0; i < nspecial; i++) {
LOOP(j) {
int nreal;
compute_word(ws, special[3*i+j], gen, argv[i+5], j, 0);
nreal = eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws);
if(nreal == 3) {
special_attracting[nspecial_hyp] = 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[nspecial_hyp] = column(special_eigenvectors[3*i+j], 0);
nspecial_hyp++;
} else {
special_rotation[nspecial_rot] = column(special_eigenvectors[3*i+j], 0);
nspecial_rot++;
}
}
}
fprintf(stderr, "%d special elements, %d rotations, %d hyperbolic\n", nspecial, nspecial_rot, nspecial_hyp);
// 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 < 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));
@ -333,69 +338,80 @@ 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);
*/
gsl_matrix_free(frame_correction);
// the actual drawing
point transformed[3];
point reflection_lines[3];
point triangle_points[3];
print_svg_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) 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;");
}
char stylestring[100];
char colors[3][20] = {"red", "blue", "green"};
// draw special elements
for(int k = 0; k < elements; k++) {
// if(group[k].length % 2)
// continue;
for(int i = 0; i < nspecial_hyp; i+=3) {
// draw_dot(apply(matrices[k], special_repelling[i]), frame, "fill:red;stroke-width:1;");
// draw_dot(apply(matrices[k], special_attracting[i]), frame, "fill:blue;stroke-width:1;");
snprintf(stylestring, sizeof(stylestring), "fill:none;stroke:%s;stroke-width:1;", colors[i%3]);
draw_line(apply(matrices[k], special_repelling[i]),
apply(matrices[k], special_attracting[i]),
frame, stylestring);
}
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;");
for(int i = 0; i < nspecial_rot; i+=3) {
snprintf(stylestring, sizeof(stylestring), "fill:%s;stroke:%s;stroke-width:1;", colors[i%3], colors[i%3]);
// fprintf(stderr, "%f %f\n", special_rotation[i].x[0], special_rotation[i].x[1]);
draw_dot(apply(matrices[k], special_rotation[i]), frame, stylestring);
}
// 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");
}
/*
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);
free(special_rotation);
}

BIN
limit_curve.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 15 KiB

BIN
limit_curve2.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 20 KiB

BIN
limit_curve3.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 198 KiB

BIN
limit_curve4.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 139 KiB

28
linalg.c
View File

@ -145,11 +145,11 @@ double determinant(gsl_matrix *g, workspace_t *ws)
return gsl_linalg_LU_det(ws->tmp, s);
}
int jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws)
void 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_ABS_DESC);
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_VAL_DESC);
int real = 1;
for(int i = 0; i < ws->n; i++)
@ -165,11 +165,9 @@ int 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;
}
int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, 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);
@ -178,16 +176,20 @@ int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
int real = 0;
for(int j = 0; j < ws->n; j++) {
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0) {
for(int i = 0; i < ws->n; i++)
gsl_matrix_set(evec_real, i, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
real++;
}
int real = 1;
for(int i = 0; i < ws->n; i++)
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
real = 0;
if(!real) {
printf("We have non-real eigenvalues!\n");
exit(1);
}
return real;
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)

4
linalg.h
View File

@ -38,10 +38,10 @@ 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);
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
double trace(gsl_matrix *g);
double determinant(gsl_matrix *g, workspace_t *ws);
int eigenvectors(gsl_matrix *g, gsl_matrix *evec, 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);

12
maketiling.sh
View File

@ -1,12 +0,0 @@
#!/bin/bash
filename="output/tiling_$1$2$3"
for i in $(seq 5 $#); do filename+="_${!i}"; done
filename_svg="${filename}.svg"
filename_pdf="${filename}.pdf"
if [ $# -lt 4 ]; then
./hyperbolic
else
./hyperbolic $* > "$filename_svg" && rsvg-convert --format=pdf "$filename_svg" > "$filename_pdf"
fi