generate beautiful pictures of the complex Anosov set

.gitignore

@@ -1,17 +1,7 @@
 *.o
-triangle_group/singular_values
-.#*
-singular_values
-singular_values_mpi
-singular_values_barbot
-output/
-special_element
-max_slope_picture/generate
-convert
-billiard_words
-*.pnm
-*.png
-*.hi
+.\#*
+\#*\#
 gmon.out
-restart
 core
+output/
+complex_anosov

Makefile

@@ -1,50 +1,27 @@
-HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h qext.h
+HEADERS=mat.h coxeter.h enumerate.h generators.h qext.h
 
-SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
+#SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
 #SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline
-#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
+SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
 #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL
 #SPECIAL_OPTIONS=
 
 OPTIONS=-I../mps/include -L../mps/lib -pthread -m64 -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
+CC=gcc
 
-all: singular_values special_element convert billiard_words
+all: complex_anosov
 
-convert: convert.hs
-	ghc --make -dynamic convert.hs
+complex_anosov: complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o
+	$(CC) $(OPTIONS) -o complex_anosov -lm complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o -lgmp -lmps
 
-billiard_words: billiard_words.hs
-	ghc --make -dynamic billiard_words.hs
+complex_anosov.o: complex_anosov.c $(HEADERS)
+	gcc $(OPTIONS) -c complex_anosov.c
 
-singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o parallel.o
-	mpicc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps
+enumerate.o: enumerate.c $(HEADERS)
+	gcc $(OPTIONS) -c enumerate.c
 
-singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o qext.o
-	mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o qext.o -lm -lgmp -lmps
-
-#singular_values_mpi: singular_values_mpi.o coxeter.o mat.o
-#	mpicc $(OPTIONS) -o singular_values_mpi coxeter.o singular_values_mpi.o mat.o -lm -lgmp -lmps
-
-special_element: special_element.o coxeter.o linalg.o mat.o enumerate_triangle_group.o
-	gcc $(OPTIONS) -o special_element coxeter.o linalg.o special_element.o mat.o enumerate_triangle_group.o -lm -lgmp -lmps -lgsl -lcblas
-
-singular_values.o: singular_values.c $(HEADERS)
-	gcc $(OPTIONS) -c singular_values.c
-
-singular_values_barbot.o: singular_values_barbot.c $(HEADERS)
-	gcc $(OPTIONS) -c singular_values_barbot.c
-
-#singular_values_mpi.o: singular_values_mpi.c $(HEADERS)
-#	mpicc $(OPTIONS) -c singular_values_mpi.c
-
-special_element.o: special_element.c $(HEADERS)
-	gcc $(OPTIONS) -c special_element.c
-
-enumerate_triangle_group.o: enumerate_triangle_group.c $(HEADERS)
-	gcc $(OPTIONS) -c enumerate_triangle_group.c
-
-linalg.o: linalg.c $(HEADERS)
-	gcc $(OPTIONS) -c linalg.c
+generators.o: generators.c $(HEADERS)
+	gcc $(OPTIONS) -c generators.c
 
 coxeter.o: coxeter.c $(HEADERS)
 	gcc $(OPTIONS) -c coxeter.c
@@ -52,11 +29,8 @@ coxeter.o: coxeter.c $(HEADERS)
 mat.o: mat.c $(HEADERS)
 	gcc $(OPTIONS) -c mat.c
 
-parallel.o: parallel.c $(HEADERS)
-	gcc $(OPTIONS) -c parallel.c
-
 qext.o: qext.c $(HEADERS)
 	gcc $(OPTIONS) -c qext.c
 
 clean:
-	rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.o singular_values_mpi.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words enumerate_triangle_group.o parallel.o qext.o
+	rm -f complex_anosov complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o

complex_anosov.c

@@ -0,0 +1,304 @@
+#include "coxeter.h"
+#include "enumerate.h"
+#include "generators.h"
+#include "mat.h"
+#include "qext.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
+
+#define LOOP(i,n) for(int i = 0; i < (n); i++)
+#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
+//#define INFO(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__)
+#define INFO(msg, ...)
+
+/*
+  Number of elements
+  up to length 0:        1
+  up to length 1:        4
+  up to length 2:       10
+  up to length 3:       22
+  up to length 4:       46
+  up to length 5:       91
+  up to length 6:      175
+  up to length 7:      334
+  up to length 8:      634
+  up to length 9:     1198
+  up to length 10:    2260
+  up to length 11:    4261
+  up to length 12:    8029
+  up to length 13:   15124
+  up to length 14:   28486
+  up to length 15:   53650
+  up to length 16:  101038
+  up to length 17:  190279
+  up to length 18:  358339
+  up to length 19:  674830
+  up to length 20: 1270846
+  up to length 21: 2393266
+  up to length 22: 4507012
+  up to length 23: 8487625
+*/
+
+static double gaussian_sqrt5_real(NUMBER x)
+{
+	double result = 0.0;
+
+	mpq_t tmp;
+	mpq_init(tmp);
+
+	// a_0 + sqrt(5)a_1 + 4a_2 + 2sqrt(5)a_3
+	mpq_set_si(tmp, 4, 1);
+	mpq_mul(tmp, tmp, x->a[2]);
+	mpq_add(tmp, tmp, x->a[0]);
+	result = mpq_get_d(tmp);
+	mpq_set_si(tmp, 2, 1);
+	mpq_mul(tmp, tmp, x->a[3]);
+	mpq_add(tmp, tmp, x->a[1]);
+	result += mpq_get_d(tmp)*sqrt(5);
+
+	mpq_clear(tmp);
+
+	return result;
+}
+
+static double gaussian_sqrt5_imag(NUMBER x)
+{
+	double result = 0.0;
+
+	mpq_t tmp;
+	mpq_init(tmp);
+
+	// a_1 + 2sqrt(5)a_2 + 14a_3
+	mpq_set_si(tmp, 14, 1);
+	mpq_mul(tmp, tmp, x->a[3]);
+	mpq_add(tmp, tmp, x->a[1]);
+	result = mpq_get_d(tmp);
+	mpq_set_si(tmp, 2, 1);
+	mpq_mul(tmp, tmp, x->a[2]);
+	result += mpq_get_d(tmp)*sqrt(5);
+
+	mpq_clear(tmp);
+
+	return result;
+}
+
+static int read_idlist(const char *filename, int *list)
+{
+	FILE *f = fopen(filename, "r");
+	if(f == NULL) {
+		fprintf(stderr, "Could not open %s\n", filename);
+		exit(1);
+	}
+
+    char *line = NULL;
+    size_t len = 0;
+    ssize_t read;
+    int n = 0;
+
+    while ((read = getline(&line, &len, f)) != -1)
+	    list[n++] = atoi(line);
+
+    if (line)
+        free(line);
+
+	fclose(f);
+
+	return n;
+}
+
+static int compare_int(const void *a, const void *b)
+{
+    const int *aa = a;
+    const int *bb = b;
+    return (*aa > *bb) - (*aa < *bb);
+}
+
+static double runtime()
+{
+	static struct timespec starttime;
+	static int started = 0;
+
+	if(!started) {
+		clock_gettime(CLOCK_MONOTONIC, &starttime);
+		started = 1;
+	}
+
+	struct timespec curtime;
+	double diff;
+	clock_gettime(CLOCK_MONOTONIC, &curtime);
+	return (curtime.tv_sec - starttime.tv_sec) + (curtime.tv_nsec - starttime.tv_nsec) / 1e9;
+}
+
+enum mode {
+	MODE_HELP,
+	MODE_EIGENVALUES,
+	MODE_SUMMARY,
+	MODE_TRACE_IDS
+};
+
+int main(int argc, char *argv[])
+{
+	char buf[100];
+	int mode;
+
+	runtime(); // start timer
+
+	// parse arguments
+	if(argc < 2 || strcmp(argv[1], "help") == 0)
+		mode = MODE_HELP;
+	else if(strcmp(argv[1], "evs") == 0)
+		mode = MODE_EIGENVALUES;
+	else if(strcmp(argv[1], "summary") == 0)
+		mode = MODE_SUMMARY;
+	else if(strcmp(argv[1], "trace_ids") == 0)
+		mode = MODE_TRACE_IDS;
+	else
+		mode = MODE_HELP;
+
+	if(mode == MODE_HELP) {
+		fprintf(stderr, "Usage: %s <help|evs|summary|trace_ids> [arguments]\n", argv[0]);
+		fprintf(stderr, "%s help                           display this page\n", argv[0]);
+		fprintf(stderr, "%s evs <n> <q1> <q2> <q3> <treal> <timag>        enumerate group and output unique (log) eigenvalue triples\n", argv[0]);
+		fprintf(stderr, "%s summary <n> <q1> <q2> <q3> <treal> <timag>    only output max slope etc.\n", argv[0]);
+		fprintf(stderr, "%s trace_ids <n> <q1> <q2> <q3> <treal> <timag>  list of ids of unique traces\n", argv[0]);
+		return 0;
+	}
+
+	if(argc < 8) {
+		fprintf(stderr, "Not enough arguments!\n");
+		return 0;
+	}
+
+	int n = atoi(argv[2]);
+	int nlist, nuniq;
+
+	int q[3];
+	mpq_t treal, timag;
+	q[0] = atoi(argv[3]);
+	q[1] = atoi(argv[4]);
+	q[2] = atoi(argv[5]);
+	mpq_inits(treal, timag, NULL);
+	mpq_set_str(treal, argv[6], 10);
+	mpq_set_str(timag, argv[7], 10);
+	mpq_canonicalize(treal);
+	mpq_canonicalize(timag);
+
+	// enumerate group
+	INFO("generate group\n");
+	group_t *group = coxeter_init_triangle(5, 5, 5, n);
+
+	// read list of elements we need to compute, or just fill it with all elements
+	INFO("prepare list\n");
+	int *idlist = malloc(n*sizeof(int));
+	char *id_file = getenv("IDLIST");
+	if(id_file != NULL) {
+		nlist = read_idlist(id_file, idlist);
+		// sort and symmetrize the list
+		qsort(idlist, nlist, sizeof(int), compare_int);
+		int ninverses = 0;
+		LOOP(i, nlist) {
+			int id = idlist[i];
+			int invid = group->elements[id].inverse->id;
+			if(!bsearch(&invid, idlist, nlist, sizeof(int), compare_int)) {
+				idlist[nlist+ninverses] = invid;
+				ninverses++;
+			}
+		}
+		nlist += ninverses;
+		qsort(idlist, nlist, sizeof(int), compare_int);
+	} else {
+		// just list all elements which have inverses
+		nlist = 0;
+		LOOP(i, n) {
+			if(group->elements[i].inverse)
+				idlist[nlist++] = i;
+		}
+	}
+
+	// get generator matrices
+	INFO("make generators\n");
+	mat gen[3];
+	LOOP(i, 3) mat_init(gen[i], 3, QT_GAUSS_SQRT5);
+	generators_triangle_reflection_group_555_complex(gen, q[0], q[1], q[2], treal, timag);
+
+	// compute the traces of all elements in idlist
+	INFO("enumerate traces\n");
+	struct tracedata *traces;
+	nuniq = enumerate_coxeter_group_traces(group, gen, &traces, idlist, nlist, 1);
+
+	// compute eigenvalues out of traces
+	INFO("compute eigenvalues\n");
+	mps_context *solver = mps_context_new();
+	mps_monomial_poly *poly = mps_monomial_poly_new(solver, 3);
+	mps_context_set_output_prec(solver, 20); // relative precision
+	mps_context_set_output_goal(solver, MPS_OUTPUT_GOAL_APPROXIMATE);
+
+	double ev_real[3], ev_imag[3], ev_abs2[3];
+	double max_slope = 0;
+	int max_slope_id = 0;
+
+	LOOP(i, nuniq) {
+		solve_characteristic_polynomial_d(solver, poly,
+		                                  gaussian_sqrt5_real(traces[i].tr),
+		                                  gaussian_sqrt5_imag(traces[i].tr),
+		                                  gaussian_sqrt5_real(traces[i].trinv),
+		                                  gaussian_sqrt5_imag(traces[i].trinv),
+		                                  ev_real, ev_imag);
+
+		LOOP(j, 3) ev_abs2[j] = ev_real[j]*ev_real[j] + ev_imag[j]*ev_imag[j];
+
+		if(fabs(ev_abs2[0]) < fabs(ev_abs2[1]))
+			SWAP(double, ev_abs2[0], ev_abs2[1]);
+		if(fabs(ev_abs2[1]) < fabs(ev_abs2[2]))
+			SWAP(double, ev_abs2[1], ev_abs2[2]);
+		if(fabs(ev_abs2[0]) < fabs(ev_abs2[1]))
+			SWAP(double, ev_abs2[0], ev_abs2[1]);
+
+		if(mode == MODE_TRACE_IDS) {
+			// we only want to record unordered pairs here
+			if(CMP(traces[i].tr, traces[i].trinv) >= 0)
+				printf("%d\n", traces[i].id);
+			continue;
+		}
+
+		if(log(ev_abs2[0]) < 1e-6) // we regard this as a finite order element
+			continue;
+
+		double slope = - log(ev_abs2[0]) / log(ev_abs2[2]);
+		if(slope > max_slope) {
+			max_slope = slope;
+			max_slope_id = traces[i].id;
+		}
+
+		if(mode == MODE_EIGENVALUES) {
+			printf("%d %f %f %f\n",
+			       traces[i].id, log(ev_abs2[0])/2, log(ev_abs2[1])/2, log(ev_abs2[2])/2);
+		}
+	}
+
+	// output summary
+	coxeter_snprint(buf, sizeof(buf), &group->elements[max_slope_id]);
+	if(mode == MODE_SUMMARY) {
+		gmp_fprintf(stdout, "%f %f %d %d %d %f %s\n", mpq_get_d(treal), mpq_get_d(timag), n, nlist, nuniq, max_slope, buf);
+	} else if(mode == MODE_EIGENVALUES) {
+		gmp_fprintf(stderr, "q = %.2f + i*%.2f\tn = %d\tnlist = %d\tnuniq = %d\tMaximal slope: %f at %s\n", mpq_get_d(treal), mpq_get_d(timag), n, nlist, nuniq, max_slope, buf);
+	} else if(mode == MODE_TRACE_IDS) {
+		gmp_fprintf(stderr, "q = %.2f + i*%.2f\tElements: %d\tTraces: %d\n", mpq_get_d(treal), mpq_get_d(timag), n, nuniq);
+	}
+
+	// clean up
+	INFO("clean up\n");
+//	mps_monomial_poly_free(solver, MPS_POLYNOMIAL(poly));
+	mps_context_free(solver);
+	enumerate_tracedata_clear(traces, nuniq);
+	free(idlist);
+	LOOP(i, 3) mat_clear(gen[i]);
+	coxeter_clear(group);
+	mpq_clears(treal, timag, NULL);
+
+	return 0;
+}

compute_picture.sh

@@ -0,0 +1,9 @@
+#!/bin/bash
+
+for i in $(seq -50 50); do
+	for j in $(seq 0 50); do
+		if [ $i -ne 0 ] || [ $j -ne 0 ]; then
+			IDLIST=output/idlist_13 ./complex_anosov summary 15124 1 1 1 $i/50 $j/50
+		fi
+	done
+done

coxeter.c

@@ -5,6 +5,8 @@
 #include "coxeter.h"
 
 #define LOOP(i,n) for(int i = 0; i < (n); i++)
+#define MIN(x,y) ((x) > (y) ? (y) : (x))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
 
 group_t *coxeter_init_triangle(int p, int q, int r, int nmax)
 {
@@ -185,3 +187,17 @@ void coxeter_clear(group_t *g)
 	free(g->lists);
 	free(g);
 }
+
+int coxeter_snprint(char *str, int buflen, groupelement_t *g)
+{
+	int n = MIN(g->length, buflen-1); // number of characters without null byte
+
+	str[n] = 0;
+	int i = n-1;
+	while(g->parent) {
+		str[i--] = 'a' + g->letter;
+		g = g->parent;
+	}
+
+	return n;
+}

coxeter.h

@@ -13,6 +13,7 @@ typedef struct _groupelement {
 	struct _groupelement *inverse;
 	struct _groupelement **left;
 	struct _groupelement **right;
+	struct _groupelement *conjugacy_class;
 } groupelement_t;
 
 typedef struct _group {
@@ -26,5 +27,6 @@ typedef struct _group {
 group_t *coxeter_init(int rank, int *coxeter_matrix, int nmax);
 group_t *coxeter_init_triangle(int p, int q, int r, int nmax);
 void coxeter_clear(group_t *g);
+int coxeter_snprint(char *str, int buflen, groupelement_t *g);
 
 #endif

enumerate.c

@@ -0,0 +1,231 @@
+#include "enumerate.h"
+
+#include <stdlib.h>
+
+#define LOOP(i,n) for(int i = 0; i < (n); i++)
+// #define INFO(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__)
+#define INFO(msg, ...)
+
+static int compare_int(const void *a, const void *b)
+{
+    const int *aa = a;
+    const int *bb = b;
+    return (*aa > *bb) - (*aa < *bb);
+}
+
+static int index_in_list(const int *list, int n, int key)
+{
+	int *ptr = bsearch(&key, list, n, sizeof(int), compare_int);
+	if(ptr)
+		return ptr - list;
+	else
+		return -1;
+}
+
+// idlist is assumed to be sorted and symmetric
+void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices, const int *idlist, int nlist)
+{
+	TYPE type = GETTYPE(gen[0]->x[0]);
+	mat_workspace *ws;
+
+	// mark elements which we need to compute
+	int ncompute = 0;
+	LOOP(i, group->size) group->elements[i].need_to_compute = 0;
+	LOOP(i, nlist) {
+		groupelement_t *cur = &group->elements[idlist[i]];
+		while(cur && cur->need_to_compute == 0) {
+			cur->need_to_compute = 1;
+			ncompute++;
+			cur = cur->parent;
+		}
+	}
+	INFO("need to compute %d elements for %d elements in list\n", ncompute, nlist);
+
+	// compute them depth first
+	int maxlen = group->elements[group->size-1].length;
+
+	groupelement_t **stack = malloc((maxlen+1)*sizeof(groupelement_t));
+	int *letter_stack = malloc((maxlen+1)*sizeof(int));
+	mat *matrix_stack = malloc((maxlen+1)*sizeof(mat));
+	int level = 0;
+
+	LOOP(i, maxlen+1) mat_init(matrix_stack[i], 3, type);
+	ws = mat_workspace_init(3, type);
+
+	stack[0] = &group->elements[0];
+	mat_identity(matrix_stack[0]);
+	letter_stack[0] = 0;
+
+	while(level >= 0) {
+		int letter = letter_stack[level];
+		groupelement_t *next = stack[level]->left[letter];
+		if(next && next->length > level && next->need_to_compute) {
+			mat_multiply(ws, matrix_stack[level+1], gen[letter], matrix_stack[level]);
+			int id = next->id;
+			int listid = index_in_list(idlist, nlist, id);
+			if(listid != -1)
+				mat_copy(matrices[listid], matrix_stack[level+1]);
+			next->need_to_compute = 0;
+			stack[level+1] = next;
+			letter_stack[level+1] = 0;
+			level++;
+		} else {
+			letter = ++letter_stack[level];
+			while(letter >= group->rank) { // done with this level, go back down
+				level--;
+				if(level < 0)
+					break;
+				letter = ++letter_stack[level];
+			}
+		}
+	}
+
+	LOOP(i, maxlen+1) mat_clear(matrix_stack[i]);
+	mat_workspace_clear(ws);
+}
+
+static int compare_tracedata(const void *a_, const void *b_)
+{
+        int d = 0;
+
+        struct tracedata **a = (struct tracedata **)a_;
+        struct tracedata **b = (struct tracedata **)b_;
+
+        d = CMP((*a)->tr,(*b)->tr);
+        if(d == 0) {
+                d = CMP((*a)->trinv, (*b)->trinv);
+        }
+
+        return d;
+}
+
+static int compare_tracedata_id(const void *a, const void *b)
+{
+        int ida = (*(struct tracedata **)a)->id;
+        int idb = (*(struct tracedata **)b)->id;
+
+        return ida > idb ? 1 : ida < idb ? -1 : 0;
+}
+
+int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out, const int *idlist, int nlist, int unique)
+{
+	TYPE t = GETTYPE(gen[0]->x[0]);
+	int nuniq;
+
+	mat *matrices = malloc(nlist*sizeof(mat));
+	struct tracedata *traces = malloc(nlist*sizeof(struct tracedata));
+	struct tracedata **distinct_traces = malloc(nlist*sizeof(struct tracedata*));
+
+	LOOP(i, nlist) mat_init(matrices[i], 3, t);
+	enumerate_coxeter_group(group, gen, matrices, idlist, nlist);
+
+	LOOP(i, nlist) {
+		INIT(traces[i].tr, t);
+		INIT(traces[i].trinv, t);
+		int inv_id_in_list = index_in_list(idlist, nlist, group->elements[idlist[i]].inverse->id);
+		mat_trace(traces[i].tr, matrices[i]);
+		mat_trace(traces[i].trinv, matrices[inv_id_in_list]);
+		traces[i].id = idlist[i];
+		distinct_traces[i] = &traces[i];
+	}
+
+	if(unique) {
+		qsort(distinct_traces, nlist, sizeof(struct tracedata*), compare_tracedata);
+
+		nuniq = 0;
+		LOOP(i, nlist) {
+			if(i == 0 || compare_tracedata(&distinct_traces[i], &distinct_traces[nuniq-1]) != 0) {
+				distinct_traces[nuniq] = distinct_traces[i];
+				nuniq++;
+			} else {
+				int oldlength = group->elements[distinct_traces[nuniq-1]->id].length;
+				int newlength = group->elements[distinct_traces[i]->id].length;
+				if(newlength < oldlength)
+					distinct_traces[nuniq-1]->id = distinct_traces[i]->id;
+			}
+		}
+	} else {
+		nuniq = nlist;
+	}
+
+	qsort(distinct_traces, nuniq, sizeof(struct tracedata*), compare_tracedata_id);
+
+	struct tracedata *unique_traces = malloc(nuniq*sizeof(struct tracedata));
+	LOOP(i, nuniq) {
+		INIT(unique_traces[i].tr, t);
+		INIT(unique_traces[i].trinv, t);
+		SET(unique_traces[i].tr, distinct_traces[i]->tr);
+		SET(unique_traces[i].trinv, distinct_traces[i]->trinv);
+		unique_traces[i].id = distinct_traces[i]->id;
+	}
+
+	LOOP(i, nlist) mat_clear(matrices[i]);
+	LOOP(i, nlist) {
+		CLEAR(traces[i].tr);
+		CLEAR(traces[i].trinv);
+	}
+	free(matrices);
+	free(traces);
+	free(distinct_traces);
+
+	*traces_out = unique_traces;
+	return nuniq;
+}
+
+void enumerate_tracedata_clear(struct tracedata *traces, int n)
+{
+	LOOP(i, n) {
+		CLEAR(traces[i].tr);
+		CLEAR(traces[i].trinv);
+	}
+	free(traces);
+}
+
+// returns squares of absolute values
+int solve_characteristic_polynomial_d(mps_context *solv, mps_monomial_poly *poly, double tr_real, double tr_imag, double trinv_real, double trinv_imag, double *eigenvalues_real, double *eigenvalues_imag)
+{
+//        mpq_t coeff1, coeff2, zero;
+        cplx_t *roots;
+        double radii[3];
+        double *radii_p[3];
+        mps_boolean errors;
+        int result = 0;
+
+        /*
+        mpq_inits(coeff1, coeff2, zero, NULL);
+        mpq_set(coeff1, trinv);
+        mpq_sub(coeff2, zero, tr);
+        */
+
+        mps_monomial_poly_set_coefficient_d(solv, poly, 0, -1, 0);
+        mps_monomial_poly_set_coefficient_d(solv, poly, 1, trinv_real, trinv_imag);
+        mps_monomial_poly_set_coefficient_d(solv, poly, 2, -tr_real, -tr_imag);
+//        mps_monomial_poly_set_coefficient_q(solv, poly, 1, coeff1, zero);
+//        mps_monomial_poly_set_coefficient_q(solv, poly, 2, coeff2, zero);
+        mps_monomial_poly_set_coefficient_d(solv, poly, 3, 1, 0);
+
+        mps_context_set_input_poly(solv, (mps_polynomial*)poly);
+        mps_mpsolve(solv);
+
+        roots = cplx_valloc(3);
+        for(int i = 0; i < 3; i++)
+                radii_p[i] = &(radii[i]);
+        mps_context_get_roots_d(solv, &roots, radii_p);
+        errors = mps_context_has_errors(solv);
+
+        if(errors) {
+                result = 1;
+        } else {
+                for(int i = 0; i < 3; i++) {
+                        eigenvalues_real[i] = cplx_Re(roots[i]);
+                        eigenvalues_imag[i] = cplx_Im(roots[i]);
+//                        if(fabs(cplx_Im(roots[i])) > radii[i]) // non-real root
+//                                result = 2;
+                }
+        }
+
+        cplx_vfree(roots);
+//        mpq_clears(coeff1, coeff2, zero, NULL);
+
+        return result;
+}

enumerate.h

@@ -0,0 +1,22 @@
+#ifndef ENUMERATE_H
+#define ENUMERATE_H
+
+#include "mat.h"
+#include "coxeter.h"
+
+#include <mps/mps.h>
+
+struct tracedata {
+	int id;
+	NUMBER tr;
+	NUMBER trinv;
+};
+
+void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices, const int *idlist, int nlist);
+int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out, const int *idlist, int nlist, int unique);
+
+void enumerate_tracedata_clear(struct tracedata *traces, int n);
+
+int solve_characteristic_polynomial_d(mps_context *solv, mps_monomial_poly *poly, double tr_real, double tr_imag, double trinv_real, double trinv_imag, double *eigenvalues_real, double *eigenvalues_imag);
+
+#endif

generators.c

@@ -0,0 +1,148 @@
+#include "generators.h"
+#include "mat.h"
+
+#define LOOP(i,n) for(int i = 0; i < (n); i++)
+
+void generators_triangle_reflection_generic(mat *gen, NUMBER sqrho1, NUMBER sqrho2, NUMBER sqrho3, NUMBER t)
+{
+	NUMBER tmp;
+	NUMBER tinv;
+	TYPE type = GETTYPE(sqrho1);
+
+	INIT(tmp, type);
+	INIT(tinv, type);
+	INV(tinv, t);
+
+	LOOP(i, 3) {
+		mat_init(gen[i], 3, type);
+		mat_zero(gen[i]);
+		LOOP(j, 3) {
+			SET_INT(*mat_ref(gen[i], j, j), i == j ? 1 : -1);
+		}
+	}
+
+	MUL(tmp, t, sqrho1);
+	NEG(*mat_ref(gen[2], 1, 2), tmp);
+	NEG(*mat_ref(gen[0], 2, 0), sqrho2);
+	NEG(*mat_ref(gen[1], 0, 1), sqrho3);
+
+	MUL(tmp, tinv, sqrho1);
+	NEG(*mat_ref(gen[1], 2, 1), tmp);
+	NEG(*mat_ref(gen[2], 0, 2), sqrho2);
+	NEG(*mat_ref(gen[0], 1, 0), sqrho3);
+
+	CLEAR(tmp);
+	CLEAR(tinv);
+}
+
+// warning: this is not compatible anymore!
+int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t t)
+{
+	int p[3] = {p1, p2, p3};
+	int q[3] = {q1, q2, q3};
+	TYPE type;
+
+	LOOP(i, 3)
+		if(p[i] < 2 || p[i] > 6)
+			return 0;
+
+#ifdef QEXT
+	type = QT_TRIVIAL;
+	LOOP(i, 3)
+		if(p[i] == 5)
+			type = QT_SQRT5;
+#else
+	LOOP(i, 3)
+		if(p[i] == 5)
+			return 0;
+#endif
+
+
+	NUMBER rho[3];
+	LOOP(i, 3) INIT(rho[i], type);
+
+	// compute rho1, rho2, rho3
+	LOOP(i, 3) {
+		if(p[i] == 2 && q[i] == 1) {
+			SET_INT(rho[i], 0);
+		} else if(p[i] == 3 && q[i] == 1) {
+			SET_INT(rho[i], 1);
+		} else if(p[i] == 4 && q[i] == 1) {
+			SET_INT(rho[i], 2);
+		} else if(p[i] == 4 && q[i] == 2) {
+			SET_INT(rho[i], 0);
+		} else if(p[i] == 6 && q[i] == 1) {
+			SET_INT(rho[i], 3);
+		} else if(p[i] == 6 && q[i] == 2) {
+			SET_INT(rho[i], 1);
+		} else if(p[i] == 6 && q[i] == 3) {
+			SET_INT(rho[i], 0);
+#ifdef QEXT
+		} else if(p[i] == 5 && q[i] == 1) {
+			mpq_set_si(rho[i]->a[0], 3, 2);
+			mpq_set_si(rho[i]->a[1], 1, 2);
+		} else if(p[i] == 5 && q[i] == 2) {
+			mpq_set_si(rho[i]->a[0], 3, 2);
+			mpq_set_si(rho[i]->a[1], -1, 2);
+#endif
+		} else {
+			return 0;
+		}
+	}
+
+	NUMBER param;
+	INIT(param, type);
+	SET_Q(param, t);
+
+	generators_triangle_reflection_generic(gen, rho[0], rho[1], rho[2], param);
+
+	LOOP(i, 3) CLEAR(rho[i]);
+	CLEAR(param);
+
+	return 1;
+}
+
+int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, int q3, mpq_t treal, mpq_t timag)
+{
+
+	int q[3] = {q1, q2, q3};
+
+	NUMBER sqrho[3];
+	LOOP(i, 3) INIT(sqrho[i], QT_GAUSS_SQRT5);
+
+	// compute sqrho1, sqrho2, sqrho3
+	LOOP(i, 3) {
+		if(q[i] == 1) {
+			// 2cos(pi/5) = 1/2 + 1/2*sqrt(5) = 1/2 + 7/12 a - 1/24 a^3
+			mpq_set_si(sqrho[i]->a[0], 1, 2);
+			mpq_set_si(sqrho[i]->a[1], 7, 12);
+			mpq_set_si(sqrho[i]->a[2], 0, 1);
+			mpq_set_si(sqrho[i]->a[3], -1, 24);
+		} else if(q[i] == 2) {
+			// 2cos(pi/5) = -1/2 + 1/2*sqrt(5)
+			mpq_set_si(sqrho[i]->a[0], -1, 2);
+			mpq_set_si(sqrho[i]->a[1], 7, 12);
+			mpq_set_si(sqrho[i]->a[2], 0, 1);
+			mpq_set_si(sqrho[i]->a[3], -1, 24);
+		} else {
+			return 0;
+		}
+	}
+
+	NUMBER param;
+	INIT(param, QT_GAUSS_SQRT5);
+	mpq_set(param->a[0], treal);
+	mpq_set_si(param->a[1], -1, 6);
+	mpq_mul(param->a[1], param->a[1], timag);
+	mpq_set_si(param->a[2], 0, 1);
+	mpq_set_si(param->a[3], 1, 12);
+	mpq_mul(param->a[3], param->a[3], timag);
+
+	generators_triangle_reflection_generic(gen, sqrho[0], sqrho[1], sqrho[2], param);
+
+	LOOP(i, 3) CLEAR(sqrho[i]);
+	CLEAR(param);
+
+	return 1;
+
+}

generators.h

@@ -0,0 +1,13 @@
+#ifndef GENERATORS_H
+#define GENERATORS_H
+
+#include "mat.h"
+
+void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER q);
+int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t q);
+
+#ifdef QEXT
+int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, int q3, mpq_t real, mpq_t imag);
+#endif
+
+#endif

make_picture.py

@@ -0,0 +1,67 @@
+#!/usr/bin/python
+
+import sys
+import math
+
+def sharpen(x):
+	alpha = 1
+	y = abs(2*x-1)**alpha # between 0 and 1
+	if x > 0.499 and x < 0.501:
+		return 0.5
+	elif x > 0.5:
+		return (y+1)/2
+	else:
+		return (1-y)/2
+
+f = open(sys.argv[1])
+lines = [x.split() for x in f]
+f.close()
+
+res1 = int(sys.argv[2]) # 400 pixel per unit
+res2 = int(sys.argv[3]) # 50 pixel in picture (acutally one quadrant)
+
+data = {(round(float(l[0])*res1), round(float(l[1])*res1)) : float(l[5]) for l in lines}
+data[(0,0)] = 2.0
+
+print("P3")
+print("1000 1000")
+print("255")
+
+for i in range (-500,500):
+	for j in range(-500,500):
+		x = j/500*res2
+		y = i/500*res2
+
+		if y < 0:
+			y = -y
+
+		x0 = math.floor(x)
+		y0 = math.floor(y)
+		x1 = math.ceil(x)
+		y1 = math.ceil(y)
+		tx = sharpen(x-x0)
+		ty = sharpen(y-y0)
+
+		if not (x0,y0) in data or not (x0,y1) in data or not (x1,y0) in data or not (x1,y1) in data:
+			value = 0
+		else:
+			value = 0.0
+			value += data[(x0,y0)]*(1-tx)*(1-ty)
+			value += data[(x0,y1)]*(1-tx)*ty
+			value += data[(x1,y0)]*tx*(1-ty)
+			value += data[(x1,y1)]*tx*ty
+			value -= 1
+
+		value = 0 if value < 0 else 1 if value > 1 else value
+
+		r = round(255*math.sqrt(value))
+		g = round(255*(value)**3)
+		b = round(255*math.sin(2*math.pi*(value)))
+
+		r = 0 if r < 0 else 255 if r > 255 else r
+		g = 0 if g < 0 else 255 if g > 255 else g
+		b = 0 if b < 0 else 255 if b > 255 else b
+
+		print("%d %d %d" % (r,g,b))
+
+#print(data)

mat.h

@@ -24,7 +24,8 @@
 #define SUB qext_sub
 #define NEG qext_neg
 #define MUL qext_mul
-// #define DIV qext_div
+#define DIV qext_div
+#define INV qext_inv
 #define CMP qext_cmp
 #define PRINT qext_print
 #define SNPRINT qext_snprint
@@ -46,6 +47,7 @@
 #define NEG mpq_neg
 #define MUL mpq_mul
 #define DIV mpq_div
+#define INV mpq_inv
 #define PRINT(x) gmp_printf("%Qd", x)
 #define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x)
 #define TYPE int

qext.c

@@ -1,5 +1,7 @@
 #include "qext.h"
 
+#include <assert.h>
+
 int qext_rank;
 mpq_t *qext_coefficient;
 
@@ -9,13 +11,20 @@ mpq_t *qext_coefficient;
 #define LOOP(i,n) for(int i = 0; i < (n); i++)
 #define RANK(x) ((x)->type->rank)
 
-const int qext_type_trivial_coeffs[] = {-1, 1};
-struct qext_type qext_type_trivial_v = {1, qext_type_trivial_coeffs, NULL};
-struct qext_type *QT_TRIVIAL = &qext_type_trivial_v;
+struct qext_type *QT_TRIVIAL = &(struct qext_type) {
+	.rank = 1,
+	.integer_coeffs = (const int[]) {-1, 1}
+};
 
-const int qext_type_sqrt5_coeffs[] = {-5, 0, 1};
-struct qext_type qext_type_sqrt5_v = {2, qext_type_sqrt5_coeffs, NULL};
-struct qext_type *QT_SQRT5 = &qext_type_sqrt5_v;
+struct qext_type *QT_SQRT5 = &(struct qext_type) {
+	.rank = 2,
+	.integer_coeffs = (const int[]) {-5, 0, 1}
+};
+
+struct qext_type *QT_GAUSS_SQRT5 = &(struct qext_type) {
+	.rank = 4,
+	.integer_coeffs = (const int[]) {36, 0, -8, 0, 1}
+};
 
 static void qext_init_type(struct qext_type *type)
 {
@@ -161,7 +170,86 @@ void qext_mul(qext_number out, qext_number x, qext_number y)
 	mpq_clear(tmp);
 }
 
+void qext_inv(qext_number out, qext_number in)
+{
+	qext_number one;
+	qext_init(one, in->type);
+	qext_set_int(one, 1);
+	qext_div(out, one, in);
+	qext_clear(one);
+}
+
 void qext_div(qext_number out, qext_number x, qext_number y)
 {
-	// todo: implement this by solving a linear system
+	int rk = RANK(x);
+	struct qext_type *type = x->type;
+
+	mpq_t tmp;
+	mpq_t result[rk];
+	mpq_t matrix[rk*rk];
+
+	mpq_init(tmp);
+	LOOP(i, rk) mpq_init(result[i]);
+	LOOP(i, rk*rk) mpq_init(matrix[i]);
+
+	// initialize a matrix which represents multiplication by y
+	// that means the i-th column is y*a^i
+	LOOP(i, rk)
+		mpq_set(matrix[i*rk], y->a[i]);
+	LOOP(j, rk-1) { // columns
+		mpq_mul(matrix[j+1], matrix[(rk-1)*rk+j], type->coeffs[0]);
+		LOOP(i, rk-1) { // rows
+			mpq_mul(matrix[(i+1)*rk+(j+1)], matrix[(rk-1)*rk+j], type->coeffs[i+1]);
+			mpq_add(matrix[(i+1)*rk+(j+1)], matrix[(i+1)*rk+(j+1)], matrix[i*rk+j]);
+		}
+	}
+
+	// initialize result as x
+	LOOP(i, rk) mpq_set(result[i], x->a[i]);
+
+	// use Gaussian elimination to solve the system matrix * out = result
+	LOOP(j, rk)
+	{
+		// find nonzero entry
+		int row = j;
+		while(row < rk && mpq_sgn(matrix[row*rk+j]) == 0) row++;
+
+		// if the entire column is zero, the matrix was not invertible
+		// this could happen if the polynomial is not irreducible / we're not in a field
+		assert(row != rk);
+
+		// permute rows
+		if(row != j) {
+			mpq_swap(result[row], result[j]);
+			LOOP(k, rk) {
+				mpq_swap(matrix[row*rk+k], matrix[j*rk+k]);
+			}
+		}
+
+		// normalize
+		LOOP(k, rk)
+			if(k != j)
+				mpq_div(matrix[j*rk+k], matrix[j*rk+k], matrix[j*rk+j]);
+		mpq_div(result[j], result[j], matrix[j*rk+j]);
+		mpq_set_ui(matrix[j*rk+j], 1, 1);
+
+		// subtract
+		LOOP(i, rk) {
+			if(i == j)
+				continue;
+			mpq_mul(tmp, matrix[i*rk+j], result[j]);
+			mpq_sub(result[i], result[i], tmp);
+			for(int k = j+1; k < rk; k++) {
+				mpq_mul(tmp, matrix[i*rk+j], matrix[j*rk+k]);
+				mpq_sub(matrix[i*rk+k], matrix[i*rk+k], tmp);
+			}
+			mpq_set_ui(matrix[i*rk+j], 0, 1); // it isn't strictly necessary to do this as we won't use this entry again, but helps debugging
+		}
+	}
+
+	LOOP(i, rk) mpq_set(out->a[i], result[i]);
+
+	mpq_clear(tmp);
+	LOOP(i, rk) mpq_clear(result[i]);
+	LOOP(i, rk*rk) mpq_clear(matrix[i]);
 }

qext.h

@@ -5,8 +5,8 @@
 
 struct qext_type {
 	int rank;
-	const int *integer_coeffs;
-	mpq_t *coeffs;
+	const int *integer_coeffs;  // actually the rank+1 coefficients of the polynomial
+	mpq_t *coeffs;  // components of a^rank
 };
 
 struct qext_number_internal {
@@ -18,6 +18,7 @@ typedef struct qext_number_internal qext_number[1];
 
 extern struct qext_type *QT_TRIVIAL;
 extern struct qext_type *QT_SQRT5;
+extern struct qext_type *QT_GAUSS_SQRT5;
 
 struct qext_type *qext_newtype(int rank, const int *coeffs);
 void qext_init(qext_number x, struct qext_type *type);
@@ -29,6 +30,7 @@ void qext_add(qext_number result, qext_number x, qext_number y);
 void qext_sub(qext_number result, qext_number x, qext_number y);
 void qext_neg(qext_number result, qext_number x);
 void qext_mul(qext_number out, qext_number x, qext_number y);
+void qext_inv(qext_number out, qext_number in);
 void qext_div(qext_number out, qext_number x, qext_number y);
 int qext_cmp(qext_number x, qext_number y);
 void qext_print(qext_number x);