add documentation

.gitignore

@@ -0,0 +1,4 @@
+*.o
+hyperbolic
+output/
+README.html

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

README.md

@@ -0,0 +1,49 @@
+# hyperbolic tilings #
+
+A program to draw regular tilings of the hyperbolic plane by triangles as SVG files.
+
+## Setup ##
+
+GCC and GNU make are needed. To compile, just run `make`.
+
+## Usage ##
+
+The executable is called `hyperbolic` and is invoked as follows:
+
+    $ ./hyperbolic
+    Usage: ./hyperbolic <p> <q> <r> <n_elements> <word1> <word2> ...
+
+The arguments `p,q,r` specify the orders of the rotations `bc`, `ca` and `ab`, where `a,b,c` are the reflections along the three sides of the triangle. The parameter `n_elements` limits how many translates of the triangle are being drawn (1000 is a good default value, although for some values of `p,q,r` more are needed).
+
+The program can not only draw the tiling by triangles, but also the axes of specific hyperbolic elements in the triangle group or the fixed points of elliptic elements. To use this feature, just specify the group elements as additional arguments, written as words in the three reflections `a,b,c`. This will draw its axis/fixed point and that of all of its conjugates.
+
+## Examples ##
+
+    ./hyperbolic 5 5 5 5000 abc
+
+![(5,5,5) tiling with axes of abc and its conjugates](https://florianstecker.net/git_files/hyperbolic_tilings/tiling_555_abc.png)
+
+    ./hyperbolic 2 3 7 20000
+
+![(2,3,7) tiling](https://florianstecker.net/git_files/hyperbolic_tilings/tiling_237.png)
+
+    ./hyperbolic 2 3 7 20000 ab
+
+![(2,3,7) tiling](https://florianstecker.net/git_files/hyperbolic_tilings/tiling_237_ab.png)
+
+
+## Conversion to PDF or PNG ##
+
+If a different image format is required, the SVG files can easily be converted. Using [rsvg-convert] gives good results.
+
+To convert to PNG (used for the images above):
+
+    rsvg-convert --format=png -w 600 -h 600 tiling_555_abc.svg -o tiling_555_abc.png
+
+To convert to PDF:
+
+    rsvg-convert --format=pdf -w tiling_555_abc.svg -o tiling_555_abc.pdf
+
+`maketiling.sh` is a convenience script which takes the same arguments as `hyperbolic`, but saves the SVG file in the folder `output/` (which has to be present) with a descriptive name and also converts it to PDF.
+
+[rsvg-convert]: https://gitlab.gnome.org/GNOME/librsvg

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,56 @@ 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++) {
+	    int j = 0;
+	    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 +332,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