+#!/usr/bin/env -S nix -L shell .#ghc --command runghc
+-- #!/usr/bin/env -S nix -L shell .#ghc .#ghcid --command ghcid --test :main
+-- SPDX-FileCopyrightText: 2022 Julien Mouyinho <julm+most-perfect-magic-squares@sourcephile.fr>
+-- SPDX-License-Identifier: CC0-1.0
+{-# OPTIONS_GHC -Wno-missing-signatures #-}
+{-# LANGUAGE BlockArguments #-}
+{-# LANGUAGE NoImplicitPrelude #-}
+{-# LANGUAGE OverloadedStrings #-}
+module Main where
+
+-- import Blaze.ByteString.Builder (Builder)
+import Data.ByteString qualified as BS
+import Relude
+import System.IO qualified as IO
+import Text.Blaze
+import qualified Text.Blaze.Html5 as H
+import qualified Text.Blaze.Html5.Attributes as HA
+import qualified Text.Blaze.Renderer.Utf8 as Blaze
+import Data.List ((\\))
+import Data.List qualified as List
+
+data I = D | M | C | Y
+ deriving (Eq, Ord, Show)
+data Cell = I I | H I
+ deriving (Eq, Ord, Show)
+
+mpms4 :: a -> a -> a -> a -> a -> a -> a -> a -> [[a]]
+mpms4 iD iM iC iY hD hM hC hY =
+ [ [iD, iM, iC, iY]
+ , [iC, iY, iD, iM]
+ , [hC, hY, hD, hM]
+ , [hD, hM, hC, hY]
+ ]
+
+isValid :: [Cell] -> Bool
+isValid [I D, I M, I C, I Y] = True
+isValid [H D, H M, H C, H Y] = True
+isValid [I a, I b, H aa, H bb] | a == aa && b == bb = True
+isValid _ = False
+
+squary :: [[((Integer, Integer), Cell)]]
+squary =
+ zipWith (\r -> zipWith (\c -> ((r,c),)) [0..]) [0..] $
+ mpms4 (I D) (I M) (I C) (I Y) (H D) (H M) (H C) (H Y)
+allPaths = combins 4 (concat squary)
+validPaths = filter (isValid . sort . (snd <$>)) $ allPaths
+validSortedPaths = sortBy (compare `on` (snd <$>)) validPaths
+validSortedCoords = sort $ (fst <$>) <$> validPaths
+
+mostPerfectMagicSquare4 :: Fractional a => a -> a -> a -> a -> [[a]]
+mostPerfectMagicSquare4 d m c y = mpms4 d m c y (h d) (h m) (h c) (h y)
+ where h x = ((d + m + c + y) / 2) - x
+
+square :: [[Int]]
+square = (round @Rational <$>) <$> mostPerfectMagicSquare4 24 03 19 74
+
+pathsByCombin :: [[[(Integer, Integer)]]]
+pathsByCombin =
+ [ [ [ (0,c) | c <- [0..3] ] ]
+ , [ [(0,c0), (0,c1), (0,c2), (1,c3)]
+ | [c0, c1, c2] <- combins 3 [0..3]
+ , let c3 = (List.head ([0..3] \\ [c0, c1, c2]) - 2) `mod` 4
+ ]
+ , [ [(0,c0), (0,c1), (1,c2), (1,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = ([0..3] \\ [c0, c1]) <&> (\c -> (c - 2) `mod` 4)
+ ]
+ , [ [(0,c0), (0,c1), (2,c2), (2,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = [c0, c1] <&> (\c -> (c - 2) `mod` 4)
+ ]
+ , [ [(0,c0), (0,c1), (3,c2), (3,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = [c0, c1]
+ ]
+ , [ [(0,c0), (1,c1), (1,c2), (1,c3)]
+ | [c0] <- combins 1 [0..3]
+ , let [c1, c2, c3] = ([0..3] \\ [c0]) <&> (\c -> (c - 2) `mod` 4)
+ ]
+ , [ [(0,c0), (1,c1), (2,c2), (2,c3)]
+ | [c0] <- combins 1 [0..3]
+ , [c1] <- combins 1 ([0..3] \\ [(c0 + 2) `mod` 4])
+ , let [c2, c3] = [(c0 + 2) `mod` 4, c1]
+ ]
+ , [ [(0,c0), (1,c1), (2,c2), (3,c3)]
+ | [c0] <- combins 1 [0..3]
+ , [c1] <- combins 1 ([0..3] \\ [(c0 + 2) `mod` 4])
+ , [c2, c3] <- [[(c0+2)`mod`4, (c1+2)`mod`4], [c1, c0]]
+ ]
+ , [ [(0,c0), (1,c1), (3,c2), (3,c3)]
+ | [c0] <- combins 1 [0..3]
+ , [c1] <- combins 1 ([0..3] \\ [(c0 + 2) `mod` 4])
+ , let [c2, c3] = [c0, (c1+2)`mod`4]
+ ]
+ , [ [(1,c0), (1,c1), (1,c2), (1,c3)]
+ | [c0, c1, c2, c3] <- combins 4 [0..3]
+ ]
+ , [ [(1,c0), (1,c1), (2,c2), (2,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = [c0, c1]
+ ]
+ , [ [(1,c0), (1,c1), (2,c2), (3,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = [c0, (c1 + 2) `mod` 4]
+ ]
+ , [ [(1,c0), (1,c1), (3,c2), (3,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = [(c0 + 2) `mod` 4, (c1 + 2) `mod` 4]
+ ]
+ , [ [(2,c0), (2,c1), (2,c2), (2,c3)]
+ | [c0, c1, c2, c3] <- combins 4 [0..3]
+ ]
+ , [ [(2,c0), (2,c1), (2,c2), (3,c3)]
+ | [c0, c1, c2] <- combins 3 [0..3]
+ , let [c3] = [0..3] \\ [c0, c1, c2]
+ ]
+ , [ [(2,c0), (2,c1), (3,c2), (3,c3)]
+ | [c0, c1] <- combins 2 [0..3]
+ , let [c2, c3] = [c0, c1]
+ ]
+ , [ [(2,c0), (3,c1), (3,c2), (3,c3)]
+ | [c0] <- combins 1 [0..3]
+ , let [c1, c2, c3] = ([0..3] \\ [(c0 + 2) `mod` 4])
+ ]
+ , [ [(3,c0), (3,c1), (3,c2), (3,c3)]
+ | [c0, c1, c2, c3] <- combins 4 [0..3]
+ ]
+ ]
+ol :: [b] -> [(Integer, b)]
+ol = zip [0..]
+
+main :: IO ()
+main = do
+ IO.withFile "mpms.html" IO.WriteMode $ \h ->
+ Blaze.renderMarkupToByteStringIO (BS.hPutStr h) do
+ H.docTypeHtml do
+ H.head do
+ H.title "Most-Perfect Magic Square"
+ H.link ! HA.rel "stylesheet"
+ ! HA.type_ "text/css"
+ ! HA.href "mpms.css"
+ H.body do
+ H.ul do
+ H.li $ "#squares = " <> show (length validSortedCoords)
+ H.li $ "sum = " <> show (sum (List.head square))
+ H.div ! HA.class_ "squares" $ do
+ forM_ validSortedCoords $ \path ->
+ H.div ! HA.class_ "square" $ do
+ forM_ (ol square) $ \(rowCoord, row) ->
+ forM_ (ol row) $ \(colCoord, num) ->
+ H.span ! HA.class_ ("square-num"<>""<>(if (rowCoord,colCoord) `elem` path then " num-path" else "")) $
+ fromString $ show num
+
+-- | @'nCk' n k@ retourne le nombre de combinaisons
+-- de longueur 'k' d’un ensemble de longueur 'n'.
+--
+-- Computed using the formula:
+-- @'nCk' n (k+1) == 'nCk' n (k-1) * (n-k+1) / k@
+nCk :: Integral i => i -> i -> i
+n`nCk`k | n<0||k<0||n<k = error "nCk"
+ | otherwise = go 1 1
+ where
+ go i acc = if k' < i then acc else go (i+1) (acc * (n-i+1) `div` i)
+ -- Use a symmetry to compute over smaller numbers,
+ -- which is more efficient and safer
+ k' = if n`div`2 < k then n-k else k
+
+-- | @combinOfRank n k r@ retourne les indices de permutation
+-- de la combinaison de 'k' entiers parmi @[1..n]@
+-- au rang lexicographique 'r' dans @[0..'nCk' n k - 1]@.
+--
+-- Construit chaque choix de la combinaison en prenant le prochain plus grand
+-- dont le successeur engendre un nombre de combinaisons
+-- qui dépasse le rang restant à atteindre.
+--
+-- DOC: <http://www.site.uottawa.ca/~lucia/courses/5165-09/GenCombObj.pdf>, p.26
+combinOfRank :: Integral i => i -> i -> i -> [i]
+combinOfRank n k rk | rk<0||n`nCk`k<rk = error "combinOfRank"
+ | otherwise = for1K 1 1 rk
+ where
+ for1K i j r | i < k = uptoRank i j r
+ | i == k = [j+r] -- because when i == k, nbCombs is always 1
+ | otherwise = []
+ uptoRank i j r | nbCombs <- (n-j)`nCk`(k-i)
+ , nbCombs <= r = uptoRank i (j+1) (r-nbCombs)
+ | otherwise = j : for1K (i+1) (j+1) r
+
+-- | @combins k xs@ retourne les combinaisons
+-- de longueur 'k' des éléments de la liste 'xs'.
+combins :: Int -> [a] -> [[a]]
+combins k xs = combinsK xs List.!! (length xs - k)
+
+-- | @combinsK xs@ retourne toutes les combinaisons
+-- de longueur allant de @length xs@ à 0 de la liste 'xs',
+--
+-- Algorithme dynamique permettant un calcul de 'combins'
+-- relativement rapide du fait du partage de 'combinsKmoins1'.
+combinsK :: [a] -> [[[a]]]
+combinsK [] = [[[]]]
+combinsK (x : xs) =
+ zipWith (++) ([] : combinsKmoins1)
+ (map (map (x :)) combinsKmoins1 ++ [[]])
+ where combinsKmoins1 = combinsK xs
sourcephile / git / haskell / mpms.git / commitdiff