{-# LANGUAGE GeneralizedNewtypeDeriving,ScopedTypeVariables #-}
module Weyl where


import           Control.Applicative  ((<$>))
import           Control.Monad        (replicateM)
import           Control.Monad.Random (Random(random),evalRandIO,RandomGen,MonadRandom,Rand,StdGen,runRand,getRandom,getRandomR)
import           Data.List            (sort,foldl')
import qualified Data.Set       as Set
import           Prelude hiding ((*),Rational,signum)

import qualified Algebra.Additive
import qualified Algebra.Module
import qualified Algebra.Ring
import qualified MathObj.Matrix as Matrix
import qualified Number.Ratio
import           NumericPrelude ((*),(*>),Rational,(%),signum,numerator,denominator)


-- A value of type 'G' is an element of the Weyl group.
newtype G
  = G { unG :: Matrix.T (Number.Ratio.T Int) }
  deriving (Eq,Ord,Algebra.Additive.C,Algebra.Ring.C)

instance Show G where
  show = Matrix.format . unG

identity :: G
identity = G $ Matrix.one 8

hadamard :: G
hadamard = G $ (1 % (2 :: Int)) *> Matrix.fromRows 8 8
  [
    [ 1, 1, 1, 1] ++ [ 0, 0, 0, 0]
  , [ 1,-1, 1,-1] ++ [ 0, 0, 0, 0]
  , [ 1, 1,-1,-1] ++ [ 0, 0, 0, 0]
  , [ 1,-1,-1, 1] ++ [ 0, 0, 0, 0]
  , [ 0, 0, 0, 0] ++ [ 1, 1, 1, 1]
  , [ 0, 0, 0, 0] ++ [ 1,-1, 1,-1]
  , [ 0, 0, 0, 0] ++ [ 1, 1,-1,-1]
  , [ 0, 0, 0, 0] ++ [ 1,-1,-1, 1]
  ]

instance Random G where
  random = runRand $ iter identity where
    iter x = do
      again <- getRandom
      if again
        then do
          r1 : rs <- permute . Matrix.rows . unG $ x
          signs <- replicateM 7 $ getRandomR (0,1)
          let s1 = sum signs :: Int
          iter $ (hadamard *) . G . Matrix.fromRows 8 8 $
            (((-1) ^ s1 % (1 :: Int)) *> r1)
           :
            zipWith (\ s r -> ((-1) ^ s % (1 :: Int)) *> r) signs rs
        else return x

-- H acts from the left on the elements of a (right) coset, by permuting and
-- sign-changing the rows of the matrix. To normalise a matrix, first permute
-- the rows, such that, after taking the lex-absolute value of each, they
-- appear lex-nondecreasing in the matrix; after that flip the sign of any row
-- but the first that is lex-negative, and flip the sign of the first row, if
-- one would otherwise have changed an odd number of signs.
normalise :: G -> G
normalise = G . Matrix.fromRows 8 8 . changeSigns . sort' . Matrix.rows . unG where
  sort' = map snd . sort . map (\ x -> (sign x *> x,x))
  changeSigns (r1 : rs) = case unzip . map (\ r -> let s = sign r in (s *> r,s)) $ rs of
    (rs',signs) -> foldl' (*) 1 signs *> r1 : rs'
  sign = signum . head . dropWhile (== 0)

representatives :: forall m. (MonadRandom m,Functor m) => m [G]
representatives = Set.toList <$> go (Set.singleton $ normalise identity) where
  go :: Set.Set G -> m (Set.Set G)
  go s = if Set.size s == 135 then return s else do
    g <- getRandom
    go $ Set.insert (normalise g) s

main :: IO ()
main = do
  rs <- evalRandIO representatives
  mapM_ (putStrLn . Matrix.format . fmap (\ x -> 4 * numerator x `div` denominator x) . unG) rs

-- Random permutation

type RandomList g a = Rand g [a]

empty :: RandomList g a
empty = return []

singleton :: a -> RandomList g a
singleton x = return [x]

merge :: (RandomGen g) => RandomList g a -> RandomList g a -> RandomList g a
merge rxs rys = do
  xs <- rxs
  ys <- rys
  merge' (length xs, xs) (length ys, ys) where
    merge' (0 , [])   (_ , ys)   = return ys
    merge' (_ , xs)   (0 , [])   = return xs
    merge' (nx, x:xs) (ny, y:ys) = do
      k <- getRandomR (1,nx+ny)
      if k <= nx
        then (x:) <$> ((nx-1, xs) `merge'` (ny, y:ys))
        else (y:) <$> ((nx, x:xs) `merge'` (ny-1, ys))

permute :: (RandomGen g) => [a] -> RandomList g a
permute = fromList where
  fromList []  = empty
  fromList [x] = singleton x
  fromList xs  = (fromList l) `merge` (fromList r) where
    (l,r) = splitAt (length xs `div` 2) xs