From a8b4bb7c2cdcc95fc8d7b163c7cd1f75fce1c424 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Fri, 29 Jul 2022 14:38:48 +0900 Subject: [PATCH] implement rotation group and play around with it a bit --- draw.c | 396 ++++++++++++++++++++++++++-------------------------- limit_set.c | 107 +++++++++++--- main.c | 52 ++++--- main.h | 15 +- triangle.h | 1 + 5 files changed, 332 insertions(+), 239 deletions(-) diff --git a/draw.c b/draw.c index 17cb4b0..4844638 100644 --- a/draw.c +++ b/draw.c @@ -46,21 +46,22 @@ int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) { gsl_matrix *tmp = getTempMatrix(ctx->ws); gsl_matrix *ev = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); gsl_matrix_set_identity(tmp); for(int i = 0; i < strlen(word); i++) { - if(word[i] == ' ') - continue; - multiply_right(tmp, gen[word[i]-'a'], ctx->ws); + if(word[i] >= 'a' && word[i] <= 'c') + multiply_right(tmp, gen[word[i]-'a'], ctx->ws); + else if(word[i] >= 'A' && word[i] <= 'C') + multiply_right(tmp, gen[word[i]-'A'+3], ctx->ws); } int count = real_eigenvectors(tmp, ev, ctx->ws); LOOP(i) LOOP(j) out[i].x[j] = gsl_matrix_get(ev, j, i); - releaseTempMatrices(ctx->ws, 5); + releaseTempMatrices(ctx->ws, 8); return count; } @@ -69,9 +70,9 @@ int wordEigenvalues(DrawingContext *ctx, const char *word, double *out) { gsl_matrix *tmp = getTempMatrix(ctx->ws); gsl_vector *ev = getTempVector(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); gsl_matrix_set_identity(tmp); for(int i = 0; i < strlen(word); i++) { @@ -83,7 +84,7 @@ int wordEigenvalues(DrawingContext *ctx, const char *word, double *out) LOOP(i) out[i] = gsl_vector_get(ev, i); - releaseTempMatrices(ctx->ws, 4); + releaseTempMatrices(ctx->ws, 7); releaseTempVectors(ctx->ws, 1); return count; @@ -298,8 +299,6 @@ 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]; @@ -329,13 +328,13 @@ void drawRotationOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, vector_t sta vector_t v[3], w; point_t p; double parameter, startangle; - int iterations = 200; + int iterations = 2000; gsl_matrix *inverse = getTempMatrix(ctx->ws); gsl_vector *start_v = getTempVector(ctx->ws); gsl_vector *start_in_frame = getTempVector(ctx->ws); cairo_t *C = ctx->cairo; -// computeRotationMatrix(ctx, frame, word); +// computeRotationMatrixFrame(ctx, frame, word); LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(frame, j, i); LOOP(i) gsl_vector_set(start_v, i, start.x[i]); @@ -371,7 +370,7 @@ void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) { gsl_matrix *frame = getTempMatrix(ctx->ws); - computeRotationMatrix(ctx, frame, word); + computeRotationMatrixFrame(ctx, frame, word); drawRotationOrbitFrame(ctx, frame, start); releaseTempMatrices(ctx->ws, 1); @@ -389,7 +388,7 @@ void drawDualRotationOrbit(DrawingContext *ctx, const char *word, vector_t start gsl_vector *start_in_frame = getTempVector(ctx->ws); cairo_t *C = ctx->cairo; - computeRotationMatrix(ctx, frame, word); + computeRotationMatrixFrame(ctx, frame, word); LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(frame, j, i); LOOP(i) gsl_vector_set(start_v, i, start.x[i]); @@ -434,7 +433,7 @@ void drawArcWithOutput(DrawingContext *ctx, const char *word, vector_t start, ve gsl_vector *vector_in_frame = getTempVector(ctx->ws); cairo_t *C = ctx->cairo; - computeRotationMatrix(ctx, frame, word); + computeRotationMatrixFrame(ctx, frame, word); LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(frame, j, i); LOOP(i) gsl_vector_set(vector, i, start.x[i]); @@ -576,9 +575,9 @@ void drawAttractors(DrawingContext *ctx) vector_t p[6][3]; vector_t l[6][3]; - fixedPoints(ctx, "abc", p[0]); - fixedPoints(ctx, "bca", p[1]); - fixedPoints(ctx, "cab", p[2]); + fixedPoints(ctx, "cba", p[0]); + fixedPoints(ctx, "bac", p[1]); + fixedPoints(ctx, "acb", p[2]); fixedPoints(ctx, "a cab a", p[3]); fixedPoints(ctx, "b abc b", p[4]); fixedPoints(ctx, "c bca c", p[5]); @@ -705,200 +704,180 @@ void drawCurvedBox(DrawingContext *ctx, int base, const char *conj, int style) } } +groupelement_t *left(const char *word, groupelement_t *g) +{ + int n = strlen(word); + + for(int i = n-1; i >= 0; i--) { + if(word[i] == ' ') + continue; + + g = g->adj[word[i]-'a']; + + if(!g) + break; + } + + return g; +} void drawBoxes(DrawingContext *ctx) { - gsl_matrix *rot = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); gsl_matrix *frame = getTempMatrix(ctx->ws); + gsl_matrix *frame2 = getTempMatrix(ctx->ws); + gsl_vector *startpoint_drawbasis = getTempVector(ctx->ws); + gsl_vector *startpoint_globalbasis = getTempVector(ctx->ws); + gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements); + cairo_t *C = ctx->cairo; cairo_save(C); - - vector_t p[22][3]; - vector_t l[22][3]; - vector_t alpha[6]; - vector_t ptmp[3]; - char word[100], word2[100]; - - fixedPoints(ctx, "abc", p[0]); - fixedPoints(ctx, "bca", p[1]); - fixedPoints(ctx, "cab", p[2]); - fixedPoints(ctx, "bacabab", p[3]); - fixedPoints(ctx, "bcacabacb", p[4]); - 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]); + vector_t p[22][3]; + vector_t fp[3]; + vector_t l[22][3]; + vector_t alpha[6]; + vector_t ptmp[3]; + vector_t start; + vector_t start2; + char word[100], word2[100]; - if(ctx->mode >= 2) { - cairo_set_source_rgb(C, 0.6, 0.6, 1); - drawRotationOrbit(ctx, "bcabcb", p[1][0]); // bcC - drawRotationOrbit(ctx, "abcabcba", p[0][0]); // abcC -// drawRotationOrbit(ctx, "bcabcabacb", p[1][0]); // bcabC'' -// drawRotationOrbit(ctx, "bacabcacab", p[3][0]); // bacaC' -// drawRotationOrbit(ctx, "bcacabcacacb", p[4][0]); // bcacaC' bcacabacb + fixedPoints(ctx, "cba", p[0]); + fixedPoints(ctx, "acb", p[1]); + fixedPoints(ctx, "bac", p[2]); + + initializeTriangleGeneratorsCurrent(gen, ctx); + gsl_matrix_set_identity(elements[0]); + for(int i = 1; i < ctx->n_group_elements; i++) { + if(ctx->group[i].length % 2) + continue; + int letter = ROTATION_LETTER(ctx->group[i].letter, ctx->group[i].parent->letter); + multiply(gen[letter], elements[ctx->group[i].parent->parent->id], elements[i]); } - cairo_set_source_rgb(C, 1, 0, 1); -// drawRotationOrbit(ctx, "bacabcacab", p[3][0]); // ababcba - cairo_set_source_rgb(C, 0, 0, 0); - fixedPoints(ctx, "abababcbaba", p[3]); -// drawRotationOrbit(ctx, "abababcabcbababa", p[3][0]); // bab abc bab - - cairo_set_source_rgb(C, 0, 0, 1); -// fixedPoints(ctx, "cab", p[3]); -// drawRotationOrbit(ctx, "cabc", p[3][0]); -// fixedPoints(ctx, "bca", p[3]); -// drawRotationOrbit(ctx, "bcabcb", p[3][0]); -// fixedPoints(ctx, "bc abc cb", p[3]); -// drawRotationOrbit(ctx, "bcabcb", p[3][0]); - - fixedPoints(ctx, "bc bca cb", p[3]); -// drawRotationOrbit(ctx, "bcbcacbc", p[3][0]); + gsl_vector_set(startpoint_drawbasis, 0, ctx->marking.x); + gsl_vector_set(startpoint_drawbasis, 1, ctx->marking.y); + gsl_vector_set(startpoint_drawbasis, 2, 1); + solve(ctx->cob, startpoint_drawbasis, startpoint_globalbasis, ctx->ws); + LOOP(i) start.x[i] = gsl_vector_get(startpoint_globalbasis, i); /* - cairo_set_source_rgb(C, 0, 0, 0); - strncpy(word,"abc",100); - for(int i = 0; i < 9; i++) { - conjugate_word(word, 0, "ab", word2); - strncpy(word, word2, 100); - fixedPoints(ctx, word, ptmp); - drawVector(ctx, ptmp[0]); -// drawVector(ctx, ptmp[2]); + fixedPoints(ctx, "ABAB", fp); + drawCovector(ctx, cross(fp[0], fp[2])); + drawVector(ctx, fp[1]); + computeRotationMatrixFrame(ctx, frame, "c"); + drawRotationOrbit(ctx, "c", start); // ba cb ab ac = C A c B + drawVector(ctx, start); + + computeMatrix(ctx, tmp, "ABABc"); + start2 = apply(tmp, start); + computeRotationMatrixFrame(ctx, frame, "c"); + multiply(tmp, frame, frame2); + drawVector(ctx, start2); + drawRotationOrbitFrame(ctx, frame2, apply(tmp, start)); + + computeMatrix(ctx, tmp, "ABABcABABC"); + start2 = apply(tmp, start); + computeRotationMatrixFrame(ctx, frame, "c"); + multiply(tmp, frame, frame2); + drawVector(ctx, start2); + drawRotationOrbitFrame(ctx, frame2, apply(tmp, start)); + */ + + for(groupelement_t *cur = &ctx->group[180]; cur->parent; cur = cur->parent) + fputc('a'+cur->letter, stdout); + fputc('\n', stdout); + + queue_t queue; + queue_init(&queue); + queue_put(&queue, 0); + int current; + groupelement_t *cur, *next; + + for(int i = 0; i < ctx->n_group_elements_combinatorial; i++) + ctx->group[i].visited = 0; + + cur = &ctx->group[0]; + + computeRotationMatrixFrame(ctx, frame, "c"); + + for(int i = 0; i < 1000; i++) { + if(ctx->group[i].length % 2 != 0) + continue; + + multiply(elements[i], frame, frame2); + drawRotationOrbitFrame(ctx, frame2, apply(elements[i], start)); } - strncpy(word,"bca",100); - for(int i = 0; i < 9; i++) { - conjugate_word(word, 0, "ab", word2); - strncpy(word, word2, 100); - fixedPoints(ctx, word, ptmp); - drawVector(ctx, ptmp[0]); -// drawVector(ctx, ptmp[2]); - } + /* + while((current = queue_get(&queue)) != -1) { + cur = &ctx->group[current]; + if(cur->visited > 4) + continue; - strncpy(word,"abc",100); - for(int i = 0; i < 9; i++) { - conjugate_word(word, 0, "bc", word2); - strncpy(word, word2, 100); - fixedPoints(ctx, word, ptmp); - drawVector(ctx, ptmp[0]); - } - - strncpy(word,"cab",100); - for(int i = 0; i < 9; i++) { - conjugate_word(word, 0, "bc", word2); - strncpy(word, word2, 100); - fixedPoints(ctx, word, ptmp); - drawVector(ctx, ptmp[0]); - } - - strncpy(word,"cab",100); - for(int i = 0; i < 9; i++) { - conjugate_word(word, 0, "ca", word2); - strncpy(word, word2, 100); - fixedPoints(ctx, word, ptmp); - drawVector(ctx, ptmp[0]); - } - - strncpy(word,"abc",100); - for(int i = 0; i < 9; i++) { - conjugate_word(word, 0, "ca", word2); - strncpy(word, word2, 100); - fixedPoints(ctx, word, ptmp); - drawVector(ctx, ptmp[0]); + if(cur->id < ctx->n_group_elements) { + multiply(elements[cur->id], frame, frame2); + drawRotationOrbitFrame(ctx, frame2, apply(elements[cur->id], start)); } + + next = left("ab ab", cur); + if(next && next->visited == 0) { + queue_put(&queue, next->id); + next->visited = cur->visited+1; + } + next = left("cbac cbac", cur); + if(next && next->visited == 0) { + queue_put(&queue, next->id); + next->visited = cur->visited+1; + } + next = left("cacbca cacbca", cur); + if(next && next->visited == 0) { + queue_put(&queue, next->id); + next->visited = cur->visited+1; + } + next = left("cabcacbcac cabcacbcac", cur); + if(next && next->visited == 0) { + queue_put(&queue, next->id); + next->visited = cur->visited+1; + } + next = left("acbcacba acbcacba", cur); + if(next && next->visited == 0) { + queue_put(&queue, next->id); + next->visited = cur->visited+1; + } + next = left("bcacbc bcacbc", cur); + if(next && next->visited == 0) { + queue_put(&queue, next->id); + next->visited = cur->visited+1; + } + } */ - /* - cairo_set_source_rgb(C, 1, 0, 0); - drawVector(ctx, p[0][0]); - cairo_set_source_rgb(C, 0, 0.6, 0); - drawVector(ctx, p[1][0]); - cairo_set_source_rgb(C, 0, 0, 1); - drawVector(ctx, p[2][0]); - */ - - /* - fixedPoints(ctx, "ab abc ba", p[4]); - fixedPoints(ctx, "abab abc baba", p[5]); - fixedPoints(ctx, "ababab abc bababa", p[6]); - fixedPoints(ctx, "abababab abc babababa", p[7]); - fixedPoints(ctx, "babababa abc abababab", p[8]); - fixedPoints(ctx, "bababa abc ababab", p[9]); - fixedPoints(ctx, "baba abc abab", p[10]); - fixedPoints(ctx, "ba abc ab", p[11]); - fixedPoints(ctx, "bca", p[12]); - fixedPoints(ctx, "b abc b", p[13]); - fixedPoints(ctx, "bab abc bab", p[14]); - fixedPoints(ctx, "babab abc babab", p[15]); - fixedPoints(ctx, "bababab abc bababab", p[16]); - fixedPoints(ctx, "abababab bca babababa", p[17]); - fixedPoints(ctx, "ababab bca bababa", p[18]); - fixedPoints(ctx, "abab bca baba", p[19]); - fixedPoints(ctx, "ab bca ba", p[20]); - */ - -// initializeTriangleGenerators(gen, ctx->cartan); - -// for(int i = 0; i < 22; i++) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); -// LOOP(i) LOOP(j) alpha[i].x[j] = gsl_matrix_get(ctx->cartan, i, j); - - /* - gsl_matrix_set(frame, 0, 0, 2.0); - gsl_matrix_set(frame, 0, 1, 0.0); - gsl_matrix_set(frame, 0, 2, 1.0); - gsl_matrix_set(frame, 1, 0, -1.0); - gsl_matrix_set(frame, 1, 1, sqrt(3)); - gsl_matrix_set(frame, 1, 2, 1.0); - gsl_matrix_set(frame, 2, 0, -1.0); - gsl_matrix_set(frame, 2, 1, -sqrt(3)); - gsl_matrix_set(frame, 2, 2, 1.0);*/ -// drawRotationOrbitFrame(ctx, frame, p[0][0]); - -// drawRotationOrbit(ctx, "bc", p[0][0]); - -// drawRotationOrbit(ctx, "ca", p[0][0]); - -/* - for(int i = 0; i < 18; i++) { - if(i == 0) - cairo_set_source_rgb(C, 1, 0, 0); - else if(i == 8) - cairo_set_source_rgb(C, 0, 0, 1); - else if(i == 9) - cairo_set_source_rgb(C, 0, 0.6, 0); - else - cairo_set_source_rgb(C, 0, 0, 0); - drawVector(ctx, p[3+i][0]); -// drawCovector(ctx, l[3+i][0]); -} -*/ - -// drawRotationOrbit(ctx, "ab", cross(l[0][0], l[2][1])); - -// drawRotationOrbit(ctx, "abca", p[0][0]); +// drawRotationOrbit(ctx, "a", p[1][0]); +// drawRotationOrbit(ctx, "b", p[2][0]); cairo_restore(C); - releaseTempMatrices(ctx->ws, 5); + releaseTempMatrices(ctx->ws, 9 + ctx->n_group_elements); + releaseTempVectors(ctx->ws, 2); } void drawBoxes2(DrawingContext *ctx) { gsl_matrix *rot = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); cairo_t *C = ctx->cairo; cairo_save(C); - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); vector_t p[4][3]; - fixedPoints(ctx, "abc", p[0]); - fixedPoints(ctx, "bca", p[1]); - fixedPoints(ctx, "cab", p[2]); + fixedPoints(ctx, "cba", p[0]); + fixedPoints(ctx, "acb", p[1]); + fixedPoints(ctx, "bac", p[2]); cairo_set_line_width(C, 2.5/ctx->dim->scalefactor); @@ -1262,13 +1241,13 @@ void drawBoxes2(DrawingContext *ctx) */ cairo_restore(C); - releaseTempMatrices(ctx->ws, 4); + releaseTempMatrices(ctx->ws, 7); } void drawRotatedReflectors(DrawingContext *ctx) { gsl_matrix *rot = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); cairo_t *C = ctx->cairo; vector_t fp[3], fp2[3]; vector_t w; @@ -1276,7 +1255,7 @@ void drawRotatedReflectors(DrawingContext *ctx) cairo_set_source_rgb(C, 0.7, 0.7, 0.7); - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(ctx->cartan, i, j); multiply(gen[0], gen[1], rot); @@ -1297,7 +1276,7 @@ void drawRotatedReflectors(DrawingContext *ctx) fixedPoints(ctx, "cacabac", fp2); drawRotationOrbit(ctx, "ac", fp[0]); - releaseTempMatrices(ctx->ws, 4); + releaseTempMatrices(ctx->ws, 7); } void drawDualLimitCurve(DrawingContext *ctx) @@ -1307,12 +1286,36 @@ void drawDualLimitCurve(DrawingContext *ctx) cairo_save(C); cairo_set_source_rgb(C, 0.5, 0.5, 1); - int n = 18; - vector_t p[n][3]; - vector_t l[n][3]; - vector_t ptmp[3], ltmp[3]; + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); + gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements); - fixedPoints(ctx, "abc", p[0]); +// wordEigenvalues(ctx, "abc", ev); +// LOOP(i) LOOP(j) gsl_matrix_set(coxeter_fixedpoints, j, i, cox[0][i].x[j]); + + initializeTriangleGeneratorsCurrent(gen, ctx); + gsl_matrix_set_identity(elements[0]); + for(int i = 1; i < ctx->n_group_elements; i++) + multiply(gen[ctx->group[i].letter], elements[ctx->group[i].parent->id], elements[i]); + + vector_t p[3], l[3], v; + + fixedPoints(ctx, "cba", p); + drawVector(ctx, p[0]); + drawVector(ctx, p[1]); + drawVector(ctx, p[2]); + + LOOP(i) l[i] = cross(p[(i+1)%3], p[(i+2)%3]); + + for(int i = 0; i < ctx->n_group_elements; i++) { + v = apply_transpose(elements[i], l[0]); + drawCovector(ctx, v); + } + + releaseTempMatrices(ctx->ws, 3 + ctx->n_group_elements); +// releaseTempVectors(ctx->ws, 4); + + +/* fixedPoints(ctx, "ab abc ba", p[1]); fixedPoints(ctx, "abab abc baba", p[2]); fixedPoints(ctx, "ababab abc bababa", p[3]); @@ -1331,6 +1334,7 @@ void drawDualLimitCurve(DrawingContext *ctx) fixedPoints(ctx, "bababa bca ababab", p[15]); fixedPoints(ctx, "baba bca abab", p[16]); fixedPoints(ctx, "ba bca ab", p[17]); +*/ /* fixedPoints(ctx, "abc", p[0]); @@ -1368,6 +1372,7 @@ void drawDualLimitCurve(DrawingContext *ctx) // drawCovector(ctx, l[i][2]); }*/ + /* fixedPoints(ctx, "abc", ptmp); drawCovector(ctx, cross(ptmp[0], ptmp[1])); fixedPoints(ctx, "bca", ptmp); @@ -1380,6 +1385,7 @@ void drawDualLimitCurve(DrawingContext *ctx) drawCovector(ctx, cross(ptmp[0], ptmp[1])); fixedPoints(ctx, "acaba", ptmp); drawCovector(ctx, cross(ptmp[0], ptmp[1])); + */ cairo_restore(C); } @@ -1435,7 +1441,7 @@ void drawLimitCurve(DrawingContext *ctx) void drawCoxeterOrbit(DrawingContext *ctx) { gsl_matrix *rot = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); gsl_vector *eval = getTempVector(ctx->ws); gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws); gsl_vector *startpoint_coxeterbasis = getTempVector(ctx->ws); @@ -1452,7 +1458,7 @@ void drawCoxeterOrbit(DrawingContext *ctx) int first = 1; cairo_save(C); - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); cairo_set_source_rgb(C, 0, 0, 1); @@ -1463,7 +1469,7 @@ void drawCoxeterOrbit(DrawingContext *ctx) wordEigenvalues(ctx, "abc", ev); LOOP(i) LOOP(j) gsl_matrix_set(coxeter_fixedpoints, j, i, cox[0][i].x[j]); - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); gsl_matrix_set_identity(elements[0]); for(int i = 1; i < ctx->n_group_elements; i++) multiply(gen[ctx->group[i].letter], elements[ctx->group[i].parent->id], elements[i]); @@ -1477,7 +1483,7 @@ void drawCoxeterOrbit(DrawingContext *ctx) gsl_vector_set(startpoint_drawbasis, 2, 1); solve(ctx->cob, startpoint_drawbasis, startpoint_globalbasis, ctx->ws); - solve(coxeter_fixedpoints, startpoint_globalbasis, startpoint_coxeterbasis, ctx->ws); +// solve(coxeter_fixedpoints, startpoint_globalbasis, startpoint_coxeterbasis, ctx->ws); // LOOP(i) start.x[i] = gsl_vector_get(startpoint_coxeterbasis, i); LOOP(i) start.x[i] = gsl_vector_get(startpoint_globalbasis, i); @@ -1517,7 +1523,7 @@ void drawCoxeterOrbit(DrawingContext *ctx) // LOOP(i) drawVector(ctx, abcb[i]); cairo_restore(C); - releaseTempMatrices(ctx->ws, 5 + ctx->n_group_elements); + releaseTempMatrices(ctx->ws, 8 + ctx->n_group_elements); releaseTempVectors(ctx->ws, 4); } @@ -1526,7 +1532,7 @@ 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); + sprintf(buf, "t = exp(%.8f) = %.8f, s = exp(%.8f) = %.8f, marking = (%.5f, %.5f)", log(ctx->parameter), ctx->parameter, log(ctx->parameter2), ctx->parameter2, ctx->marking.x, ctx->marking.y); cairo_show_text(ctx->cairo, buf); } diff --git a/limit_set.c b/limit_set.c index 91537b3..ef5233a 100644 --- a/limit_set.c +++ b/limit_set.c @@ -5,6 +5,7 @@ static int compareAngle(const void *x, const void *y) return ((double*)x)[2] > ((double*)y)[2] ? 1 : -1; } +// might need a rewrite void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s) { gsl_matrix_set(cartan, 0, 0, 2); @@ -20,11 +21,58 @@ void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s) gsl_matrix_set(cartan, 2, 2, 2); } -void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan) +void initializeTriangleGenerators(gsl_matrix **gen, double a1, double a2, double a3, double s, double t, workspace_t *ws) { - LOOP(i) gsl_matrix_set_identity(gen[i]); - LOOP(i) LOOP(j) *gsl_matrix_ptr(gen[i], j, j) = -1.0; - LOOP(i) LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j); + gsl_matrix *reflection_gen[3]; + + LOOP(i) { + reflection_gen[i] = gsl_matrix_alloc(3, 3); + gsl_matrix_set_identity(reflection_gen[i]); + } + + double rho[3]; + rho[0] = sqrt(s*s + 2*s*cos(a1) + 1); + rho[1] = sqrt(s*s + 2*s*cos(a2) + 1); + rho[2] = sqrt(s*s + 2*s*cos(a3) + 1); + + gsl_matrix_set(reflection_gen[0], 0, 0, -1.0); + gsl_matrix_set(reflection_gen[0], 0, 1, rho[2]*t); + gsl_matrix_set(reflection_gen[0], 0, 2, rho[1]/t); + + gsl_matrix_set(reflection_gen[1], 1, 0, rho[2]/t); + gsl_matrix_set(reflection_gen[1], 1, 1, -1.0); + gsl_matrix_set(reflection_gen[1], 1, 2, rho[0]*t); + + gsl_matrix_set(reflection_gen[2], 2, 0, rho[1]*t); + gsl_matrix_set(reflection_gen[2], 2, 1, rho[0]/t); + gsl_matrix_set(reflection_gen[2], 2, 2, -1.0); + + LOOP(i) { + gsl_matrix_set_identity(gen[i]); + gsl_matrix_set(gen[i], (i+1)%3, (i+1)%3, s); + gsl_matrix_set(gen[i], (i+2)%3, (i+2)%3, 1/s); + + gsl_matrix_set_identity(gen[i+3]); + gsl_matrix_set(gen[i+3], (i+1)%3, (i+1)%3, 1/s); + gsl_matrix_set(gen[i+3], (i+2)%3, (i+2)%3, s); + } + + LOOP(i) { + multiply_left(reflection_gen[i], gen[(i+2)%3], ws); + multiply_right(gen[(i+2)%3], reflection_gen[(i+1)%3], ws); + + multiply_left(reflection_gen[(i+1)%3], gen[(i+2)%3+3], ws); + multiply_right(gen[(i+2)%3+3], reflection_gen[i], ws); + } + + LOOP(i) gsl_matrix_free(reflection_gen[i]); +} + +void initializeTriangleGeneratorsCurrent(gsl_matrix **gen, DrawingContext *ctx) +{ + double angle[3]; + LOOP(i) angle[i] = 2*M_PI*ctx->k[i]/ctx->p[i]; + initializeTriangleGenerators(gen, angle[0], angle[1], angle[2], ctx->parameter2, ctx->parameter, ctx->ws); } int computeLimitCurve(DrawingContext *ctx) @@ -38,7 +86,7 @@ int computeLimitCurve(DrawingContext *ctx) gsl_matrix *coxeter = getTempMatrix(ctx->ws); gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws); gsl_matrix *fixedpoints = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements); groupelement_t *group = ctx->group; int success = 0; @@ -50,31 +98,43 @@ int computeLimitCurve(DrawingContext *ctx) // do first in the Fuchsian positive case to get the angles cartanMatrix(cartan_pos, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0); - initializeTriangleGenerators(gen, cartan_pos); + initializeTriangleGenerators(gen, 2*M_PI/ctx->p[0], 2*M_PI/ctx->p[1], 2*M_PI/ctx->p[2], 1.0, 1.0, ctx->ws); gsl_matrix_set_identity(elements[0]); - for(int i = 1; i < ctx->n_group_elements; i++) - multiply(gen[group[i].letter], elements[group[i].parent->id], elements[i]); + for(int i = 1; i < ctx->n_group_elements; i++) { + if(group[i].length % 2) + continue; + int letter = ROTATION_LETTER(group[i].letter, group[i].parent->letter); + multiply(gen[letter], elements[group[i].parent->parent->id], elements[i]); + } diagonalize_symmetric_form(cartan_pos, cob_pos, ws); - multiply_many(ws, coxeter_pos, 3, gen[0], gen[1], gen[2]); + multiply_many(ws, coxeter_pos, 3, gen[2], gen[1], gen[0]); int ev_count_pos = real_eigenvectors(coxeter_pos, coxeter_fixedpoints_pos, ws); if(ev_count_pos != 3) goto error_out; + int n = 0; for(int i = 0; i < ctx->n_group_elements; i++) { + if(group[i].length % 2) + continue; multiply_many(ws, fixedpoints_pos, 3, cob_pos, elements[i], coxeter_fixedpoints_pos); - ctx->limit_curve[3*i+2] = atan2( + ctx->limit_curve[3*n+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)); + n++; } // now do it again to calculate x and y coordinates - initializeTriangleGenerators(gen, ctx->cartan); + initializeTriangleGeneratorsCurrent(gen, ctx); gsl_matrix_set_identity(elements[0]); - for(int i = 1; i < ctx->n_group_elements; i++) - multiply(gen[group[i].letter], elements[group[i].parent->id], elements[i]); + for(int i = 1; i < ctx->n_group_elements; i++) { + if(group[i].length % 2) + continue; + int letter = ROTATION_LETTER(group[i].letter, group[i].parent->letter); + multiply(gen[letter], elements[group[i].parent->parent->id], elements[i]); + } - multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]); + multiply_many(ws, coxeter, 3, gen[2], gen[1], gen[0]); int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws); if(ev_count == 1) @@ -82,11 +142,17 @@ int computeLimitCurve(DrawingContext *ctx) if(ev_count == 0) goto error_out; + ctx->limit_curve_count = 0; for(int i = 0; i < ctx->n_group_elements; i++) { + if(group[i].length % 2) + continue; + 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[3*ctx->limit_curve_count ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column); + y = ctx->limit_curve[3*ctx->limit_curve_count+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column); + + ctx->limit_curve_count++; if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10) { @@ -95,19 +161,16 @@ int computeLimitCurve(DrawingContext *ctx) fputc('a' + cur->letter, stdout); // bcbcbca, bacbcacab, bc bca cb fputc('\n',stdout); } - - // bca abc acb = abc - } - qsort(ctx->limit_curve, ctx->n_group_elements, 3*sizeof(double), compareAngle); + qsort(ctx->limit_curve, ctx->limit_curve_count, 3*sizeof(double), compareAngle); - ctx->limit_curve_count = ctx->n_group_elements; +// ctx->limit_curve_count = ctx->n_group_elements; success = 1; error_out: - releaseTempMatrices(ctx->ws, 11+ctx->n_group_elements); + releaseTempMatrices(ctx->ws, 14+ctx->n_group_elements); return success; } diff --git a/main.c b/main.c index 3803ced..03e1c0b 100644 --- a/main.c +++ b/main.c @@ -18,6 +18,7 @@ DrawingContext *screen_context; void setupContext(DrawingContext *ctx, int argc, char *argv[]) { ctx->n_group_elements = NUM_GROUP_ELEMENTS; + ctx->n_group_elements_combinatorial = NUM_GROUP_ELEMENTS_COMBINATORIAL; ctx->p[0] = atoi(argv[1]); ctx->p[1] = atoi(argv[2]); ctx->p[2] = atoi(argv[3]); @@ -28,6 +29,10 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]) ctx->parameter = atof(argv[7]); else ctx->parameter = 1.0; + if(argc > 8) + ctx->parameter2 = atof(argv[8]); + else + ctx->parameter2 = 1.0; // ctx->parameter = 2.77; // ctx->parameter = 0.1; ctx->show_boxes = 0; @@ -35,12 +40,12 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]) ctx->show_attractors = 0; ctx->show_reflectors = 0; ctx->show_rotated_reflectors = 0; - ctx->show_limit= 0; - ctx->show_dual_limit= 0; + ctx->show_limit = 0; + ctx->show_dual_limit = 0; ctx->show_text = 1; ctx->mode = 0; - ctx->use_rotation_basis = 0; - ctx->limit_with_lines = 1; + ctx->use_rotation_basis = 2; + ctx->limit_with_lines = 0; ctx->use_repelling = 0; ctx->show_marking = 1; ctx->marking.x = -0.73679; @@ -50,8 +55,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]) 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)); - generate_triangle_group(ctx->group, ctx->n_group_elements, ctx->p[0], ctx->p[1], ctx->p[2]); + ctx->group = malloc(ctx->n_group_elements_combinatorial*sizeof(groupelement_t)); + generate_triangle_group(ctx->group, ctx->n_group_elements_combinatorial, ctx->p[0], ctx->p[1], ctx->p[2]); // the temporary stuff ctx->cartan = gsl_matrix_alloc(3, 3); @@ -70,21 +75,30 @@ void destroyContext(DrawingContext *ctx) workspace_free(ctx->ws); } -void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type) +void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type) { - gsl_matrix *tmp = getTempMatrix(ctx->ws); - gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **gen = getTempMatrices(ctx->ws, 6); // 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); + initializeTriangleGeneratorsCurrent(gen, ctx); + gsl_matrix_set_identity(result); + for(int i = 0; i < strlen(type); i++) { + if(type[i] >= 'a' && type[i] <= 'c') + multiply_right(result, gen[type[i]-'a'], ctx->ws); + else if(type[i] >= 'A' && type[i] <= 'C') + multiply_right(result, gen[type[i]-'A'+3], ctx->ws); + } + releaseTempMatrices(ctx->ws, 6); +} + +void computeRotationMatrixFrame(DrawingContext *ctx, gsl_matrix *result, const char *type) +{ + gsl_matrix *tmp = getTempMatrix(ctx->ws); + computeMatrix(ctx, tmp, type); rotation_frame(tmp, result, ctx->ws); - - releaseTempMatrices(ctx->ws, 4); + releaseTempMatrices(ctx->ws, 1); } void computeBoxTransform(DrawingContext *ctx, char *word1, char *word2, gsl_matrix *result) @@ -159,10 +173,10 @@ void updateMatrices(DrawingContext *ctx) } else if(ctx->use_rotation_basis % 5 == 1) { gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason? } else if(ctx->use_rotation_basis % 5 == 2) { - computeRotationMatrix(ctx, tmp, "ba"); + computeRotationMatrixFrame(ctx, tmp, "C"); invert(tmp, ctx->cob, ctx->ws); } else if(ctx->use_rotation_basis % 5 == 3) { - computeBoxTransform(ctx, "bca", "abc", ctx->cob); + computeBoxTransform(ctx, "acb", "cba", ctx->cob); // computeBoxTransform(ctx, "cab", "bca", ctx->cob); // computeBoxTransform(ctx, "acb", "cba", ctx->cob); } else { @@ -258,12 +272,12 @@ int processEvent(GraphicsInfo *info, XEvent *ev) computeLimitCurve(screen_context); break; case XK_Left: - screen_context->parameter /= exp(0.00002); + screen_context->parameter2 /= exp(0.002); updateMatrices(screen_context); computeLimitCurve(screen_context); break; case XK_Right: - screen_context->parameter *= exp(0.00002); + screen_context->parameter2 *= exp(0.002); updateMatrices(screen_context); computeLimitCurve(screen_context); break; diff --git a/main.h b/main.h index 9829f45..f365728 100644 --- a/main.h +++ b/main.h @@ -11,7 +11,12 @@ #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 50000 +#define NUM_GROUP_ELEMENTS 10000 +#define NUM_GROUP_ELEMENTS_COMBINATORIAL 100000 + +// (0,1) -> 2, (1,2) -> 0, (2,0) -> 1 +// (1,0) -> 5, (2,1) -> 3, (0,2) -> 4 +#define ROTATION_LETTER(x,y) (((y)-(x)+3)%3 == 1 ? ((y)+1)%3 : ((x)+1)%3+3) typedef struct { double x[3]; @@ -30,7 +35,9 @@ typedef struct { // a priori parameter int p[3],k[3]; double parameter; + double parameter2; int n_group_elements; + int n_group_elements_combinatorial; int show_boxes; int show_boxes2; int show_attractors; @@ -66,7 +73,8 @@ typedef enum { // implemented in limit_set.c void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s); -void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan); +void initializeTriangleGenerators(gsl_matrix **gen, double a1, double a2, double a3, double s, double t, workspace_t *ws); +void initializeTriangleGeneratorsCurrent(gsl_matrix **gen, DrawingContext *ctx); int computeLimitCurve(DrawingContext *ctx); // implemented in draw.c @@ -96,7 +104,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]); void destroyContext(DrawingContext *ctx); void print(DrawingContext *screen); int processEvent(GraphicsInfo *info, XEvent *ev); -void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type); +void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type); +void computeRotationMatrixFrame(DrawingContext *ctx, gsl_matrix *result, const char *type); void updateMatrices(DrawingContext *ctx); static vector_t vectorFromGsl(gsl_vector *v) diff --git a/triangle.h b/triangle.h index 076f6d4..65f2221 100644 --- a/triangle.h +++ b/triangle.h @@ -14,6 +14,7 @@ typedef struct _groupelement { struct _groupelement *parent; struct _groupelement *inverse; int letter; + int visited; } groupelement_t; int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);