From 67e5b30808881924cbb1a5377beb34567b05768b Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sun, 3 Oct 2021 19:18:10 -0500 Subject: [PATCH] compute character variety picture based only on billiard words --- .gitignore | 12 ++ Makefile | 13 +- billiard_words.hs | 31 ++++ .../max_slope_highres.plt | 15 ++ max_slope_picture_billiards/words | 132 ++++++++++++++++++ special_element.c | 116 +++++++++------ 6 files changed, 268 insertions(+), 51 deletions(-) create mode 100644 .gitignore create mode 100644 billiard_words.hs create mode 100644 max_slope_picture_billiards/max_slope_highres.plt create mode 100644 max_slope_picture_billiards/words diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c6e4cb2 --- /dev/null +++ b/.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 diff --git a/Makefile b/Makefile index b045eeb..0ff4f43 100644 --- a/Makefile +++ b/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 -flto -funroll-loops -Winline +SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline +#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 diff --git a/billiard_words.hs b/billiard_words.hs new file mode 100644 index 0000000..92ef72c --- /dev/null +++ b/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" diff --git a/max_slope_picture_billiards/max_slope_highres.plt b/max_slope_picture_billiards/max_slope_highres.plt new file mode 100644 index 0000000..7a3d237 --- /dev/null +++ b/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 diff --git a/max_slope_picture_billiards/words b/max_slope_picture_billiards/words new file mode 100644 index 0000000..52d73c5 --- /dev/null +++ b/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 + diff --git a/special_element.c b/special_element.c index e4ed229..84c382c 100644 --- a/special_element.c +++ b/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 j = 1; j <= 100; j++) { - mpq_set_ui(t, j, 100); - mpq_set_ui(m, i, 100); // 414/1000 ~ sqrt(2)-1 <-> s=1 - s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m)); + for(int i = STARTX; i <= WIDTH; i++) { + for(int j = 1; j <= HEIGHT; j++) { + for(int w = 1; w < argc; w++) { + mpq_set_ui(t, j, DENOMINATOR); + 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"); - initialize_triangle_generators(ws, gen, m, t); + DEBUG("Compute matrix\n"); + initialize_triangle_generators(ws, gen, m, t); - mat_identity(element); - mat_identity(inverse); - for(int i = 0; i < strlen(argv[1]); i+=2) { - letter1 = argv[1][i] - 'a'; - letter2 = argv[1][i+1] - 'a'; + mat_identity(element); + mat_identity(inverse); + for(int k = 0; k < strlen(argv[w]); k+=2) { + letter1 = argv[w][k] - 'a'; + letter2 = argv[w][k+1] - 'a'; - if(letter1 == 1 && letter2 == 2) - letter = 0; // p = bc - else if(letter1 == 2 && letter2 == 0) - letter = 1; // q = ca - else if(letter1 == 0 && letter2 == 1) - letter = 2; // r = ab - else if(letter1 == 2 && letter2 == 1) - letter = 3; // p^{-1} = cb - else if(letter1 == 0 && letter2 == 2) - letter = 4; // q^{-1} = ac - else if(letter1 == 1 && letter2 == 0) - letter = 5; // r^{-1} = ba + if(letter1 == 1 && letter2 == 2) + letter = 0; // p = bc + else if(letter1 == 2 && letter2 == 0) + letter = 1; // q = ca + else if(letter1 == 0 && letter2 == 1) + letter = 2; // r = ab + else if(letter1 == 2 && letter2 == 1) + letter = 3; // p^{-1} = cb + else if(letter1 == 0 && letter2 == 2) + letter = 4; // q^{-1} = ac + else if(letter1 == 1 && letter2 == 0) + letter = 5; // r^{-1} = ba - mat_multiply(ws, element, element, gen[letter]); - mat_multiply(ws, inverse, gen[(letter+3)%6], inverse); - } + mat_multiply(ws, element, element, gen[letter]); + mat_multiply(ws, inverse, gen[(letter+3)%6], inverse); + } - DEBUG("Compute traces\n"); + DEBUG("Compute traces\n"); - mat_trace(tr, element); - mat_trace(trinv, inverse); + mat_trace(tr, element); + mat_trace(trinv, inverse); - DEBUG("Solve characteristic polynomials\n"); - retval = solve_characteristic_polynomial(solver, tr, trinv, evs); - if(retval == 1) { - fprintf(stderr, "Error! Could not solve polynomial.\n"); - return 1; - } + DEBUG("Solve characteristic polynomials\n"); + retval = solve_characteristic_polynomial(solver, tr, trinv, evs); + if(retval == 1) { + fprintf(stderr, "Error! Could not solve polynomial.\n"); + return 1; + } - if(fabs(evs[0]) < fabs(evs[1])) - SWAP(double, evs[0], evs[1]); - if(fabs(evs[1]) < fabs(evs[2])) - SWAP(double, evs[1], evs[2]); - if(fabs(evs[0]) < fabs(evs[1])) - SWAP(double, evs[0], evs[1]); + if(fabs(evs[0]) < fabs(evs[1])) + SWAP(double, evs[0], evs[1]); + if(fabs(evs[1]) < fabs(evs[2])) + SWAP(double, evs[1], evs[2]); + if(fabs(evs[0]) < fabs(evs[1])) + SWAP(double, evs[0], evs[1]); - x = log(fabs(evs[0])); - y = -log(fabs(evs[2])); + x = log(fabs(evs[0])); + 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); }