orthogonal words and some comments

billiard_picture.sh

@@ -10,11 +10,20 @@ qdenom=100
 qstart=1
 qend=100  # 1/sqrt(2) = 0.7071...
 
-words="$(./billiard_words $wordlength | awk '{print $1}')"
+#words="$(./billiard_words $wordlength | awk '{print $1}')"
+#words="cbabcabacabcacbcab cabacabcacbcabcbab cabcacbcabcbabcaba"
+words="abcb acbc baca"
 
 for s in $(seq $sstart $send); do
 	for q in $(seq $qstart $qend); do
+		i=0
 		echo -n "$s/$sdenom $q/$qdenom "
-		MAXIMUM=only ./special_element $s/$sdenom $q/$qdenom $words
+#		MAXIMUM=only ./special_element $s/$sdenom $q/$qdenom $words
+#		MAXIMUM=no ./special_element $s/$sdenom $q/$qdenom abcb
+		MAXIMUM=no ./special_element $s/$sdenom $q/$qdenom $words | while read line; do
+			echo -n "$line "
+			((i=i+1))
+		done
+		echo
 	done
 done

billiard_words.hs

@@ -9,20 +9,21 @@ main = do
 
 listWordsUpToLength :: Int -> IO ()
 listWordsUpToLength n = do
-  putStrLn $ unlines [printf "%s %d/%d %f"
+  putStrLn $ unlines [printf "%s %d/%d %f %s"
                              w
-                             (x `div` gcd x y)
+                             (p `div` gcd p q)
-                             (y `div` gcd x y)
+                             (q `div` gcd p q)
-                             (fromIntegral x / fromIntegral y :: Double) |
+                             (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double)
-                      ((p,q),w) <- wordlist (n `div` 2) (n `div` 2),
+                             (slopeWord "bca" (orthogonalSlope (p,q))) |
+                      ((p,q),w) <- wordlist (n `div` 2, n `div` 2),
                       length w <= n,
                       let x = 2*q + p,
                       let y = 2*p + q]
 
-wordlist :: Int -> Int -> [((Int,Int),String)]
+wordlist :: (Int,Int) -> [((Int,Int),String)]
-wordlist pmax qmax = nub $
+wordlist (pmax,qmax) = nub $
                      sortBy (comparing sl)
-                     [((p `div` gcd p q, q `div` gcd p q), slopeWord "bca" p q) |
+                     [((p `div` gcd p q, q `div` gcd p q), slopeWord "bca" (p,q)) |
                       p <- [0..pmax],
                       q <- [0..qmax],
                       q /= 0]   -- use p /= 0 || q /= 0 for more symmetric output
@@ -30,9 +31,11 @@ wordlist pmax qmax = nub $
       sl ((p,q),_) = fromIntegral p / fromIntegral q
 
 -- letters: reflection along e_1, reflection along e_2, other one; p,q >= 0
-slopeWord :: [Char] -> Int -> Int -> String
+-- the "slope" (p,q) means the Euclidean vector q*e_1 + p*e_2, where e_1,e_2 are at a 60 degree angle
-slopeWord [x,y,z] p q
+-- in Euclidean coordinates this is (q + p/2, sqrt(3)/2 * p)
-    | p > q = slopeWord [y,x,z] q p
+slopeWord :: [Char] -> (Int,Int) -> String
+slopeWord [x,y,z] (p,q)
+    | p > q = slopeWord [y,x,z] (q,p)
     | otherwise = concat $ map word $ zipWith step list (tail list)
     where
       p_ = p `div` gcd p q
@@ -46,3 +49,10 @@ slopeWord [x,y,z] p q
       word (0,1) = [z,y,x,y]
       word (1,1) = [y,x,z,x]
       word (2,1) = [x,z,y,z]
+
+-- assuming p, q >= 0
+orthogonalSlope :: (Int, Int) -> (Int, Int)
+orthogonalSlope (p,q)
+    | p > q = (p-q, p+2*q)
+    | p < q = (q+2*p, q-p)
+    | otherwise = (1,0)

enumerate_triangle_group.c

@@ -107,14 +107,17 @@ void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2,
 	mat_init(r2, 3);
 	mat_init(r3, 3);
 
+	// sinv = s^{-1}
 	mpq_set_ui(sinv, 1, 1);
 	mpq_div(sinv, sinv, s);
 
+	// rho_i = s^2 + 2s cos(2 pi / p_i) + 1
 	// coefficient 2 is the value for p=infinity, not sure if that would even work
 	quartic(rho1, s, 0, 0, 1, p1 == 2 ? -2 : p1 == 3 ? -1 : p1 == 4 ? 0 : p1 == 6 ? 1 : 2, 1);
 	quartic(rho2, s, 0, 0, 1, p2 == 2 ? -2 : p2 == 3 ? -1 : p2 == 4 ? 0 : p2 == 6 ? 1 : 2, 1);
 	quartic(rho3, s, 0, 0, 1, p3 == 2 ? -2 : p3 == 3 ? -1 : p3 == 4 ? 0 : p3 == 6 ? 1 : 2, 1);
 
+	// c1 = rho2 q, a2 = rho3 q, b3 = rho1 q, b1 = c2 = a3 = 1/q
 	mpq_mul(c1, rho2, q);
 	mpq_mul(a2, rho3, q);
 	mpq_mul(b3, rho1, q);
@@ -125,7 +128,6 @@ void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2,
 	mpq_div(c2, c2, q);
 	mpq_div(a3, a3, q);
 
-	// actually, we want minus everything
 	mat_zero(r1);
 	mat_zero(r2);
 	mat_zero(r3);
@@ -150,18 +152,24 @@ void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2,
 	mat_zero(gen[1]);
 	mat_zero(gen[2]);
 
+	// gen[0] = diag(1,s^{-1},s)
 	mpq_set_ui(*mat_ref(gen[0], 0, 0), 1, 1);
 	mat_set(gen[0], 1, 1, sinv);
 	mat_set(gen[0], 2, 2, s);
 
+	// gen[1] = diag(s,1,s^{-1})
 	mat_set(gen[1], 0, 0, s);
 	mpq_set_ui(*mat_ref(gen[1], 1, 1), 1, 1);
 	mat_set(gen[1], 2, 2, sinv);
 
+	// gen[3] = diag(s^{-1},s,1)
 	mat_set(gen[2], 0, 0, sinv);
 	mat_set(gen[2], 1, 1, s);
 	mpq_set_ui(*mat_ref(gen[2], 2, 2), 1, 1);
 
+	// gen[0] = r2 * gen[0] * r3
+	// gen[1] = r3 * gen[1] * r1
+	// gen[2] = r1 * gen[2] * r2
 	mat_multiply(ws, gen[0], r2, gen[0]);
 	mat_multiply(ws, gen[0], gen[0], r3);
 	mat_multiply(ws, gen[1], r3, gen[1]);
@@ -169,6 +177,9 @@ void initialize_triangle_generators(mat_workspace *ws, mat *gen, int p1, int p2,
 	mat_multiply(ws, gen[2], r1, gen[2]);
 	mat_multiply(ws, gen[2], gen[2], r2);
 
+	// gen[3] = gen[0]^{-1}
+	// gen[4] = gen[1]^{-1}
+	// gen[5] = gen[2]^{-1}
 	mat_pseudoinverse(ws, gen[3], gen[0]);
 	mat_pseudoinverse(ws, gen[4], gen[1]);
 	mat_pseudoinverse(ws, gen[5], gen[2]);

special_element.c

@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
 	}
 
 	for(int w = 0; w < argc - 3; w++) {
-		initialize_triangle_generators(ws, gen, 4, 4, 4, s, q);
+		initialize_triangle_generators(ws, gen, 6, 4, 3, s, q);
 
 		mat_identity(element);
 		mat_identity(inverse);