Main.hs

   1 #!/usr/bin/env -S nix -L shell .#ghc --command runghc
   2 -- #!/usr/bin/env -S nix -L shell .#ghc .#ghcid --command ghcid --test :main
   3 -- SPDX-FileCopyrightText: 2022 Julien Mouyinho <julm+most-perfect-magic-squares@sourcephile.fr>
   4 -- SPDX-License-Identifier: CC0-1.0
   5 {-# OPTIONS_GHC -Wno-missing-signatures #-}
   6 {-# LANGUAGE BlockArguments #-}
   7 {-# LANGUAGE NoImplicitPrelude #-}
   8 {-# LANGUAGE OverloadedStrings #-}
   9 module Main where
  10
  11 -- import Blaze.ByteString.Builder (Builder)
  12 import Data.ByteString qualified as BS
  13 import Relude
  14 import System.IO qualified as IO
  15 import Text.Blaze
  16 import qualified Text.Blaze.Html5 as H
  17 import qualified Text.Blaze.Html5.Attributes as HA
  18 import qualified Text.Blaze.Renderer.Utf8 as Blaze
  19 import Data.List ((\\))
  20 import Data.List qualified as List
  21
  22 data I = D | M | C | Y
  23   deriving (Eq, Ord, Show)
  24 data Cell = I I | H I
  25   deriving (Eq, Ord, Show)
  26
  27 mpms4 :: a -> a -> a -> a -> a -> a -> a -> a -> [[a]]
  28 mpms4 iD iM iC iY hD hM hC hY =
  29   [ [iD, iM, iC, iY]
  30   , [iC, iY, iD, iM]
  31   , [hC, hY, hD, hM]
  32   , [hD, hM, hC, hY]
  33   ]
  34
  35 isValid :: [Cell] -> Bool
  36 isValid [I D, I M, I C, I Y] = True
  37 isValid [H D, H M, H C, H Y] = True
  38 isValid [I a, I b, H aa, H bb] | a == aa && b == bb = True
  39 isValid _ = False
  40
  41 squary :: [[((Integer, Integer), Cell)]]
  42 squary =
  43   zipWith (\r -> zipWith (\c -> ((r,c),)) [0..]) [0..] $
  44   mpms4 (I D) (I M) (I C) (I Y) (H D) (H M) (H C) (H Y)
  45 allPaths = combins 4 (concat squary)
  46 validPaths = filter (isValid . sort . (snd <$>)) $ allPaths
  47 validSortedPaths = sortBy (compare `on` (snd <$>)) validPaths
  48 validSortedCoords = sort $ (fst <$>) <$> validPaths
  49
  50 mostPerfectMagicSquare4 :: Fractional a => a -> a -> a -> a -> [[a]]
  51 mostPerfectMagicSquare4 d m c y = mpms4 d m c y (h d) (h m) (h c) (h y)
  52   where h x = ((d + m + c + y) / 2) - x
  53
  54 square :: [[Int]]
  55 square = (round @Rational <$>) <$> mostPerfectMagicSquare4 24 03 19 74
  56
  57 pathsByCombin :: [[[(Integer, Integer)]]]
  58 pathsByCombin =
  59   [ [ [ (0,c) | c <- [0..3] ] ]
  60   , [ [(0,c0), (0,c1), (0,c2), (1,c3)]
  61       | [c0, c1, c2] <- combins 3 [0..3]
  62       , let c3 = (List.head ([0..3] \\ [c0, c1, c2]) - 2) `mod` 4
  63     ]
  64   , [ [(0,c0), (0,c1), (1,c2), (1,c3)]
  65       | [c0, c1] <- combins 2 [0..3]
  66       , let [c2, c3] = ([0..3] \\ [c0, c1]) <&> (\c -> (c - 2) `mod` 4)
  67     ]
  68   , [ [(0,c0), (0,c1), (2,c2), (2,c3)]
  69       | [c0, c1] <- combins 2 [0..3]
  70       , let [c2, c3] = [c0, c1] <&> (\c -> (c - 2) `mod` 4)
  71     ]
  72   , [ [(0,c0), (0,c1), (3,c2), (3,c3)]
  73       | [c0, c1] <- combins 2 [0..3]
  74       , let [c2, c3] = [c0, c1]
  75     ]
  76   , [ [(0,c0), (1,c1), (1,c2), (1,c3)]
  77       | [c0] <- combins 1 [0..3]
  78       , let [c1, c2, c3] = ([0..3] \\ [c0]) <&> (\c -> (c - 2) `mod` 4)
  79     ]
  80   , [ [(0,c0), (1,c1), (2,c2), (2,c3)]
  81       | [c0] <- combins 1 [0..3]
  82       , [c1] <- combins 1 ([0..3] \\ [(c0 + 2) `mod` 4])
  83       , let [c2, c3] = [(c0 + 2) `mod` 4, c1]
  84     ]
  85   , [ [(0,c0), (1,c1), (2,c2), (3,c3)]
  86       | [c0] <- combins 1 [0..3]
  87       , [c1] <- combins 1 ([0..3] \\ [(c0 + 2) `mod` 4])
  88       , [c2, c3] <- [[(c0+2)`mod`4, (c1+2)`mod`4], [c1, c0]]
  89     ]
  90   , [ [(0,c0), (1,c1), (3,c2), (3,c3)]
  91       | [c0] <- combins 1 [0..3]
  92       , [c1] <- combins 1 ([0..3] \\ [(c0 + 2) `mod` 4])
  93       , let [c2, c3] = [c0, (c1+2)`mod`4]
  94     ]
  95   , [ [(1,c0), (1,c1), (1,c2), (1,c3)]
  96       | [c0, c1, c2, c3] <- combins 4 [0..3]
  97     ]
  98   , [ [(1,c0), (1,c1), (2,c2), (2,c3)]
  99       | [c0, c1] <- combins 2 [0..3]
 100       , let [c2, c3] = [c0, c1]
 101     ]
 102   , [ [(1,c0), (1,c1), (2,c2), (3,c3)]
 103       | [c0, c1] <- combins 2 [0..3]
 104       , let [c2, c3] = [c0, (c1 + 2) `mod` 4]
 105     ]
 106   , [ [(1,c0), (1,c1), (3,c2), (3,c3)]
 107       | [c0, c1] <- combins 2 [0..3]
 108       , let [c2, c3] = [(c0 + 2) `mod` 4, (c1 + 2) `mod` 4]
 109     ]
 110   , [ [(2,c0), (2,c1), (2,c2), (2,c3)]
 111       | [c0, c1, c2, c3] <- combins 4 [0..3]
 112     ]
 113   , [ [(2,c0), (2,c1), (2,c2), (3,c3)]
 114       | [c0, c1, c2] <- combins 3 [0..3]
 115       , let [c3] = [0..3] \\ [c0, c1, c2]
 116     ]
 117   , [ [(2,c0), (2,c1), (3,c2), (3,c3)]
 118       | [c0, c1] <- combins 2 [0..3]
 119       , let [c2, c3] = [c0, c1]
 120     ]
 121   , [ [(2,c0), (3,c1), (3,c2), (3,c3)]
 122       | [c0] <- combins 1 [0..3]
 123       , let [c1, c2, c3] = ([0..3] \\ [(c0 + 2) `mod` 4])
 124     ]
 125   , [ [(3,c0), (3,c1), (3,c2), (3,c3)]
 126       | [c0, c1, c2, c3] <- combins 4 [0..3]
 127     ]
 128   ]
 129 ol :: [b] -> [(Integer, b)]
 130 ol = zip [0..]
 131
 132 main :: IO ()
 133 main = do
 134   IO.withFile "mpms.html" IO.WriteMode $ \h ->
 135     Blaze.renderMarkupToByteStringIO (BS.hPutStr h) do
 136       H.docTypeHtml do
 137         H.head do
 138           H.title "Most-Perfect Magic Square"
 139           H.link ! HA.rel "stylesheet"
 140                  ! HA.type_ "text/css"
 141                  ! HA.href "mpms.css"
 142         H.body do
 143           H.ul do
 144             H.li $ "#squares = " <> show (length validSortedCoords)
 145             H.li $ "sum = " <> show (sum (List.head square))
 146           H.div ! HA.class_ "squares" $ do
 147             forM_ validSortedCoords $ \path ->
 148               H.div ! HA.class_ "square" $ do
 149                 forM_ (ol square) $ \(rowCoord, row) ->
 150                     forM_ (ol row) $ \(colCoord, num) ->
 151                       H.span ! HA.class_ ("square-num"<>""<>(if (rowCoord,colCoord) `elem` path then " num-path" else "")) $
 152                         fromString $ show num
 153
 154 -- | @'nCk' n k@ retourne le nombre de combinaisons
 155 -- de longueur 'k' d’un ensemble de longueur 'n'.
 156 --
 157 -- Computed using the formula:
 158 -- @'nCk' n (k+1) == 'nCk' n (k-1) * (n-k+1) / k@
 159 nCk :: Integral i => i -> i -> i
 160 n`nCk`k | n<0||k<0||n<k = error "nCk"
 161         | otherwise     = go 1 1
 162         where
 163         go i acc = if k' < i then acc else go (i+1) (acc * (n-i+1) `div` i)
 164         -- Use a symmetry to compute over smaller numbers,
 165         -- which is more efficient and safer
 166         k' = if n`div`2 < k then n-k else k
 167
 168 -- | @combinOfRank n k r@ retourne les indices de permutation
 169 -- de la combinaison de 'k' entiers parmi @[1..n]@
 170 -- au rang lexicographique 'r' dans @[0..'nCk' n k - 1]@.
 171 --
 172 -- Construit chaque choix de la combinaison en prenant le prochain plus grand
 173 -- dont le successeur engendre un nombre de combinaisons
 174 -- qui dépasse le rang restant à atteindre.
 175 --
 176 -- DOC: <http://www.site.uottawa.ca/~lucia/courses/5165-09/GenCombObj.pdf>, p.26
 177 combinOfRank :: Integral i => i -> i -> i -> [i]
 178 combinOfRank n k rk | rk<0||n`nCk`k<rk = error "combinOfRank"
 179                     | otherwise = for1K 1 1 rk
 180  where
 181  for1K i j r | i <  k    = uptoRank i j r
 182              | i == k    = [j+r] -- because when i == k, nbCombs is always 1
 183              | otherwise = []
 184  uptoRank i j r | nbCombs <- (n-j)`nCk`(k-i)
 185                 , nbCombs <= r = uptoRank i (j+1) (r-nbCombs)
 186                 | otherwise    = j : for1K (i+1) (j+1) r
 187
 188 -- | @combins k xs@ retourne les combinaisons
 189 -- de longueur 'k' des éléments de la liste 'xs'.
 190 combins :: Int -> [a] -> [[a]]
 191 combins k xs = combinsK xs List.!! (length xs - k)
 192
 193 -- | @combinsK xs@ retourne toutes les combinaisons
 194 -- de longueur allant de @length xs@ à 0 de la liste 'xs',
 195 --
 196 -- Algorithme dynamique permettant un calcul de 'combins'
 197 -- relativement rapide du fait du partage de 'combinsKmoins1'.
 198 combinsK :: [a] -> [[[a]]]
 199 combinsK [] = [[[]]]
 200 combinsK (x : xs) =
 201        zipWith (++) ([] : combinsKmoins1)
 202                     (map (map (x :)) combinsKmoins1 ++ [[]])
 203        where combinsKmoins1 = combinsK xs