module Math.Combinat.Permutations
(
Permutation
, fromPermutation
, toPermutationUnsafe
, isPermutation
, toPermutation
, permutationSize
, permute
, permuteList
, multiply
, inverse
, permutations
, _permutations
, permutationsNaive
, _permutationsNaive
, countPermutations
, randomPermutation
, _randomPermutation
, randomCyclicPermutation
, _randomCyclicPermutation
, randomPermutationDurstenfeld
, randomCyclicPermutationSattolo
, permuteMultiset
, countPermuteMultiset
, fasc2B_algorithm_L
)
where
import Control.Monad
import Control.Monad.ST
import Data.List hiding (permutations)
import Data.Array
import Data.Array.ST
import Math.Combinat.Helper
import System.Random
newtype Permutation = Permutation (Array Int Int) deriving (Eq,Ord,Show,Read)
fromPermutation :: Permutation -> [Int]
fromPermutation (Permutation ar) = elems ar
toPermutationUnsafe :: [Int] -> Permutation
toPermutationUnsafe xs = Permutation perm where
n = length xs
perm = listArray (1,n) xs
isPermutation :: [Int] -> Bool
isPermutation xs = (ar!0 == 0) && and [ ar!j == 1 | j<-[1..n] ] where
n = length xs
ar = accumArray (+) 0 (0,n) $ map f xs
f :: Int -> (Int,Int)
f j = if j<1 || j>n then (0,1) else (j,1)
toPermutation :: [Int] -> Permutation
toPermutation xs = if isPermutation xs
then toPermutationUnsafe xs
else error "toPermutation: not a permutation"
permutationSize :: Permutation -> Int
permutationSize (Permutation ar) = snd $ bounds ar
permute :: Permutation -> Array Int a -> Array Int a
permute pi@(Permutation perm) ar =
if (a==1) && (b==n)
then listArray (1,n) [ ar!(perm!i) | i <- [1..n] ]
else error "permute: array bounds do not match"
where
(_,n) = bounds perm
(a,b) = bounds ar
permuteList :: Permutation -> [a] -> [a]
permuteList perm xs = elems $ permute perm $ listArray (1,n) xs where
n = permutationSize perm
multiply :: Permutation -> Permutation -> Permutation
multiply pi1@(Permutation perm1) (Permutation perm2) =
if (n==m)
then Permutation result
else error "multiply: permutations of different sets"
where
(_,n) = bounds perm1
(_,m) = bounds perm2
result = permute pi1 perm2
infixr 7 `multiply`
inverse :: Permutation -> Permutation
inverse (Permutation perm1) = Permutation result
where
result = array (1,n) $ map swap $ assocs perm1
(_,n) = bounds perm1
permutations :: Int -> [Permutation]
permutations = permutationsNaive
_permutations :: Int -> [[Int]]
_permutations = _permutationsNaive
permutationsNaive :: Int -> [Permutation]
permutationsNaive n = map toPermutationUnsafe $ _permutations n
_permutationsNaive :: Int -> [[Int]]
_permutationsNaive 0 = [[]]
_permutationsNaive 1 = [[1]]
_permutationsNaive n = helper [1..n] where
helper [] = [[]]
helper xs = [ i : ys | i <- xs , ys <- helper (xs `minus` i) ]
minus [] _ = []
minus (x:xs) i = if x < i then x : minus xs i else xs
countPermutations :: Int -> Integer
countPermutations = factorial
randomPermutation :: RandomGen g => Int -> g -> (Permutation,g)
randomPermutation = randomPermutationDurstenfeld
_randomPermutation :: RandomGen g => Int -> g -> ([Int],g)
_randomPermutation n rndgen = (fromPermutation perm, rndgen') where
(perm, rndgen') = randomPermutationDurstenfeld n rndgen
randomCyclicPermutation :: RandomGen g => Int -> g -> (Permutation,g)
randomCyclicPermutation = randomCyclicPermutationSattolo
_randomCyclicPermutation :: RandomGen g => Int -> g -> ([Int],g)
_randomCyclicPermutation n rndgen = (fromPermutation perm, rndgen') where
(perm, rndgen') = randomCyclicPermutationSattolo n rndgen
randomPermutationDurstenfeld :: RandomGen g => Int -> g -> (Permutation,g)
randomPermutationDurstenfeld = randomPermutationDurstenfeldSattolo False
randomCyclicPermutationSattolo :: RandomGen g => Int -> g -> (Permutation,g)
randomCyclicPermutationSattolo = randomPermutationDurstenfeldSattolo True
randomPermutationDurstenfeldSattolo :: RandomGen g => Bool -> Int -> g -> (Permutation,g)
randomPermutationDurstenfeldSattolo isSattolo n rnd = res where
res = runST $ do
ar <- newArray_ (1,n)
forM_ [1..n] $ \i -> writeArray ar i i
rnd' <- worker n (if isSattolo then n1 else n) rnd ar
perm <- unsafeFreeze ar
return (Permutation perm, rnd')
worker :: RandomGen g => Int -> Int -> g -> STUArray s Int Int -> ST s g
worker n m rnd ar =
if n==1
then return rnd
else do
let (k,rnd') = randomR (1,m) rnd
when (k /= n) $ do
y <- readArray ar k
z <- readArray ar n
writeArray ar n y
writeArray ar k z
worker (n1) (m1) rnd' ar
permuteMultiset :: (Eq a, Ord a) => [a] -> [[a]]
permuteMultiset = fasc2B_algorithm_L
countPermuteMultiset :: (Eq a, Ord a) => [a] -> Integer
countPermuteMultiset xs = factorial n `div` product [ factorial (length z) | z <- group ys ]
where
ys = sort xs
n = length xs
fasc2B_algorithm_L :: (Eq a, Ord a) => [a] -> [[a]]
fasc2B_algorithm_L xs = unfold1 next (sort xs) where
next xs = case findj (reverse xs,[]) of
Nothing -> Nothing
Just ( (l:ls) , rs) -> Just $ inc l ls (reverse rs,[])
Just ( [] , _ ) -> error "permute: should not happen"
findj ( xxs@(x:xs) , yys@(y:_) ) = if x >= y
then findj ( xs , x : yys )
else Just ( xxs , yys )
findj ( x:xs , [] ) = findj ( xs , [x] )
findj ( [] , _ ) = Nothing
inc u us ( (x:xs) , yys ) = if u >= x
then inc u us ( xs , x : yys )
else reverse (x:us) ++ reverse (u:yys) ++ xs
inc _ _ ( [] , _ ) = error "permute: should not happen"