allow rotations

.gitignore

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

Makefile

@@ -7,61 +7,19 @@ SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 
 OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
 
-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
+all: hyperbolic
 
 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 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
+	rm -f hyperbolic triangle.o linalg.o hyperbolic.o

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("<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];
@@ -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, *special_rotation;
 	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>\n", argv[0]);
+		fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements> <word1> <word2> ...\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,57 @@ 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);
-	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);
+	nspecial_hyp = nspecial_rot = 0;
+	for(int i = 0; i < nspecial; i++) {
+		LOOP(j) {
+			int nreal;
 
-	/*
-	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);
+			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++;
+			}
+		}
 	}
 
-	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();
+	fprintf(stderr, "%d special elements, %d rotations, %d hyperbolic\n", nspecial, nspecial_rot, nspecial_hyp);
 
+	/*
 	// 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 +333,69 @@ 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;");
 	}
 
+	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;
 
-		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_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);
+		}
 
-//		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 < 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(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);
 }

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,9 +165,11 @@ 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)
+int 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);
@@ -176,20 +178,16 @@ void 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 = 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);
+	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++;
+		}
 	}
 
-	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)));
-
+	return real;
 }
 
 void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)

linalg.h

@@ -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);
-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);
+int 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);
 

maketiling.sh

@@ -0,0 +1,12 @@
+#!/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