Verifying https://github.com/TokusiN/AtanMagic with Haskell

atanmagic.hs

#!/usr/bin/env cabal
{- cabal:
build-depends: base, vector, process
-}
module Main where

{-

検算

https://twitter.com/toku51n/status/1533401499983749121

-}

import System.IO (readFile')
import System.Process

import Data.List (scanl')
import Data.Foldable (foldl', for_)
import Data.Ratio
import qualified Data.Vector as V

import Debug.Trace

size :: Int
size = 16

rowIndices, colIndices :: Int -> [Int]
rowIndices r = [ r * size + c | c <- [0 .. size - 1] ]
colIndices c = [ r * size + c | r <- [0 .. size - 1] ]

diag, antiDiag :: [Int]
diag     = [i * size + i | i <- [0 .. size - 1] ]
antiDiag = [i * size + (size - 1 - i) | i <- [0 .. size - 1] ]

-- 主張されていること(atanして総和すれば2πになる)ことを確認するには、
--   x ⊕ y = z
-- のとき
--   atan x + atan y = atan z
-- となるような⊕をうまく定義できれば便利そうである。
-- 
-- いま、 x = a/b, y = c/d（a,b,c,dは整数）と書かれていたとき、
--
-- z = tan (atan z) = tan (atan x + atan y)
--   = (tan (atan x) + tan (atan y)) / (1 - tan (atan x) * tan (atan y))
--   = ((a/b) + (c/d)) / (1 - (a/b) * (c/d))
--   = (a * d + b * c) / (b * d - a * c)
--
-- なので、浮動小数点演算のような精度の落ちるものを用いなくても、
-- 整数演算だけで検算ができる。
-- 
-- tan(2π) = 0 = 0/1
-- 
-- より、x0 ⊕ x1 ⊕  ・・・ ⊕  xn = 0
-- が各列、行、対角線について成り立っている。
-- しかし、この等式を確認しただけでは、atan(x_i)の総和がπの整数倍であったことしか
-- わからない。ちょうど2πになったことはどのように確認できるだろうか？
--
-- ここでtan(θ)のグラフを思い浮かべると、tanは周期πをもち、
-- 各(π/2+kπ)にある極を除いて単調増加だった。いま、魔法陣の各マスの値xは正の数なので、
--   0 < atan(x_i) < π/2
-- である。ならば、数列
-- 
-- tan(atan(x0)), tan(atan(x0) + atan(x1)), tan(atan(x0) + atan(x1) + atan(x2)), ...
-- 
-- すなわち
--   x0, x0 ⊕ x1, x0 ⊕ x1 ⊕ x2, ...
-- は、atan(x_i)の総和がkπになるならば、ちょうどk回減少するはずである。

-- Haskellには任意精度の有理数型Rationalがあるのだが、
-- 上記の計算中に tan(π/2) に相当する値が出現したときにゼロ除算となってしまう。
-- しかし、上記の計算はtan(π/2)が"分数"(1/0)に対応すると考えても全く問題がなく行える。
-- そのため、"ゼロ除算"の結果を(1/0)で表す、有理数のような型を新たに作って対応する。
data PQ = PQ !Integer !Integer
  deriving (Show, Eq)

instance Ord PQ where
  PQ a b <= PQ c d = a * d <= c * b

rationalToPQ :: Rational -> PQ
rationalToPQ x = PQ (numerator x) (denominator x)

pqToDouble :: PQ -> Double
pqToDouble (PQ a b) = fromInteger a / fromInteger b

pqZero :: PQ
pqZero = PQ 0 1

atanPlus :: PQ -> PQ -> PQ
atanPlus (PQ a b) (PQ c d) = PQ x y
  where
    numer = a * d + b * c
    denom = b * d - a * c
    e = gcd numer denom * (if denom < 0 then -1 else 1)
    x = numer `quot` e
    y = denom `quot` e

atanSum :: [PQ] -> PQ
atanSum = foldl' atanPlus pqZero

atanCumsum :: [PQ] -> [PQ]
atanCumsum = scanl' atanPlus pqZero

countDecreases :: Ord a => [a] -> Int
countDecreases as = length $ filter (\(x,y) -> x > y) $ zip as (tail as)

dataFileName :: String
dataFileName = "atan-magic-data.txt"

dataUrl :: String
dataUrl = "https://raw.githubusercontent.com/TokusiN/AtanMagic/a6c4b4acd40e5869b7c983e5455ecef3668d1e89/data.txt"

parseDataFileText :: String -> [Rational]
parseDataFileText = concatMap (\s -> read (brackets s) :: [Rational]) . lines . map (replace '/' '%')
  where
    replace :: Eq a => a -> a -> a -> a
    replace x y z | x == z    = y
                  | otherwise = z
    
    brackets :: String -> String
    brackets str = "[" ++ str ++ "]"

main :: IO ()
main = do
  -- Fetch the magic square data from GitHub
  callProcess "wget" ["-O", dataFileName, dataUrl]
  -- Read the downloaded file
  fileContent <- readFile' dataFileName
  let -- Parse the read file into list of Rationals
      magicSquareDataList = parseDataFileText fileContent
      -- Convert Rationals to PQ; also store in Vector instead of list
      square = V.fromList (map rationalToPQ magicSquareDataList)
      
      -- Define a utility function to check for the specific part of 
      -- (e.g. row, column, ...) of the loaded data
      validation ix =
        let ps = (square V.!) <$> ix
            total = atanSum ps
            sumPs = atanCumsum ps
            wc = countDecreases sumPs
        in (total, wc, total == pqZero && wc == 2)
  
  for_ [0 .. size - 1] $ \r ->
     putStrLn $ "行" ++ show r ++ ":" ++ show (validation (rowIndices r))
  for_ [0 .. size - 1] $ \c ->
     putStrLn $ "列" ++ show c ++ ":" ++ show (validation (colIndices c))
  putStrLn $ "斜＼:" ++ show (validation diag)
  putStrLn $ "斜／:" ++ show (validation antiDiag)