diff --git a/draw.c b/draw.c index 169c74e..6092222 100644 --- a/draw.c +++ b/draw.c @@ -29,6 +29,43 @@ double scal(vector_t a, vector_t b) return result; } +// degenerate conic whose zero set is the union of two lines +// Q = a^T b + b^T a +void degenerate_conic(vector_t a, vector_t b, gsl_matrix *conic) +{ + LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, a.x[i]*b.x[j] + a.x[j]*b.x[i]); +} + +double conic_value(gsl_matrix *conic, vector_t v) +{ + double result = 0.0; + LOOP(i) LOOP(j) result += gsl_matrix_get(conic, i, j) * v.x[i] * v.x[j]; + return result; +} + +// finds conic in the pencil through (l1,l2) and (L1,L2) going through p +void conic_from_lines(workspace_t *ws, vector_t l1, vector_t l2, vector_t L1, vector_t L2, vector_t p, gsl_matrix *conic) +{ + gsl_matrix *degenerate1 = gsl_matrix_alloc(3, 3); + gsl_matrix *degenerate2 = gsl_matrix_alloc(3, 3); + gsl_matrix *tmp = gsl_matrix_alloc(3, 3); + + degenerate_conic(l1, l2, degenerate1); + degenerate_conic(L1, L2, degenerate2); + + double t = - conic_value(degenerate1, p) / conic_value(degenerate2, p); + LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, + gsl_matrix_get(degenerate1, i, j) + + gsl_matrix_get(degenerate2, i, j)*t); + int positives = diagonalize_symmetric_form(conic, tmp, ws); + if(positives == 2) + LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; + + gsl_matrix_free(degenerate1); + gsl_matrix_free(degenerate2); + gsl_matrix_free(tmp); +} + vector_t apply(gsl_matrix *m, vector_t x) { vector_t out; @@ -49,7 +86,7 @@ vector_t apply_transpose(gsl_matrix *m, vector_t x) return out; } -vector_t apply_inverse(gsl_matrix *m, vector_t in) +vector_t apply_pseudoinverse_transpose(gsl_matrix *m, vector_t in) { vector_t out; double cofactor; @@ -64,6 +101,80 @@ vector_t apply_inverse(gsl_matrix *m, vector_t in) return out; } +vector_t apply_pseudoinverse(gsl_matrix *m, vector_t in) +{ + vector_t out; + double cofactor; + + LOOP(i) out.x[i] = 0.0; + LOOP(i) LOOP(j) { + cofactor = gsl_matrix_get(m, (j+1)%3, (i+1)%3) * gsl_matrix_get(m, (j+2)%3, (i+2)%3) + - gsl_matrix_get(m, (j+1)%3, (i+2)%3) * gsl_matrix_get(m, (j+2)%3, (i+1)%3); + out.x[i] += cofactor * in.x[j]; + } + + return out; +} + +int intersect_line_and_conic(DrawingContext *ctx, vector_t line, gsl_matrix *conic, vector_t *intersection1, vector_t *intersection2) +{ + gsl_matrix *frame = getTempMatrix(ctx->ws); + + int pos = diagonalize_symmetric_form(conic, frame, ctx->ws); + + vector_t a = apply_pseudoinverse_transpose(frame, line); + + double disc = a.x[0]*a.x[0] + a.x[1]*a.x[1] - a.x[2]*a.x[2]; + if(disc < 0) + { + return 0; + } + + vector_t intersection1_in_frame; + intersection1_in_frame.x[0] = -a.x[0]*a.x[2] + a.x[1]*sqrt(disc); + intersection1_in_frame.x[1] = -a.x[1]*a.x[2] - a.x[0]*sqrt(disc); + intersection1_in_frame.x[2] = a.x[0]*a.x[0] + a.x[1]*a.x[1]; + *intersection1 = apply_pseudoinverse(frame, intersection1_in_frame); + + vector_t intersection2_in_frame; + intersection2_in_frame.x[0] = -a.x[0]*a.x[2] - a.x[1]*sqrt(disc); + intersection2_in_frame.x[1] = -a.x[1]*a.x[2] + a.x[0]*sqrt(disc); + intersection2_in_frame.x[2] = a.x[0]*a.x[0] + a.x[1]*a.x[1]; + *intersection2 = apply_pseudoinverse(frame, intersection2_in_frame); + + releaseTempMatrices(ctx->ws, 1); +} + +// should be three collinear vectors! +double halfCR(vector_t x, vector_t y, vector_t z) +{ + double xy = scal(x,y); + double yz = scal(y,z); + double zx = scal(z,x); + double yy = scal(y,y); + double xx = scal(x,x); + + return fabs((yz*xy-zx*yy)/(yz*xx-zx*xy)); +} + +void getMinMaxHalfCrossRatio(double *limit_curve, int n, vector_t v1, vector_t v2, int max, int *min_index, double* min_cr) +{ + vector_t fp[3]; + vector_t tangent; + double cr; + *min_cr = max ? 0 : INFINITY; + for(int i = 0; i < n; i++) { + LOOP(j) LOOP(k) fp[j].x[k] = limit_curve[12*i+3+3*j+k]; + + tangent = cross(fp[0],fp[1]); + cr = fabs(scal(v2, tangent) / scal(v1, tangent)); + if(!max && cr < *min_cr || max && cr > *min_cr) { + *min_cr = cr; + *min_index = i; + } + } +} + int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) { gsl_matrix *tmp = getTempMatrix(ctx->ws); @@ -389,6 +500,28 @@ void drawRotationOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, vector_t sta releaseTempVectors(ctx->ws, 2); } + +void drawConic(DrawingContext *ctx, gsl_matrix *form) +{ + gsl_matrix *orthogonal_frame = getTempMatrix(ctx->ws); + gsl_matrix *frame = getTempMatrix(ctx->ws); + double eval[3]; + vector_t start; + + diagonalize_symmetric_matrix(form, orthogonal_frame, eval, ctx->ws); + + LOOP(i) LOOP(j) + gsl_matrix_set(frame, j, i, + gsl_matrix_get(orthogonal_frame, i, j)/sqrt(fabs(eval[i]))); + + // find any timelike vector; a simple method is adding the first and third column of frame + LOOP(i) start.x[i] = gsl_matrix_get(frame, i, 0) + gsl_matrix_get(frame, i, 2); + + drawRotationOrbitFrame(ctx, frame, start); + + releaseTempMatrices(ctx->ws, 2); +} + void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) { gsl_matrix *frame = getTempMatrix(ctx->ws); @@ -609,14 +742,18 @@ void drawAttractors(DrawingContext *ctx) for(int i = 0; i < n; i++) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); - for(int i = 0; i < n; i++) LOOP(j) { - cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); - drawVector(ctx, p[i][j]); + for(int i = 0; i < n; i++) { + for(int j = 0; j < 3; j += 1) { + cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); + drawVector(ctx, p[i][j]); + } } - for(int i = 0; i < n; i++) LOOP(j) { - cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); - drawCovector(ctx, l[i][j]); + for(int i = 0; i < n; i++) { + for(int j = 0; j < 3; j += 1) { + cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); + drawCovector(ctx, l[i][j]); + } } } @@ -732,6 +869,10 @@ void drawBoxes(DrawingContext *ctx) { gsl_matrix *rot = getTempMatrix(ctx->ws); gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix *order3 = getTempMatrix(ctx->ws); + gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_matrix *tmp2 = getTempMatrix(ctx->ws); + gsl_matrix *conic = getTempMatrix(ctx->ws); gsl_matrix *frame = getTempMatrix(ctx->ws); cairo_t *C = ctx->cairo; cairo_save(C); @@ -751,180 +892,95 @@ void drawBoxes(DrawingContext *ctx) 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]); - - 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 - } - - 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]); - - /* - 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]); - } - - 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]); - } - - - 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]); - } - */ - - /* - 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, "ab", 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", p[0][2]); +// drawRotationOrbit(ctx, "bc", p[0][2]); +// drawRotationOrbit(ctx, "ca", p[0][2]); -// drawRotationOrbit(ctx, "ab", cross(l[0][0], l[2][1])); +// gsl_matrix_set_zero(order3); +// LOOP(i) gsl_matrix_set(order3, (i+1)%3, i, 1); +// rotation_frame(order3, frame, ctx->ws); +// drawRotationOrbitFrame(ctx, frame, p[0][0]); +// drawRotationOrbitFrame(ctx, frame, p[0][2]); -// drawRotationOrbit(ctx, "abca", p[0][0]); + vector_t line1; + vector_t line2; + double t; + int positives; + + /* + // conic 1 + line1 = cross(p[0][0],p[0][1]); + line2 = cross(p[1][0],p[1][1]); + degenerate_conic(line1, line2, tmp); + line1 = cross(p[0][0], p[1][0]); + line2 = cross(p[0][0], p[1][0]); + degenerate_conic(line1, line2, tmp2); + t = - conic_value(tmp, p[2][0]) / conic_value(tmp2, p[2][0]); + LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, + gsl_matrix_get(tmp, i, j) + + gsl_matrix_get(tmp2, i, j)*t); + positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); + if(positives == 2) + LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; + drawConic(ctx, conic); + */ + + // conic 2 + line1 = cross(p[0][0],p[0][1]); + line2 = cross(p[2][0],p[2][1]); + degenerate_conic(line1, line2, tmp); + line1 = cross(p[0][0], p[2][0]); + line2 = cross(p[0][0], p[2][0]); + degenerate_conic(line1, line2, tmp2); + t = - conic_value(tmp, p[1][0]) / conic_value(tmp2, p[1][0]); + LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, + gsl_matrix_get(tmp, i, j) + + gsl_matrix_get(tmp2, i, j)*t); + positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); + if(positives == 2) + LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; + drawConic(ctx, conic); + + /* + // conic 3 + line1 = cross(p[0][2],p[0][1]); + line2 = cross(p[1][2],p[1][1]); + degenerate_conic(line1, line2, tmp); + line1 = cross(p[0][2], p[1][2]); + line2 = cross(p[0][2], p[1][2]); + degenerate_conic(line1, line2, tmp2); + t = - conic_value(tmp, p[2][2]) / conic_value(tmp2, p[2][2]); + LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, + gsl_matrix_get(tmp, i, j) + + gsl_matrix_get(tmp2, i, j)*t); + positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); + if(positives == 2) + LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; + drawConic(ctx, conic); + */ + + // conic 4 + line1 = cross(p[0][2],p[0][1]); + line2 = cross(p[2][2],p[2][1]); + degenerate_conic(line1, line2, tmp); + line1 = cross(p[0][2], p[2][2]); + line2 = cross(p[0][2], p[2][2]); + degenerate_conic(line1, line2, tmp2); + t = - conic_value(tmp, p[1][2]) / conic_value(tmp2, p[1][2]); + LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, + gsl_matrix_get(tmp, i, j) + + gsl_matrix_get(tmp2, i, j)*t); + positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); + if(positives == 2) + LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; + drawConic(ctx, conic); cairo_restore(C); - releaseTempMatrices(ctx->ws, 5); -} - -void getMinMaxHalfCrossRatio(double *limit_curve, int n, vector_t v1, vector_t v2, int max, int *min_index, double* min_cr) -{ - vector_t fp[3]; - vector_t tangent; - double cr; - *min_cr = max ? 0 : INFINITY; - for(int i = 0; i < n; i++) { - LOOP(j) LOOP(k) fp[j].x[k] = limit_curve[12*i+3+3*j+k]; - tangent = cross(fp[0],fp[1]); - - cr = fabs(scal(v1, tangent) / scal(v2, tangent)); - if(!max && cr < *min_cr || max && cr > *min_cr) { - *min_cr = cr; - *min_index = i; - } - } - + releaseTempMatrices(ctx->ws, 9); } void drawBoxes2(DrawingContext *ctx) @@ -933,104 +989,173 @@ void drawBoxes2(DrawingContext *ctx) gsl_matrix **gen = getTempMatrices(ctx->ws, 3); gsl_vector *marking_globalbasis = getTempVector(ctx->ws); gsl_vector *marking_drawbasis = getTempVector(ctx->ws); + gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_matrix *tmp2 = getTempMatrix(ctx->ws); + gsl_matrix *conic = getTempMatrix(ctx->ws); cairo_t *C = ctx->cairo; cairo_save(C); initializeTriangleGenerators(gen, ctx->cartan); - vector_t p[4][3]; + vector_t p[3][3]; + cairo_set_line_width(C, 1/ctx->dim->scalefactor); + cairo_set_source_rgb(C, 0, 0, 0); + + // find coxeter fixed points and triangle vertices fixedPoints(ctx, "abc", p[0]); fixedPoints(ctx, "bca", p[1]); fixedPoints(ctx, "cab", p[2]); - cairo_set_line_width(C, 1/ctx->dim->scalefactor); - - cairo_set_source_rgb(C, 0, 0, 0); - - vector_t att_line[3]; - vector_t rep_line[3]; - - LOOP(i) att_line[i] = cross(p[i][0], p[i][1]); - LOOP(i) rep_line[i] = cross(p[i][2], p[i][1]); - vector_t vertex[3]; LOOP(i) vertex[i] = cross( cross(p[(i+1)%3][0],p[(i+1)%3][2]), cross(p[(i+2)%3][0],p[(i+2)%3][2])); + vector_t a,b; + double t = ctx->distance_parameter1, s = ctx->distance_parameter2; + + vector_t phiplus, phiminus, eta; + double dvert, param; + + // want to choose a so that d(vertex[0],a) = t*d(vertex[0],vertex[1]) + phiplus = cross(p[2][0], p[2][1]); + phiminus = cross(p[2][2], p[2][1]); + dvert = log( + (scal(phiplus, vertex[0])*scal(phiminus, vertex[1])) / + (scal(phiplus, vertex[1])*scal(phiminus, vertex[0]))); + LOOP(i) eta.x[i] = + scal(phiminus, vertex[0]) / scal(phiplus, vertex[0]) + * exp(t*dvert) * phiplus.x[i] + - phiminus.x[i]; + param = 1/(1-scal(eta,vertex[0])/scal(eta,vertex[1])); + LOOP(i) a.x[i] = param*vertex[0].x[i] + (1-param)*vertex[1].x[i]; + + // want to choose b so that d(vertex[0],b) = s*d(vertex[0],vertex[2]) + phiplus = cross(p[1][0], p[1][1]); + phiminus = cross(p[1][2], p[1][1]); + dvert = log( + (scal(phiplus, vertex[0])*scal(phiminus, vertex[2])) / + (scal(phiplus, vertex[2])*scal(phiminus, vertex[0]))); + LOOP(i) eta.x[i] = + scal(phiminus, vertex[0]) / scal(phiplus, vertex[0]) + * exp(s*dvert) * phiplus.x[i] + - phiminus.x[i]; + param = 1/(1-scal(eta,vertex[0])/scal(eta,vertex[2])); + LOOP(i) b.x[i] = param*vertex[0].x[i] + (1-param)*vertex[2].x[i]; + + // draw and output stuff LOOP(i) drawVector(ctx, vertex[i]); -// LOOP(i) drawCovector(ctx, att_line[i]); + drawVector(ctx, a); + drawVector(ctx, b); + drawCovector(ctx, cross(a,b)); + drawCovector(ctx, cross(b,vertex[0])); + drawCovector(ctx, cross(vertex[0],a)); - double half_cross_ratio[3]; - LOOP(i) half_cross_ratio[i] = scal(vertex[(i+1)%3],att_line[i])/scal(vertex[(i+2)%3],att_line[i]); - double triple_ratio = half_cross_ratio[0]*half_cross_ratio[1]*half_cross_ratio[2]; - LOOP(i) half_cross_ratio[i] = scal(vertex[(i+1)%3],rep_line[i])/scal(vertex[(i+2)%3],rep_line[i]); - double triple_ratio2 = half_cross_ratio[0]*half_cross_ratio[1]*half_cross_ratio[2]; + // find the intersections with the limit set and the asymmetric distances + int boundary_index; + double value; + double perimeter1 = 0.0, perimeter2 = 0.0; + double approx_perimeter1 = 0.0, approx_perimeter2 = 0.0; + vector_t boundary_point; + cairo_set_source_rgb(C, 1, 0, 0); - - - vector_t marking, marking2, marking3; - gsl_vector_set(marking_drawbasis, 0, ctx->marking.x); - gsl_vector_set(marking_drawbasis, 1, ctx->marking.y); - gsl_vector_set(marking_drawbasis, 2, 1); - solve(ctx->cob, marking_drawbasis, marking_globalbasis, ctx->ws); - LOOP(i) marking.x[i] = gsl_vector_get(marking_globalbasis, i); - - gsl_vector_set(marking_drawbasis, 0, ctx->marking2.x); - gsl_vector_set(marking_drawbasis, 1, ctx->marking2.y); - gsl_vector_set(marking_drawbasis, 2, 1); - solve(ctx->cob, marking_drawbasis, marking_globalbasis, ctx->ws); - LOOP(i) marking2.x[i] = gsl_vector_get(marking_globalbasis, i); - - gsl_vector_set(marking_drawbasis, 0, ctx->marking3.x); - gsl_vector_set(marking_drawbasis, 1, ctx->marking3.y); - gsl_vector_set(marking_drawbasis, 2, 1); - solve(ctx->cob, marking_drawbasis, marking_globalbasis, ctx->ws); - LOOP(i) marking3.x[i] = gsl_vector_get(marking_globalbasis, i); - - /* - vector_t line_drawbasis = apply_inverse(ctx->cob, cross(marking, vertex[0])); - double last, current; - int boundary_points[2] = {0,0}; - int current_boundary_point = 0; - last = - ctx->limit_curve[12*(ctx->n_group_elements-1) ]*line_drawbasis.x[0] + - ctx->limit_curve[12*(ctx->n_group_elements-1)+1]*line_drawbasis.x[1] + - line_drawbasis.x[2]; - */ - - - double value[6]; - int boundary_point[6]; getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, - marking, marking2, 1, &boundary_point[0], &value[0]); - getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, - marking2, marking3, 1, &boundary_point[1], &value[1]); - getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, - marking3, marking, 1, &boundary_point[2], &value[2]); - getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, - marking, marking2, 0, &boundary_point[0], &value[3]); - getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, - marking2, marking3, 0, &boundary_point[1], &value[4]); - getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, - marking3, marking, 0, &boundary_point[2], &value[5]); + a, b, 1, &boundary_index, &value); + LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k]; + drawVector(ctx, boundary_point); + perimeter1 += log(value); -/* vector_t fp[3], tangent; - LOOP(j) LOOP(k) fp[j].x[k] = ctx->limit_curve[12*boundary_point[0]+3+3*j+k]; - tangent = cross(fp[0],fp[1]); - drawCovector(ctx, tangent);*/ + getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, + b, vertex[0], 1, &boundary_index, &value); + LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k]; + drawVector(ctx, boundary_point); + perimeter1 += log(value); + approx_perimeter1 += log(value); - snprintf(ctx->extra_text, 1000, "tr1 = %f, tr2 2 = %f, aa = %f", - log(value[0])+log(value[1])+log(value[2]), - -log(value[3])-log(value[4])-log(value[5]), - log(value[0])+log(value[1])+log(value[2])+log(value[3])+log(value[4])+log(value[5])); + getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, + vertex[0], a, 1, &boundary_index, &value); + LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k]; + drawVector(ctx, boundary_point); + perimeter1 += log(value); + approx_perimeter1 += log(value); - drawCovector(ctx, cross(marking, marking2)); - drawCovector(ctx, cross(marking2, marking3)); - drawCovector(ctx, cross(marking3, marking)); + cairo_set_source_rgb(C, 0, 0, 1); + + getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, + b, a, 1, &boundary_index, &value); + LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k]; + drawVector(ctx, boundary_point); + perimeter2 += log(value); + + getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, + vertex[0], b, 1, &boundary_index, &value); + LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k]; + drawVector(ctx, boundary_point); + perimeter2 += log(value); + approx_perimeter2 += log(value); + + getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements, + a, vertex[0], 1, &boundary_index, &value); + LOOP(k) boundary_point.x[k] = ctx->limit_curve[12*boundary_index+3+k]; + drawVector(ctx, boundary_point); + perimeter2 += log(value); + approx_perimeter2 += log(value); + + vector_t conic_intersection[2]; + double hC0, hC1; + + // find the outer conic + cairo_set_source_rgb(C, 0, 0.7, 0); + + conic_from_lines(ctx->ws, + cross(p[0][0], p[0][1]), cross(p[1][0], p[1][1]), + cross(p[0][0], p[1][0]), cross(p[0][0], p[1][0]), + p[2][0], conic); + drawConic(ctx, conic); + + conic_intersection[2]; + intersect_line_and_conic(ctx, cross(a, b), conic, + &conic_intersection[0], &conic_intersection[1]); + + hC0 = halfCR(b, a, conic_intersection[0]); + hC1 = halfCR(b, a, conic_intersection[1]); + if(hC0 > hC1) + drawVector(ctx, conic_intersection[0]); + else + drawVector(ctx, conic_intersection[1]); + approx_perimeter2 += log(hC0 > hC1 ? hC0 : hC1); + + // find the inner conic + cairo_set_source_rgb(C, 0, 0.5, 0.5); + + conic_from_lines(ctx->ws, + cross(p[0][2], p[0][1]), cross(p[1][2], p[1][1]), + cross(p[0][2], p[1][2]), cross(p[0][2], p[1][2]), + p[2][2], conic); + drawConic(ctx, conic); + + conic_intersection[2]; + intersect_line_and_conic(ctx, cross(a, b), conic, + &conic_intersection[0], &conic_intersection[1]); + + hC0 = halfCR(a, b, conic_intersection[0]); + hC1 = halfCR(a, b, conic_intersection[1]); + if(hC0 > hC1) + drawVector(ctx, conic_intersection[0]); + else + drawVector(ctx, conic_intersection[1]); + approx_perimeter1 += log(hC0 > hC1 ? hC0 : hC1); + + // output + snprintf(ctx->extra_text, 1000, "perimeter1 = %f (%f), perimeter2 = %f (%f)", + perimeter1, approx_perimeter1, + perimeter2, approx_perimeter2); + + printf("%f %f %f %f %f %f\n", + t, s, perimeter1, perimeter2, approx_perimeter1, approx_perimeter2); cairo_restore(C); - releaseTempMatrices(ctx->ws, 4); + releaseTempMatrices(ctx->ws, 7); releaseTempVectors(ctx->ws, 2); } @@ -1341,9 +1466,9 @@ void draw(DrawingContext *ctx) if(ctx->show_marking) { cairo_set_source_rgb(C, 0, 0, 1); - drawPoint(ctx, ctx->marking); - drawPoint(ctx, ctx->marking2); - drawPoint(ctx, ctx->marking3); +// drawPoint(ctx, ctx->marking); +// drawPoint(ctx, ctx->marking2); +// drawPoint(ctx, ctx->marking3); } if(ctx->show_coxeter_orbit) diff --git a/linalg.c b/linalg.c index 4327edb..4e28bbc 100644 --- a/linalg.c +++ b/linalg.c @@ -315,7 +315,6 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws) ERROR(r, "gsl_eigen_symmv failed!\n"); gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC); - gsl_matrix_transpose(cob); int positive = 0; @@ -332,6 +331,32 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws) return positive; } +int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, workspace_t *ws) +{ + gsl_matrix *A_ = getTempMatrix(ws); + + gsl_matrix_memcpy(A_, A); + int r = gsl_eigen_symmv (A_, ws->eval_real, cob, ws->work_symmv); + ERROR(r, "gsl_eigen_symmv failed!\n"); + + gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC); + + gsl_matrix_transpose(cob); + + int positive = 0; + + for(int i = 0; i < ws->n; i++) { + if(eval) + eval[i] = gsl_vector_get(ws->eval_real, i); + + if(gsl_vector_get(ws->eval_real, i) > 0) + positive++; + } + + releaseTempMatrices(ws, 1); + return positive; +} + // computes a matrix in SL(3, R) which projectively transforms (e1, e2, e3, e1+e2+e3) to the 4 given vectors void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws) { diff --git a/linalg.h b/linalg.h index f2d9bbe..e788c3e 100644 --- a/linalg.h +++ b/linalg.h @@ -53,6 +53,7 @@ int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws); int real_eigenvalues(gsl_matrix *g, gsl_vector *evec, workspace_t *ws); void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws); int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws); +int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, workspace_t *ws); void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws); void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws); diff --git a/main.c b/main.c index 2c0e9ce..4a3f0df 100644 --- a/main.c +++ b/main.c @@ -49,6 +49,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]) ctx->marking2.y = -0.11873; ctx->marking3.x = -0.73679; ctx->marking3.y = -0.21873; + ctx->distance_parameter1 = 0.95; + ctx->distance_parameter2 = 0.95; ctx->show_coxeter_orbit = 0; ctx->extra_text = malloc(1000*sizeof(char)); memset(ctx->extra_text, 0, 1000*sizeof(char)); @@ -77,17 +79,15 @@ void destroyContext(DrawingContext *ctx) workspace_free(ctx->ws); } -void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type) +void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *word) { gsl_matrix *tmp = getTempMatrix(ctx->ws); gsl_matrix **gen = getTempMatrices(ctx->ws, 3); -// ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n"); - initializeTriangleGenerators(gen, ctx->cartan); gsl_matrix_set_identity(tmp); - for(int i = 0; i < strlen(type); i++) - multiply_right(tmp, gen[type[i]-'a'], ctx->ws); + for(int i = 0; i < strlen(word); i++) + multiply_right(tmp, gen[word[i]-'a'], ctx->ws); rotation_frame(tmp, result, ctx->ws); @@ -268,7 +268,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev) case KeyPress: state = ev->xkey.state & (ShiftMask | LockMask | ControlMask); key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask)); - printf("Key pressed: %ld\n", key); +// printf("Key pressed: %ld\n", key); switch(key) { case XK_Down: @@ -355,6 +355,16 @@ int processEvent(GraphicsInfo *info, XEvent *ev) case 'p': print(screen_context); break; + case 'P': + for(int i = 0; i <= 100; i++) { + for(int j = 0; j <= 100; j++) { + screen_context->distance_parameter1 = 0.02*i; + screen_context->distance_parameter2 = 0.02*j; + draw(screen_context); + } + fflush(stdout); + } + break; case 'M': /* screen_context->limit_with_lines = 0; @@ -468,7 +478,7 @@ int main(int argc, char *argv[]) gettimeofday(¤t_time, 0); end_time = current_time.tv_sec + current_time.tv_usec*1e-6; - printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000); +// printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000); } waitUpdateTimer(info); } diff --git a/main.h b/main.h index 014a763..5889879 100644 --- a/main.h +++ b/main.h @@ -46,6 +46,8 @@ typedef struct { point_t marking; point_t marking2; point_t marking3; + double distance_parameter1; + double distance_parameter2; int show_coxeter_orbit; int use_repelling; gsl_matrix *cartan, *cob;