triangle_reflection_complex/billiard_words.hs

63 lines
2.3 KiB
Haskell
Raw Normal View History

import Data.List
import Data.Ord
import Text.Printf
import System.Environment
main = do
argv <- getArgs
listWordsUpToLength $ read $ argv !! 0
listWordsUpToLength :: Int -> IO ()
listWordsUpToLength n = do
putStr $ unlines [printf "%s %d/%d %f"
w
2022-02-03 17:33:20 +00:00
(p `div` gcd p q)
(q `div` gcd p q)
(atan (sqrt 3 / (2*q_/p_ + 1))) |
2022-02-03 17:33:20 +00:00
((p,q),w) <- wordlist (n `div` 2, n `div` 2),
let p_ = fromIntegral p :: Double,
let q_ = fromIntegral q :: Double,
length w <= n,
let x = 2*q + p,
let y = 2*p + q]
-- (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double) |
-- (slopeWord "bca" (orthogonalSlope (p,q))) |
2022-02-03 17:33:20 +00:00
wordlist :: (Int,Int) -> [((Int,Int),String)]
wordlist (pmax,qmax) = nub $
sortBy (comparing sl)
2022-02-03 17:33:20 +00:00
[((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
where
sl ((p,q),_) = fromIntegral p / fromIntegral q
-- letters: reflection along e_1, reflection along e_2, other one; p,q >= 0
2022-02-03 17:33:20 +00:00
-- 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
-- in Euclidean coordinates this is (q + p/2, sqrt(3)/2 * 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
q_ = q `div` gcd p q
2021-10-05 01:12:47 +00:00
xmax = if (p_-q_) `mod` 3 == 0 then q_ else 3*q_ :: Int
list = [(x,(x*p) `div` q) | x <- [0..xmax]]
step (x1,y1) (x2,y2) = ((x1-y1) `mod` 3, y2-y1)
word (0,0) = [z,x]
word (1,0) = [y,z]
word (2,0) = [x,y]
word (0,1) = [z,y,x,y]
word (1,1) = [y,x,z,x]
word (2,1) = [x,z,y,z]
2022-02-03 17:33:20 +00:00
-- 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)