compute character variety picture based only on billiard words

.gitignore

@@ -0,0 +1,12 @@
+*.o
+triangle_group/singular_values
+.#*
+singular_values
+output/
+special_element
+max_slope_picture/generate
+convert
+billiard_words
+*.pnm
+*.png
+*.hi

Makefile

@@ -1,23 +1,26 @@
 HEADERS=linalg.h mat.h coxeter.h
 
 #SPECIAL_OPTIONS=-O0 -g -D_DEBUG
-#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
+SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
-SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
+#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 #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)
 
-all: singular_values special_element convert
+all: singular_values special_element convert billiard_words
 
 convert: convert.hs
 	ghc --make -dynamic convert.hs
 
+billiard_words: billiard_words.hs
+	ghc --make -dynamic billiard_words.hs
+
 singular_values: singular_values.o coxeter.o mat.o
 	mpicc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o -lm -lgmp -lmps
 
 special_element: special_element.o coxeter.o linalg.o mat.o
-	gcc $(OPTIONS) -o special_element coxeter.o linalg.o special_element.o mat.o -lm -lgmp -lmps
+	gcc $(OPTIONS) -o special_element coxeter.o linalg.o special_element.o mat.o -lm -lgmp -lmps -lgsl -lcblas
 
 singular_values.o: singular_values.c $(HEADERS)
 	mpicc $(OPTIONS) -c singular_values.c
@@ -35,4 +38,4 @@ mat.o: mat.c $(HEADERS)
 	gcc $(OPTIONS) -c mat.c
 
 clean:
-	rm -f singular_values special_element coxeter.o linalg.o singular_values.o mat.o special_element.o convert.hi convert.o convert
+	rm -f singular_values special_element coxeter.o linalg.o singular_values.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words

billiard_words.hs

@@ -0,0 +1,31 @@
+import Data.List
+import Data.Ord
+import Text.Printf
+
+main = do
+  putStrLn $ unlines [printf "%d/%d\t%.5f\t%.5f\t%d\t%s" p q (fromIntegral p / fromIntegral q :: Double) (sqrt 3 / (1 + 2*fromIntegral q/fromIntegral p) :: Double) (length w) w | ((p,q),w) <- wordlist, length w <= 100]
+
+wordlist = nub $ sortBy (comparing sl) [((p `div` gcd p q, q `div` gcd p q), slopeWord p q) | p <- [0..50], q <- [1..50], p <= q]
+    where
+      sl ((p,q),_) = fromIntegral p / fromIntegral q
+
+-- only allows slopes 0 <= p/q <= 1
+slopeWord :: Int -> Int -> String
+slopeWord p q =
+    concat $ map word $ zipWith step list (tail list)
+    where
+      p_ = p `div` gcd p q
+      q_ = q `div` gcd p q
+      xmax = if (p_ - q_) `mod` 3 == 0 then q_ else 3*q_
+      list = [(x,y) | x <- ([0..xmax] :: [Int]), let y = floor (fromIntegral x*fromIntegral p/fromIntegral q)]
+
+      step :: (Int,Int) -> (Int,Int) -> (Int,Int)
+      step (x1,y1) (x2,y2) = ((x1-y1) `mod` 3, y2 - y1)
+
+      word :: (Int,Int) -> String
+      word (0,0) = "bc"
+      word (1,0) = "ab"
+      word (2,0) = "ca"
+      word (0,1) = "baca"
+      word (1,1) = "acbc"
+      word (2,1) = "cbab"

max_slope_picture_billiards/max_slope_highres.plt

@@ -0,0 +1,15 @@
+if(!exists("index")) index = 50
+
+set xrange [0:0.45]
+set yrange [0:1]
+set xyplane at 0
+#plot "max_slopes_billiard" using 1:2:($3 - 10*floor($3/10)) w p pt 7 ps 1.1 lc palette t ""
+plot "max_slopes_billiard" using 1:2:3 w p pt 7 ps 1.1 lc palette t ""
+#splot "max_slopes_billiard" using 1:2:3 w p pt 7 ps 0.3 t ""
+#plot "max_slopes_billiard" using 1:2:4 w p pt 7 ps 1.1 lc palette t ""
+
+pause 10
+#if(MOUSE_KEY == 60) index=index-1
+#if(MOUSE_KEY == 62) index=index+1
+#print MOUSE_KEY
+reread

max_slope_picture_billiards/words

@@ -0,0 +1,132 @@
+bcabca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcabcabcacbcabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabcacbcabcabcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabcacbcabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcacbcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabcacbcabcabcabcabcabcabcabcbabcabcabcabcabcabcabcabacabcabcabcabcabcabcabaca
+bcabcabcabcabcabcabaca
+bcabcabcabcabcabcabacabcabcabcabcabcabcbabcabcabcabcabcabcacbcabcabcabcabcabcabaca
+bcabcabcabcabcabcbabcabcabcabcabcabcacbcabcabcabcabcabcabaca
+bcabcabcabcabcabcbabcabcabcabcabcabcacbcabcabcabcabcabcabacabcabcabcabcabcabcbabcabcabcabcabcabaca
+bcabcabcabcabcabcbabcabcabcabcabcabaca
+bcabcabcabcabcabcbabcabcabcabcabcabacabcabcabcabcabcacbcabcabcabcabcabcbabcabcabcabcabcabaca
+bcabcabcabcabcacbcabcabcabcabcabcbabcabcabcabcabcabaca
+bcabcabcabcabcacbcabcabcabcabcabcbabcabcabcabcabcabacabcabcabcabcabaca
+bcabcabcabcabcacbcabcabcabcabcabcbabcabcabcabcabcbabcabcabcabcabcabacabcabcabcabcabaca
+bcabcabcabcabaca
+bcabcabcabcabacabcabcabcabcbabcabcabcabcabcbabcabcabcabcacbcabcabcabcabcacbcabcabcabcabaca
+bcabcabcabcabacabcabcabcabcbabcabcabcabcabcbabcabcabcabcacbcabcabcabcabaca
+bcabcabcabcabacabcabcabcabcbabcabcabcabcacbcabcabcabcabaca
+bcabcabcabcabacabcabcabcabcbabcabcabcabcacbcabcabcabcabacabcabcabcabcbabcabcabcabcacbcabcabcabcabaca
+bcabcabcabcbabcabcabcabcacbcabcabcabcabaca
+bcabcabcabcbabcabcabcabcacbcabcabcabcabacabcabcabcabcbabcabcabcabaca
+bcabcabcabcbabcabcabcabcacbcabcabcabcabacabcabcabcacbcabcabcabcabacabcabcabcabcbabcabcabcabaca
+bcabcabcabcbabcabcabcabaca
+bcabcabcabcbabcabcabcabacabcabcabcacbcabcabcabcabacabcabcabcacbcabcabcabcbabcabcabcabaca
+bcabcabcabcbabcabcabcabacabcabcabcacbcabcabcabcbabcabcabcabaca
+bcabcabcabcbabcabcabcabacabcabcabcacbcabcabcabcbabcabcabcabacabcabcabcacbcabcabcabcbabcabcabcabaca
+bcabcabcacbcabcabcabcbabcabcabcabaca
+bcabcabcacbcabcabcabcbabcabcabcabacabcabcabcacbcabcabcabcbabcabcabcabacabcabcabaca
+bcabcabcacbcabcabcabcbabcabcabcabacabcabcabaca
+bcabcabcacbcabcabcabcbabcabcabcbabcabcabcabacabcabcabaca
+bcabcabcacbcabcabcacbcabcabcabcbabcabcabcbabcabcabcabacabcabcabaca
+bcabcabcacbcabcabcacbcabcabcabcbabcabcabcbabcabcabcabacabcabcabacabcabcabaca
+bcabcabcacbcabcabcacbcabcabcabcbabcabcabcbabcabcabcbabcabcabcabacabcabcabacabcabcabaca
+bcabcabcacbcabcabcacbcabcabcacbcabcabcabcbabcabcabcbabcabcabcbabcabcabcabacabcabcabacabcabcabaca
+bcabcabaca
+bcabcabacabcabcabacabcabcabacabcabcbabcabcabcbabcabcabcbabcabcacbcabcabcacbcabcabcacbcabcabaca
+bcabcabacabcabcabacabcabcbabcabcabcbabcabcabcbabcabcacbcabcabcacbcabcabcacbcabcabaca
+bcabcabacabcabcabacabcabcbabcabcabcbabcabcabcbabcabcacbcabcabcacbcabcabaca
+bcabcabacabcabcabacabcabcbabcabcabcbabcabcacbcabcabcacbcabcabaca
+bcabcabacabcabcbabcabcabcbabcabcacbcabcabcacbcabcabaca
+bcabcabacabcabcbabcabcabcbabcabcacbcabcabcacbcabcabacabcabcabacabcabcbabcabcabcbabcabcacbcabcabaca
+bcabcabacabcabcbabcabcabcbabcabcacbcabcabaca
+bcabcabacabcabcbabcabcacbcabcabcacbcabcabacabcabcbabcabcabcbabcabcacbcabcabaca
+bcabcabacabcabcbabcabcacbcabcabaca
+bcabcabacabcabcbabcabcacbcabcabacabcabcbabcabcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabaca
+bcabcabacabcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabaca
+bcabcabacabcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabaca
+bcabcbabcabcacbcabcabaca
+bcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcabacabcabcbabcabcacbcabcabacabcacbcabcabacabcabcbabcabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcabacabcacbcabcabacabcabcbabcabcacbcabcbabcabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcabacabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcbabcabcacbcabcabacabcacbcabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcbabcabcacbcabcabacabcacbcabcabacabcabcbabcabacabcabcbabcabaca
+bcabcbabcabcacbcabcbabcabcacbcabcabacabcacbcabcabacabcacbcabcabacabcabcbabcabacabcabcbabcabaca
+bcabcbabcabaca
+bcabcbabcabacabcabcbabcabacabcacbcabcabacabcacbcabcabacabcacbcabcbabcabcacbcabcbabcabaca
+bcabcbabcabacabcabcbabcabacabcacbcabcabacabcacbcabcbabcabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcabacabcacbcabcbabcabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcabacabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcbabcabcacbcabcbabcabacabcacbcabcabacabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcbabcabacabcacbcabcabacabcacbcabcbabcabacabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabaca
+bcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabaca
+bcacbcabcbabcabaca
+bcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabcacbcabcbabcabacabcacbcabcbabcbabcabacabcacbcabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabcacbcabcbabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabcacbcacbcabcbabcabacabcacbcabcbabcbabcabacabcacbcabcbabcabacabaca
+bcacbcabcbabcabacabaca
+bcacbcabcbabcabacabacabcacbcabcbabcbabcabacabcacbcabcbabcbabcabacabcacbcacbcabcbabcabacabaca
+bcacbcabcbabcabacabacabcacbcabcbabcbabcabacabcacbcacbcabcbabcabacabaca
+bcacbcabcbabcbabcabacabcacbcacbcabcbabcabacabaca
+bcacbcabcbabcbabcabacabcacbcacbcabcbabcabacabacabcacbcabcbabcbabcabacabaca
+bcacbcabcbabcbabcabacabcacbcacbcabcbabcabacabacabcacbcacbcabcbabcabacabacabcacbcabcbabcbabcabacabaca
+bcacbcabcbabcbabcabacabaca
+bcacbcabcbabcbabcabacabacabcacbcacbcabcbabcabacabacabcacbcacbcabcbabcbabcabacabaca
+bcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabaca
+bcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabaca
+bcacbcacbcabcbabcbabcabacabaca
+bcacbcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabacabaca
+bcacbcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabacabaca
+bcacbcacbcabcbabcbabcabacabacabcacbcacbcabcbabcbabcbabcabacabacabcacbcacbcabcbabcbabcabacabacabaca
+bcacbcacbcabcbabcbabcabacabacabaca
+bcacbcacbcabcbabcbabcbabcabacabacabcacbcacbcacbcabcbabcbabcabacabacabaca
+bcacbcacbcabcbabcbabcbabcabacabacabaca
+bcacbcacbcabcbabcbabcbabcabacabacabacabcacbcacbcacbcabcbabcbabcbabcabacabacabaca
+bcacbcacbcacbcabcbabcbabcbabcabacabacabaca
+bcacbcacbcacbcabcbabcbabcbabcabacabacabacabcacbcacbcacbcabcbabcbabcbabcabacabacabacabaca
+bcacbcacbcacbcabcbabcbabcbabcabacabacabacabaca
+bcacbcacbcacbcabcbabcbabcbabcbabcabacabacabacabcacbcacbcacbcacbcabcbabcbabcbabcabacabacabacabaca
+bcacbcacbcacbcabcbabcbabcbabcbabcabacabacabacabaca
+bcacbcacbcacbcacbcabcbabcbabcbabcbabcabacabacabacabaca
+bcacbcacbcacbcacbcabcbabcbabcbabcbabcabacabacabacabacabaca
+bcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabacabacabaca
+bcacbcacbcacbcacbcacbcacbcacbcabcbabcbabcbabcbabcbabcbabcbabcbabcabacabacabacabacabacabacabacabaca
+baca
+

special_element.c

@@ -8,6 +8,11 @@
 
 #define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
 #define DEBUG(msg, ...)
+#define DENOMINATOR 300
+#define WIDTH 135
+#define STARTX 121
+#define HEIGHT 300
+#define IDX(i,j) (((i)-1)*HEIGHT + ((j)-1))
 
 int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues)
 {
@@ -189,10 +194,12 @@ int main(int argc, char *argv[])
 	mat element, inverse;
 	int letter1, letter2, letter;
 	mpq_t tr, trinv;
-	double x, y;
+	double x, y, slope;
 	int retval;
 	double evs[3];
 	char buf[100];
+	double *max_slope;
+	int *max_slope_index;
 
 	DEBUG("Allocate\n");
 
@@ -202,68 +209,83 @@ int main(int argc, char *argv[])
 		mat_init(gen[i], 3);
 	mat_init(element, 3);
 	mat_init(inverse, 3);
+	max_slope = malloc(sizeof(double)*WIDTH*HEIGHT);
+	max_slope_index = malloc(sizeof(int)*WIDTH*HEIGHT);
+	memset(max_slope_index, 0, sizeof(int)*WIDTH*HEIGHT);
+	memset(max_slope, 0, sizeof(int)*WIDTH*HEIGHT);
 
 	solver = mps_context_new();
 	mps_context_set_output_prec(solver, 20); // relative precision
 	mps_context_set_output_goal(solver, MPS_OUTPUT_GOAL_APPROXIMATE);
 
-	for(int i = 1; i <= 99; i++) {
+	for(int i = STARTX; i <= WIDTH; i++) {
-		for(int j = 1; j <= 100; j++) {
+		for(int j = 1; j <= HEIGHT; j++) {
-			mpq_set_ui(t, j, 100);
+			for(int w = 1; w < argc; w++) {
-			mpq_set_ui(m, i, 100); // 414/1000 ~ sqrt(2)-1 <-> s=1
+				mpq_set_ui(t, j, DENOMINATOR);
-			s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m));
+				mpq_set_ui(m, i, DENOMINATOR); // 414/1000 ~ sqrt(2)-1 <-> s=1
+				s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m));
 
-			DEBUG("Compute matrix\n");
+				DEBUG("Compute matrix\n");
-			initialize_triangle_generators(ws, gen, m, t);
+				initialize_triangle_generators(ws, gen, m, t);
 
-			mat_identity(element);
+				mat_identity(element);
-			mat_identity(inverse);
+				mat_identity(inverse);
-			for(int i = 0; i < strlen(argv[1]); i+=2) {
+				for(int k = 0; k < strlen(argv[w]); k+=2) {
-				letter1 = argv[1][i] - 'a';
+					letter1 = argv[w][k] - 'a';
-				letter2 = argv[1][i+1] - 'a';
+					letter2 = argv[w][k+1] - 'a';
 
-				if(letter1 == 1 && letter2 == 2)
+					if(letter1 == 1 && letter2 == 2)
-					letter = 0; // p = bc
+						letter = 0; // p = bc
-				else if(letter1 == 2 && letter2 == 0)
+					else if(letter1 == 2 && letter2 == 0)
-					letter = 1; // q = ca
+						letter = 1; // q = ca
-				else if(letter1 == 0 && letter2 == 1)
+					else if(letter1 == 0 && letter2 == 1)
-					letter = 2; // r = ab
+						letter = 2; // r = ab
-				else if(letter1 == 2 && letter2 == 1)
+					else if(letter1 == 2 && letter2 == 1)
-					letter = 3; // p^{-1} = cb
+						letter = 3; // p^{-1} = cb
-				else if(letter1 == 0 && letter2 == 2)
+					else if(letter1 == 0 && letter2 == 2)
-					letter = 4; // q^{-1} = ac
+						letter = 4; // q^{-1} = ac
-				else if(letter1 == 1 && letter2 == 0)
+					else if(letter1 == 1 && letter2 == 0)
-					letter = 5; // r^{-1} = ba
+						letter = 5; // r^{-1} = ba
 
-				mat_multiply(ws, element, element, gen[letter]);
+					mat_multiply(ws, element, element, gen[letter]);
-				mat_multiply(ws, inverse, gen[(letter+3)%6], inverse);
+					mat_multiply(ws, inverse, gen[(letter+3)%6], inverse);
-			}
+				}
 
-			DEBUG("Compute traces\n");
+				DEBUG("Compute traces\n");
 
-			mat_trace(tr, element);
+				mat_trace(tr, element);
-			mat_trace(trinv, inverse);
+				mat_trace(trinv, inverse);
 
-			DEBUG("Solve characteristic polynomials\n");
+				DEBUG("Solve characteristic polynomials\n");
-			retval = solve_characteristic_polynomial(solver, tr, trinv, evs);
+				retval = solve_characteristic_polynomial(solver, tr, trinv, evs);
-			if(retval == 1) {
+				if(retval == 1) {
-				fprintf(stderr, "Error! Could not solve polynomial.\n");
+					fprintf(stderr, "Error! Could not solve polynomial.\n");
-				return 1;
+					return 1;
-			}
+				}
 
-			if(fabs(evs[0]) < fabs(evs[1]))
+				if(fabs(evs[0]) < fabs(evs[1]))
-				SWAP(double, evs[0], evs[1]);
+					SWAP(double, evs[0], evs[1]);
-			if(fabs(evs[1]) < fabs(evs[2]))
+				if(fabs(evs[1]) < fabs(evs[2]))
-				SWAP(double, evs[1], evs[2]);
+					SWAP(double, evs[1], evs[2]);
-			if(fabs(evs[0]) < fabs(evs[1]))
+				if(fabs(evs[0]) < fabs(evs[1]))
-				SWAP(double, evs[0], evs[1]);
+					SWAP(double, evs[0], evs[1]);
 
-			x = log(fabs(evs[0]));
+				x = log(fabs(evs[0]));
-			y = -log(fabs(evs[2]));
+				y = -log(fabs(evs[2]));
+				slope = y/x > 1 ? y/x : x/y;
+
+				if(slope > max_slope[IDX(i,j)]) {
+					max_slope[IDX(i,j)] = slope;
+					max_slope_index[IDX(i,j)] = w;
+				}
 
 //	gmp_printf("%Qd %Qd %f %f %f\n", tr, trinv, x, y, y/x);
-			gmp_printf("%.5f %.5f %.7f %.9f\n", mpq_get_d(t), mpq_get_d(m), s, y/x);
+//				gmp_printf("%.5f %.5f %.7f %.9f\n", mpq_get_d(t), mpq_get_d(m), s, slope);
 
+			}
+
+			printf("%.5f %.5f %d %.9f\n", (double)i/DENOMINATOR, (double)j/DENOMINATOR, max_slope_index[IDX(i,j)], max_slope[IDX(i,j)]);
+			fflush(stdout);
 		}
 	}
 
@@ -275,4 +297,6 @@ int main(int argc, char *argv[])
 	mat_clear(element);
 	mat_clear(inverse);
 	mps_context_free(solver);
+	free(max_slope);
+	free(max_slope_index);
 }