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

Compare commits

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

2 Commits
015b391cc0 ... 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
10 changed files with 1062 additions and 1170 deletions
Expand all files Collapse all files

11
Makefile
View File

@@ -1,4 +1,4 @@
HEADERS=triangle.h linalg.h queue.h initcairo.h main.h exp_equation.h HEADERS=triangle.h linalg.h queue.h initcairo.h main.h
SPECIAL_OPTIONS=-O0 -g -D_DEBUG SPECIAL_OPTIONS=-O0 -g -D_DEBUG
#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
@@ -11,8 +11,8 @@ OPTIONS=$(GENERAL_OPTIONS) $(CAIRO_OPTIONS) $(SPECIAL_OPTIONS)
all: limit_set all: limit_set
limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o
gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o -lm -lgsl -lcblas -lcairo -lX11 gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o -lm -lgsl -lcblas -lcairo -lX11
linalg.o: linalg.c $(HEADERS) linalg.o: linalg.c $(HEADERS)
gcc $(OPTIONS) -c linalg.c gcc $(OPTIONS) -c linalg.c
@@ -32,8 +32,5 @@ draw.o: draw.c $(HEADERS)
main.o: main.c $(HEADERS) main.o: main.c $(HEADERS)
gcc $(OPTIONS) -c main.c gcc $(OPTIONS) -c main.c
exp_equation.o: exp_equation.c $(HEADERS)
gcc $(OPTIONS) -c exp_equation.c
clean: clean:
rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o exp_equation.o rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o

1086
draw.c
View File

File diff suppressed because it is too large Load Diff

183
exp_equation.c
View File

@@ -1,183 +0,0 @@
#include "main.h"
#include "exp_equation.h"
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#define EPSILON 1e-9
#define LOOP(i) for(int i = 0; i < 3; i++)
struct solve_exp_plus_exp_params
{
double alpha;
double beta;
double x;
};
static double solve_exp_plus_exp_f(double t, void *_params)
{
struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;
// return exp(params->alpha * t) + exp(params->beta * t) - params->x;
return log(exp(params->alpha * t) + exp(params->beta * t)) - log(params->x);
}
static double solve_exp_plus_exp_df(double t, void *_params)
{
struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;
// return params->alpha * exp(params->alpha * t) + params->beta * exp(params->beta * t);
return params->alpha + (params->beta - params->alpha) / (1 + exp((params->alpha-params->beta)*t));
}
static void solve_exp_plus_exp_fdf(double t, void *params, double *f, double *df)
{
*f = solve_exp_plus_exp_f(t, params);
*df = solve_exp_plus_exp_df(t, params);
}
// solve the equation exp(alpha t) + exp(beta t) = x for t
int solve_exp_plus_exp(double alpha, double beta, double x, double *t)
{
if(alpha <= 0 && beta >= 0 || alpha >= 0 && beta <= 0) {
double critical =
pow(-beta/alpha, alpha/(alpha-beta)) +
pow(-beta/alpha, beta/(alpha-beta));
if(x < critical)
return 0;
else if (x < critical + EPSILON) {
t[0] = log(-beta/alpha)/(alpha-beta);
return 1;
}
// Newton this
gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
struct solve_exp_plus_exp_params params;
params.alpha = alpha;
params.beta = beta;
params.x = x;
gsl_function_fdf FDF;
FDF.f = &solve_exp_plus_exp_f;
FDF.df = &solve_exp_plus_exp_df;
FDF.fdf = &solve_exp_plus_exp_fdf;
FDF.params = (void *)&params;
int status;
double root, lastroot;
for(int r = 0; r < 2; r++) {
root = r == 0 ? log(x)/beta : log(x)/alpha;
gsl_root_fdfsolver_set(solver, &FDF, root);
for(int i = 0; i < 100; i++) {
gsl_root_fdfsolver_iterate(solver);
lastroot = root;
root = gsl_root_fdfsolver_root(solver);
status = gsl_root_test_delta(root, lastroot, 0, 1e-9);
if(status == GSL_SUCCESS)
break;
// printf("iteration %d, root %f\n", i, root);
}
t[r] = root;
}
gsl_root_fdfsolver_free(solver);
return 2;
} else {
// Newton with start value 0
gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
struct solve_exp_plus_exp_params params;
params.alpha = alpha;
params.beta = beta;
params.x = x;
gsl_function_fdf FDF;
FDF.f = &solve_exp_plus_exp_f;
FDF.df = &solve_exp_plus_exp_df;
FDF.fdf = &solve_exp_plus_exp_fdf;
FDF.params = (void *)&params;
int status;
double root, lastroot;
root = 0;
gsl_root_fdfsolver_set(solver, &FDF, root);
for(int i = 0; i < 100; i++) {
gsl_root_fdfsolver_iterate(solver);
lastroot = root;
root = gsl_root_fdfsolver_root(solver);
status = gsl_root_test_delta(root, lastroot, 0, 1e-9);
// printf("iteration %d, root %f\n", i, root);
if(status == GSL_SUCCESS)
break;
}
t[0] = root;
gsl_root_fdfsolver_free(solver);
return 1;
}
}
// solve the equation x1 exp(a1 t) + x2 exp(a2 t) + x3 exp(a3 t) = 0
int solve_linear_exp(vector_t a, vector_t x, double *t)
{
if(x.x[0] > 0 && x.x[1] > 0 && x.x[2] > 0 ||
x.x[0] < 0 && x.x[1] < 0 && x.x[2] < 0)
return 0;
// ensure that y[0] < 0 and y[1], y[2] > 0
int j;
vector_t y, b;
for(j = 0; j < 3; j++)
if(x.x[(j+1)%3] * x.x[(j+2)%3] > 0)
break;
LOOP(i) y.x[i] = x.x[(i+j)%3];
if(y.x[0] > 0)
LOOP(i) y.x[i] *= -1;
LOOP(i) b.x[i] = a.x[(i+j)%3];
double T = (log(y.x[1]) - log(y.x[2])) / (b.x[2] - b.x[1]);
double rhs = - y.x[0] *
pow(y.x[1], (b.x[0]-b.x[2])/(b.x[2]-b.x[1])) *
pow(y.x[2], (b.x[0]-b.x[1])/(b.x[1]-b.x[2]));
int n = solve_exp_plus_exp(b.x[1] - b.x[0], b.x[2] - b.x[0], rhs, t);
for(int i = 0; i < n; i++)
t[i] += T;
return n;
}
/*
int main(int argc, char *argv[])
{
// int n = solve_exp_plus_exp(atof(argv[1]), atof(argv[2]), atof(argv[3]), result);
double result[2];
vector_t a, x;
a.x[0] = atof(argv[1]);
a.x[1] = atof(argv[2]);
a.x[2] = atof(argv[3]);
x.x[0] = atof(argv[4]);
x.x[1] = atof(argv[5]);
x.x[2] = atof(argv[6]);
int n = solve_linear_exp(a, x, result);
if(n == 0)
printf("0 results found\n");
else if(n == 1)
printf("1 result found: %.9f\n", result[0]);
else if(n == 2)
printf("2 results found: %.9f and %.9f\n", result[0], result[1]);
else
printf("%d results found\n", n);
return 0;
}
*/

14
exp_equation.h
View File

@@ -1,14 +0,0 @@
#ifndef EXP_EQUATION_H
#define EXP_EQUATION_H
#include "main.h"
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
int solve_exp_plus_exp(double alpha, double beta, double x, double *t);
int solve_linear_exp(vector_t a, vector_t x, double *t);
#endif

110
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; 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) void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
{ {
gsl_matrix_set(cartan, 0, 0, 2); 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); 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]); gsl_matrix *reflection_gen[3];
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); 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) int computeLimitCurve(DrawingContext *ctx)
@@ -38,7 +86,7 @@ int computeLimitCurve(DrawingContext *ctx)
gsl_matrix *coxeter = getTempMatrix(ctx->ws); gsl_matrix *coxeter = getTempMatrix(ctx->ws);
gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws); gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws);
gsl_matrix *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); gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);
groupelement_t *group = ctx->group; groupelement_t *group = ctx->group;
int success = 0; int success = 0;
@@ -50,31 +98,43 @@ int computeLimitCurve(DrawingContext *ctx)
// do first in the Fuchsian positive case to get the angles // 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); 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]); gsl_matrix_set_identity(elements[0]);
for(int i = 1; i < ctx->n_group_elements; i++) for(int i = 1; i < ctx->n_group_elements; i++) {
multiply(gen[group[i].letter], elements[group[i].parent->id], 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); 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); int ev_count_pos = real_eigenvectors(coxeter_pos, coxeter_fixedpoints_pos, ws);
if(ev_count_pos != 3) if(ev_count_pos != 3)
goto error_out; goto error_out;
int n = 0;
for(int i = 0; i < ctx->n_group_elements; i++) { 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); multiply_many(ws, fixedpoints_pos, 3, cob_pos, elements[i], coxeter_fixedpoints_pos);
ctx->limit_curve[12*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, 2, column)/gsl_matrix_get(fixedpoints_pos, 0, column),
gsl_matrix_get(fixedpoints_pos, 1, 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 // now do it again to calculate x and y coordinates
initializeTriangleGenerators(gen, ctx->cartan); initializeTriangleGeneratorsCurrent(gen, ctx);
gsl_matrix_set_identity(elements[0]); gsl_matrix_set_identity(elements[0]);
for(int i = 1; i < ctx->n_group_elements; i++) for(int i = 1; i < ctx->n_group_elements; i++) {
multiply(gen[group[i].letter], elements[group[i].parent->id], 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); int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);
if(ev_count == 1) if(ev_count == 1)
@@ -82,11 +142,17 @@ int computeLimitCurve(DrawingContext *ctx)
if(ev_count == 0) if(ev_count == 0)
goto error_out; goto error_out;
ctx->limit_curve_count = 0;
for(int i = 0; i < ctx->n_group_elements; i++) { 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); multiply_many(ws, fixedpoints, 3, ctx->cob, elements[i], coxeter_fixedpoints);
x = ctx->limit_curve[12*i ] = gsl_matrix_get(fixedpoints, 0, 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[12*i+1] = gsl_matrix_get(fixedpoints, 1, 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) if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10)
{ {
@@ -95,22 +161,16 @@ int computeLimitCurve(DrawingContext *ctx)
fputc('a' + cur->letter, stdout); // bcbcbca, bacbcacab, bc bca cb fputc('a' + cur->letter, stdout); // bcbcbca, bacbcacab, bc bca cb
fputc('\n',stdout); fputc('\n',stdout);
} }
multiply_many(ws, fixedpoints, 2, elements[i], coxeter_fixedpoints);
LOOP(j) LOOP(k) ctx->limit_curve[12*i+3+3*j+k] = gsl_matrix_get(fixedpoints, k, j);
// bca abc acb = abc
} }
qsort(ctx->limit_curve, ctx->n_group_elements, 12*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; success = 1;
error_out: error_out:
releaseTempMatrices(ctx->ws, 11+ctx->n_group_elements); releaseTempMatrices(ctx->ws, 14+ctx->n_group_elements);
return success; return success;
} }

27
linalg.c
View File

@@ -315,6 +315,7 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
ERROR(r, "gsl_eigen_symmv failed!\n"); ERROR(r, "gsl_eigen_symmv failed!\n");
gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC); gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC);
gsl_matrix_transpose(cob); gsl_matrix_transpose(cob);
int positive = 0; int positive = 0;
@@ -331,32 +332,6 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
return positive; 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 // 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) void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws)
{ {

1
linalg.h
View File

@@ -53,7 +53,6 @@ int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
int real_eigenvalues(gsl_matrix *g, gsl_vector *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); 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_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 projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws);
void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws); void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws);

767
main.c
View File

@@ -11,481 +11,474 @@
#include "linalg.h" #include "linalg.h"
#define TOGGLE(a) do { (a) = !(a); } while(0) #define TOGGLE(a) do { (a) = !(a); } while(0)
#define SIGN(x) ((x) > 0 ? 1.0 : -1.0)
DrawingContext *screen_context; DrawingContext *screen_context;
// setup everything except cairo and dim, which will be provided by the graphics system // setup everything except cairo and dim, which will be provided by the graphics system
void setupContext(DrawingContext *ctx, int argc, char *argv[]) void setupContext(DrawingContext *ctx, int argc, char *argv[])
{ {
ctx->n_group_elements = NUM_GROUP_ELEMENTS; ctx->n_group_elements = NUM_GROUP_ELEMENTS;
ctx->p[0] = atoi(argv[1]); ctx->n_group_elements_combinatorial = NUM_GROUP_ELEMENTS_COMBINATORIAL;
ctx->p[1] = atoi(argv[2]); ctx->p[0] = atoi(argv[1]);
ctx->p[2] = atoi(argv[3]); ctx->p[1] = atoi(argv[2]);
ctx->k[0] = atoi(argv[4]); ctx->p[2] = atoi(argv[3]);
ctx->k[1] = atoi(argv[5]); ctx->k[0] = atoi(argv[4]);
ctx->k[2] = atoi(argv[6]); ctx->k[1] = atoi(argv[5]);
if(argc > 7) ctx->k[2] = atoi(argv[6]);
ctx->parameter = atof(argv[7]); if(argc > 7)
else ctx->parameter = atof(argv[7]);
ctx->parameter = 1.0; else
// ctx->parameter = 2.77; ctx->parameter = 1.0;
// ctx->parameter = 0.1; if(argc > 8)
ctx->show_boxes = 0; ctx->parameter2 = atof(argv[8]);
ctx->show_boxes2 = 0; else
ctx->show_attractors = 0; ctx->parameter2 = 1.0;
ctx->show_reflectors = 0; if(argc > 12) {
ctx->show_rotated_reflectors = 0; ctx->movie_filename = argv[9];
ctx->show_limit= 0; ctx->movie_parameter_duration = atof(argv[10]);
ctx->show_dual_limit= 0; ctx->movie_parameter2_duration = atof(argv[11]);
ctx->show_text = 1; ctx->movie_n_frames = atoi(argv[12]);
ctx->mode = 0; } else {
ctx->use_rotation_basis = 0; ctx->movie_n_frames = 0;
ctx->limit_with_lines = 1; }
ctx->use_repelling = 0; // ctx->parameter = 2.77;
ctx->show_marking = 1; // ctx->parameter = 0.1;
ctx->marking.x = -0.73679; ctx->show_boxes = 0;
ctx->marking.y = -0.01873; ctx->show_boxes2 = 0;
ctx->marking2.x = -0.73679; ctx->show_attractors = 0;
ctx->marking2.y = -0.11873; ctx->show_reflectors = 0;
ctx->marking3.x = -0.73679; ctx->show_rotated_reflectors = 0;
ctx->marking3.y = -0.21873; ctx->show_limit = 0;
ctx->distance_parameter1 = 0.45; ctx->show_dual_limit = 0;
ctx->distance_parameter2 = 0.2; ctx->show_text = 1;
ctx->show_coxeter_orbit = 0; ctx->mode = 0;
ctx->extra_text = malloc(1000*sizeof(char)); ctx->use_rotation_basis = 1;
memset(ctx->extra_text, 0, 1000*sizeof(char)); ctx->limit_with_lines = 0;
ctx->use_repelling = 0;
ctx->show_marking = 0;
ctx->marking.x = -0.73679;
ctx->marking.y = -0.01873;
ctx->show_coxeter_orbit = 0;
ctx->limit_curve = malloc(12*ctx->n_group_elements*sizeof(double)); ctx->limit_curve = malloc(3*ctx->n_group_elements*sizeof(double));
ctx->limit_curve_count = -1; ctx->limit_curve_count = -1;
ctx->group = malloc(ctx->n_group_elements*sizeof(groupelement_t)); ctx->group = malloc(ctx->n_group_elements_combinatorial*sizeof(groupelement_t));
generate_triangle_group(ctx->group, ctx->n_group_elements, ctx->p[0], ctx->p[1], ctx->p[2]); generate_triangle_group(ctx->group, ctx->n_group_elements_combinatorial, ctx->p[0], ctx->p[1], ctx->p[2]);
// the temporary stuff // the temporary stuff
ctx->cartan = gsl_matrix_alloc(3, 3); ctx->cartan = gsl_matrix_alloc(3, 3);
ctx->cob = gsl_matrix_alloc(3, 3); ctx->cob = gsl_matrix_alloc(3, 3);
ctx->ws = workspace_alloc(3); ctx->ws = workspace_alloc(3);
} }
void destroyContext(DrawingContext *ctx) void destroyContext(DrawingContext *ctx)
{ {
free(ctx->limit_curve); free(ctx->limit_curve);
free(ctx->group); free(ctx->group);
free(ctx->extra_text);
gsl_matrix_free(ctx->cartan); gsl_matrix_free(ctx->cartan);
gsl_matrix_free(ctx->cob); gsl_matrix_free(ctx->cob);
workspace_free(ctx->ws); workspace_free(ctx->ws);
} }
void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *word) void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
{ {
gsl_matrix *tmp = getTempMatrix(ctx->ws); gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
initializeTriangleGenerators(gen, ctx->cartan); // ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");
gsl_matrix_set_identity(tmp);
for(int i = 0; i < strlen(word); i++)
multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
rotation_frame(tmp, result, 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, 4); 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, 1);
} }
void computeBoxTransform(DrawingContext *ctx, char *word1, char *word2, gsl_matrix *result) void computeBoxTransform(DrawingContext *ctx, char *word1, char *word2, gsl_matrix *result)
{ {
vector_t p[2][3],i[2]; vector_t p[2][3],i[2];
vector_t std[4] = { vector_t std[4] = {
{-1, -1, 1}, {-1, -1, 1},
{-1, 1, 1}, {-1, 1, 1},
{1, 1, 1}, {1, 1, 1},
{1, -1, 1} {1, -1, 1}
}; };
gsl_vector **vertices = getTempVectors(ctx->ws, 4); gsl_vector **vertices = getTempVectors(ctx->ws, 4);
gsl_vector **std_vertices = getTempVectors(ctx->ws, 4); gsl_vector **std_vertices = getTempVectors(ctx->ws, 4);
gsl_matrix *tmp = getTempMatrix(ctx->ws); gsl_matrix *tmp = getTempMatrix(ctx->ws);
gsl_matrix *to_frame = getTempMatrix(ctx->ws); gsl_matrix *to_frame = getTempMatrix(ctx->ws);
gsl_matrix *to_std_frame = getTempMatrix(ctx->ws); gsl_matrix *to_std_frame = getTempMatrix(ctx->ws);
fixedPoints(ctx, word1, p[0]); fixedPoints(ctx, word1, p[0]);
fixedPoints(ctx, word2, p[1]); fixedPoints(ctx, word2, p[1]);
// intersect attracting line with neutral line of the other element // intersect attracting line with neutral line of the other element
for(int j = 0; j < 2; j++) 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])); i[j] = cross(cross(p[j%2][0],p[j%2][1]),cross(p[(j+1)%2][0],p[(j+1)%2][2]));
// box consists of p[0][0], i[0], p[1][0], i[1] // box consists of p[0][0], i[0], p[1][0], i[1]
for(int i = 0; i < 4; i++) for(int i = 0; i < 4; i++)
vectorToGsl(std[i], std_vertices[i]); vectorToGsl(std[i], std_vertices[i]);
vectorToGsl(p[0][0], vertices[0]); vectorToGsl(p[0][0], vertices[0]);
vectorToGsl(i[0], vertices[1]); vectorToGsl(i[0], vertices[1]);
vectorToGsl(p[1][0], vertices[2]); vectorToGsl(p[1][0], vertices[2]);
vectorToGsl(i[1], vertices[3]); vectorToGsl(i[1], vertices[3]);
projective_frame(std_vertices, to_std_frame, ctx->ws); projective_frame(std_vertices, to_std_frame, ctx->ws);
projective_frame(vertices, to_frame, ctx->ws); projective_frame(vertices, to_frame, ctx->ws);
invert(to_frame, tmp, ctx->ws); invert(to_frame, tmp, ctx->ws);
multiply(to_std_frame, tmp, result); multiply(to_std_frame, tmp, result);
/* /*
LOOP(i) { LOOP(i) {
LOOP(j) { LOOP(j) {
printf("%.4f ", gsl_matrix_get(result, i, j)); printf("%.4f ", gsl_matrix_get(result, i, j));
} }
printf("\n"); printf("\n");
}*/ }*/
releaseTempVectors(ctx->ws, 8); releaseTempVectors(ctx->ws, 8);
releaseTempMatrices(ctx->ws, 3); releaseTempMatrices(ctx->ws, 3);
} }
void updateMatrices(DrawingContext *ctx) void updateMatrices(DrawingContext *ctx)
{ {
double angle[3]; double angle[3];
LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i]; LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i];
cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter); cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);
gsl_matrix *tmp = getTempMatrix(ctx->ws); 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, 0, 0.0);
gsl_matrix_set(tmp, 0, 1, sqrt(3.0)/2.0); gsl_matrix_set(tmp, 0, 1, sqrt(3.0)/2.0);
gsl_matrix_set(tmp, 0, 2, -sqrt(3.0)/2.0); gsl_matrix_set(tmp, 0, 2, -sqrt(3.0)/2.0);
gsl_matrix_set(tmp, 1, 0, 1.0); gsl_matrix_set(tmp, 1, 0, 1.0);
gsl_matrix_set(tmp, 1, 1, -0.5); gsl_matrix_set(tmp, 1, 1, -0.5);
gsl_matrix_set(tmp, 1, 2, -0.5); gsl_matrix_set(tmp, 1, 2, -0.5);
gsl_matrix_set(tmp, 2, 0, 1.0); gsl_matrix_set(tmp, 2, 0, 1.0);
gsl_matrix_set(tmp, 2, 1, 1.0); gsl_matrix_set(tmp, 2, 1, 1.0);
gsl_matrix_set(tmp, 2, 2, 1.0); gsl_matrix_set(tmp, 2, 2, 1.0);
gsl_matrix_memcpy(ctx->cob, tmp); gsl_matrix_memcpy(ctx->cob, tmp);
} else if(ctx->use_rotation_basis % 5 == 1) { } else if(ctx->use_rotation_basis % nmodes == 1) {
gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason? gsl_matrix_set(tmp, 0, 0, 1.0);
} else if(ctx->use_rotation_basis % 5 == 2) { gsl_matrix_set(tmp, 0, 1, -1.0);
computeRotationMatrix(ctx, tmp, "ba"); gsl_matrix_set(tmp, 0, 2, 0.0);
invert(tmp, ctx->cob, ctx->ws); gsl_matrix_set(tmp, 1, 0, 1.0);
} else if(ctx->use_rotation_basis % 5 == 3) { gsl_matrix_set(tmp, 1, 1, 1.0);
computeBoxTransform(ctx, "bca", "abc", ctx->cob); gsl_matrix_set(tmp, 1, 2, 0.0);
// computeBoxTransform(ctx, "cab", "bca", ctx->cob); gsl_matrix_set(tmp, 2, 0, 0.0);
// computeBoxTransform(ctx, "acb", "cba", ctx->cob); gsl_matrix_set(tmp, 2, 1, 0.0);
} else { gsl_matrix_set(tmp, 2, 2, 1.0);
cartanMatrix(tmp, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0); gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason?
diagonalize_symmetric_form(tmp, ctx->cob, ctx->ws); 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 % nmodes == 3) {
computeBoxTransform(ctx, "acb", "cba", ctx->cob);
// computeBoxTransform(ctx, "cab", "bca", ctx->cob);
// computeBoxTransform(ctx, "acb", "cba", ctx->cob);
} else {
cartanMatrix(tmp, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0);
diagonalize_symmetric_form(tmp, ctx->cob, ctx->ws);
}
releaseTempMatrices(ctx->ws, 1); releaseTempMatrices(ctx->ws, 1);
} }
void output_info(DrawingContext *ctx) void output_info(DrawingContext *ctx)
{ {
vector_t p[4][3]; vector_t p[4][3];
point_t pt; point_t pt;
fixedPoints(ctx, "abc", p[0]); fixedPoints(ctx, "abc", p[0]);
fixedPoints(ctx, "bca", p[1]); fixedPoints(ctx, "bca", p[1]);
fixedPoints(ctx, "cab", p[2]); fixedPoints(ctx, "cab", p[2]);
pt = vectorToPoint(ctx, p[0][0]); pt = vectorToPoint(ctx, p[0][0]);
printf("(abc)-+ = (%f %f)\n", pt.x, pt.y); printf("(abc)-+ = (%f %f)\n", pt.x, pt.y);
pt = vectorToPoint(ctx, p[1][0]); pt = vectorToPoint(ctx, p[1][0]);
printf("(bca)-+ = (%f %f)\n", pt.x, pt.y); printf("(bca)-+ = (%f %f)\n", pt.x, pt.y);
} }
void print(DrawingContext *screen) void print(DrawingContext *screen)
{ {
DrawingContext file; DrawingContext file;
DimensionsInfo dim; DimensionsInfo dim;
cairo_surface_t *surface; cairo_surface_t *surface;
char filename[100]; char filename[100];
time_t t = time(NULL); time_t t = time(NULL);
strftime(filename, sizeof(filename), "screenshot_%Y%m%d_%H%M%S.pdf", localtime(&t)); strftime(filename, sizeof(filename), "screenshot_%Y%m%d_%H%M%S.pdf", localtime(&t));
memcpy(&file, screen, sizeof(file)); memcpy(&file, screen, sizeof(file));
dim.width = screen->dim->width; dim.width = screen->dim->width;
dim.height = screen->dim->width / sqrt(2.0); dim.height = screen->dim->width / sqrt(2.0);
dim.matrix = screen->dim->matrix; dim.matrix = screen->dim->matrix;
dim.matrix.y0 += ((double)dim.height - (double)screen->dim->height) / 2.0; // recenter vertically dim.matrix.y0 += ((double)dim.height - (double)screen->dim->height) / 2.0; // recenter vertically
updateDimensions(&dim); updateDimensions(&dim);
file.dim = &dim; file.dim = &dim;
surface = cairo_pdf_surface_create(filename, (double)dim.width, (double)dim.height); surface = cairo_pdf_surface_create(filename, (double)dim.width, (double)dim.height);
file.cairo = cairo_create(surface); file.cairo = cairo_create(surface);
draw(&file); draw(&file);
cairo_destroy(file.cairo); cairo_destroy(file.cairo);
cairo_surface_destroy(surface); cairo_surface_destroy(surface);
printf("Wrote sceenshot to file: %s\n", filename); printf("Wrote sceenshot to file: %s\n", filename);
} }
int processEvent(GraphicsInfo *info, XEvent *ev) int processEvent(GraphicsInfo *info, XEvent *ev)
{ {
int state; int state;
unsigned long key; unsigned long key;
char filename[100]; char filename[100];
// fprintf(stderr, "Event: %d\n", ev->type); // fprintf(stderr, "Event: %d\n", ev->type);
switch(ev->type) { switch(ev->type) {
case ButtonPress: case ButtonPress:
state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask); state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask);
if(ev->xbutton.button == 1 && state & ShiftMask && state & ControlMask) { if(ev->xbutton.button == 1 && state & ShiftMask) {
screen_context->marking.x = (double)ev->xbutton.x; screen_context->marking.x = (double)ev->xbutton.x;
screen_context->marking.y = (double)ev->xbutton.y; screen_context->marking.y = (double)ev->xbutton.y;
printf("mouse button pressed: %f, %f\n", screen_context->marking.x, screen_context->marking.y); printf("mouse button pressed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix); cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
cairo_device_to_user(screen_context->cairo, &screen_context->marking.x, &screen_context->marking.y); cairo_device_to_user(screen_context->cairo, &screen_context->marking.x, &screen_context->marking.y);
printf("mouse button pressed transformed: %f, %f\n", screen_context->marking.x, screen_context->marking.y); printf("mouse button pressed transformed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
return STATUS_REDRAW; return STATUS_REDRAW;
} else if(ev->xbutton.button == 1 && state & ControlMask) { }
screen_context->marking2.x = (double)ev->xbutton.x; break;
screen_context->marking2.y = (double)ev->xbutton.y;
printf("mouse button pressed: %f, %f\n", screen_context->marking2.x, screen_context->marking2.y);
cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
cairo_device_to_user(screen_context->cairo, &screen_context->marking2.x, &screen_context->marking2.y);
printf("mouse button pressed transformed: %f, %f\n", screen_context->marking2.x, screen_context->marking2.y);
return STATUS_REDRAW;
} else if(ev->xbutton.button == 1 && state & ShiftMask) {
screen_context->marking3.x = (double)ev->xbutton.x;
screen_context->marking3.y = (double)ev->xbutton.y;
printf("mouse button pressed: %f, %f\n", screen_context->marking3.x, screen_context->marking3.y);
cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
cairo_device_to_user(screen_context->cairo, &screen_context->marking3.x, &screen_context->marking3.y);
printf("mouse button pressed transformed: %f, %f\n", screen_context->marking3.x, screen_context->marking3.y);
return STATUS_REDRAW;
}
break; 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);
case KeyPress: switch(key) {
state = ev->xkey.state & (ShiftMask | LockMask | ControlMask); case XK_Down:
key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask)); if(ev->xkey.state & ShiftMask)
// printf("Key pressed: %ld\n", key); 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:
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:
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_Page_Down:
screen_context->parameter /= exp(0.02);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Page_Up:
screen_context->parameter *= exp(0.02);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case ' ':
screen_context->parameter = 5.57959706;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Return:
// screen_context->parameter = 2.76375163;
screen_context->parameter = 5.29063366;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case 'm':
printf("matrix.xx = %f;\n", info->dim->matrix.xx);
printf("matrix.xy = %f;\n", info->dim->matrix.xy);
printf("matrix.x0 = %f;\n", info->dim->matrix.x0);
printf("matrix.yx = %f;\n", info->dim->matrix.yx);
printf("matrix.yy = %f;\n", info->dim->matrix.yy);
printf("matrix.y0 = %f;\n", info->dim->matrix.y0);
break;
case 'i':
output_info(screen_context);
break;
case 'b':
TOGGLE(screen_context->show_boxes);
break;
case 'B':
TOGGLE(screen_context->show_boxes2);
break;
case 'a':
TOGGLE(screen_context->show_attractors);
break;
case 'r':
TOGGLE(screen_context->show_reflectors);
break;
case 'x':
TOGGLE(screen_context->show_rotated_reflectors);
break;
case 'L':
TOGGLE(screen_context->limit_with_lines);
break;
case 'l':
TOGGLE(screen_context->show_limit);
break;
case 'd':
TOGGLE(screen_context->show_dual_limit);
break;
case 'R':
screen_context->use_rotation_basis++;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case 'p':
print(screen_context);
break;
case 'M':
screen_context->limit_with_lines = 0;
double parameter_start = screen_context->parameter;
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, "output/%s%03d.png", screen_context->movie_filename, i);
cairo_surface_write_to_png(info->buffer_surface, filename);
printf("Finished drawing %s\n", filename);
}
switch(key) { case 'f':
case XK_Down: TOGGLE(screen_context->use_repelling);
screen_context->parameter /= exp(0.002); computeLimitCurve(screen_context);
updateMatrices(screen_context); break;
computeLimitCurve(screen_context); case 't':
break; TOGGLE(screen_context->show_text);
case XK_Up: break;
screen_context->parameter *= exp(0.002); case 'c':
updateMatrices(screen_context); TOGGLE(screen_context->show_coxeter_orbit);
computeLimitCurve(screen_context); break;
break; case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': case '0':
case XK_Left: screen_context->mode = key - '0';
screen_context->parameter /= exp(0.00002); break;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Right:
screen_context->parameter *= exp(0.00002);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Page_Down:
screen_context->parameter /= exp(0.02);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Page_Up:
screen_context->parameter *= exp(0.02);
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case ' ':
screen_context->parameter = 5.57959706;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case XK_Return:
// screen_context->parameter = 2.76375163;
screen_context->parameter = 5.29063366;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
case 'm':
printf("matrix.xx = %f;\n", info->dim->matrix.xx);
printf("matrix.xy = %f;\n", info->dim->matrix.xy);
printf("matrix.x0 = %f;\n", info->dim->matrix.x0);
printf("matrix.yx = %f;\n", info->dim->matrix.yx);
printf("matrix.yy = %f;\n", info->dim->matrix.yy);
printf("matrix.y0 = %f;\n", info->dim->matrix.y0);
break;
case 'i':
output_info(screen_context);
break;
case 'b':
TOGGLE(screen_context->show_boxes);
break;
case 'B':
TOGGLE(screen_context->show_boxes2);
break;
case 'a':
TOGGLE(screen_context->show_attractors);
break;
case 'r':
TOGGLE(screen_context->show_reflectors);
break;
case 'x':
TOGGLE(screen_context->show_rotated_reflectors);
break;
case 'L':
TOGGLE(screen_context->limit_with_lines);
break;
case 'l':
TOGGLE(screen_context->show_limit);
break;
case 'd':
TOGGLE(screen_context->show_dual_limit);
break;
case 'R':
screen_context->use_rotation_basis++;
updateMatrices(screen_context);
computeLimitCurve(screen_context);
break;
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;
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));
updateMatrices(screen_context);
computeLimitCurve(screen_context);
draw(screen_context);
sprintf(filename, "movie3/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;
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);
}
case 'f':
TOGGLE(screen_context->use_repelling);
computeLimitCurve(screen_context);
break;
case 't':
TOGGLE(screen_context->show_text);
break;
case 'c':
TOGGLE(screen_context->show_coxeter_orbit);
break;
case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': case '0':
screen_context->mode = key - '0';
break;
}
return STATUS_REDRAW;
} }
return STATUS_NOTHING; return STATUS_REDRAW;
}
return STATUS_NOTHING;
} }
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
GraphicsInfo *info; GraphicsInfo *info;
screen_context = malloc(sizeof(DrawingContext)); screen_context = malloc(sizeof(DrawingContext));
setupContext(screen_context, argc, argv); setupContext(screen_context, argc, argv);
updateMatrices(screen_context); updateMatrices(screen_context);
computeLimitCurve(screen_context); computeLimitCurve(screen_context);
info = initCairo(0, KeyPressMask, 200, 200, "Triangle group"); info = initCairo(0, KeyPressMask, 200, 200, "Triangle group");
if(!info) if(!info)
return 1; return 1;
/* /*
info->dim->matrix.xx = 274.573171; info->dim->matrix.xx = 274.573171;
info->dim->matrix.xy = 0.000000; info->dim->matrix.xy = 0.000000;
info->dim->matrix.x0 = 583.073462; info->dim->matrix.x0 = 583.073462;
info->dim->matrix.yx = 0.000000; info->dim->matrix.yx = 0.000000;
info->dim->matrix.yy = 274.573171; info->dim->matrix.yy = 274.573171;
info->dim->matrix.y0 = 777.225293; info->dim->matrix.y0 = 777.225293;
*/ */
info->dim->matrix.xx = 274.573171; info->dim->matrix.xx = 274.573171;
info->dim->matrix.xy = 0.000000; info->dim->matrix.xy = 0.000000;
info->dim->matrix.x0 = 910.073462; info->dim->matrix.x0 = 910.073462;
info->dim->matrix.yx = 0.000000; info->dim->matrix.yx = 0.000000;
info->dim->matrix.yy = 274.573171; info->dim->matrix.yy = 274.573171;
info->dim->matrix.y0 = 509.225293; info->dim->matrix.y0 = 509.225293;
updateDimensions(info->dim); updateDimensions(info->dim);
screen_context->dim = info->dim; screen_context->dim = info->dim;
screen_context->cairo = info->buffer_context; screen_context->cairo = info->buffer_context;
startTimer(info); startTimer(info);
while(1) { while(1) {
int result = checkEvents(info, processEvent, NULL); int result = checkEvents(info, processEvent, NULL);
if(result == STATUS_QUIT) if(result == STATUS_QUIT)
return 0; return 0;
else if(result == STATUS_REDRAW) { else if(result == STATUS_REDRAW) {
struct timeval current_time; struct timeval current_time;
double start_time, intermediate_time, end_time; double start_time, intermediate_time, end_time;
gettimeofday(&current_time, 0); gettimeofday(&current_time, 0);
start_time = current_time.tv_sec + current_time.tv_usec*1e-6; start_time = current_time.tv_sec + current_time.tv_usec*1e-6;
draw(screen_context); draw(screen_context);
gettimeofday(&current_time, 0); gettimeofday(&current_time, 0);
intermediate_time = current_time.tv_sec + current_time.tv_usec*1e-6; intermediate_time = current_time.tv_sec + current_time.tv_usec*1e-6;
cairo_set_source_surface(info->front_context, info->buffer_surface, 0, 0); cairo_set_source_surface(info->front_context, info->buffer_surface, 0, 0);
cairo_paint(info->front_context); cairo_paint(info->front_context);
gettimeofday(&current_time, 0); gettimeofday(&current_time, 0);
end_time = current_time.tv_sec + current_time.tv_usec*1e-6; 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);
} }
waitUpdateTimer(info);
}
free(screen_context); free(screen_context);
destroyCairo(info); destroyCairo(info);
destroyContext(screen_context); destroyContext(screen_context);
return 0; return 0;
} }

32
main.h
View File

@@ -11,7 +11,12 @@
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);} #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#define LOOP(i) for(int i = 0; i < 3; i++) #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 { typedef struct {
double x[3]; double x[3];
@@ -22,14 +27,6 @@ typedef struct {
double y; double y;
} point_t; } point_t;
typdef struct {
double x;
double y;
double angle;
vector_t fp[3];
groupelement_t *g;
} limit_curve_t;
typedef struct { typedef struct {
// infos about the screen to draw on // infos about the screen to draw on
cairo_t *cairo; cairo_t *cairo;
@@ -38,7 +35,13 @@ typedef struct {
// a priori parameter // a priori parameter
int p[3],k[3]; int p[3],k[3];
double parameter; 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;
int n_group_elements_combinatorial;
int show_boxes; int show_boxes;
int show_boxes2; int show_boxes2;
int show_attractors; int show_attractors;
@@ -52,14 +55,9 @@ typedef struct {
int limit_with_lines; int limit_with_lines;
int show_marking; int show_marking;
point_t marking; point_t marking;
point_t marking2;
point_t marking3;
double distance_parameter1;
double distance_parameter2;
int show_coxeter_orbit; int show_coxeter_orbit;
int use_repelling; int use_repelling;
gsl_matrix *cartan, *cob; gsl_matrix *cartan, *cob;
char *extra_text;
// precomputed and constant // precomputed and constant
groupelement_t *group; groupelement_t *group;
@@ -79,7 +77,8 @@ typedef enum {
// implemented in limit_set.c // implemented in limit_set.c
void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s); 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); int computeLimitCurve(DrawingContext *ctx);
// implemented in draw.c // implemented in draw.c
@@ -109,7 +108,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]);
void destroyContext(DrawingContext *ctx); void destroyContext(DrawingContext *ctx);
void print(DrawingContext *screen); void print(DrawingContext *screen);
int processEvent(GraphicsInfo *info, XEvent *ev); 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); void updateMatrices(DrawingContext *ctx);
static vector_t vectorFromGsl(gsl_vector *v) static vector_t vectorFromGsl(gsl_vector *v)

1
triangle.h
View File

@@ -14,6 +14,7 @@ typedef struct _groupelement {
struct _groupelement *parent; struct _groupelement *parent;
struct _groupelement *inverse; struct _groupelement *inverse;
int letter; int letter;
int visited;
} groupelement_t; } groupelement_t;
int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3); int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);