#!/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