some clarification

.gitignore

@@ -1,3 +0,0 @@
-*.o
-limit_set
-*.pdf

Makefile

@@ -1,4 +1,4 @@
-HEADERS=triangle.h linalg.h queue.h initcairo.h main.h exp_equation.h
+HEADERS=triangle.h linalg.h queue.h initcairo.h main.h
 
 SPECIAL_OPTIONS=-O0 -g -D_DEBUG
 #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
@@ -11,8 +11,8 @@ OPTIONS=$(GENERAL_OPTIONS) $(CAIRO_OPTIONS) $(SPECIAL_OPTIONS)
 
 all: limit_set
 
-limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o
+limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o
-	gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o -lm -lgsl -lcblas -lcairo -lX11
+	gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o -lm -lgsl -lcblas -lcairo -lX11
 
 linalg.o: linalg.c $(HEADERS)
 	gcc $(OPTIONS) -c linalg.c
@@ -32,8 +32,5 @@ draw.o: draw.c $(HEADERS)
 main.o: main.c $(HEADERS)
 	gcc $(OPTIONS) -c main.c
 
-exp_equation.o: exp_equation.c $(HEADERS)
-	gcc $(OPTIONS) -c exp_equation.c
-
 clean:
-	rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o exp_equation.o
+	rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o

README.md

@@ -0,0 +1,50 @@
+# triangle_group_limit_set - visualizer for triangle group limit sets in the projective plane
+
+This program visualizes the fractal limit set of triangle groups in SL(3,R) in the projective plane. For the mathematical details, see [our paper].
+
+## Installation & Prerequisites ##
+
+This program is written C and only runs on Linux with the X window system (because that's what I happen to be using). It uses [Cairo] for drawing and X11 for input. We also require the [Gnu Scientific Library (GSL)] including CBLAS.
+
+To build, just make sure GCC as well as these libraries are installed, and run
+
+    make
+
+Then we can run it for example like this:
+
+    ./limit_set 5 5 5 1 1 1 1.0
+
+The arguments are 6 integers (p1, p2, p3, q1, q2, q3) and a floating point numbers t. Together they discribe the triangle group under consideration (t can be changed with the arrow keys).
+
+## Key bindings
+
+Drag the mouse to move the image, drag with Shift pressed to rotate.
+
+| Key      | function                                                                                            |
+|----------|-----------------------------------------------------------------------------------------------------|
+| PageUp   | increase t by 2%                                                                                    |
+| PageDown | decrease t by 2%                                                                                    |
+| Right    | increase t by 0.2%                                                                                  |
+| Left     | decreaes t by 0.2%                                                                                  |
+| Up       | increase t by 0.002%                                                                                |
+| Down     | decrease t by 0.002%                                                                                |
+| Space    | (reset t to original value)                                                                         |
+| R        | cycle through affine charts                                                                         |
+| l        | show limit curve                                                                                    |
+| L        | show limit curve as line or as points                                                               |
+| f        | generate limit curve using attracting / repelling fixed points of conjugates of the Coxeter element |
+| d        | show dual limit curve                                                                               |
+| r        | show fixed points and lines of generating reflections                                               |
+| a        | show fixed points and lines of Coxeter elements                                                     |
+| c        | show the orbit of an arbitrary point (point can be changed with Shift+Click)                        |
+| b        | show certain conics touching the limit curve                                                        |
+| B        | show the "boxes" used to approximate the limit curve                                                |
+| p        | save a screenshot of the current image (in PDF format)                                              |
+| t        | toggle info text                                                                                    |
+| i        | (print some info to terminal?)                                                                      |
+| x        | (show rotated reflectors?)                                                                          |
+| M        | (make a movie)                                                                                      |
+
+[Cairo]: https://cairographics.org/
+[Gnu Scientific Library (GSL)]: https://www.gnu.org/software/gsl/
+[our paper]: https://arxiv.org/abs/2106.11349

draw.c

@@ -1,5 +1,4 @@
 #include "main.h"
-#include "exp_equation.h"
 
 #define FMOD(x,y) (fmod(x,y) < 0 ? fmod(x,y) + y : fmod(x,y))
 #define ANGLE_DIFF(x,y) (FMOD((x)-(y), 2*M_PI))
@@ -23,50 +22,6 @@ vector_t cross(vector_t a, vector_t b)
 	return result;
 }
 
-double scal(vector_t a, vector_t b)
-{
-	double result = 0.0;
-	LOOP(i) result += a.x[i]*b.x[i];
-	return result;
-}
-
-// degenerate conic whose zero set is the union of two lines
-// Q = a^T b + b^T a
-void degenerate_conic(vector_t a, vector_t b, gsl_matrix *conic)
-{
-	LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, a.x[i]*b.x[j] + a.x[j]*b.x[i]);
-}
-
-double conic_value(gsl_matrix *conic, vector_t v)
-{
-	double result = 0.0;
-	LOOP(i) LOOP(j) result += gsl_matrix_get(conic, i, j) * v.x[i] * v.x[j];
-	return result;
-}
-
-// finds conic in the pencil through (l1,l2) and (L1,L2) going through p
-void conic_from_lines(workspace_t *ws, vector_t l1, vector_t l2, vector_t L1, vector_t L2, vector_t p, gsl_matrix *conic)
-{
-	gsl_matrix *degenerate1 = gsl_matrix_alloc(3, 3);
-	gsl_matrix *degenerate2 = gsl_matrix_alloc(3, 3);
-	gsl_matrix *tmp = gsl_matrix_alloc(3, 3);
-
-	degenerate_conic(l1, l2, degenerate1);
-	degenerate_conic(L1, L2, degenerate2);
-
-	double t = - conic_value(degenerate1, p) / conic_value(degenerate2, p);
-	LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j,
-	                               gsl_matrix_get(degenerate1, i, j) +
-	                               gsl_matrix_get(degenerate2, i, j)*t);
-	int positives = diagonalize_symmetric_form(conic, tmp, ws);
-	if(positives == 2)
-		LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1;
-
-	gsl_matrix_free(degenerate1);
-	gsl_matrix_free(degenerate2);
-	gsl_matrix_free(tmp);
-}
-
 vector_t apply(gsl_matrix *m, vector_t x)
 {
 	vector_t out;
@@ -87,133 +42,6 @@ vector_t apply_transpose(gsl_matrix *m, vector_t x)
 	return out;
 }
 
-vector_t apply_pseudoinverse_transpose(gsl_matrix *m, vector_t in)
-{
-	vector_t out;
-	double cofactor;
-
-	LOOP(i) out.x[i] = 0.0;
-	LOOP(i) LOOP(j) {
-		cofactor = gsl_matrix_get(m, (i+1)%3, (j+1)%3) * gsl_matrix_get(m, (i+2)%3, (j+2)%3)
-			- gsl_matrix_get(m, (i+1)%3, (j+2)%3) * gsl_matrix_get(m, (i+2)%3, (j+1)%3);
-		out.x[i] += cofactor * in.x[j];
-	}
-
-	return out;
-}
-
-vector_t apply_pseudoinverse(gsl_matrix *m, vector_t in)
-{
-	vector_t out;
-	double cofactor;
-
-	LOOP(i) out.x[i] = 0.0;
-	LOOP(i) LOOP(j) {
-		cofactor = gsl_matrix_get(m, (j+1)%3, (i+1)%3) * gsl_matrix_get(m, (j+2)%3, (i+2)%3)
-			- gsl_matrix_get(m, (j+1)%3, (i+2)%3) * gsl_matrix_get(m, (j+2)%3, (i+1)%3);
-		out.x[i] += cofactor * in.x[j];
-	}
-
-	return out;
-}
-
-int intersect_line_and_conic(DrawingContext *ctx, vector_t line, gsl_matrix *conic, vector_t *intersection1, vector_t *intersection2)
-{
-	gsl_matrix *frame = getTempMatrix(ctx->ws);
-
-	int pos = diagonalize_symmetric_form(conic, frame, ctx->ws);
-
-	vector_t a = apply_pseudoinverse_transpose(frame, line);
-
-	double disc = a.x[0]*a.x[0] + a.x[1]*a.x[1] - a.x[2]*a.x[2];
-	if(disc < 0)
-	{
-		return 0;
-	}
-
-	vector_t intersection1_in_frame;
-	intersection1_in_frame.x[0] = -a.x[0]*a.x[2] + a.x[1]*sqrt(disc);
-	intersection1_in_frame.x[1] = -a.x[1]*a.x[2] - a.x[0]*sqrt(disc);
-	intersection1_in_frame.x[2] = a.x[0]*a.x[0] + a.x[1]*a.x[1];
-	*intersection1 = apply_pseudoinverse(frame, intersection1_in_frame);
-
-	vector_t intersection2_in_frame;
-	intersection2_in_frame.x[0] = -a.x[0]*a.x[2] - a.x[1]*sqrt(disc);
-	intersection2_in_frame.x[1] = -a.x[1]*a.x[2] + a.x[0]*sqrt(disc);
-	intersection2_in_frame.x[2] = a.x[0]*a.x[0] + a.x[1]*a.x[1];
-	*intersection2 = apply_pseudoinverse(frame, intersection2_in_frame);
-
-	releaseTempMatrices(ctx->ws, 1);
-}
-
-// intersect the line given by the covector "line" with the orbit of "orbit_point"
-// by the one--parameter subgroup of SL(3,R) which contains the element "loxodromic"
-// in an eigenbasis of "loxodromic", this corresponds
-int intersect_line_and_loxodromic_orbit(DrawingContext *ctx, vector_t line, gsl_matrix *frame, double *logeigenvalues, vector_t start, vector_t *out)
-{
-	vector_t line_in_frame = apply_transpose(frame, line);
-	vector_t start_in_frame = apply_pseudoinverse(frame, start);
-
-	vector_t a, x;
-	LOOP(i) a.x[i] = logeigenvalues[i];
-	LOOP(i) x.x[i] = line_in_frame.x[i]*start_in_frame.x[i];
-
-	double t[2];
-	vector_t v[2];
-	int n1, n2;
-
-	n1 = solve_linear_exp(a, x, t);
-	for(int i = 0; i < n1; i++) {
-		LOOP(j) v[i].x[j] = exp(a.x[j]*t[i]) * start_in_frame.x[j];
-		out[i] = apply(frame, v[i]);
-	}
-
-	x.x[1] *= -1;
-	n2 = solve_linear_exp(a, x, t);
-	for(int i = 0; i < n2; i++) {
-		LOOP(j) v[i].x[j] = exp(a.x[j]*t[i]) * start_in_frame.x[j];
-		v[i].x[1] *= -1;
-		out[i+n1] = apply(frame, v[i]);
-	}
-
-	if(n1+n2 > 2) {
-		fprintf(stderr, "more than 2 solutions in intersect_line_and_loxodromic_orbit()!\n");
-		exit(1);
-	}
-
-	return n1+n2;
-}
-
-// should be three collinear vectors!
-double halfCR(vector_t x, vector_t y, vector_t z)
-{
-	double xy = scal(x,y);
-	double yz = scal(y,z);
-	double zx = scal(z,x);
-	double yy = scal(y,y);
-	double xx = scal(x,x);
-
-	return fabs((yz*xy-zx*yy)/(yz*xx-zx*xy));
-}
-
-void getMinMaxHalfCrossRatio(double *limit_curve, int n, vector_t v1, vector_t v2, int max, int *min_index, double* min_cr)
-{
-	vector_t fp[3];
-	vector_t tangent;
-	double cr;
-	*min_cr = max ? 0 : INFINITY;
-	for(int i = 0; i < n; i++) {
-		LOOP(j) LOOP(k) fp[j].x[k] = limit_curve[12*i+3+3*j+k];
-
-		tangent = cross(fp[0],fp[1]);
-		cr = fabs(scal(v2, tangent) / scal(v1, tangent));
-		if(!max && cr < *min_cr || max && cr > *min_cr) {
-			*min_cr = cr;
-			*min_index = i;
-		}
-	}
-}
-
 int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out)
 {
 	gsl_matrix *tmp = getTempMatrix(ctx->ws);
@@ -268,20 +96,10 @@ void drawPoint(DrawingContext *ctx, point_t p)
 	cairo_t *C = ctx->cairo;
 
 	cairo_save(C);
-	cairo_arc(C, p.x, p.y, 3.0/ctx->dim->scalefactor, 0, 2*M_PI);
+	cairo_arc(C, p.x, p.y, 4.0/ctx->dim->scalefactor, 0, 2*M_PI);
 	cairo_close_path(C);
 	cairo_fill(C);
 	cairo_restore(C);
-
-	/*
-	cairo_save(C);
-	cairo_move_to(C, p.x, p.y);
-	cairo_close_path(C);
-	cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);
-	cairo_set_line_width(C, 10.0/ctx->dim->scalefactor);
-	cairo_stroke(C);
-	cairo_restore(C);
-	*/
 }
 
 void drawSegment2d(DrawingContext *ctx, point_t a, point_t b)
@@ -365,22 +183,13 @@ void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line,
 	LOOP(i) x[i] = 0.0;
 	LOOP(i) LOOP(j) {
 		cofactor = gsl_matrix_get(ctx->cob, (i+1)%3, (j+1)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+2)%3)
-			- gsl_matrix_get(ctx->cob, (i+1)%3, (j+2)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+1)%3);
+		    - gsl_matrix_get(ctx->cob, (i+1)%3, (j+2)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+1)%3);
 		x[i] += cofactor * line.x[j];
 	}
 
 	// t = parameter on segment of intersection with line, s(t) = a + (b-a)*t
 	tline = -(a_.x*x[0] + a_.y*x[1] + x[2])/((b_.x - a_.x)*x[0] + (b_.y - a_.y)*x[1]);
 
-	/*
-	printf("tline: %f\n", tline);
-
-	point_t s;
-	s.x = a_.x - (b_.x-a_.x)*tline;
-	s.y = a_.y - (b_.y-a_.y)*tline;
-	drawPoint(ctx, s);
-	*/
-
 	if((tline < 0 || tline > 1) != contains) {
 		drawSegment2d(ctx, a_, b_);
 	} else { // need to draw complementary segment
@@ -388,6 +197,7 @@ void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line,
 		m.x = ctx->dim->center_x;
 		m.y = ctx->dim->center_y;
 		r = ctx->dim->radius;
+
 		// equation is coeff2 t^2 + 2 coeff1 t + coeff0 = 0
 		coeff0 = (a_.x - m.x)*(a_.x - m.x) + (a_.y - m.y)*(a_.y - m.y) - r*r;
 		coeff1 = (a_.x - m.x)*(b_.x - a_.x) + (a_.y - m.y)*(b_.y - a_.y);
@@ -408,7 +218,6 @@ void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line,
 }
 
 // level 3: boxes and polygons
-
 void drawPolygon(DrawingContext *ctx, int segments, int sides, ...)
 {
 	va_list args;
@@ -470,8 +279,6 @@ void drawBoxLines(DrawingContext *ctx, const char *word1, const char *word2)
 	drawPolygon(ctx, 0, 4, p[0][0], i[0], p[1][0], i[1]);
 }
 
-
-
 void drawBoxStd(DrawingContext *ctx, const char *word, char base)
 {
 	char word1[100];
@@ -539,70 +346,6 @@ void drawRotationOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, vector_t sta
 	releaseTempVectors(ctx->ws, 2);
 }
 
-void drawLoxodromicOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, double *logeigenvalues, vector_t start)
-{
-	vector_t start_in_frame = apply_pseudoinverse(frame, start);
-	int iterations = 500;
-	double stepsize = 0.02;
-	vector_t x, w;
-	point_t p;
-	cairo_t *C = ctx->cairo;
-	double t;
-
-	int previous_inside = 0;
-	for(int k = 0; k <= iterations; k++) {
-		// 0 = repelling fixed point, iterations/2 = attracting fixed point
-		if(k == 0 || k == iterations) {
-			w.x[0] = 0.0; w.x[1] = 0.0; w.x[2] = 1.0;
-		} else if(k == iterations/2) {
-			w.x[0] = 1.0; w.x[1] = 0.0; w.x[2] = 0.0;
-		} else if(k < iterations/2) {
-			t = (k-(double)iterations/4.0)*stepsize;
-			LOOP(i) w.x[i] = start_in_frame.x[i] * exp(logeigenvalues[i]*t);
-			w.x[1] *= -1;
-		} else {
-			t = ((double)iterations*3.0/4.0-k)*stepsize;
-			LOOP(i) w.x[i] = start_in_frame.x[i] * exp(logeigenvalues[i]*t);
-		}
-
-		x = apply(frame, w);
-		p = vectorToPoint(ctx, x);
-
-		if(isInsideBB(ctx, p)) {
-			if(!previous_inside)
-				cairo_move_to(C, p.x, p.y);
-			else
-				cairo_line_to(C, p.x, p.y);
-			previous_inside = 1;
-		} else {
-			previous_inside = 0;
-		}
-	}
-
-	cairo_stroke(C);
-}
-
-void drawConic(DrawingContext *ctx, gsl_matrix *form)
-{
-	gsl_matrix *orthogonal_frame = getTempMatrix(ctx->ws);
-	gsl_matrix *frame = getTempMatrix(ctx->ws);
-	double eval[3];
-	vector_t start;
-
-	diagonalize_symmetric_matrix(form, orthogonal_frame, eval, ctx->ws);
-
-	LOOP(i) LOOP(j)
-		gsl_matrix_set(frame, j, i,
-		               gsl_matrix_get(orthogonal_frame, i, j)/sqrt(fabs(eval[i])));
-
-	// find any timelike vector; a simple method is adding the first and third column of frame
-	LOOP(i) start.x[i] = gsl_matrix_get(frame, i, 0) + gsl_matrix_get(frame, i, 2);
-
-	drawRotationOrbitFrame(ctx, frame, start);
-
-	releaseTempMatrices(ctx->ws, 2);
-}
-
 void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start)
 {
 	gsl_matrix *frame = getTempMatrix(ctx->ws);
@@ -815,7 +558,8 @@ void drawAttractors(DrawingContext *ctx)
 	fixedPoints(ctx, "abc", p[0]);
 	fixedPoints(ctx, "bca", p[1]);
 	fixedPoints(ctx, "cab", p[2]);
-	fixedPoints(ctx, "a cab a", p[3]);
+	//	fixedPoints(ctx, "abacacbabc", p[3]);
+	//	fixedPoints(ctx, "abcabcabcabcab", p[3]);
 	fixedPoints(ctx, "b abc b", p[4]);
 	fixedPoints(ctx, "c bca c", p[5]);
 
@@ -823,18 +567,14 @@ void drawAttractors(DrawingContext *ctx)
 
 	for(int i = 0; i < n; i++) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]);
 
-	for(int i = 0; i < n; i++) {
+	for(int i = 0; i < n; i++)  LOOP(j) {
-		for(int j = 0; j < 3; j += 1) {
+		cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
-			cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
+		drawVector(ctx, p[i][j]);
-			drawVector(ctx, p[i][j]);
-		}
 	}
 
-	for(int i = 0; i < n; i++) {
+	for(int i = 0; i < n; i++)  LOOP(j) {
-		for(int j = 0; j < 3; j += 1) {
+		cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
-			cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
+		drawCovector(ctx, l[i][j]);
-			drawCovector(ctx, l[i][j]);
-		}
 	}
 }
 
@@ -893,12 +633,6 @@ void drawCurvedBox(DrawingContext *ctx, int base, const char *conj, int style)
 	fixedPoints(ctx, word, p[10]);
 
 
-//	conjugate_word("bca b abc b acb", modifier, conj, word);
-//	fixedPoints(ctx, word, p[6]);
-//	conjugate_word("bca baca cab acab acb", modifier, conj, word);
-//	fixedPoints(ctx, word, p[7]);
-
-
 	LOOP(j) l[0][j] = cross(p[0][(3-j)%3], p[0][(4-j)%3]);
 	LOOP(j) l[1][j] = cross(p[1][(3-j)%3], p[1][(4-j)%3]);
 	LOOP(j) l[10][j] = cross(p[10][(3-j)%3], p[10][(4-j)%3]);
@@ -950,10 +684,6 @@ void drawBoxes(DrawingContext *ctx)
 {
 	gsl_matrix *rot = getTempMatrix(ctx->ws);
 	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
-	gsl_matrix *order3 = getTempMatrix(ctx->ws);
-	gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	gsl_matrix *tmp2 = getTempMatrix(ctx->ws);
-	gsl_matrix *conic = getTempMatrix(ctx->ws);
 	gsl_matrix *frame = getTempMatrix(ctx->ws);
 	cairo_t *C = ctx->cairo;
 	cairo_save(C);
@@ -973,236 +703,105 @@ void drawBoxes(DrawingContext *ctx)
 	cairo_set_line_width(C, 2.0/ctx->dim->scalefactor);
 	cairo_set_source_rgb(C, 0.6, 0.6, 0.6);
 
-	LOOP(i) LOOP(j) gsl_matrix_set(tmp, i, j, p[0][j].x[i]);
+	drawRotationOrbit(ctx, "ab", p[0][0]);
-	double evs[3];
+	drawRotationOrbit(ctx, "bc", p[0][0]);
-	wordEigenvalues(ctx, "abc", evs);
+	drawRotationOrbit(ctx, "ca", p[0][0]);
-	LOOP(i) evs[i] = log(fabs(evs[i]));
-	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][0]);
-	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][2]);
 
-	vector_t x;
+	if(ctx->mode >= 2) {
-	for(int i = 0; i < ctx->n_group_elements; i++) {
+		drawRotationOrbit(ctx, "bcabcb", p[1][0]); // bcC
-		LOOP(j) x.x[j] = ctx->limit_curve[12*i+3+j];
+		drawRotationOrbit(ctx, "abcabcba", p[0][0]); // abcC
-		x = apply_pseudoinverse(tmp, x);
-		printf("%f\n",
-		       pow(fabs(x.x[0]), evs[1]-evs[2]) *
-		       pow(fabs(x.x[1]), evs[2]-evs[0]) *
-		       pow(fabs(x.x[2]), evs[0]-evs[1]));
 	}
 
+	cairo_set_source_rgb(C, 0, 0, 0);
+	fixedPoints(ctx, "abababcbaba", p[3]);
+
+	cairo_set_source_rgb(C, 0, 0, 1);
+	fixedPoints(ctx, "bc bca cb", p[3]);
+
 	cairo_restore(C);
-	releaseTempMatrices(ctx->ws, 9);
+	releaseTempMatrices(ctx->ws, 5);
 }
 
 void drawBoxes2(DrawingContext *ctx)
 {
 	gsl_matrix *rot = getTempMatrix(ctx->ws);
 	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
-	gsl_vector *marking_globalbasis = getTempVector(ctx->ws);
-	gsl_vector *marking_drawbasis = getTempVector(ctx->ws);
-	gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	gsl_matrix *tmp2 = getTempMatrix(ctx->ws);
-	gsl_matrix *conic = getTempMatrix(ctx->ws);
 	cairo_t *C = ctx->cairo;
 	cairo_save(C);
 	initializeTriangleGenerators(gen, ctx->cartan);
 
-	vector_t p[3][3];
+	vector_t p[4][3];
 
-	cairo_set_line_width(C, 1/ctx->dim->scalefactor);
-	cairo_set_source_rgb(C, 0, 0, 0);
-
-	// find coxeter fixed points and triangle vertices
 	fixedPoints(ctx, "abc", p[0]);
 	fixedPoints(ctx, "bca", p[1]);
 	fixedPoints(ctx, "cab", p[2]);
 
-	vector_t vertex[3];
+	cairo_set_line_width(C, 2.5/ctx->dim->scalefactor);
-	LOOP(i) vertex[i] = cross(
-		cross(p[(i+1)%3][0],p[(i+1)%3][2]),
-		cross(p[(i+2)%3][0],p[(i+2)%3][2]));
 
-	vector_t a,b;
+	cairo_set_source_rgb(C, 0, 0, 0);
-	double t = ctx->distance_parameter1, s = ctx->distance_parameter2;
+	drawCurvedBox(ctx, 'A', "", 2);
 
-	vector_t phiplus, phiminus, eta;
+	if(ctx->mode >= 3 && ctx->mode != 4) {
-	double dvert, param;
+		cairo_set_source_rgb(C, 0, 0, 0);
+		drawCurvedBox(ctx, 'B', "", 2);
+		drawCurvedBox(ctx, 'C', "", 2);
+	}
 
-	// want to choose a so that d(vertex[0],a) = t*d(vertex[0],vertex[1])
+	if(ctx->mode == 4 || ctx->mode == 5) {
-	phiplus = cross(p[2][0], p[2][1]);
+		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-	phiminus = cross(p[2][2], p[2][1]);
+		drawCurvedBox(ctx, 'A', "", 2);
-	dvert = log(
+		drawCurvedBox(ctx, 'A', "bc", 2);
-		(scal(phiplus, vertex[0])*scal(phiminus, vertex[1])) /
+		drawCurvedBox(ctx, 'A', "bcbc", 2);
-		(scal(phiplus, vertex[1])*scal(phiminus, vertex[0])));
+		drawCurvedBox(ctx, 'A', "bcbcbc", 2);
-	LOOP(i) eta.x[i] =
+		drawCurvedBox(ctx, 'A', "bcbcbcbc", 2);
-		scal(phiminus, vertex[0]) / scal(phiplus, vertex[0])
+		drawCurvedBox(ctx, 'A', "bcbcbcbcbc", 2);
-		* exp(t*dvert) * phiplus.x[i]
+		drawCurvedBox(ctx, 'A', "bcbcbcbcbcbc", 2);
-		- phiminus.x[i];
+	}
-	param = 1/(1-scal(eta,vertex[0])/scal(eta,vertex[1]));
-	LOOP(i) a.x[i] = param*vertex[0].x[i] + (1-param)*vertex[1].x[i];
 
-	// want to choose b so that d(vertex[0],b) = s*d(vertex[0],vertex[2])
+	if(ctx->mode == 6) {
-	phiplus = cross(p[1][0], p[1][1]);
+		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-	phiminus = cross(p[1][2], p[1][1]);
+		drawCurvedBox(ctx, 'C', "ababab", 2);
-	dvert = log(
+		drawCurvedBox(ctx, 'C', "abababab", 2);
-		(scal(phiplus, vertex[0])*scal(phiminus, vertex[2])) /
+		drawCurvedBox(ctx, 'C', "ababababab", 2);
-		(scal(phiplus, vertex[2])*scal(phiminus, vertex[0])));
+	}
-	LOOP(i) eta.x[i] =
-		scal(phiminus, vertex[0]) / scal(phiplus, vertex[0])
-		* exp(s*dvert) * phiplus.x[i]
-		- phiminus.x[i];
-	param = 1/(1-scal(eta,vertex[0])/scal(eta,vertex[2]));
-	LOOP(i) b.x[i] = param*vertex[0].x[i] + (1-param)*vertex[2].x[i];
 
-	// draw and output stuff
+	if(ctx->mode >= 6) {
-	LOOP(i) drawVector(ctx, vertex[i]);
+		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-	drawVector(ctx, a);
+		drawCurvedBox(ctx, 'C', "ab", 2);
-	drawVector(ctx, b);
+		drawCurvedBox(ctx, 'C', "abab", 2);
-	drawCovector(ctx, cross(a,b));
+	}
-	drawCovector(ctx, cross(b,vertex[0]));
-	drawCovector(ctx, cross(vertex[0],a));
 
-	// find the intersections with the limit set and the asymmetric distances
+	if(ctx->mode >= 8) {
-	int boundary_index;
+		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-	double value;
+		drawCurvedBox(ctx, 'C', "b", 2);
-	double perimeter1 = 0.0, perimeter2 = 0.0;
+		drawCurvedBox(ctx, 'C', "bab", 2);
-	double approx_perimeter1 = 0.0, approx_perimeter2 = 0.0;
+	}
-	vector_t boundary_point;
-	cairo_set_source_rgb(C, 1, 0, 0);
 
-	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	if(ctx->mode >= 9) {
-	                        a, b, 1, &boundary_index, &value);
+	    cairo_set_source_rgb(C, 0.8, 0, 0.5);
-	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	    drawCurvedBox(ctx, 'B', "abca", 2);
-	drawVector(ctx, boundary_point);
+	    drawCurvedBox(ctx, 'B', "abcaca", 2);
-	perimeter1 += log(value);
+	    drawCurvedBox(ctx, 'B', "aba", 2);
+	    drawCurvedBox(ctx, 'B', "abaca", 2);
 
-	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	    drawCurvedBox(ctx, 'B', "ababca", 2);
-	                        b, vertex[0], 1, &boundary_index, &value);
+	    drawCurvedBox(ctx, 'B', "ababcaca", 2);
-	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	    drawCurvedBox(ctx, 'B', "ababa", 2);
-	drawVector(ctx, boundary_point);
+	    drawCurvedBox(ctx, 'B', "ababaca", 2);
-	perimeter1 += log(value);
-	approx_perimeter1 += log(value);
 
-	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	    drawCurvedBox(ctx, 'B', "bca", 2);
-	                        vertex[0], a, 1, &boundary_index, &value);
+	    drawCurvedBox(ctx, 'B', "bcaca", 2);
-	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	    drawCurvedBox(ctx, 'B', "ba", 2);
-	drawVector(ctx, boundary_point);
+	    drawCurvedBox(ctx, 'B', "baca", 2);
-	perimeter1 += log(value);
-	approx_perimeter1 += log(value);
 
-	cairo_set_source_rgb(C, 0, 0, 1);
+	    drawCurvedBox(ctx, 'B', "babca", 2);
+	    drawCurvedBox(ctx, 'B', "babcaca", 2);
+	    drawCurvedBox(ctx, 'B', "baba", 2);
+	    drawCurvedBox(ctx, 'B', "babaca", 2);
 
-	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	    cairo_restore(C);
-	                        b, a, 1, &boundary_index, &value);
+	    releaseTempMatrices(ctx->ws, 4);
-	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	}
-	drawVector(ctx, boundary_point);
-	perimeter2 += log(value);
-
-	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
-	                        vertex[0], b, 1, &boundary_index, &value);
-	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
-	drawVector(ctx, boundary_point);
-	perimeter2 += log(value);
-	approx_perimeter2 += log(value);
-
-	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
-	                        a, vertex[0], 1, &boundary_index, &value);
-	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
-	drawVector(ctx, boundary_point);
-	perimeter2 += log(value);
-	approx_perimeter2 += log(value);
-
-	cairo_set_source_rgb(C, 0, 0.6, 0.1);
-
-	LOOP(i) LOOP(j) gsl_matrix_set(tmp, i, j, p[0][j].x[i]);
-	double evs[3];
-	wordEigenvalues(ctx, "abc", evs);
-	LOOP(i) evs[i] = log(fabs(evs[i]));
-	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][0]);
-	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][2]);
-
-	vector_t intersection[2];
-	double hC0, hC1;
-
-	cairo_set_source_rgb(C, 1, 0, 0);
-	intersect_line_and_loxodromic_orbit(ctx, cross(a,b), tmp, evs, p[1][2], intersection);
-	hC0 = halfCR(a, b, intersection[0]);
-	hC1 = halfCR(a, b, intersection[1]);
-	if(hC0 > hC1)
-		drawVector(ctx, intersection[0]);
-	else
-		drawVector(ctx, intersection[1]);
-	approx_perimeter1 += log(hC0 > hC1 ? hC0 : hC1);
-
-	cairo_set_source_rgb(C, 0, 0, 1);
-	intersect_line_and_loxodromic_orbit(ctx, cross(a,b), tmp, evs, p[1][0], intersection);
-	hC0 = halfCR(b, a, intersection[0]);
-	hC1 = halfCR(b, a, intersection[1]);
-	if(hC0 > hC1)
-		drawVector(ctx, intersection[0]);
-	else
-	drawVector(ctx, intersection[1]);
-	approx_perimeter2 += log(hC0 > hC1 ? hC0 : hC1);
-
-	/*
-	vector_t conic_intersection[2];
-	double hC0, hC1;
-
-	// find the outer conic
-	cairo_set_source_rgb(C, 0, 0.7, 0);
-
-	conic_from_lines(ctx->ws,
-		cross(p[0][0], p[0][1]), cross(p[1][0], p[1][1]),
-		cross(p[0][0], p[1][0]), cross(p[0][0], p[1][0]),
-		p[2][0], conic);
-	drawConic(ctx, conic);
-
-	conic_intersection[2];
-	intersect_line_and_conic(ctx, cross(a, b), conic,
-	                         &conic_intersection[0], &conic_intersection[1]);
-
-	hC0 = halfCR(b, a, conic_intersection[0]);
-	hC1 = halfCR(b, a, conic_intersection[1]);
-	if(hC0 > hC1)
-		drawVector(ctx, conic_intersection[0]);
-	else
-		drawVector(ctx, conic_intersection[1]);
-	approx_perimeter2 += log(hC0 > hC1 ? hC0 : hC1);
-
-	// find the inner conic
-	cairo_set_source_rgb(C, 0, 0.5, 0.5);
-
-	conic_from_lines(ctx->ws,
-		cross(p[0][2], p[0][1]), cross(p[1][2], p[1][1]),
-		cross(p[0][2], p[1][2]), cross(p[0][2], p[1][2]),
-		p[2][2], conic);
-	drawConic(ctx, conic);
-
-	conic_intersection[2];
-	intersect_line_and_conic(ctx, cross(a, b), conic,
-	                         &conic_intersection[0], &conic_intersection[1]);
-
-	hC0 = halfCR(a, b, conic_intersection[0]);
-	hC1 = halfCR(a, b, conic_intersection[1]);
-	if(hC0 > hC1)
-		drawVector(ctx, conic_intersection[0]);
-	else
-		drawVector(ctx, conic_intersection[1]);
-	approx_perimeter1 += log(hC0 > hC1 ? hC0 : hC1);
-
-	// output
-	snprintf(ctx->extra_text, 1000, "perimeter1 = %f (%f), perimeter2 = %f (%f)",
-	         perimeter1, approx_perimeter1,
-	         perimeter2, approx_perimeter2);
-	*/
-
-	printf("%f %f %f %f %f %f\n",
-	       t, s, perimeter1, perimeter2, approx_perimeter1, approx_perimeter2);
-
-
-	cairo_restore(C);
-	releaseTempMatrices(ctx->ws, 7);
-	releaseTempVectors(ctx->ws, 2);
 }
 
 void drawRotatedReflectors(DrawingContext *ctx)
@@ -1340,8 +939,8 @@ void drawLimitCurve(DrawingContext *ctx)
 
 		for(int i = 0; i < ctx->limit_curve_count; i++) {
 			point_t p;
-			p.x = ctx->limit_curve[12*i];
+			p.x = ctx->limit_curve[3*i];
-			p.y = ctx->limit_curve[12*i+1];
+			p.y = ctx->limit_curve[3*i+1];
 
 			if(isInsideBB(ctx, p)) {
 				if(!previous_inside)
@@ -1358,8 +957,8 @@ void drawLimitCurve(DrawingContext *ctx)
 	} else {
 		for(int i = 0; i < ctx->limit_curve_count; i++) {
 			point_t p;
-			p.x = ctx->limit_curve[12*i];
+			p.x = ctx->limit_curve[3*i];
-			p.y = ctx->limit_curve[12*i+1];
+			p.y = ctx->limit_curve[3*i+1];
 
 			if(isInsideBB(ctx, p)) {
 				cairo_arc(C, p.x, p.y, 2.0/ctx->dim->scalefactor, 0, 2*M_PI);
@@ -1465,8 +1064,8 @@ void drawText(DrawingContext *ctx)
 {
 	cairo_move_to(ctx->cairo, 15, 30);
 	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
-	char buf[1000];
+	char buf[100];
-	snprintf(buf, 1000, "t = exp(%.8f) = %.8f, %s", log(ctx->parameter), ctx->parameter, ctx->extra_text);
+	sprintf(buf, "t = exp(%.8f) = %.8f, marking = (%.5f, %.5f)", log(ctx->parameter), ctx->parameter, ctx->marking.x, ctx->marking.y);
 	cairo_show_text(ctx->cairo, buf);
 }
 
@@ -1512,9 +1111,6 @@ void draw(DrawingContext *ctx)
 		if(ctx->show_marking)
 		{
 			cairo_set_source_rgb(C, 0, 0, 1);
-//			drawPoint(ctx, ctx->marking);
-//			drawPoint(ctx, ctx->marking2);
-//			drawPoint(ctx, ctx->marking3);
 		}
 
 		if(ctx->show_coxeter_orbit)
@@ -1524,7 +1120,7 @@ void draw(DrawingContext *ctx)
 	cairo_identity_matrix(C); // text is in screen coordinates
 
 	if(ctx->show_text)
-		drawText(ctx);
+	    drawText(ctx);
 
 	cairo_surface_flush(cairo_get_target(C));
 }

exp_equation.c

@@ -1,183 +0,0 @@
-#include "main.h"
-#include "exp_equation.h"
-
-#include <stdio.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_roots.h>
-
-#define EPSILON 1e-9
-#define LOOP(i) for(int i = 0; i < 3; i++)
-
-struct solve_exp_plus_exp_params
-{
-	double alpha;
-	double beta;
-	double x;
-};
-
-static double solve_exp_plus_exp_f(double t, void *_params)
-{
-	struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;
-
-//	return exp(params->alpha * t) + exp(params->beta * t) - params->x;
-	return log(exp(params->alpha * t) + exp(params->beta * t)) - log(params->x);
-}
-
-static double solve_exp_plus_exp_df(double t, void *_params)
-{
-	struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;
-
-//	return params->alpha * exp(params->alpha * t) + params->beta * exp(params->beta * t);
-	return params->alpha + (params->beta - params->alpha) / (1 + exp((params->alpha-params->beta)*t));
-}
-
-static void solve_exp_plus_exp_fdf(double t, void *params, double *f, double *df)
-{
-	*f = solve_exp_plus_exp_f(t, params);
-	*df = solve_exp_plus_exp_df(t, params);
-}
-
-// solve the equation exp(alpha t) + exp(beta t) = x for t
-int solve_exp_plus_exp(double alpha, double beta, double x, double *t)
-{
-	if(alpha <= 0 && beta >= 0 || alpha >= 0 && beta <= 0) {
-		double critical =
-			pow(-beta/alpha, alpha/(alpha-beta)) +
-			pow(-beta/alpha, beta/(alpha-beta));
-
-		if(x < critical)
-			return 0;
-		else if (x < critical + EPSILON) {
-			t[0] = log(-beta/alpha)/(alpha-beta);
-			return 1;
-		}
-
-		// Newton this
-		gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
-		struct solve_exp_plus_exp_params params;
-		params.alpha = alpha;
-		params.beta = beta;
-		params.x = x;
-		gsl_function_fdf FDF;
-		FDF.f = &solve_exp_plus_exp_f;
-		FDF.df = &solve_exp_plus_exp_df;
-		FDF.fdf = &solve_exp_plus_exp_fdf;
-		FDF.params = (void *)&params;
-		int status;
-		double root, lastroot;
-
-		for(int r = 0; r < 2; r++) {
-			root = r == 0 ? log(x)/beta : log(x)/alpha;
-			gsl_root_fdfsolver_set(solver, &FDF, root);
-			for(int i = 0; i < 100; i++) {
-				gsl_root_fdfsolver_iterate(solver);
-				lastroot = root;
-				root = gsl_root_fdfsolver_root(solver);
-				status = gsl_root_test_delta(root, lastroot, 0, 1e-9);
-
-				if(status == GSL_SUCCESS)
-					break;
-
-//				printf("iteration %d, root %f\n", i, root);
-			}
-			t[r] = root;
-		}
-
-		gsl_root_fdfsolver_free(solver);
-
-		return 2;
-	} else {
-		// Newton with start value 0
-		gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
-		struct solve_exp_plus_exp_params params;
-		params.alpha = alpha;
-		params.beta = beta;
-		params.x = x;
-		gsl_function_fdf FDF;
-		FDF.f = &solve_exp_plus_exp_f;
-		FDF.df = &solve_exp_plus_exp_df;
-		FDF.fdf = &solve_exp_plus_exp_fdf;
-		FDF.params = (void *)&params;
-		int status;
-		double root, lastroot;
-
-		root = 0;
-		gsl_root_fdfsolver_set(solver, &FDF, root);
-		for(int i = 0; i < 100; i++) {
-			gsl_root_fdfsolver_iterate(solver);
-			lastroot = root;
-			root = gsl_root_fdfsolver_root(solver);
-			status = gsl_root_test_delta(root, lastroot, 0, 1e-9);
-
-//			printf("iteration %d, root %f\n", i, root);
-
-			if(status == GSL_SUCCESS)
-				break;
-		}
-		t[0] = root;
-
-		gsl_root_fdfsolver_free(solver);
-
-		return 1;
-	}
-}
-
-// solve the equation x1 exp(a1 t) + x2 exp(a2 t) + x3 exp(a3 t) = 0
-int solve_linear_exp(vector_t a, vector_t x, double *t)
-{
-	if(x.x[0] > 0 && x.x[1] > 0 && x.x[2] > 0 ||
-	   x.x[0] < 0 && x.x[1] < 0 && x.x[2] < 0)
-		return 0;
-
-	// ensure that y[0] < 0 and y[1], y[2] > 0
-	int j;
-	vector_t y, b;
-	for(j = 0; j < 3; j++)
-		if(x.x[(j+1)%3] * x.x[(j+2)%3] > 0)
-			break;
-	LOOP(i) y.x[i] = x.x[(i+j)%3];
-	if(y.x[0] > 0)
-		LOOP(i) y.x[i] *= -1;
-	LOOP(i) b.x[i] = a.x[(i+j)%3];
-
-	double T = (log(y.x[1]) - log(y.x[2])) / (b.x[2] - b.x[1]);
-	double rhs = - y.x[0] *
-		pow(y.x[1], (b.x[0]-b.x[2])/(b.x[2]-b.x[1])) *
-		pow(y.x[2], (b.x[0]-b.x[1])/(b.x[1]-b.x[2]));
-
-	int n = solve_exp_plus_exp(b.x[1] - b.x[0], b.x[2] - b.x[0], rhs, t);
-
-	for(int i = 0; i < n; i++)
-		t[i] += T;
-
-	return n;
-}
-
-/*
-int main(int argc, char *argv[])
-{
-
-//	int n = solve_exp_plus_exp(atof(argv[1]), atof(argv[2]), atof(argv[3]), result);
-
-	double result[2];
-	vector_t a, x;
-	a.x[0] = atof(argv[1]);
-	a.x[1] = atof(argv[2]);
-	a.x[2] = atof(argv[3]);
-	x.x[0] = atof(argv[4]);
-	x.x[1] = atof(argv[5]);
-	x.x[2] = atof(argv[6]);
-	int n = solve_linear_exp(a, x, result);
-
-	if(n == 0)
-		printf("0 results found\n");
-	else if(n == 1)
-		printf("1 result found: %.9f\n", result[0]);
-	else if(n == 2)
-		printf("2 results found: %.9f and %.9f\n", result[0], result[1]);
-	else
-		printf("%d results found\n", n);
-	return 0;
-}
-*/

exp_equation.h

@@ -1,14 +0,0 @@
-#ifndef EXP_EQUATION_H
-#define EXP_EQUATION_H
-
-#include "main.h"
-
-#include <stdio.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_roots.h>
-
-int solve_exp_plus_exp(double alpha, double beta, double x, double *t);
-int solve_linear_exp(vector_t a, vector_t x, double *t);
-
-#endif

limit_set.c

@@ -63,7 +63,7 @@ int computeLimitCurve(DrawingContext *ctx)
 
 	for(int i = 0; i < ctx->n_group_elements; i++) {
 		multiply_many(ws, fixedpoints_pos, 3, cob_pos, elements[i], coxeter_fixedpoints_pos);
-		ctx->limit_curve[12*i+2] = atan2(
+		ctx->limit_curve[3*i+2] = atan2(
 			gsl_matrix_get(fixedpoints_pos, 2, column)/gsl_matrix_get(fixedpoints_pos, 0, column),
 			gsl_matrix_get(fixedpoints_pos, 1, column)/gsl_matrix_get(fixedpoints_pos, 0, column));
 	}
@@ -85,8 +85,8 @@ int computeLimitCurve(DrawingContext *ctx)
 	for(int i = 0; i < ctx->n_group_elements; i++) {
 		multiply_many(ws, fixedpoints, 3, ctx->cob, elements[i], coxeter_fixedpoints);
 
-		x = ctx->limit_curve[12*i  ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
+		x = ctx->limit_curve[3*i  ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
-		y = ctx->limit_curve[12*i+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
+		y = ctx->limit_curve[3*i+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
 
 		if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10)
 		{
@@ -96,14 +96,11 @@ int computeLimitCurve(DrawingContext *ctx)
 			fputc('\n',stdout);
 		}
 
-		multiply_many(ws, fixedpoints, 2, elements[i], coxeter_fixedpoints);
-		LOOP(j) LOOP(k) ctx->limit_curve[12*i+3+3*j+k] = gsl_matrix_get(fixedpoints, k, j);
-
 		// bca abc acb = abc
 
 	}
 
-	qsort(ctx->limit_curve, ctx->n_group_elements, 12*sizeof(double), compareAngle);
+	qsort(ctx->limit_curve, ctx->n_group_elements, 3*sizeof(double), compareAngle);
 
 	ctx->limit_curve_count = ctx->n_group_elements;
 

linalg.c

@@ -315,6 +315,7 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
 	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;
@@ -331,32 +332,6 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
 	return positive;
 }
 
-int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, workspace_t *ws)
-{
-	gsl_matrix *A_ = getTempMatrix(ws);
-
-	gsl_matrix_memcpy(A_, A);
-	int r = gsl_eigen_symmv (A_, 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(eval)
-			eval[i] = gsl_vector_get(ws->eval_real, i);
-
-		if(gsl_vector_get(ws->eval_real, i) > 0)
-			positive++;
-	}
-
-	releaseTempMatrices(ws, 1);
-	return positive;
-}
-
 // computes a matrix in SL(3, R) which projectively transforms (e1, e2, e3, e1+e2+e3) to the 4 given vectors
 void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws)
 {

linalg.h

@@ -13,7 +13,7 @@
 #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
 #define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
 
-#define MAX_TEMP_MATRICES 600000
+#define MAX_TEMP_MATRICES 1200000
 #define MAX_TEMP_VECTORS 100
 
 typedef struct _workspace {
@@ -53,7 +53,6 @@ int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
 int real_eigenvalues(gsl_matrix *g, gsl_vector *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);
-int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, workspace_t *ws);
 void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws);
 void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws);
 

main.c

@@ -45,17 +45,9 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[])
 	ctx->show_marking = 1;
 	ctx->marking.x = -0.73679;
 	ctx->marking.y = -0.01873;
-	ctx->marking2.x = -0.73679;
-	ctx->marking2.y = -0.11873;
-	ctx->marking3.x = -0.73679;
-	ctx->marking3.y = -0.21873;
-	ctx->distance_parameter1 = 0.45;
-	ctx->distance_parameter2 = 0.2;
 	ctx->show_coxeter_orbit = 0;
-	ctx->extra_text = malloc(1000*sizeof(char));
-	memset(ctx->extra_text, 0, 1000*sizeof(char));
 
-	ctx->limit_curve = malloc(12*ctx->n_group_elements*sizeof(double));
+	ctx->limit_curve = malloc(3*ctx->n_group_elements*sizeof(double));
 	ctx->limit_curve_count = -1;
 
 	ctx->group = malloc(ctx->n_group_elements*sizeof(groupelement_t));
@@ -71,7 +63,6 @@ void destroyContext(DrawingContext *ctx)
 {
 	free(ctx->limit_curve);
 	free(ctx->group);
-	free(ctx->extra_text);
 
 	gsl_matrix_free(ctx->cartan);
 	gsl_matrix_free(ctx->cob);
@@ -79,15 +70,17 @@ void destroyContext(DrawingContext *ctx)
 	workspace_free(ctx->ws);
 }
 
-void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *word)
+void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
 {
 	gsl_matrix *tmp = getTempMatrix(ctx->ws);
 	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
 
+//	ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");
+
 	initializeTriangleGenerators(gen, ctx->cartan);
 	gsl_matrix_set_identity(tmp);
-	for(int i = 0; i < strlen(word); i++)
+	for(int i = 0; i < strlen(type); i++)
-		multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
+		multiply_right(tmp, gen[type[i]-'a'], ctx->ws);
 
 	rotation_frame(tmp, result, ctx->ws);
 
@@ -237,7 +230,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 	case ButtonPress:
 		state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask);
 
-		if(ev->xbutton.button == 1 && state & ShiftMask && state & ControlMask) {
+		if(ev->xbutton.button == 1 && state & ShiftMask) {
 			screen_context->marking.x = (double)ev->xbutton.x;
 			screen_context->marking.y = (double)ev->xbutton.y;
 			printf("mouse button pressed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
@@ -245,30 +238,13 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 			cairo_device_to_user(screen_context->cairo, &screen_context->marking.x, &screen_context->marking.y);
 			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
 			return STATUS_REDRAW;
-		} else if(ev->xbutton.button == 1 && state & ControlMask) {
-			screen_context->marking2.x = (double)ev->xbutton.x;
-			screen_context->marking2.y = (double)ev->xbutton.y;
-			printf("mouse button pressed: %f, %f\n", screen_context->marking2.x, screen_context->marking2.y);
-			cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
-			cairo_device_to_user(screen_context->cairo, &screen_context->marking2.x, &screen_context->marking2.y);
-			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking2.x, screen_context->marking2.y);
-			return STATUS_REDRAW;
-		} else if(ev->xbutton.button == 1 && state & ShiftMask) {
-			screen_context->marking3.x = (double)ev->xbutton.x;
-			screen_context->marking3.y = (double)ev->xbutton.y;
-			printf("mouse button pressed: %f, %f\n", screen_context->marking3.x, screen_context->marking3.y);
-			cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
-			cairo_device_to_user(screen_context->cairo, &screen_context->marking3.x, &screen_context->marking3.y);
-			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking3.x, screen_context->marking3.y);
-			return STATUS_REDRAW;
 		}
-
 		break;
 
 	case KeyPress:
 		state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
 		key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask));
-//		printf("Key pressed: %ld\n", key);
+		printf("Key pressed: %ld\n", key);
 
 		switch(key) {
 		case XK_Down:
@@ -355,16 +331,6 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 		case 'p':
 			print(screen_context);
 			break;
-		case 'P':
-			for(int i = 0; i <= 100; i++) {
-				for(int j = 0; j <= 100; j++) {
-					screen_context->distance_parameter1 = 0.02*i;
-					screen_context->distance_parameter2 = 0.02*j;
-					draw(screen_context);
-				}
-				fflush(stdout);
-			}
-			break;
 		case 'M':
 			/*
 			screen_context->limit_with_lines = 0;
@@ -384,6 +350,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 				printf("Finished drawing %s\n", filename);
 			}
 			*/
+		    /*
 			screen_context->limit_with_lines = 0;
 			double parameter_start = screen_context->parameter;
 			for(int i = 0; i <= 1300; i++) {
@@ -400,7 +367,17 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 				cairo_surface_write_to_png(info->buffer_surface, filename);
 				printf("Finished drawing %s\n", filename);
 			}
-
+		    */
+		    screen_context->limit_with_lines = 0;
+		    for(int i = 0; i <= 1000; i++) {
+			screen_context->parameter = exp(log(3.0)*((double)i/500-1));
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			draw(screen_context);
+			sprintf(filename, "output/movie8/test%04d.png", i);
+			//			cairo_surface_write_to_png(info->buffer_surface, filename);
+			printf("Finished drawing %s\n", filename);
+		    }
 		case 'f':
 			TOGGLE(screen_context->use_repelling);
 			computeLimitCurve(screen_context);
@@ -478,7 +455,7 @@ int main(int argc, char *argv[])
 
 			gettimeofday(&current_time, 0);
 			end_time = current_time.tv_sec + current_time.tv_usec*1e-6;
-//			printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000);
+			printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000);
 		}
 		waitUpdateTimer(info);
 	}

main.h

@@ -11,7 +11,7 @@
 #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
 #define LOOP(i) for(int i = 0; i < 3; i++)
 
-#define NUM_GROUP_ELEMENTS 50000
+#define NUM_GROUP_ELEMENTS 10000
 
 typedef struct {
 	double x[3];
@@ -22,14 +22,6 @@ typedef struct {
 	double y;
 } point_t;
 
-typdef struct {
-	double x;
-	double y;
-	double angle;
-	vector_t fp[3];
-	groupelement_t *g;
-} limit_curve_t;
-
 typedef struct {
 	// infos about the screen to draw on
 	cairo_t *cairo;
@@ -52,14 +44,9 @@ typedef struct {
 	int limit_with_lines;
 	int show_marking;
 	point_t marking;
-	point_t marking2;
-	point_t marking3;
-	double distance_parameter1;
-	double distance_parameter2;
 	int show_coxeter_orbit;
 	int use_repelling;
 	gsl_matrix *cartan, *cob;
-	char *extra_text;
 
     // precomputed and constant
 	groupelement_t *group;