From e631dee6613d33b1d99c6b2b73efa9e12e0f8e94 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sat, 27 Oct 2018 14:56:46 -0700 Subject: [PATCH] continue limit curve past degeneration point --- limit_set.c | 105 ++++++++++++++++++++++++++++++++++++++++------------ linalg.c | 23 ++++++++++++ linalg.h | 1 + 3 files changed, 106 insertions(+), 23 deletions(-) diff --git a/limit_set.c b/limit_set.c index b5edf73..3606f24 100644 --- a/limit_set.c +++ b/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, 1, 0)/gsl_matrix_get(fixedpoints, 2, 0)); + gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column), + 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+1] = gsl_matrix_get(fixedpoints, 1, 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, 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) { - cairo_move_to(C, limit_curve[3*i], limit_curve[3*i+1]); + if(limit_with_lines) { + 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; diff --git a/linalg.c b/linalg.c index 472e8ff..998496f 100644 --- a/linalg.c +++ b/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); diff --git a/linalg.h b/linalg.h index 880d160..5e5a946 100644 --- a/linalg.h +++ b/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);