draw x^alpha curves; no want to refactor limit curve

.gitignore

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

Makefile

@@ -1,4 +1,4 @@
-HEADERS=triangle.h linalg.h queue.h initcairo.h main.h
+HEADERS=triangle.h linalg.h queue.h initcairo.h main.h exp_equation.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
-	gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o -lm -lgsl -lcblas -lcairo -lX11
+limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.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
 
 linalg.o: linalg.c $(HEADERS)
 	gcc $(OPTIONS) -c linalg.c
@@ -32,5 +32,8 @@ 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
+	rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o exp_equation.o

README.md

@@ -1,50 +0,0 @@
-# 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.
-
-## Installation & Prerequisites ##
-
-This program is written C and exclusively for Linux. It uses [Cairo] for drawing and X11 for input, which means it will only run in the X window system. Besides these, the only required library is the [Gnu Scientific Library (GSL)] including CBLAS.
-
-To build, 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 1.0
-
-The arguments are 6 integers (p1, p2, p3, q1, q2, q3) and two floating point numbers (t,s) describing the triangle group under consideration (initially, t and s can be changed with the arrow keys).
-The last argument is optional; if it is not given, it defaults to s = 1.
-
-## 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/

draw.c

@@ -1,4 +1,5 @@
 #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))
@@ -22,6 +23,50 @@ 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;
@@ -42,6 +87,133 @@ 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);
@@ -96,10 +268,20 @@ void drawPoint(DrawingContext *ctx, point_t p)
 	cairo_t *C = ctx->cairo;
 
 	cairo_save(C);
-	cairo_arc(C, p.x, p.y, 4.0/ctx->dim->scalefactor, 0, 2*M_PI);
+	cairo_arc(C, p.x, p.y, 3.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)
@@ -183,13 +365,22 @@ 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
@@ -197,7 +388,6 @@ 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);
@@ -218,6 +408,7 @@ 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;
@@ -279,6 +470,8 @@ 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];
@@ -346,6 +539,70 @@ 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);
@@ -558,8 +815,7 @@ void drawAttractors(DrawingContext *ctx)
 	fixedPoints(ctx, "abc", p[0]);
 	fixedPoints(ctx, "bca", p[1]);
 	fixedPoints(ctx, "cab", p[2]);
-	//	fixedPoints(ctx, "abacacbabc", p[3]);
-	//	fixedPoints(ctx, "abcabcabcabcab", p[3]);
+	fixedPoints(ctx, "a cab a", p[3]);
 	fixedPoints(ctx, "b abc b", p[4]);
 	fixedPoints(ctx, "c bca c", p[5]);
 
@@ -567,14 +823,18 @@ 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++)  LOOP(j) {
-		cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
-		drawVector(ctx, p[i][j]);
+	for(int i = 0; i < n; i++) {
+		for(int j = 0; j < 3; j += 1) {
+			cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
+			drawVector(ctx, p[i][j]);
+		}
 	}
 
-	for(int i = 0; i < n; i++)  LOOP(j) {
-		cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
-		drawCovector(ctx, l[i][j]);
+	for(int i = 0; i < n; i++) {
+		for(int j = 0; j < 3; j += 1) {
+			cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]);
+			drawCovector(ctx, l[i][j]);
+		}
 	}
 }
 
@@ -633,6 +893,12 @@ 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]);
@@ -684,6 +950,10 @@ 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);
@@ -703,105 +973,236 @@ 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);
 
-	drawRotationOrbit(ctx, "ab", p[0][0]);
-	drawRotationOrbit(ctx, "bc", p[0][0]);
-	drawRotationOrbit(ctx, "ca", p[0][0]);
+	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]);
 
-	if(ctx->mode >= 2) {
-		drawRotationOrbit(ctx, "bcabcb", p[1][0]); // bcC
-		drawRotationOrbit(ctx, "abcabcba", p[0][0]); // abcC
+	vector_t x;
+	for(int i = 0; i < ctx->n_group_elements; i++) {
+		LOOP(j) x.x[j] = ctx->limit_curve[12*i+3+j];
+		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, 5);
+	releaseTempMatrices(ctx->ws, 9);
 }
 
 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[4][3];
+	vector_t p[3][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]);
 
-	cairo_set_line_width(C, 2.5/ctx->dim->scalefactor);
+	vector_t vertex[3];
+	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]));
 
-	cairo_set_source_rgb(C, 0, 0, 0);
-	drawCurvedBox(ctx, 'A', "", 2);
+	vector_t a,b;
+	double t = ctx->distance_parameter1, s = ctx->distance_parameter2;
 
-	if(ctx->mode >= 3 && ctx->mode != 4) {
-		cairo_set_source_rgb(C, 0, 0, 0);
-		drawCurvedBox(ctx, 'B', "", 2);
-		drawCurvedBox(ctx, 'C', "", 2);
-	}
+	vector_t phiplus, phiminus, eta;
+	double dvert, param;
 
-	if(ctx->mode == 4 || ctx->mode == 5) {
-		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-		drawCurvedBox(ctx, 'A', "", 2);
-		drawCurvedBox(ctx, 'A', "bc", 2);
-		drawCurvedBox(ctx, 'A', "bcbc", 2);
-		drawCurvedBox(ctx, 'A', "bcbcbc", 2);
-		drawCurvedBox(ctx, 'A', "bcbcbcbc", 2);
-		drawCurvedBox(ctx, 'A', "bcbcbcbcbc", 2);
-		drawCurvedBox(ctx, 'A', "bcbcbcbcbcbc", 2);
-	}
+	// want to choose a so that d(vertex[0],a) = t*d(vertex[0],vertex[1])
+	phiplus = cross(p[2][0], p[2][1]);
+	phiminus = cross(p[2][2], p[2][1]);
+	dvert = log(
+		(scal(phiplus, vertex[0])*scal(phiminus, vertex[1])) /
+		(scal(phiplus, vertex[1])*scal(phiminus, vertex[0])));
+	LOOP(i) eta.x[i] =
+		scal(phiminus, vertex[0]) / scal(phiplus, vertex[0])
+		* exp(t*dvert) * phiplus.x[i]
+		- 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];
 
-	if(ctx->mode == 6) {
-		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-		drawCurvedBox(ctx, 'C', "ababab", 2);
-		drawCurvedBox(ctx, 'C', "abababab", 2);
-		drawCurvedBox(ctx, 'C', "ababababab", 2);
-	}
+	// want to choose b so that d(vertex[0],b) = s*d(vertex[0],vertex[2])
+	phiplus = cross(p[1][0], p[1][1]);
+	phiminus = cross(p[1][2], p[1][1]);
+	dvert = log(
+		(scal(phiplus, vertex[0])*scal(phiminus, vertex[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];
 
-	if(ctx->mode >= 6) {
-		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-		drawCurvedBox(ctx, 'C', "ab", 2);
-		drawCurvedBox(ctx, 'C', "abab", 2);
-	}
+	// draw and output stuff
+	LOOP(i) drawVector(ctx, vertex[i]);
+	drawVector(ctx, a);
+	drawVector(ctx, b);
+	drawCovector(ctx, cross(a,b));
+	drawCovector(ctx, cross(b,vertex[0]));
+	drawCovector(ctx, cross(vertex[0],a));
 
-	if(ctx->mode >= 8) {
-		cairo_set_source_rgb(C, 0, 0.8, 0.5);
-		drawCurvedBox(ctx, 'C', "b", 2);
-		drawCurvedBox(ctx, 'C', "bab", 2);
-	}
+	// find the intersections with the limit set and the asymmetric distances
+	int boundary_index;
+	double value;
+	double perimeter1 = 0.0, perimeter2 = 0.0;
+	double approx_perimeter1 = 0.0, approx_perimeter2 = 0.0;
+	vector_t boundary_point;
+	cairo_set_source_rgb(C, 1, 0, 0);
 
-	if(ctx->mode >= 9) {
-	    cairo_set_source_rgb(C, 0.8, 0, 0.5);
-	    drawCurvedBox(ctx, 'B', "abca", 2);
-	    drawCurvedBox(ctx, 'B', "abcaca", 2);
-	    drawCurvedBox(ctx, 'B', "aba", 2);
-	    drawCurvedBox(ctx, 'B', "abaca", 2);
+	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	                        a, b, 1, &boundary_index, &value);
+	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	drawVector(ctx, boundary_point);
+	perimeter1 += log(value);
 
-	    drawCurvedBox(ctx, 'B', "ababca", 2);
-	    drawCurvedBox(ctx, 'B', "ababcaca", 2);
-	    drawCurvedBox(ctx, 'B', "ababa", 2);
-	    drawCurvedBox(ctx, 'B', "ababaca", 2);
+	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	                        b, vertex[0], 1, &boundary_index, &value);
+	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	drawVector(ctx, boundary_point);
+	perimeter1 += log(value);
+	approx_perimeter1 += log(value);
 
-	    drawCurvedBox(ctx, 'B', "bca", 2);
-	    drawCurvedBox(ctx, 'B', "bcaca", 2);
-	    drawCurvedBox(ctx, 'B', "ba", 2);
-	    drawCurvedBox(ctx, 'B', "baca", 2);
+	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	                        vertex[0], a, 1, &boundary_index, &value);
+	LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k];
+	drawVector(ctx, boundary_point);
+	perimeter1 += log(value);
+	approx_perimeter1 += log(value);
 
-	    drawCurvedBox(ctx, 'B', "babca", 2);
-	    drawCurvedBox(ctx, 'B', "babcaca", 2);
-	    drawCurvedBox(ctx, 'B', "baba", 2);
-	    drawCurvedBox(ctx, 'B', "babaca", 2);
+	cairo_set_source_rgb(C, 0, 0, 1);
 
-	    cairo_restore(C);
-	    releaseTempMatrices(ctx->ws, 4);
-	}
+	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
+	                        b, a, 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);
+
+	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)
@@ -939,8 +1340,8 @@ void drawLimitCurve(DrawingContext *ctx)
 
 		for(int i = 0; i < ctx->limit_curve_count; i++) {
 			point_t p;
-			p.x = ctx->limit_curve[3*i];
-			p.y = ctx->limit_curve[3*i+1];
+			p.x = ctx->limit_curve[12*i];
+			p.y = ctx->limit_curve[12*i+1];
 
 			if(isInsideBB(ctx, p)) {
 				if(!previous_inside)
@@ -957,8 +1358,8 @@ void drawLimitCurve(DrawingContext *ctx)
 	} else {
 		for(int i = 0; i < ctx->limit_curve_count; i++) {
 			point_t p;
-			p.x = ctx->limit_curve[3*i];
-			p.y = ctx->limit_curve[3*i+1];
+			p.x = ctx->limit_curve[12*i];
+			p.y = ctx->limit_curve[12*i+1];
 
 			if(isInsideBB(ctx, p)) {
 				cairo_arc(C, p.x, p.y, 2.0/ctx->dim->scalefactor, 0, 2*M_PI);
@@ -1064,8 +1465,8 @@ void drawText(DrawingContext *ctx)
 {
 	cairo_move_to(ctx->cairo, 15, 30);
 	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
-	char buf[100];
-	sprintf(buf, "t = exp(%.8f) = %.8f, marking = (%.5f, %.5f)", log(ctx->parameter), ctx->parameter, ctx->marking.x, ctx->marking.y);
+	char buf[1000];
+	snprintf(buf, 1000, "t = exp(%.8f) = %.8f, %s", log(ctx->parameter), ctx->parameter, ctx->extra_text);
 	cairo_show_text(ctx->cairo, buf);
 }
 
@@ -1111,6 +1512,9 @@ 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)
@@ -1120,7 +1524,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

@@ -0,0 +1,183 @@
+#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

@@ -0,0 +1,14 @@
+#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[3*i+2] = atan2(
+		ctx->limit_curve[12*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[3*i  ] = gsl_matrix_get(fixedpoints, 0, 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);
+		x = ctx->limit_curve[12*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);
 
 		if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10)
 		{
@@ -96,11 +96,14 @@ 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, 3*sizeof(double), compareAngle);
+	qsort(ctx->limit_curve, ctx->n_group_elements, 12*sizeof(double), compareAngle);
 
 	ctx->limit_curve_count = ctx->n_group_elements;
 

linalg.c

@@ -315,7 +315,6 @@ 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;
@@ -332,6 +331,32 @@ 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 1200000
+#define MAX_TEMP_MATRICES 600000
 #define MAX_TEMP_VECTORS 100
 
 typedef struct _workspace {
@@ -53,6 +53,7 @@ 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,9 +45,17 @@ 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(3*ctx->n_group_elements*sizeof(double));
+	ctx->limit_curve = malloc(12*ctx->n_group_elements*sizeof(double));
 	ctx->limit_curve_count = -1;
 
 	ctx->group = malloc(ctx->n_group_elements*sizeof(groupelement_t));
@@ -63,6 +71,7 @@ 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);
@@ -70,17 +79,15 @@ void destroyContext(DrawingContext *ctx)
 	workspace_free(ctx->ws);
 }
 
-void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
+void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *word)
 {
 	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(type); i++)
-		multiply_right(tmp, gen[type[i]-'a'], ctx->ws);
+	for(int i = 0; i < strlen(word); i++)
+		multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
 
 	rotation_frame(tmp, result, ctx->ws);
 
@@ -230,7 +237,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 	case ButtonPress:
 		state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask);
 
-		if(ev->xbutton.button == 1 && state & ShiftMask) {
+		if(ev->xbutton.button == 1 && state & ShiftMask && state & ControlMask) {
 			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);
@@ -238,13 +245,30 @@ 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:
@@ -331,6 +355,16 @@ 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;
@@ -350,7 +384,6 @@ 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++) {
@@ -367,17 +400,7 @@ 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);
@@ -455,7 +478,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 10000
+#define NUM_GROUP_ELEMENTS 50000
 
 typedef struct {
 	double x[3];
@@ -22,6 +22,14 @@ 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;
@@ -44,9 +52,14 @@ 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;