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

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

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

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

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