#include "main.h"
#include "exp_equation.h"

#define FMOD(x,y) (fmod(x,y) < 0 ? fmod(x,y) + y : fmod(x,y))
#define ANGLE_DIFF(x,y) (FMOD((x)-(y), 2*M_PI))
#define ANGLE_IN_INTERVAL(a,b,x) (ANGLE_DIFF(x,a) < ANGLE_DIFF(b,a))
#define FLIP(x,y) do {double tmp = x; x = y; y = tmp;} while(0)

// level 0: helper functions

int isInsideBB(DrawingContext *ctx, point_t p)
{
	cairo_user_to_device(ctx->cairo, &p.x, &p.y);
	return -p.x < ctx->dim->width && p.x < 3*ctx->dim->width && -p.y < ctx->dim->height && p.y < 3*ctx->dim->height;
}

vector_t cross(vector_t a, vector_t b)
{
	vector_t result;
	result.x[0] = a.x[1]*b.x[2] - a.x[2]*b.x[1];
	result.x[1] = a.x[2]*b.x[0] - a.x[0]*b.x[2];
	result.x[2] = a.x[0]*b.x[1] - a.x[1]*b.x[0];
	return result;
}

double scal(vector_t a, vector_t b)
{
	double result = 0.0;
	LOOP(i) result += a.x[i]*b.x[i];
	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;

	LOOP(i) out.x[i] = 0.0;
	LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(m, i, j) * x.x[j];

	return out;
}

vector_t apply_transpose(gsl_matrix *m, vector_t x)
{
	vector_t out;

	LOOP(i) out.x[i] = 0.0;
	LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(m, j, i) * x.x[j];

	return out;
}

vector_t apply_pseudoinverse_transpose(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, (i+1)%3, (j+1)%3) * gsl_matrix_get(m, (i+2)%3, (j+2)%3)
			- gsl_matrix_get(m, (i+1)%3, (j+2)%3) * gsl_matrix_get(m, (i+2)%3, (j+1)%3);
		out.x[i] += cofactor * in.x[j];
	}

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

// intersect the line given by the covector "line" with the orbit of "orbit_point"
// by the one--parameter subgroup of SL(3,R) which contains the element "loxodromic"
// in an eigenbasis of "loxodromic", this corresponds
int intersect_line_and_loxodromic_orbit(DrawingContext *ctx, vector_t line, gsl_matrix *frame, double *logeigenvalues, vector_t start, vector_t *out)
{
	vector_t line_in_frame = apply_transpose(frame, line);
	vector_t start_in_frame = apply_pseudoinverse(frame, start);

	vector_t a, x;
	LOOP(i) a.x[i] = logeigenvalues[i];
	LOOP(i) x.x[i] = line_in_frame.x[i]*start_in_frame.x[i];

	double t[2];
	vector_t v[2];
	int n1, n2;

	n1 = solve_linear_exp(a, x, t);
	for(int i = 0; i < n1; i++) {
		LOOP(j) v[i].x[j] = exp(a.x[j]*t[i]) * start_in_frame.x[j];
		out[i] = apply(frame, v[i]);
	}

	x.x[1] *= -1;
	n2 = solve_linear_exp(a, x, t);
	for(int i = 0; i < n2; i++) {
		LOOP(j) v[i].x[j] = exp(a.x[j]*t[i]) * start_in_frame.x[j];
		v[i].x[1] *= -1;
		out[i+n1] = apply(frame, v[i]);
	}

	if(n1+n2 > 2) {
		fprintf(stderr, "more than 2 solutions in intersect_line_and_loxodromic_orbit()!\n");
		exit(1);
	}

	return n1+n2;
}

// 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);
	gsl_matrix *ev = getTempMatrix(ctx->ws);
	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);

	initializeTriangleGenerators(gen, ctx->cartan);

	gsl_matrix_set_identity(tmp);
	for(int i = 0; i < strlen(word); i++) {
		if(word[i] == ' ')
			continue;
		multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
	}
	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);

	return count;
}

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

	initializeTriangleGenerators(gen, ctx->cartan);

	gsl_matrix_set_identity(tmp);
	for(int i = 0; i < strlen(word); i++) {
		if(word[i] == ' ')
			continue;
		multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
	}
	int count = real_eigenvalues(tmp, ev, ctx->ws);

	LOOP(i) out[i] = gsl_vector_get(ev, i);

	releaseTempMatrices(ctx->ws, 4);
	releaseTempVectors(ctx->ws, 1);

	return count;
}

// level 1: the elementary drawing functions, drawPoint, drawSegment2d

void drawPoint(DrawingContext *ctx, point_t p)
{
	cairo_t *C = ctx->cairo;

	cairo_save(C);
	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)
{
	cairo_t *C = ctx->cairo;
	cairo_move_to(C, a.x, a.y);
	cairo_line_to(C, b.x, b.y);
	cairo_stroke(C);
}

// level 2: drawVector, drawCovector, drawSegment

point_t vectorToPoint(DrawingContext *ctx, vector_t in)
{
	double x[3];
	point_t out;

	LOOP(i) x[i] = 0.0;
	LOOP(i) LOOP(j) x[i] += gsl_matrix_get(ctx->cob, i, j) * in.x[j];

	out.x = x[0] / x[2];
	out.y = x[1] / x[2];

	return out;
}

void drawVector(DrawingContext *ctx, vector_t v)
{
	drawPoint(ctx, vectorToPoint(ctx, v));
}

static void drawImplicitLine(DrawingContext *ctx, double a, double b, double c)
{
	double norm, lambda;
	point_t m, s, xminus, xplus;
	m.x = ctx->dim->center_x;
	m.y = ctx->dim->center_y;
	lambda = (a*m.x + b*m.y + c)/(a*a + b*b);
	s.x = m.x - lambda*a;
	s.y = m.y - lambda*b;
	norm = sqrt(a*a + b*b);
	xminus.x = s.x - ctx->dim->radius * b / norm;
	xminus.y = s.y + ctx->dim->radius * a / norm;
	xplus.x = s.x + ctx->dim->radius * b / norm;
	xplus.y = s.y - ctx->dim->radius * a / norm;

	drawSegment2d(ctx, xminus, xplus);
}

void drawCovector(DrawingContext *ctx, vector_t v)
{
	double x[3];
	double cofactor;

	LOOP(i) x[i] = 0.0;
	LOOP(i) LOOP(j) {
		cofactor = gsl_matrix_get(ctx->cob, (i+1)%3, (j+1)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+2)%3)
			- gsl_matrix_get(ctx->cob, (i+1)%3, (j+2)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+1)%3);
		x[i] += cofactor * v.x[j];
	}

	drawImplicitLine(ctx, x[0], x[1], x[2]);
}

void drawSegment(DrawingContext *ctx, vector_t a, vector_t b)
{
	drawSegment2d(ctx, vectorToPoint(ctx, a), vectorToPoint(ctx, b));
}

void drawSegmentWith(DrawingContext *ctx, vector_t a, vector_t b, vector_t line, int contains)
{
	point_t a_ = vectorToPoint(ctx,a);
	point_t b_ = vectorToPoint(ctx,b);
	double x[3];
	double cofactor;
	double r, tline, tminus, tplus;
	double coeff0, coeff1, coeff2;
	point_t m, xminus, xplus;

	// multiply line with inverse of cob to get it as implicit line x in chart
	LOOP(i) x[i] = 0.0;
	LOOP(i) LOOP(j) {
		cofactor = gsl_matrix_get(ctx->cob, (i+1)%3, (j+1)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+2)%3)
			- gsl_matrix_get(ctx->cob, (i+1)%3, (j+2)%3) * gsl_matrix_get(ctx->cob, (i+2)%3, (j+1)%3);
		x[i] += cofactor * line.x[j];
	}

	// 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
		// find t so that s(t) is at radius r from center, |s(t)-m| = r
		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);
		coeff2 = (b_.x - a_.x)*(b_.x - a_.x) + (b_.y - a_.y)*(b_.y - a_.y);
		if(coeff1*coeff1 - coeff0*coeff2 <= 0)
			return;
		tplus  = (- coeff1 + sqrt(coeff1*coeff1 - coeff0*coeff2))/coeff2;
		tminus = (- coeff1 - sqrt(coeff1*coeff1 - coeff0*coeff2))/coeff2;
		xplus.x = a_.x + tplus * (b_.x - a_.x);
		xplus.y = a_.y + tplus * (b_.y - a_.y);
		xminus.x = a_.x + tminus * (b_.x - a_.x);
		xminus.y = a_.y + tminus * (b_.y - a_.y);
		if(tplus > 1)
			drawSegment2d(ctx, b_, xplus);
		if(tminus < 0)
			drawSegment2d(ctx, a_, xminus);
	}
}

// level 3: boxes and polygons

void drawPolygon(DrawingContext *ctx, int segments, int sides, ...)
{
	va_list args;
	vector_t first, prev, current;

	va_start(args, sides);

	first = va_arg(args, vector_t);
	current = first;
	for(int i = 0; i < sides-1; i++) {
		prev = current;
		current = va_arg(args, vector_t);
		if(segments)
			drawSegment(ctx, prev, current);
		else
			drawCovector(ctx, cross(prev, current));
	}
	if(segments)
		drawSegment(ctx, current, first);
	else
		drawCovector(ctx, cross(current, first));

	va_end(args);
}

void drawTriangle(DrawingContext *ctx, const char *word)
{
	vector_t p[3];

	fixedPoints(ctx, word, p);
	drawPolygon(ctx, 1, 3, p[0], p[1], p[2]);
}

void drawBox(DrawingContext *ctx, const char *word1, const char *word2)
{
	vector_t p[2][3],i[2];

	fixedPoints(ctx, word1, p[0]);
	fixedPoints(ctx, word2, p[1]);

	// 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]);
}

void drawBoxLines(DrawingContext *ctx, const char *word1, const char *word2)
{
	vector_t p[2][3],i[2];

	fixedPoints(ctx, word1, p[0]);
	fixedPoints(ctx, word2, p[1]);

	// 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, 0, 4, p[0][0], i[0], p[1][0], i[1]);
}



void drawBoxStd(DrawingContext *ctx, const char *word, char base)
{
	char word1[100];
	char word2[100];

	int len = strlen(word);
	if(len*2 + 4 > 100)
		return;

	for(int i = 0; i < len; i++) {
		word1[i] = word1[2*len+2-i] = word[i];
		word2[i] = word2[2*len+2-i] = word[i];
	}
	word1[2*len+3] = 0;
	word2[2*len+3] = 0;

	LOOP(i) word1[len+i] = (base-'A'+6+i+1)%3+'a';
	LOOP(i) word2[len+i] = (base-'A'+6-i-1)%3+'a';

//	printf("Words: %s %s\n", word1, word2);

	drawBox(ctx, word1, word2);
}

void drawRotationOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, vector_t start)
{
	vector_t v[3], w;
	point_t p;
	double parameter, startangle;
	int iterations = 200;
	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);
	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]);
	solve(frame, start_v, start_in_frame, ctx->ws);
	parameter = gsl_vector_get(start_in_frame, 2);
	parameter /= sqrt(gsl_vector_get(start_in_frame, 0)*gsl_vector_get(start_in_frame, 0) +
	                  gsl_vector_get(start_in_frame, 1)*gsl_vector_get(start_in_frame, 1));
	startangle = atan2(gsl_vector_get(start_in_frame, 1), gsl_vector_get(start_in_frame, 0));

	int previous_inside = 0;
	for(int k = 0; k <= iterations; k++) {
		LOOP(i) w.x[i] = parameter * v[2].x[i] + cos(2*k*M_PI/iterations) * v[0].x[i] + sin(2*k*M_PI/iterations) * v[1].x[i];
		p = vectorToPoint(ctx, w);

		if(isInsideBB(ctx, p)) {
			if(!previous_inside)
				cairo_move_to(C, p.x, p.y);
			else
				cairo_line_to(C, p.x, p.y);
			previous_inside = 1;
		} else {
			previous_inside = 0;
		}
	}

	cairo_stroke(C);

	releaseTempMatrices(ctx->ws, 1);
	releaseTempVectors(ctx->ws, 2);
}

void drawLoxodromicOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, double *logeigenvalues, vector_t start)
{
	vector_t start_in_frame = apply_pseudoinverse(frame, start);
	int iterations = 500;
	double stepsize = 0.02;
	vector_t x, w;
	point_t p;
	cairo_t *C = ctx->cairo;
	double t;

	int previous_inside = 0;
	for(int k = 0; k <= iterations; k++) {
		// 0 = repelling fixed point, iterations/2 = attracting fixed point
		if(k == 0 || k == iterations) {
			w.x[0] = 0.0; w.x[1] = 0.0; w.x[2] = 1.0;
		} else if(k == iterations/2) {
			w.x[0] = 1.0; w.x[1] = 0.0; w.x[2] = 0.0;
		} else if(k < iterations/2) {
			t = (k-(double)iterations/4.0)*stepsize;
			LOOP(i) w.x[i] = start_in_frame.x[i] * exp(logeigenvalues[i]*t);
			w.x[1] *= -1;
		} else {
			t = ((double)iterations*3.0/4.0-k)*stepsize;
			LOOP(i) w.x[i] = start_in_frame.x[i] * exp(logeigenvalues[i]*t);
		}

		x = apply(frame, w);
		p = vectorToPoint(ctx, x);

		if(isInsideBB(ctx, p)) {
			if(!previous_inside)
				cairo_move_to(C, p.x, p.y);
			else
				cairo_line_to(C, p.x, p.y);
			previous_inside = 1;
		} else {
			previous_inside = 0;
		}
	}

	cairo_stroke(C);
}

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

	computeRotationMatrix(ctx, frame, word);
	drawRotationOrbitFrame(ctx, frame, start);

	releaseTempMatrices(ctx->ws, 1);
}

void drawDualRotationOrbit(DrawingContext *ctx, const char *word, vector_t start)
{
	vector_t v[3], w;
	point_t p;
	double parameter, startangle;
	int iterations = 200;
	gsl_matrix *frame = getTempMatrix(ctx->ws);
	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);
	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]);
//	solve(frame, start_v, start_in_frame, ctx->ws);
	gsl_blas_dgemv(CblasTrans, 1.0, frame, start_v, 0.0, start_in_frame);
	parameter = sqrt(gsl_vector_get(start_in_frame, 0)*gsl_vector_get(start_in_frame, 0) +
	                  gsl_vector_get(start_in_frame, 1)*gsl_vector_get(start_in_frame, 1));
	parameter /= gsl_vector_get(start_in_frame, 2);
	startangle = atan2(gsl_vector_get(start_in_frame, 1), gsl_vector_get(start_in_frame, 0));

	int previous_inside = 0;
	for(int k = 0; k <= iterations; k++) {
		LOOP(i) w.x[i] = parameter * v[2].x[i] + cos(2*k*M_PI/iterations) * v[0].x[i] + sin(2*k*M_PI/iterations) * v[1].x[i];
		p = vectorToPoint(ctx, w);

		if(isInsideBB(ctx, p)) {
			if(!previous_inside)
				cairo_move_to(C, p.x, p.y);
			else
				cairo_line_to(C, p.x, p.y);
			previous_inside = 1;
		} else {
			previous_inside = 0;
		}
	}

	cairo_stroke(C);

	releaseTempMatrices(ctx->ws, 2);
	releaseTempVectors(ctx->ws, 2);
}

void drawArcWithOutput(DrawingContext *ctx, const char *word, vector_t start, vector_type_t starttype, vector_t end, vector_type_t endtype, vector_t third, int contain, vector_t *start_vector_out, vector_t *end_vector_out, int dontdraw)
{
	vector_t v[3], w, w_;
	point_t p, p_;
	double radius, angle_start, angle_end, angle_third, angle, angle_end_delta, sign, angle_start_final, angle_end_final, angle_end_other;
	int iterations = 200;
	gsl_matrix *frame = getTempMatrix(ctx->ws);
	gsl_matrix *inverse = getTempMatrix(ctx->ws);
	gsl_vector *vector = getTempVector(ctx->ws);
	gsl_vector *vector_in_frame = getTempVector(ctx->ws);
	cairo_t *C = ctx->cairo;

	computeRotationMatrix(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]);
	if(starttype == VT_POINT) {
		solve(frame, vector, vector_in_frame, ctx->ws);
		radius = sqrt(gsl_vector_get(vector_in_frame, 0)*gsl_vector_get(vector_in_frame, 0) +
		                 gsl_vector_get(vector_in_frame, 1)*gsl_vector_get(vector_in_frame, 1));
		radius /= fabs(gsl_vector_get(vector_in_frame, 2));
		angle_start = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2),
		                    gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2));
	} else {
		gsl_blas_dgemv(CblasTrans, 1.0, frame, vector, 0.0, vector_in_frame);
		radius = fabs(gsl_vector_get(vector_in_frame, 2));
		radius /= sqrt(gsl_vector_get(vector_in_frame, 0)*gsl_vector_get(vector_in_frame, 0) +
		               gsl_vector_get(vector_in_frame, 1)*gsl_vector_get(vector_in_frame, 1));

		angle_start = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2),
		                    gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2));
	}

	LOOP(i) gsl_vector_set(vector, i, third.x[i]);
	solve(frame, vector, vector_in_frame, ctx->ws);
	angle_third = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2),
	                    gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2));

	LOOP(i) gsl_vector_set(vector, i, end.x[i]);
	if(endtype == VT_POINT) {
		solve(frame, vector, vector_in_frame, ctx->ws);
		angle_end = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2),
		                  gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2));
	} else {
		gsl_blas_dgemv(CblasTrans, 1.0, frame, vector, 0.0, vector_in_frame);

		// this is only the average angle
		angle_end = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2),
		                  gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2));

		angle_end_delta = acos(-fabs(gsl_vector_get(vector_in_frame, 2))/radius/
		                       sqrt(gsl_vector_get(vector_in_frame, 0)*gsl_vector_get(vector_in_frame, 0) +
		                            gsl_vector_get(vector_in_frame, 1)*gsl_vector_get(vector_in_frame, 1)));
	}

	int previous_inside = 0;

	for(int i = 0; i < 4; i++) {
		angle_start_final = angle_start;

		if(endtype == VT_POINT) {
			angle_end_final = angle_end;
		} else {
			if(i >= 2) {
				angle_end_final = angle_end - angle_end_delta;
				angle_end_other = angle_end + angle_end_delta;
			} else {
				angle_end_final = angle_end + angle_end_delta;
				angle_end_other = angle_end - angle_end_delta;
			}
		}

		if(i%2)
			FLIP(angle_start_final, angle_end_final);

		if(endtype == VT_LINE && ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_end_other))
			continue;
		if(contain && !ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_third))
			continue;
		if(!contain && ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_third))
			continue;

		break;
	}

	// output the start end end point
	/*
	LOOP(i) w.x[i] = v[2].x[i] / radius + cos(angle_start_final) * v[0].x[i] + sin(angle_start_final) * v[1].x[i];
	p = vectorToPoint(ctx, w);
	LOOP(i) w.x[i] = v[2].x[i] / radius + cos(angle_end_final) * v[0].x[i] + sin(angle_end_final) * v[1].x[i];
	p_ = vectorToPoint(ctx, w);
	printf("\\draw (%f,%f) -- (%f,%f);\n", p.x, p.y, p_.x, p_.y);
	*/

	if(!dontdraw) {
		for(int k = 0; k <= iterations; k++) {
			angle = angle_start_final + (double)k/(double)iterations * ANGLE_DIFF(angle_end_final, angle_start_final);

			LOOP(i) w.x[i] = v[2].x[i] / radius + cos(angle) * v[0].x[i] + sin(angle) * v[1].x[i];
			p = vectorToPoint(ctx, w);

			if(isInsideBB(ctx, p)) {
				if(!previous_inside)
					cairo_move_to(C, p.x, p.y);
				else
					cairo_line_to(C, p.x, p.y);
				previous_inside = 1;
			} else {
				previous_inside = 0;
			}
		}
	}

	if(start_vector_out)
		LOOP(i) start_vector_out->x[i] = v[2].x[i] / radius +
			cos(angle_start_final) * v[0].x[i] +
			sin(angle_start_final) * v[1].x[i];
	if(end_vector_out)
		LOOP(i) end_vector_out->x[i] = v[2].x[i] / radius +
			cos(angle_end_final) * v[0].x[i] +
			sin(angle_end_final) * v[1].x[i];

	cairo_stroke(C);

	releaseTempMatrices(ctx->ws, 2);
	releaseTempVectors(ctx->ws, 2);
}

void drawArc(DrawingContext *ctx, const char *word, vector_t start, vector_type_t starttype, vector_t end, vector_type_t endtype, vector_t third, int contain)
{
	drawArcWithOutput(ctx, word, start, starttype, end, endtype, third, contain, 0, 0, 0);
}

// level 4: draw the actual image components

void drawReflectors(DrawingContext *ctx)
{
	vector_t v[3];

	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);

	LOOP(i) LOOP(j) { v[i].x[j] = (i==j) ? 1.0 : 0.0; }
	LOOP(i) drawVector(ctx, v[i]);

	LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(ctx->cartan, i, j);
	LOOP(i) drawCovector(ctx, v[i]);
}

void drawAttractors(DrawingContext *ctx)
{
	int n = 3;
	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, "a cab a", p[3]);
	fixedPoints(ctx, "b abc b", p[4]);
	fixedPoints(ctx, "c bca c", p[5]);

	double color[6][3] = {{1,0,0},{0,0.7,0},{0,0,1},{0,1,1},{0,1,1},{0,1,1}};

	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++) {
		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++) {
		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]);
		}
	}
}

char *conjugate_word(const char *word, int modifier, const char *conj, char *buffer)
{
	int wordlen = strlen(word);
	int conjlen = strlen(conj);

	for(int i = 0; i < conjlen; i++) {
		buffer[i] = conj[i];
		buffer[2*conjlen+wordlen-1-i] = conj[i];
	}

	for(int i = 0; i < wordlen; i++) {
		if(word[i] == ' ')
			buffer[conjlen+i] = word[i];
		else
			buffer[conjlen+i] = (word[i]+modifier-'a')%3 + 'a';
	}

	buffer[2*conjlen + wordlen] = 0;

	return buffer;
}

void drawCurvedBox(DrawingContext *ctx, int base, const char *conj, int style)
{
	vector_t p[11][3];
	vector_t l[2][3];
	vector_t corner1, corner2;
	vector_t tmp1, tmp2;
	char word[100];
	int modifier = base - 'A';

	conjugate_word("abc", modifier, conj, word);
	fixedPoints(ctx, word, p[0]);
	conjugate_word("bca", modifier, conj, word);
	fixedPoints(ctx, word, p[1]);
	conjugate_word("b abc b", modifier, conj, word);
	fixedPoints(ctx, word, p[2]);
	conjugate_word("ab abc ba", modifier, conj, word);
	fixedPoints(ctx, word, p[3]);
	conjugate_word("baca cab acab", modifier, conj, word);
	fixedPoints(ctx, word, p[4]);
	conjugate_word("abaca cab acaba", modifier, conj, word);
	fixedPoints(ctx, word, p[5]);
	conjugate_word("bca cab acb", modifier, conj, word);
	fixedPoints(ctx, word, p[6]);
	conjugate_word("abca cab acba", modifier, conj, word);
	fixedPoints(ctx, word, p[7]);
	conjugate_word("abacababa", modifier, conj, word);
	fixedPoints(ctx, word, p[8]);
	conjugate_word("bacabab", modifier, conj, word);
	fixedPoints(ctx, word, p[9]);
	conjugate_word("cab", modifier, conj, word);
	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]);

	// main conic
	conjugate_word("ab", modifier, conj, word);
	drawArcWithOutput(ctx, word, p[0][0], VT_POINT, p[2][0], VT_POINT, p[1][0], 0, &tmp1, &tmp2, style != 1);
	if(style == 2)
		drawSegmentWith(ctx, tmp1, tmp2, l[10][0], 0);
	if(style == 3) {
		drawVector(ctx, tmp1);
		drawVector(ctx, tmp2);
	}

	conjugate_word("ab", modifier, conj, word);
	drawArcWithOutput(ctx, word, p[1][0], VT_POINT, p[3][0], VT_POINT, p[0][0], 0, &tmp1, &tmp2, style != 1);
	if(style == 2)
		drawSegmentWith(ctx, tmp1, tmp2, l[10][0], 0);
	if(style == 3) {
		drawVector(ctx, tmp1);
		drawVector(ctx, tmp2);
	}

	conjugate_word("bcabcb", modifier, conj, word);
	drawArcWithOutput(ctx, word, p[2][0], VT_POINT, l[1][0], VT_LINE, p[6][0], 1, &tmp1, &tmp2, style != 1);  // only 1st cutoff
	if(style == 2)
		drawSegmentWith(ctx, tmp1, tmp2, l[10][0], 0);

	corner1 = tmp2;

	conjugate_word("abcabcba", modifier, conj, word);
	drawArcWithOutput(ctx, word, p[3][0], VT_POINT, l[0][0], VT_LINE, p[7][0], 1, &tmp1, &tmp2, style != 1);  // only 1st cutoff
	if(style == 2)
		drawSegmentWith(ctx, tmp1, tmp2, l[10][0], 0);
	corner2 = tmp2;

	if(style == 1 || style == 2) {
		drawSegmentWith(ctx, p[0][0], corner2, l[1][1], 0);
		drawSegmentWith(ctx, p[1][0], corner1, l[0][1], 0);
	} else if(style == 3)
	{
		drawVector(ctx, corner1);
		drawVector(ctx, corner2);
	}
}


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

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

	LOOP(i) LOOP(j) gsl_matrix_set(tmp, i, j, p[0][j].x[i]);
	double evs[3];
	wordEigenvalues(ctx, "abc", evs);
	LOOP(i) evs[i] = log(fabs(evs[i]));
	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][0]);
	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][2]);

	vector_t x;
	for(int i = 0; i < ctx->n_group_elements; i++) {
		LOOP(j) x.x[j] = ctx->limit_curve[12*i+3+j];
		x = apply_pseudoinverse(tmp, x);
		printf("%f\n",
		       pow(fabs(x.x[0]), evs[1]-evs[2]) *
		       pow(fabs(x.x[1]), evs[2]-evs[0]) *
		       pow(fabs(x.x[2]), evs[0]-evs[1]));
	}

	cairo_restore(C);
	releaseTempMatrices(ctx->ws, 9);
}

void drawBoxes2(DrawingContext *ctx)
{
	gsl_matrix *rot = getTempMatrix(ctx->ws);
	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[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]);

	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]);
	drawVector(ctx, a);
	drawVector(ctx, b);
	drawCovector(ctx, cross(a,b));
	drawCovector(ctx, cross(b,vertex[0]));
	drawCovector(ctx, cross(vertex[0],a));

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

	getMinMaxHalfCrossRatio(ctx->limit_curve, ctx->n_group_elements,
	                        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);

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

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

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

	cairo_set_source_rgb(C, 0, 0.6, 0.1);

	LOOP(i) LOOP(j) gsl_matrix_set(tmp, i, j, p[0][j].x[i]);
	double evs[3];
	wordEigenvalues(ctx, "abc", evs);
	LOOP(i) evs[i] = log(fabs(evs[i]));
	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][0]);
	drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][2]);

	vector_t intersection[2];
	double hC0, hC1;

	cairo_set_source_rgb(C, 1, 0, 0);
	intersect_line_and_loxodromic_orbit(ctx, cross(a,b), tmp, evs, p[1][2], intersection);
	hC0 = halfCR(a, b, intersection[0]);
	hC1 = halfCR(a, b, intersection[1]);
	if(hC0 > hC1)
		drawVector(ctx, intersection[0]);
	else
		drawVector(ctx, intersection[1]);
	approx_perimeter1 += log(hC0 > hC1 ? hC0 : hC1);

	cairo_set_source_rgb(C, 0, 0, 1);
	intersect_line_and_loxodromic_orbit(ctx, cross(a,b), tmp, evs, p[1][0], intersection);
	hC0 = halfCR(b, a, intersection[0]);
	hC1 = halfCR(b, a, intersection[1]);
	if(hC0 > hC1)
		drawVector(ctx, intersection[0]);
	else
	drawVector(ctx, intersection[1]);
	approx_perimeter2 += log(hC0 > hC1 ? hC0 : hC1);

	/*
	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, 7);
	releaseTempVectors(ctx->ws, 2);
}

void drawRotatedReflectors(DrawingContext *ctx)
{
	gsl_matrix *rot = getTempMatrix(ctx->ws);
	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
	cairo_t *C = ctx->cairo;
	vector_t fp[3], fp2[3];
	vector_t w;
	vector_t v[3];

	cairo_set_source_rgb(C, 0.7, 0.7, 0.7);

	initializeTriangleGenerators(gen, ctx->cartan);

	LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(ctx->cartan, i, j);
	multiply(gen[0], gen[1], rot);

	for(int j = 0; j < ctx->p[2]; j++) {
		drawCovector(ctx, v[0]);
		v[0] = apply_transpose(rot, v[0]);
	}

	LOOP(i) LOOP(j) { v[i].x[j] = (i==j) ? 1.0 : 0.0; }

	for(int j = 0; j < ctx->p[2]; j++) {
		drawVector(ctx, v[0]);
		v[0] = apply(rot, v[0]);
	}

	fixedPoints(ctx, "cab", fp);
	fixedPoints(ctx, "cacabac", fp2);
	drawRotationOrbit(ctx, "ac", fp[0]);

	releaseTempMatrices(ctx->ws, 4);
}

void drawDualLimitCurve(DrawingContext *ctx)
{
	cairo_t *C = ctx->cairo;

	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];

	fixedPoints(ctx, "abc", p[0]);
	fixedPoints(ctx, "ab abc ba", p[1]);
	fixedPoints(ctx, "abab abc baba", p[2]);
	fixedPoints(ctx, "ababab abc bababa", p[3]);
	fixedPoints(ctx, "abababab abc babababa", p[4]);
	fixedPoints(ctx, "babababa abc abababab", p[5]);
	fixedPoints(ctx, "bababa abc ababab", p[6]);
	fixedPoints(ctx, "baba abc abab", p[7]);
	fixedPoints(ctx, "ba abc ab", p[8]);

	fixedPoints(ctx, "bca", p[9]);
	fixedPoints(ctx, "ab bca ba", p[10]);
	fixedPoints(ctx, "abab bca baba", p[11]);
	fixedPoints(ctx, "ababab bca bababa", p[12]);
	fixedPoints(ctx, "abababab bca babababa", p[13]);
	fixedPoints(ctx, "babababa bca abababab", p[14]);
	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]);
	fixedPoints(ctx, "ac abc ca", p[1]);
	fixedPoints(ctx, "acac abc caca", p[2]);
	fixedPoints(ctx, "acacac abc cacaca", p[3]);
	fixedPoints(ctx, "acacacac abc cacacaca", p[4]);
	fixedPoints(ctx, "cacacaca abc acacacac", p[5]);
	fixedPoints(ctx, "cacaca abc acacac", p[6]);
	fixedPoints(ctx, "caca abc acac", p[7]);
	fixedPoints(ctx, "ca abc ac", p[8]);

	fixedPoints(ctx, "bca", p[9]);
	fixedPoints(ctx, "ac bca ca", p[10]);
	fixedPoints(ctx, "acac bca caca", p[11]);
	fixedPoints(ctx, "acacac bca cacaca", p[12]);
	fixedPoints(ctx, "acacacac bca cacacaca", p[13]);
	fixedPoints(ctx, "cacacaca bca acacacac", p[14]);
	fixedPoints(ctx, "cacaca bca acacac", p[15]);
	fixedPoints(ctx, "caca bca acac", p[16]);
	fixedPoints(ctx, "ca bca ac", p[17]);
	*/

	/*
	fixedPoints(ctx, "cab", p[2]);
	fixedPoints(ctx, "b abc b", p[3]);
	fixedPoints(ctx, "c bca c", p[4]);
	fixedPoints(ctx, "a cab a", p[5]);
	*/

	/*
	for(int i = 0; i < n; i++) {
		LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]);
		drawCovector(ctx, l[i][0]);
//		drawCovector(ctx, l[i][2]);
}*/

	fixedPoints(ctx, "abc", ptmp);
	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
	fixedPoints(ctx, "bca", ptmp);
	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
	fixedPoints(ctx, "cab", ptmp);
	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
	fixedPoints(ctx, "babcb", ptmp);
	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
	fixedPoints(ctx, "cbcac", ptmp);
	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
	fixedPoints(ctx, "acaba", ptmp);
	drawCovector(ctx, cross(ptmp[0], ptmp[1]));

	cairo_restore(C);
}

void drawLimitCurve(DrawingContext *ctx)
{
	int close_points = 0;
	cairo_t *C = ctx->cairo;

	cairo_save(C);

//	cairo_set_source_rgb(C, 0.6, 0.6, 0.6);
	cairo_set_source_rgb(C, 0, 0, 0);

	if(ctx->limit_with_lines) {
		int previous_inside = 0;
		cairo_new_path(C);

		for(int i = 0; i < ctx->limit_curve_count; i++) {
			point_t p;
			p.x = ctx->limit_curve[12*i];
			p.y = ctx->limit_curve[12*i+1];

			if(isInsideBB(ctx, p)) {
				if(!previous_inside)
					cairo_move_to(C, p.x, p.y);
				else
					cairo_line_to(C, p.x, p.y);
				previous_inside = 1;
			} else {
				previous_inside = 0;
			}
		}

		cairo_stroke(C);
	} else {
		for(int i = 0; i < ctx->limit_curve_count; i++) {
			point_t p;
			p.x = ctx->limit_curve[12*i];
			p.y = ctx->limit_curve[12*i+1];

			if(isInsideBB(ctx, p)) {
				cairo_arc(C, p.x, p.y, 2.0/ctx->dim->scalefactor, 0, 2*M_PI);
				cairo_close_path(C);
				cairo_fill(C);
			}
		}
	}

	cairo_restore(C);
}

void drawCoxeterOrbit(DrawingContext *ctx)
{
	gsl_matrix *rot = getTempMatrix(ctx->ws);
	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
	gsl_vector *eval = getTempVector(ctx->ws);
	gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws);
	gsl_vector *startpoint_coxeterbasis = getTempVector(ctx->ws);
	gsl_vector *startpoint_globalbasis = getTempVector(ctx->ws);
	gsl_vector *startpoint_drawbasis = getTempVector(ctx->ws);
	gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);

	cairo_t *C = ctx->cairo;
	vector_t cox[3][3];
	vector_t abcb[3];
	double ev[3];
	vector_t v, start;
	point_t p;
	int first = 1;

	cairo_save(C);
	initializeTriangleGenerators(gen, ctx->cartan);

	cairo_set_source_rgb(C, 0, 0, 1);

	fixedPoints(ctx, "abc", cox[0]);
	fixedPoints(ctx, "bca", cox[1]);
	fixedPoints(ctx, "cab", cox[2]);
//	fixedPoints(ctx, "babc", abcb);
	wordEigenvalues(ctx, "abc", ev);
	LOOP(i) LOOP(j) gsl_matrix_set(coxeter_fixedpoints, j, i, cox[0][i].x[j]);

	initializeTriangleGenerators(gen, ctx->cartan);
	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]);

	//printf("coxeter eigenvalues: %f %f %f\n", ev[0], ev[1], ev[2]);

//	LOOP(i) gsl_vector_set(startpoint_globalbasis, i, cox[1][0].x[i]);

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

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

	/*
	for(int t = -1000; t < 1000; t++) {
		LOOP(i) v.x[i] = 0;
		LOOP(i) LOOP(j) v.x[j] += pow(fabs(ev[i]),t*0.01)*start.x[i]*cox[0][i].x[j];
		p = vectorToPoint(ctx, v);

		if(first) {
			cairo_move_to(C, p.x, p.y);
			first = 0;
		} else {
			cairo_line_to(C, p.x, p.y);
		}
	}
	cairo_stroke(C);
	*/

	for(int i = 0; i < ctx->n_group_elements; i++) {
		v = apply(elements[i], start);
		drawVector(ctx, v);
//		ctx->limit_curve[3*i+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));
	}


/*	for(int t = -20; t < 20; t++) {
		LOOP(i) v.x[i] = 0;
		LOOP(i) LOOP(j) v.x[j] += (ev[i]<0&&t%2?-1:1)*pow(fabs(ev[i]),t/3.0)*start.x[i]*cox[0][i].x[j];
		drawVector(ctx, v);
		}*/

	cairo_set_source_rgb(C, 0, 0, 0);
//	LOOP(i) drawVector(ctx, abcb[i]);

	cairo_restore(C);
	releaseTempMatrices(ctx->ws, 5 + ctx->n_group_elements);
	releaseTempVectors(ctx->ws, 4);
}

void drawText(DrawingContext *ctx)
{
	cairo_move_to(ctx->cairo, 15, 30);
	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
	char buf[1000];
	snprintf(buf, 1000, "t = exp(%.8f) = %.8f, %s", log(ctx->parameter), ctx->parameter, ctx->extra_text);
	cairo_show_text(ctx->cairo, buf);
}

// level 5: put everything together

void draw(DrawingContext *ctx)
{
	cairo_t *C = ctx->cairo;

	cairo_set_source_rgb(C, 1, 1, 1);
	cairo_paint(C);

	cairo_set_matrix(C, &ctx->dim->matrix);

	// defaults; use save/restore whenever these are changed
	cairo_set_line_width(C, 1.0/ctx->dim->scalefactor);
	cairo_set_font_size(C, 16);
	cairo_set_line_join(C, CAIRO_LINE_JOIN_BEVEL);
	cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND);

	if(ctx->limit_curve_count >= 0) {
		if(ctx->show_limit)
			drawLimitCurve(ctx);

		if(ctx->show_dual_limit)
			drawDualLimitCurve(ctx);

		if(ctx->show_attractors)
			drawAttractors(ctx);

		if(ctx->show_rotated_reflectors)
			drawRotatedReflectors(ctx);

		if(ctx->show_reflectors)
			drawReflectors(ctx);

		if(ctx->show_boxes)
			drawBoxes(ctx);

		if(ctx->show_boxes2)
			drawBoxes2(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);
		}

		if(ctx->show_coxeter_orbit)
			drawCoxeterOrbit(ctx);
	}

	cairo_identity_matrix(C); // text is in screen coordinates

	if(ctx->show_text)
		drawText(ctx);

	cairo_surface_flush(cairo_get_target(C));
}