continue limit curve past degeneration point

limit_set.c

@@ -52,6 +52,8 @@ int show_boxes = 0;
 int show_attractors = 0;
 int show_reflectors = 0;
 int show_limit = 1;
+int limit_with_lines = 1;
+int use_repelling = 1;
 
 void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
 {
@@ -157,7 +159,7 @@ void drawPolygon(DrawingContext *ctx, int sides, point_t a, point_t b, point_t c
 	drawSegment(ctx, d, a);
 }
 
-void fixedPoints(const char *word, point_t *out)
+int fixedPoints(const char *word, point_t *out)
 {
 	gsl_matrix *tmp = ws->stack[10];
 	gsl_matrix *ev = ws->stack[11];
@@ -168,11 +170,13 @@ void fixedPoints(const char *word, point_t *out)
 			continue;
 		multiply_right(tmp, gen[word[i]-'a'], ws);
 	}
-	eigenvectors(tmp, ev, ws);
+	int count = real_eigenvectors(tmp, ev, ws);
 	multiply_left(cob, ev, ws);
 
 	LOOP(i) out[i].x = gsl_matrix_get(ev, 0, i) / gsl_matrix_get(ev, 2, i);
 	LOOP(i) out[i].y = gsl_matrix_get(ev, 1, i) / gsl_matrix_get(ev, 2, i);
+
+	return count;
 }
 
 void drawAttractorConnection(DrawingContext *ctx, const char *word1, const char *word2)
@@ -233,6 +237,8 @@ int computeLimitCurve()
 	int p = 5, q = 5, r = 5;
 	int k = 2, l = 2, m = 2;
 
+	int column = use_repelling ? 2 : 0;
+
 	generate_triangle_group(group, n_group_elements, p, q, r);
 
 	// do first in the Fuchsian case to get the angles
@@ -243,14 +249,16 @@ int computeLimitCurve()
 		multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
 	diagonalize_symmetric_form(cartan, cob, ws);
 	multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
-	if(!eigenvectors(coxeter, coxeter_fixedpoints, ws))
+	int ev_count_fuchsian = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);
+
+	if(ev_count_fuchsian != 3)
 		return 0;
 
 	for(int i = 0; i < n_group_elements; i++) {
 		multiply_many(ws, fixedpoints, 3, cob, matrices[i], coxeter_fixedpoints);
 		limit_curve[3*i+2] = atan2(
-			gsl_matrix_get(fixedpoints, 0, 0)/gsl_matrix_get(fixedpoints, 2, 0),
+			gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column),
-			gsl_matrix_get(fixedpoints, 1, 0)/gsl_matrix_get(fixedpoints, 2, 0));
+			gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column));
 	}
 
 	// now to it again to calculate x and y coordinates
@@ -262,14 +270,18 @@ int computeLimitCurve()
 	gsl_matrix_memcpy(cob, cartan); // is this a good choice of basis
 //	gsl_matrix_set_identity(cob);
 	multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
-	if(!eigenvectors(coxeter, coxeter_fixedpoints, ws))
+	int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);
+
+	if(ev_count == 1)
+		column = 0;
+	if(ev_count == 0)
 		return 0;
 
 	for(int i = 0; i < n_group_elements; i++) {
 		multiply_many(ws, fixedpoints, 3, cob, matrices[i], coxeter_fixedpoints);
 
-		limit_curve[3*i] = gsl_matrix_get(fixedpoints, 0, 0)/gsl_matrix_get(fixedpoints, 2, 0);
+		limit_curve[3*i] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
-		limit_curve[3*i+1] = gsl_matrix_get(fixedpoints, 1, 0)/gsl_matrix_get(fixedpoints, 2, 0);
+		limit_curve[3*i+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
 	}
 
 	qsort(limit_curve, n_group_elements, 3*sizeof(double), compareAngle);
@@ -281,14 +293,12 @@ void drawBoxes(DrawingContext *ctx)
 {
 	cairo_t *C = ctx->cairo;
 
-	/*
 	cairo_set_source_rgb(C, 1, 0, 0);
 	drawBoxStd(ctx, "c", 'C');
 	drawBoxStd(ctx, "", 'B');
 	drawBoxStd(ctx, "a", 'A');
 	drawBoxStd(ctx, "", 'C');
 	drawBoxStd(ctx, "b", 'B');
-	*/
 
 	cairo_set_source_rgb(C, 0, 0, 0);
 	drawBoxStd(ctx, "ca", 'A');
@@ -297,9 +307,6 @@ void drawBoxes(DrawingContext *ctx)
 	drawBoxStd(ctx, "acac", 'C');
 	drawBoxStd(ctx, "aca", 'A');
 	drawBoxStd(ctx, "ac", 'C');
-
-	cairo_set_source_rgb(C, 1, 0, 1);
-
 	/*
 	drawBoxStd(ctx, "aca cb", 'B');
 	drawBoxStd(ctx, "aca cbc", 'C');
@@ -307,14 +314,31 @@ void drawBoxes(DrawingContext *ctx)
 	drawBoxStd(ctx, "aca bcbc", 'C');
 	drawBoxStd(ctx, "aca bcb", 'B');
 	drawBoxStd(ctx, "aca bc", 'C');
+
+	drawBoxStd(ctx, "caca cb", 'B');
+	drawBoxStd(ctx, "caca cbc", 'C');
+	drawBoxStd(ctx, "caca cbcb", 'B');
+	drawBoxStd(ctx, "caca bcbc", 'C');
+	drawBoxStd(ctx, "caca bcb", 'B');
+	drawBoxStd(ctx, "caca bc", 'C');
 	*/
 
+	cairo_set_source_rgb(C, 1, 0, 1);
 	drawBoxStd(ctx, "ca bc", 'C');
 	drawBoxStd(ctx, "ca bcb", 'B');
 	drawBoxStd(ctx, "ca bcbc", 'C');
 	drawBoxStd(ctx, "ca cbcb", 'B');
 	drawBoxStd(ctx, "ca cbc", 'C');
 	drawBoxStd(ctx, "ca cb", 'B');
+
+	cairo_set_source_rgb(C, 0, 1, 0);
+//	drawBoxStd(ctx, "ca bc", 'C');
+	drawBoxStd(ctx, "cabc ba", 'A');
+	drawBoxStd(ctx, "cabc bab", 'B');
+	drawBoxStd(ctx, "cabc baba", 'A');
+	drawBoxStd(ctx, "cabc abab", 'B');
+	drawBoxStd(ctx, "cabc aba", 'A');
+	drawBoxStd(ctx, "cabc ab", 'B');
 }
 
 void drawReflectors(DrawingContext *ctx)
@@ -373,7 +397,9 @@ void draw(DrawingContext *ctx)
 
 	if(limit_curve_valid) {
 		if(show_limit) {
-			int last_inside = 0;
+			cairo_save(C);
+
+			int previous_inside = 0;
 			for(int i = 0; i < n_group_elements; i++) {
 				double x = limit_curve[3*i];
 				double y = limit_curve[3*i+1];
@@ -381,19 +407,30 @@ void draw(DrawingContext *ctx)
 				cairo_user_to_device(C, &x, &y);
 
 				if(-x < ctx->width && x < 3*ctx->width && -y < ctx->height && y < 3*ctx->height) {
-					if(!last_inside) {
+					if(limit_with_lines) {
-						cairo_move_to(C, limit_curve[3*i], limit_curve[3*i+1]);
+						if(!previous_inside)
+							cairo_move_to(C, limit_curve[3*i], limit_curve[3*i+1]);
+						else
+							cairo_line_to(C, limit_curve[3*i], limit_curve[3*i+1]);
 					} else {
-						cairo_line_to(C, limit_curve[3*i], limit_curve[3*i+1]);
+						cairo_move_to(C, limit_curve[3*i], limit_curve[3*i+1]);
+						cairo_close_path(C);
 					}
-					last_inside = 1;
+					previous_inside = 1;
 				} else {
-					last_inside = 0;
+					previous_inside = 0;
 				}
 			}
 
+			if(!limit_with_lines) { // draw dots instead of lines
+				cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);
+				cairo_set_line_width(C, 3.0/ctx->scalefactor);
+			}
+
 			cairo_set_source_rgb(C, 0, 0, 0);
 			cairo_stroke(C);
+
+			cairo_restore(C);
 		}
 
 		if(show_boxes)
@@ -411,7 +448,7 @@ void draw(DrawingContext *ctx)
 	cairo_move_to(C, 15, 30);
 	cairo_set_source_rgb(C, 0, 0, 0);
 	char buf[100];
-	sprintf(buf, "t = %.8f", parameter);
+	sprintf(buf, "t = exp(%.8f) = %.8f", log(parameter), parameter);
 	cairo_show_text(C, buf);
 
 	cairo_surface_flush(cairo_get_target(C));
@@ -419,7 +456,7 @@ void draw(DrawingContext *ctx)
 
 void setup()
 {
-	n_group_elements = 100000;
+	n_group_elements = 50000;
 	limit_curve = malloc(3*n_group_elements*sizeof(double));
 	group = malloc(n_group_elements*sizeof(groupelement_t));
 	ws = workspace_alloc(3);
@@ -490,14 +527,30 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 		printf("Key pressed: %ld\n", key);
 
 		switch(key) {
-		case '<':
+		case XK_Down:
 			parameter /= exp(0.002);
 			limit_curve_valid = computeLimitCurve();
 			break;
-		case '>':
+		case XK_Up:
 			parameter *= exp(0.002);
 			limit_curve_valid = computeLimitCurve();
 			break;
+		case XK_Left:
+			parameter /= exp(0.0001);
+			limit_curve_valid = computeLimitCurve();
+			break;
+		case XK_Right:
+			parameter *= exp(0.0001);
+			limit_curve_valid = computeLimitCurve();
+			break;
+		case XK_Page_Down:
+			parameter /= exp(0.2);
+			limit_curve_valid = computeLimitCurve();
+			break;
+		case XK_Page_Up:
+			parameter *= exp(0.2);
+			limit_curve_valid = computeLimitCurve();
+			break;
 		case ' ':
 			parameter = 2.890053638;
 			limit_curve_valid = computeLimitCurve();
@@ -523,12 +576,18 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 		case 'r':
 			show_reflectors = !show_reflectors;
 			break;
+		case 'L':
+			limit_with_lines = !limit_with_lines;
+			break;
 		case 'l':
 			show_limit = !show_limit;
 			break;
 		case 'p':
 			print(info->context);
 			break;
+		case 'f':
+			use_repelling = !use_repelling;
+			limit_curve_valid = computeLimitCurve();
 		}
 
 		return STATUS_REDRAW;

linalg.c

@@ -190,6 +190,29 @@ int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
 	return 1;
 }
 
+// only fills in the real eigenvectors and returns their count
+int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
+{
+	gsl_matrix_memcpy(ws->stack[0], g);
+	gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
+	int r = gsl_eigen_nonsymmv(ws->stack[0], ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
+	ERROR(r, "gsl_eigen_nonsymmv failed!\n");
+
+	gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
+
+	int real = 0;
+
+	for(int i = 0; i < ws->n; i++) {
+		if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) == 0) {// real
+			for(int j = 0; j < ws->n; j++)
+				gsl_matrix_set(evec_real, j, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, j, i)));
+			real++;
+		}
+	}
+
+	return real;
+}
+
 void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)
 {
 	gsl_matrix_memcpy(ws->stack[0], g);

linalg.h

@@ -42,6 +42,7 @@ int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
 double trace(gsl_matrix *g);
 double determinant(gsl_matrix *g, workspace_t *ws);
 int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
+int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
 void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws);
 int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws);