florian/triangle_group_limit_set
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Compare commits

base: florian:public
Branches Tags
florian:public
florian:master
florian:asymmetric_distance
..
compare: florian:master
Branches Tags
florian:public
florian:master
florian:asymmetric_distance

3 Commits
public ... master

Author SHA1 Message Date
Florian Stecker
78769593a7 movie command line arguments 2023-02-05 15:10:16 -05:00
Florian Stecker
a8b4bb7c2c implement rotation group and play around with it a bit 2022-07-29 14:38:48 +09:00
Florian Stecker
35932782e9 new .gitignore 2022-03-02 10:53:19 -06:00
8 changed files with 1040 additions and 526 deletions
Show Stats Expand all files Collapse all files

3
.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
*.o
limit_set
*.pdf

50
README.md
View File

@ -1,50 +0,0 @@
# triangle_group_limit_set - visualizer for triangle group limit sets in the projective plane
This program visualizes the fractal limit set of triangle groups in SL(3,R) in the projective plane. For the mathematical details, see [our paper].
## Installation & Prerequisites ##
This program is written C and only runs on Linux with the X window system (because that's what I happen to be using). It uses [Cairo] for drawing and X11 for input. We also require the [Gnu Scientific Library (GSL)] including CBLAS.
To build, just make sure GCC as well as these libraries are installed, and run
make
Then we can run it for example like this:
./limit_set 5 5 5 1 1 1 1.0
The arguments are 6 integers (p1, p2, p3, q1, q2, q3) and a floating point numbers t. Together they discribe the triangle group under consideration (t can be changed with the arrow keys).
## Key bindings
Drag the mouse to move the image, drag with Shift pressed to rotate.
| Key | function |
|----------|-----------------------------------------------------------------------------------------------------|
| PageUp | increase t by 2% |
| PageDown | decrease t by 2% |
| Right | increase t by 0.2% |
| Left | decreaes t by 0.2% |
| Up | increase t by 0.002% |
| Down | decrease t by 0.002% |
| Space | (reset t to original value) |
| R | cycle through affine charts |
| l | show limit curve |
| L | show limit curve as line or as points |
| f | generate limit curve using attracting / repelling fixed points of conjugates of the Coxeter element |
| d | show dual limit curve |
| r | show fixed points and lines of generating reflections |
| a | show fixed points and lines of Coxeter elements |
| c | show the orbit of an arbitrary point (point can be changed with Shift+Click) |
| b | show certain conics touching the limit curve |
| B | show the "boxes" used to approximate the limit curve |
| p | save a screenshot of the current image (in PDF format) |
| t | toggle info text |
| i | (print some info to terminal?) |
| x | (show rotated reflectors?) |
| M | (make a movie) |
[Cairo]: https://cairographics.org/
[Gnu Scientific Library (GSL)]: https://www.gnu.org/software/gsl/
[our paper]: https://arxiv.org/abs/2106.11349

604
draw.c
View File

@ -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;
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;
@ -96,10 +97,20 @@ void drawPoint(DrawingContext *ctx, point_t p)
cairo_t *C = ctx->cairo;
cairo_save(C);
cairo_arc(C, p.x, p.y, 4.0/ctx->dim->scalefactor, 0, 2*M_PI);
cairo_arc(C, p.x, p.y, 3.0/ctx->dim->scalefactor, 0, 2*M_PI);
cairo_close_path(C);
cairo_fill(C);
cairo_restore(C);
/*
cairo_save(C);
cairo_move_to(C, p.x, p.y);
cairo_close_path(C);
cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);
cairo_set_line_width(C, 10.0/ctx->dim->scalefactor);
cairo_stroke(C);
cairo_restore(C);
*/
}
void drawSegment2d(DrawingContext *ctx, point_t a, point_t b)
@ -190,6 +201,15 @@ void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line,
// t = parameter on segment of intersection with line, s(t) = a + (b-a)*t
tline = -(a_.x*x[0] + a_.y*x[1] + x[2])/((b_.x - a_.x)*x[0] + (b_.y - a_.y)*x[1]);
/*
printf("tline: %f\n", tline);
point_t s;
s.x = a_.x - (b_.x-a_.x)*tline;
s.y = a_.y - (b_.y-a_.y)*tline;
drawPoint(ctx, s);
*/
if((tline < 0 || tline > 1) != contains) {
drawSegment2d(ctx, a_, b_);
} else { // need to draw complementary segment
@ -197,7 +217,6 @@ void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line,
m.x = ctx->dim->center_x;
m.y = ctx->dim->center_y;
r = ctx->dim->radius;
// equation is coeff2 t^2 + 2 coeff1 t + coeff0 = 0
coeff0 = (a_.x - m.x)*(a_.x - m.x) + (a_.y - m.y)*(a_.y - m.y) - r*r;
coeff1 = (a_.x - m.x)*(b_.x - a_.x) + (a_.y - m.y)*(b_.y - a_.y);
@ -218,6 +237,7 @@ void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line,
}
// level 3: boxes and polygons
void drawPolygon(DrawingContext *ctx, int segments, int sides, ...)
{
va_list args;
@ -308,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]);
@ -350,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);
@ -368,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]);
@ -413,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]);
@ -555,11 +575,10 @@ 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, "abacacbabc", p[3]);
// fixedPoints(ctx, "abcabcabcabcab", p[3]);
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]);
@ -633,6 +652,12 @@ void drawCurvedBox(DrawingContext *ctx, int base, const char *conj, int style)
fixedPoints(ctx, word, p[10]);
// conjugate_word("bca b abc b acb", modifier, conj, word);
// fixedPoints(ctx, word, p[6]);
// conjugate_word("bca baca cab acab acb", modifier, conj, word);
// fixedPoints(ctx, word, p[7]);
LOOP(j) l[0][j] = cross(p[0][(3-j)%3], p[0][(4-j)%3]);
LOOP(j) l[1][j] = cross(p[1][(3-j)%3], p[1][(4-j)%3]);
LOOP(j) l[10][j] = cross(p[10][(3-j)%3], p[10][(4-j)%3]);
@ -679,68 +704,233 @@ 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) {
drawRotationOrbit(ctx, "bcabcb", p[1][0]); // bcC
drawRotationOrbit(ctx, "abcabcba", p[0][0]); // abcC
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, 0, 0, 0);
fixedPoints(ctx, "abababcbaba", p[3]);
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, 1);
fixedPoints(ctx, "bc bca cb", p[3]);
/*
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));
}
/*
while((current = queue_get(&queue)) != -1) {
cur = &ctx->group[current];
if(cur->visited > 4)
continue;
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;
}
}
*/
// 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);
/*
cairo_set_source_rgb(C, 0, 1, 1);
fixedPoints(ctx, "abc", p[0]);
drawRotationOrbit(ctx, "bc", p[0][0]); // (r^{-1}a)^2 C / cbr
fixedPoints(ctx, "bca", p[0]);
drawRotationOrbit(ctx, "abca", p[0][0]); // (r^{-1}a)^2 C / cbr
cairo_set_source_rgb(C, 1, 1, 0);
fixedPoints(ctx, "abc", p[0]);
drawRotationOrbit(ctx, "cbcabc", p[0][0]); // (r^{-1}a)^4 C / cbac r-
fixedPoints(ctx, "bca", p[0]);
drawRotationOrbit(ctx, "acbcabca", p[0][0]); // (r^{-1}a)^4 C / cbac r-
cairo_set_source_rgb(C, 1, 0, 1);
fixedPoints(ctx, "abc", p[0]);
drawRotationOrbit(ctx, "cbacabcabc", p[0][0]); // (r^{-1}a)^6 C / cbacba
fixedPoints(ctx, "bca", p[0]);
drawRotationOrbit(ctx, "acbacabcabca", p[0][0]); // (r^{-1}a)^6 C / cbacba
*/
/*
cairo_set_source_rgb(C, 1, 0, 0);
fixedPoints(ctx, "bca", p[3]);
drawRotationOrbit(ctx, "bcabcb", p[3][0]); // cb c a bc
fixedPoints(ctx, "abc", p[3]);
drawRotationOrbit(ctx, "abcabcba", p[3][0]); // cb c a bc
cairo_set_source_rgb(C, 0, 1, 0);
fixedPoints(ctx, "ababcba", p[3]);
drawRotationOrbit(ctx, "ababcababa", p[3][0]); // cb c a bc
fixedPoints(ctx, "abc", p[3]);
drawRotationOrbit(ctx, "abcabcabacba", p[3][0]); // cb c a bc
*/
// drawRotationOrbit(ctx, "bc", p[3][0]);
cairo_set_source_rgb(C, 0, 0, 0);
drawCurvedBox(ctx, 'A', "", 2);
cairo_set_source_rgb(C, 0, 0, 1);
// drawCurvedBox(ctx, 'C', "ab", 2);
// drawCurvedBox(ctx, 'C', "ab", 2);
// drawCurvedBox(ctx, 'B', "bca", 2);
// drawCurvedBox(ctx, 'B', "ba", 2);
// drawCurvedBox(ctx, 'B', "bca", 2);
// drawCurvedBox(ctx, 'B', "abca", 2);
// drawCurvedBox(ctx, 'C', "abab", 2);
// drawCurvedBox(ctx, 'C', "ababab", 2);
if(ctx->mode >= 3 && ctx->mode != 4) {
cairo_set_source_rgb(C, 0, 0, 0);
drawCurvedBox(ctx, 'B', "", 2);
@ -756,11 +946,15 @@ void drawBoxes2(DrawingContext *ctx)
drawCurvedBox(ctx, 'A', "bcbcbcbc", 2);
drawCurvedBox(ctx, 'A', "bcbcbcbcbc", 2);
drawCurvedBox(ctx, 'A', "bcbcbcbcbcbc", 2);
// drawCurvedBox(ctx, 'A', "b", 2);
// drawCurvedBox(ctx, 'A', "bcb", 2);
// drawCurvedBox(ctx, 'A', "bcbcb", 2);
// drawCurvedBox(ctx, 'A', "bcbcbcb", 2);
// drawCurvedBox(ctx, 'A', "bcbcbcbcb", 2);
}
if(ctx->mode == 6) {
cairo_set_source_rgb(C, 0, 0.8, 0.5);
drawCurvedBox(ctx, 'C', "ababab", 2);
drawCurvedBox(ctx, 'C', "abababab", 2);
drawCurvedBox(ctx, 'C', "ababababab", 2);
}
@ -769,45 +963,291 @@ void drawBoxes2(DrawingContext *ctx)
cairo_set_source_rgb(C, 0, 0.8, 0.5);
drawCurvedBox(ctx, 'C', "ab", 2);
drawCurvedBox(ctx, 'C', "abab", 2);
drawCurvedBox(ctx, 'C', "ababab", 2);
}
if(ctx->mode >= 8) {
cairo_set_source_rgb(C, 0, 0.8, 0.5);
drawCurvedBox(ctx, 'C', "b", 2);
drawCurvedBox(ctx, 'C', "bab", 2);
drawCurvedBox(ctx, 'C', "babab", 2);
}
if(ctx->mode >= 9) {
cairo_set_source_rgb(C, 0.8, 0, 0.5);
drawCurvedBox(ctx, 'B', "abca", 2);
drawCurvedBox(ctx, 'B', "abcaca", 2);
drawCurvedBox(ctx, 'B', "abcacaca", 2);
drawCurvedBox(ctx, 'B', "aba", 2);
drawCurvedBox(ctx, 'B', "abaca", 2);
drawCurvedBox(ctx, 'B', "abacaca", 2);
drawCurvedBox(ctx, 'B', "ababca", 2);
drawCurvedBox(ctx, 'B', "ababcaca", 2);
drawCurvedBox(ctx, 'B', "ababcacaca", 2);
drawCurvedBox(ctx, 'B', "ababa", 2);
drawCurvedBox(ctx, 'B', "ababaca", 2);
drawCurvedBox(ctx, 'B', "ababacaca", 2);
drawCurvedBox(ctx, 'B', "abababca", 2);
drawCurvedBox(ctx, 'B', "abababcaca", 2);
drawCurvedBox(ctx, 'B', "abababcacaca", 2);
drawCurvedBox(ctx, 'B', "abababa", 2);
drawCurvedBox(ctx, 'B', "abababaca", 2);
drawCurvedBox(ctx, 'B', "abababacaca", 2);
drawCurvedBox(ctx, 'B', "bca", 2);
drawCurvedBox(ctx, 'B', "bcaca", 2);
drawCurvedBox(ctx, 'B', "bcacaca", 2);
drawCurvedBox(ctx, 'B', "ba", 2);
drawCurvedBox(ctx, 'B', "baca", 2);
drawCurvedBox(ctx, 'B', "bacaca", 2);
drawCurvedBox(ctx, 'B', "babca", 2);
drawCurvedBox(ctx, 'B', "babcaca", 2);
drawCurvedBox(ctx, 'B', "babcacaca", 2);
drawCurvedBox(ctx, 'B', "baba", 2);
drawCurvedBox(ctx, 'B', "babaca", 2);
drawCurvedBox(ctx, 'B', "babacaca", 2);
drawCurvedBox(ctx, 'B', "bababca", 2);
drawCurvedBox(ctx, 'B', "bababcaca", 2);
drawCurvedBox(ctx, 'B', "bababcacaca", 2);
drawCurvedBox(ctx, 'B', "bababa", 2);
drawCurvedBox(ctx, 'B', "bababaca", 2);
drawCurvedBox(ctx, 'B', "bababacaca", 2);
}
// drawCurvedBox(ctx, 'B', "abab", 2);
// drawCurvedBox(ctx, 'B', "ababab", 2);
// drawCurvedBox(ctx, 'C', "b", 1);
// drawCurvedBox(ctx, 'C', "ababab");
// drawCurvedBox(ctx, 'C', "abababab");
// drawCurvedBox(ctx, 'C', "b");
// drawCurvedBox(ctx, 'C', "bab");
// drawCurvedBox(ctx, 'C', "babab");
// drawCurvedBox(ctx, 'C', "abababab");
// drawCurvedBox(ctx, 'C', "ababababab");
// drawCurvedBox(ctx, 'C', "abababababab");
cairo_set_source_rgb(C, 0, 0, 1);
// drawCurvedBox(ctx, 'B', "abca", 1);
// drawCurvedBox(ctx, 'B', "aba", 1);
// drawCurvedBox(ctx, 'B', "bca", 1);
// drawCurvedBox(ctx, 'B', "ba", 1);
// drawCurvedBox(ctx, 'B', "abaca");
// drawCurvedBox(ctx, 'B', "abcaca");
// drawCurvedBox(ctx, 'B', "abaca");
// drawCurvedBox(ctx, 'B', "aba");
// drawCurvedBox(ctx, 'B', "ababca");
// drawCurvedBox(ctx, 'B', "ababa");
// drawCurvedBox(ctx, 'B', "abababca");
// drawCurvedBox(ctx, 'B', "abababa");
// drawCurvedBox(ctx, 'B', "bca");
// drawCurvedBox(ctx, 'B', "ba");
// drawCurvedBox(ctx, 'B', "babca");
// drawCurvedBox(ctx, 'B', "baba");
// drawCurvedBox(ctx, 'B', "bababca");
// drawCurvedBox(ctx, 'B', "bababa");
cairo_set_source_rgb(C, 1, 0, 1);
// drawCurvedBox(ctx, 'A', "abcac");
// drawCurvedBox(ctx, 'A', "bcac");
// drawCurvedBox(ctx, 'A', "abcacbc");
// drawCurvedBox(ctx, 'A', "abcac");
// drawCurvedBox(ctx, 'A', "abcacbc");
// drawCurvedBox(ctx, 'A', "ababc");
// drawCurvedBox(ctx, 'A', "ababcbc");
/* drawCurvedBox(ctx, 'A', "abcac");
drawCurvedBox(ctx, 'A', "ababc");
drawCurvedBox(ctx, 'A', "abac");
drawCurvedBox(ctx, 'A', "ababcabc");
drawCurvedBox(ctx, 'A', "ababcac");
drawCurvedBox(ctx, 'A', "abababc");
drawCurvedBox(ctx, 'A', "ababac");
drawCurvedBox(ctx, 'A', "abababcabc");
drawCurvedBox(ctx, 'A', "abababcac");
drawCurvedBox(ctx, 'A', "ababababc");
drawCurvedBox(ctx, 'A', "abababac");
drawCurvedBox(ctx, 'A', "bcabc");
drawCurvedBox(ctx, 'A', "bcac");
drawCurvedBox(ctx, 'A', "babc");
drawCurvedBox(ctx, 'A', "bac");
drawCurvedBox(ctx, 'A', "babcabc");
drawCurvedBox(ctx, 'A', "babcac");
drawCurvedBox(ctx, 'A', "bababc");
drawCurvedBox(ctx, 'A', "babac");
drawCurvedBox(ctx, 'A', "bababcabc");
drawCurvedBox(ctx, 'A', "bababcac");
drawCurvedBox(ctx, 'A', "babababc");
drawCurvedBox(ctx, 'A', "bababac");
*/
// drawCurvedBox(ctx, 'A', "ababc");
cairo_set_source_rgb(C, 1, 0.5, 0);
// drawCurvedBox(ctx, 'C', "abcabcab");
// drawCurvedBox(ctx, 'B', "bca");
// drawCurvedBox(ctx, 'B', "abcacaca");
// drawCurvedBox(ctx, 'B', "abaca");
// drawCurvedBox(ctx, 'B', "aba");
// drawCurvedBox(ctx, 'B', "aba");
// drawCurvedBox(ctx, 'B', "abcaca");
// drawCurvedBox(ctx, 'B', "aba");
// drawCurvedBox(ctx, 'B', "abaca");
cairo_set_source_rgb(C, 0, 0, 1);
// drawCurvedBox(ctx, 'A', "abcabc");
// drawCurvedBox(ctx, 'A', "abcabcbc");
// drawCurvedBox(ctx, 'A', "abcacbc");
// drawCurvedBox(ctx, 'A', "abcac");
// drawCurvedBox(ctx, 'A', "bcabc");
// drawCurvedBox(ctx, 'A', "bcabcbc");
// drawCurvedBox(ctx, 'A', "bcacbc");
// drawCurvedBox(ctx, 'A', "bcac");
// drawCurvedBox(ctx, 'A', "abcac");
// drawCurvedBox(ctx, 'A', "abcac");
// drawCurvedBox(ctx, 'A', "abcacbc");
// drawCurvedBox(ctx, 'B', "abacaca");
// drawCurvedBox(ctx, 'B', "abacacaca");
cairo_set_source_rgb(C, 0, 1, 1);
fixedPoints(ctx, "ababcba", p[0]);
// drawRotationOrbit(ctx, "abcaba", p[0][0]);
/*
int p = 9;
vector_t fp[3][3],neutral_line[3],reflection_line[p],star[2*p],outer[2];
vector_t rotation_line = {0,0,1};
cairo_set_line_width(C, 1.5/ctx->dim->scalefactor);
cairo_set_source_rgb(C, 1, 0, 0);
multiply(gen[0], gen[1], rot);
fixedPoints(ctx, "abc", fp[0]);
fixedPoints(ctx, "bca", fp[1]);
fixedPoints(ctx, "cab", fp[2]);
LOOP(i) neutral_line[i] = cross(fp[i][0], fp[i][2]);
LOOP(j) reflection_line[0].x[j] = gsl_matrix_get(ctx->cartan, 0, j);
star[0] = cross(neutral_line[0],reflection_line[0]);
star[p] = cross(neutral_line[2],reflection_line[0]);
for(int j = 1; j < p; j++) {
reflection_line[j] = apply_transpose(rot, reflection_line[j-1]);
star[j] = apply(rot, star[j-1]);
star[j+p] = apply(rot, star[j+p-1]);
}
outer[0] = cross(neutral_line[0],reflection_line[5]);
outer[1] = cross(neutral_line[2],reflection_line[5]);
for(int j = 0; j < 8; j++) {
drawVector(ctx, star[j]);
drawVector(ctx, star[p+j]);
}
for(int j = 0; j < 7; j++) {
// if(j == 3)
// continue;
drawSegment(ctx, star[j%p], star[(j+1)%p+p]);
drawSegment(ctx, star[(j+1)%p+p], star[(j+1)%p]);
drawSegment(ctx, star[(j+1)%p], star[j%p+p]);
drawSegment(ctx, star[j%p+p], star[j%p]);
}
cairo_set_source_rgb(C, 0, 0, 1);
drawVector(ctx,outer[0]);
drawVector(ctx,outer[1]);
// drawCovector(ctx, neutral_line[0]);
// drawCovector(ctx, neutral_line[2]);
drawSegment(ctx, star[0], star[p]);
drawSegment(ctx, star[p], outer[1]);
drawSegment(ctx, outer[1], outer[0]);
drawSegment(ctx, outer[0], star[0]);
cairo_set_source_rgb(C, 0, 0, 0);
drawCovector(ctx, rotation_line);
*/
// for(int j = 0; j < 5; j++)
// drawCovector(ctx, reflection_line[j]);
// intersect attracting line with neutral line of the other element
// for(int j = 0; j < 2; j++)
// i[j] = cross(cross(p[j%2][0],p[j%2][1]),cross(p[(j+1)%2][0],p[(j+1)%2][2]));
// drawPolygon(ctx, 1, 4, p[0][0], i[0], p[1][0], i[1]);
// abc, ababcba, abababcbaba, ..., cab
// bca, acaba, abacababa, ..., babcb
/*
vector_t v[4][3];
vector_t i[4];
fixedPoints(ctx, "abc", v[0]);
fixedPoints(ctx, "bca", v[1]);
fixedPoints(ctx, "cab", v[2]);
fixedPoints(ctx, "acaba", v[3]);
i[0] = cross(cross(v[0][1],v[0][2]),cross(v[2][0],v[2][2]));
i[1] = cross(cross(v[1][1],v[1][2]),cross(v[3][0],v[3][2]));
i[2] = cross(cross(v[2][0],v[2][2]),cross(v[3][0],v[3][2]));
i[3] = cross(cross(v[0][0],v[0][2]),cross(v[1][0],v[1][2]));
// cairo_set_source_rgb(C, 0, 0, 1);
// drawPolygon(ctx, 1, 6, v[2][2], v[1][1], v[0][1], v[3][2], v[3][1], v[2][1]);
// drawPolygon(ctx, 1, 6, v[1][2], i[1], i[2], i[0], v[0][2], i[3]);
cairo_set_line_width(C, 1.5/ctx->dim->scalefactor);
cairo_set_source_rgb(C, 1, 0, 0);
// drawBox(ctx, "bca", "abc");
drawCurvedBox(ctx, 'A', "");
*/
/*
cairo_set_source_rgb(C, 1, 0.5, 0);
drawBox(ctx, "ab abc ba", "ab cab ba");
drawBox(ctx, "abab abc baba", "abab cab baba");
drawBox(ctx, "ababab abc bababa", "ababab cab bababa");
drawBox(ctx, "abababab abc babababa", "abababab cab babababa");
// drawBox(ctx, "ababababab abc bababababa", "ababababab cab bababababa");
drawBox(ctx, "b abc b", "b cab b");
drawBox(ctx, "bab abc bab", "bab cab bab");
drawBox(ctx, "babab abc babab", "babab cab babab");
drawBox(ctx, "bababab abc bababab", "bababab cab bababab");
// drawBox(ctx, "babababab abc babababab", "babababab cab babababab");
*/
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;
@ -815,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);
@ -836,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)
@ -846,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]);
@ -870,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]);
@ -907,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);
@ -919,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);
}
@ -974,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);
@ -991,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);
@ -1002,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]);
@ -1016,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);
@ -1056,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);
}
@ -1065,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);
}
@ -1111,6 +1578,7 @@ void draw(DrawingContext *ctx)
if(ctx->show_marking)
{
cairo_set_source_rgb(C, 0, 0, 1);
drawPoint(ctx, ctx->marking);
}
if(ctx->show_coxeter_orbit)

107
limit_set.c
View File

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

2
linalg.h
View File

@ -13,7 +13,7 @@
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
#define MAX_TEMP_MATRICES 1200000
#define MAX_TEMP_MATRICES 600000
#define MAX_TEMP_VECTORS 100
typedef struct _workspace {

136
main.c
View File

@ -11,6 +11,7 @@
#include "linalg.h"
#define TOGGLE(a) do { (a) = !(a); } while(0)
#define SIGN(x) ((x) > 0 ? 1.0 : -1.0)
DrawingContext *screen_context;
@ -18,6 +19,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 +30,18 @@ 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;
if(argc > 12) {
ctx->movie_filename = argv[9];
ctx->movie_parameter_duration = atof(argv[10]);
ctx->movie_parameter2_duration = atof(argv[11]);
ctx->movie_n_frames = atoi(argv[12]);
} else {
ctx->movie_n_frames = 0;
}
// ctx->parameter = 2.77;
// ctx->parameter = 0.1;
ctx->show_boxes = 0;
@ -39,10 +53,10 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[])
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 = 1;
ctx->limit_with_lines = 0;
ctx->use_repelling = 0;
ctx->show_marking = 1;
ctx->show_marking = 0;
ctx->marking.x = -0.73679;
ctx->marking.y = -0.01873;
ctx->show_coxeter_orbit = 0;
@ -50,8 +64,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 +84,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)
@ -144,8 +167,9 @@ void updateMatrices(DrawingContext *ctx)
cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);
gsl_matrix *tmp = getTempMatrix(ctx->ws);
int nmodes = 5;
if(ctx->use_rotation_basis % 5 == 0) {
if(ctx->use_rotation_basis % nmodes == 0) {
gsl_matrix_set(tmp, 0, 0, 0.0);
gsl_matrix_set(tmp, 0, 1, sqrt(3.0)/2.0);
gsl_matrix_set(tmp, 0, 2, -sqrt(3.0)/2.0);
@ -156,13 +180,23 @@ void updateMatrices(DrawingContext *ctx)
gsl_matrix_set(tmp, 2, 1, 1.0);
gsl_matrix_set(tmp, 2, 2, 1.0);
gsl_matrix_memcpy(ctx->cob, tmp);
} else if(ctx->use_rotation_basis % 5 == 1) {
} else if(ctx->use_rotation_basis % nmodes == 1) {
gsl_matrix_set(tmp, 0, 0, 1.0);
gsl_matrix_set(tmp, 0, 1, -1.0);
gsl_matrix_set(tmp, 0, 2, 0.0);
gsl_matrix_set(tmp, 1, 0, 1.0);
gsl_matrix_set(tmp, 1, 1, 1.0);
gsl_matrix_set(tmp, 1, 2, 0.0);
gsl_matrix_set(tmp, 2, 0, 0.0);
gsl_matrix_set(tmp, 2, 1, 0.0);
gsl_matrix_set(tmp, 2, 2, 1.0);
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");
multiply_left(tmp, ctx->cob, ctx->ws);
} else if(ctx->use_rotation_basis % nmodes == 2) {
computeRotationMatrixFrame(ctx, tmp, "C");
invert(tmp, ctx->cob, ctx->ws);
} else if(ctx->use_rotation_basis % 5 == 3) {
computeBoxTransform(ctx, "bca", "abc", ctx->cob);
} else if(ctx->use_rotation_basis % nmodes == 3) {
computeBoxTransform(ctx, "acb", "cba", ctx->cob);
// computeBoxTransform(ctx, "cab", "bca", ctx->cob);
// computeBoxTransform(ctx, "acb", "cba", ctx->cob);
} else {
@ -248,22 +282,34 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
switch(key) {
case XK_Down:
if(ev->xkey.state & ShiftMask)
screen_context->parameter /= exp(0.00005);
else
screen_context->parameter /= exp(0.002);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Up:
if(ev->xkey.state & ShiftMask)
screen_context->parameter *= exp(0.00005);
else
screen_context->parameter *= exp(0.002);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Left:
screen_context->parameter /= exp(0.00002);
if(ev->xkey.state & ShiftMask)
screen_context->parameter2 /= exp(0.00005);
else
screen_context->parameter2 /= exp(0.002);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Right:
screen_context->parameter *= exp(0.00002);
if(ev->xkey.state & ShiftMask)
screen_context->parameter2 *= exp(0.00005);
else
screen_context->parameter2 *= exp(0.002);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
@ -332,52 +378,22 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
print(screen_context);
break;
case 'M':
/*
screen_context->limit_with_lines = 0;
double parameter_start = screen_context->parameter;
for(int i = 0; i <= 1300; i++) {
if(i < 400)
screen_context->parameter = exp(log(parameter_start)+0.002*i);
else if(i < 500)
screen_context->parameter = exp(log(parameter_start)+0.002*400);
else
screen_context->parameter = exp(log(parameter_start)+0.002*(900-i));
double parameter2_start = screen_context->parameter2;
for(int i = 0; i <= screen_context->movie_n_frames; i++) {
screen_context->parameter = SIGN(parameter_start)*exp(log(fabs(parameter_start)) +
i*screen_context->movie_parameter_duration/screen_context->movie_n_frames);
screen_context->parameter2 = SIGN(parameter2_start)*exp(log(fabs(parameter2_start)) +
i*screen_context->movie_parameter2_duration/screen_context->movie_n_frames);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
draw(screen_context);
sprintf(filename, "movie3/test%03d.png", i);
sprintf(filename, "output/%s%03d.png", screen_context->movie_filename, i);
cairo_surface_write_to_png(info->buffer_surface, filename);
printf("Finished drawing %s\n", filename);
}
*/
/*
screen_context->limit_with_lines = 0;
double parameter_start = screen_context->parameter;
for(int i = 0; i <= 1300; i++) {
if(i < 400)
screen_context->parameter = exp(0.003*i);
else if(i < 500)
screen_context->parameter = exp(0.003*400);
else
screen_context->parameter = exp(0.003*(900-i));
updateMatrices(screen_context);
computeLimitCurve(screen_context);
draw(screen_context);
sprintf(filename, "movie5/test%03d.png", i);
cairo_surface_write_to_png(info->buffer_surface, filename);
printf("Finished drawing %s\n", filename);
}
*/
screen_context->limit_with_lines = 0;
for(int i = 0; i <= 1000; i++) {
screen_context->parameter = exp(log(3.0)*((double)i/500-1));
updateMatrices(screen_context);
computeLimitCurve(screen_context);
draw(screen_context);
sprintf(filename, "output/movie8/test%04d.png", i);
// cairo_surface_write_to_png(info->buffer_surface, filename);
printf("Finished drawing %s\n", filename);
}
case 'f':
TOGGLE(screen_context->use_repelling);
computeLimitCurve(screen_context);

17
main.h
View File

@ -12,6 +12,11 @@
#define LOOP(i) for(int i = 0; i < 3; i++)
#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,13 @@ typedef struct {
// a priori parameter
int p[3],k[3];
double parameter;
double parameter2;
char *movie_filename;
double movie_parameter_duration;
double movie_parameter2_duration;
int movie_n_frames;
int n_group_elements;
int n_group_elements_combinatorial;
int show_boxes;
int show_boxes2;
int show_attractors;
@ -66,7 +77,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 +108,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)

1
triangle.h
View File

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