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
-	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
+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
 
 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

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);
-		y = ctx->limit_curve[12*i+1] = gsl_matrix_get(fixedpoints, 1, 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[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

@@ -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++)
-		multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
+	for(int i = 0; i < strlen(type); i++)
+		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;
@@ -478,7 +444,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

@@ -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;