８４問目

http://projecteuler.net/index.php?section=problems&id=84
モノポリーには 40 のマスがあるが、それぞれのマスで停まる確率を求め、TOP3を答える。
Project Euler では問題文が長いのは嫌われる傾向にあるみたいですが（５４問目とか）、この問題は色々な要素があって、とても楽しかったです。
ある時点 n に、あるマスにいる確率を P_n としたとき、別のマスへ遷移する確率を求めて P_n と P_n+1 の漸化式を出します。定常状態では P_n = P_n+1 とおけて、これは40元連立1次方程式になるので、行列で表現して、ガウス・ジョルダン法で解いてます。
肝心の遷移確率ですが、サイコロの目で遷移する確率（ゾロ目を含む）を関数 p' で求めて、それにカードによる遷移を加味して関数 p を作っています。
プログラムが長くなった分、やはりバグが入りやすくなっていて、見つけるのに苦労しました。まだ数値が若干違うので（6面サイコロでJAILに停まる確率が 6.24% ではなく 6.30% になる）、もしかしたらバグが残っているかもしれません。

import Data.List
import Data.Ratio

main = print $ concat $ take 3 $ map snd $ euler084 eq

euler084 :: [[Rational]] -> [(Double, String)]
euler084 es = reverse $ sort $ zip (map (fromRational . last) $ foldl solve es [0..39]) digits
    where digits = [ [c1,c2] | c1 <- "0123456789", c2 <- "0123456789" ]

solve :: [[Rational]] -> Int -> [[Rational]]
solve es n = (\ (f, s) -> (map sub f) ++ (e' : (map sub $ tail s))) $ splitAt n es
    where e' = map (/ (es !! n !! n)) (es !! n)
          sub e = zipWith (-) e (map (* (e !! n)) e')

eq :: [[Rational]]
eq = [ [ (if m == n then 1 else 0) - (p m n) + 1 | m <- [0..39] ] ++ [1] | n <- [0..39] ]

p :: Integer -> Integer -> Rational
p m n
    | n ==  0 = (p' m n) + (cc m)/16 + (ch m)/16 -- GO
    | n == 10 = (p' m n) + (cc m)/16 + (ch m)/16 + (g2j m) -- JAIL
    | n == 11 = (p' m n) + (ch m)/16 -- C1
    | n == 24 = (p' m n) + (ch m)/16 -- E3
    | n == 39 = (p' m n) + (ch m)/16 -- H2
    | n ==  5 = (p' m n) + (ch m)/16 + (ch3 m)*2/16 -- R1
    | n == 15 = (p' m n) + (ch1 m)*2/16 -- R2
    | n == 25 = (p' m n) + (ch2 m)*2/16 -- R3
    | n == 12 = (p' m n) + (ch1 m)/16 + (ch3 m)/16 -- U1
    | n == 28 = (p' m n) + (ch2 m)/16 -- U2
    | n ==  4 = (p' m n) + (ch1 m)/16 -- 3 squares before CH1
    | n == 19 = (p' m n) + (ch2 m)/16 -- 3 squares before CH2
    | n ==  2 = (cc1 m)*14/16 -- CC1
    | n == 17 = (cc2 m)*14/16 -- CC2
    | n == 33 = (cc3 m)*14/16 -- CC3
    | n ==  7 = (ch1 m)*6/16 -- CH1
    | n == 22 = (ch2 m)*6/16 -- CH2
    | n == 36 = (ch3 m)*6/16 -- CH3
    | n == 30 = 0 -- G2J
    | otherwise = p' m n
    where cc  m = (cc1 m) + (cc2 m) + (cc3 m)
          cc1 m = p' m  2
          cc2 m = p' m 17
          cc3 m = p' m 33 + (ch3 m)/16 -- 3 squares before CC 3
          ch  m = (ch1 m) + (ch2 m) + (ch3 m)
          ch1 m = p' m  7
          ch2 m = p' m 22
          ch3 m = p' m 36
          g2j m = p' m 30

p' :: Integer -> Integer -> Rational
p' m n
    | d ==   0 = dub1
    | d <= 2*i = dub1 - dub0 + (min (d-1) (2*i-d+1))%(i^2)
    | otherwise = dub1
    where i = 4
          d = mod (n - m + 40) 40
          dub0 = if  even d then 1%(i^4) else 0
          dub1 = if n == 10 then 1%(i^3) else 0