{- |
Count and create combinatorial objects.
Also see 'combinat' package.
-}
module Combinatorics (
permute,
permuteFast,
permuteShare,
permuteMSL,
runPermuteRep,
permuteRep,
permuteRepM,
choose,
chooseMSL,
variateRep,
variateRepMSL,
variate,
variateMSL,
tuples,
tuplesMSL,
tuplesRec,
partitions,
rectifications,
setPartitions,
chooseFromIndex,
chooseFromIndexList,
chooseFromIndexMaybe,
chooseToIndex,
factorial,
binomial,
binomialSeq,
binomialGen,
binomialSeqGen,
multinomial,
factorials,
binomials,
catalanNumber,
catalanNumbers,
derangementNumber,
derangementNumbers,
derangementNumbersAlt,
derangementNumbersInclExcl,
setPartitionNumbers,
surjectiveMappingNumber,
surjectiveMappingNumbers,
surjectiveMappingNumbersStirling,
fibonacciNumber,
fibonacciNumbers,
) where
import qualified PowerSeries
import Combinatorics.Utility (scalarProduct, )
import Data.Function.HT (nest, )
import Data.Maybe.HT (toMaybe, )
import Data.Maybe (mapMaybe, catMaybes, )
import Data.Tuple.HT (mapFst, )
import qualified Data.List.Match as Match
import Data.List.HT (tails, partition, mapAdjacent, removeEach, splitEverywhere, viewL, )
import Data.List (mapAccumL, intersperse, genericIndex, genericReplicate, genericTake, )
import qualified Control.Monad.Trans.Class as MT
import qualified Control.Monad.Trans.State as MS
import Control.Monad (liftM, liftM2, replicateM, forM, guard, )
{-* Generate compositions from a list of elements. -}
-- several functions for permutation
-- cf. Equation.hs
{- |
Generate list of all permutations of the input list.
The list is sorted lexicographically.
-}
permute :: [a] -> [[a]]
permute [] = [[]]
permute x =
concatMap (\(y, ys) -> map (y:) (permute ys))
(removeEach x)
{- |
Generate list of all permutations of the input list.
It is not lexicographically sorted.
It is slightly faster and consumes less memory
than the lexicographical ordering 'permute'.
-}
permuteFast :: [a] -> [[a]]
permuteFast x = permuteFastStep [] x []
{- |
Each element of (allcycles x) has a different element at the front.
Iterate cycling on the tail elements of each element list of (allcycles x).
-}
permuteFastStep :: [a] -> [a] -> [[a]] -> [[a]]
permuteFastStep suffix [] tl = suffix:tl
permuteFastStep suffix x tl =
foldr (\c -> permuteFastStep (head c : suffix) (tail c)) tl (allCycles x)
{- |
All permutations share as much suffixes as possible.
The reversed permutations are sorted lexicographically.
-}
permuteShare :: [a] -> [[a]]
permuteShare x =
map fst $
-- map (\(y,[]) -> y) $ -- safer but inefficient
nest (length x) (concatMap permuteShareStep) [([], x)]
permuteShareStep :: ([a], [a]) -> [([a], [a])]
permuteShareStep (perm,todo) =
map
(mapFst (:perm))
(removeEach todo)
permuteMSL :: [a] -> [[a]]
permuteMSL xs =
flip MS.evalStateT xs $ replicateM (length xs) $
MS.StateT removeEach
runPermuteRep :: ([(a,Int)] -> [[a]]) -> [(a,Int)] -> [[a]]
runPermuteRep f xs =
let (ps,ns) = partition ((>0) . snd) xs
in if any ((<0) . snd) ns
then []
else f ps
permuteRep :: [(a,Int)] -> [[a]]
permuteRep = runPermuteRep permuteRepAux
permuteRepAux :: [(a,Int)] -> [[a]]
permuteRepAux [] = [[]]
permuteRepAux xs =
concatMap (\(ys,(a,n),zs) ->
let m = pred n
in map (a:) (permuteRepAux (ys ++ (m>0, (a, m)) ?: zs))) $
filter (\(_,(_,n),_) -> n>0) $
splitEverywhere xs
permuteRepM :: [(a,Int)] -> [[a]]
permuteRepM = runPermuteRep permuteRepMAux
permuteRepMAux :: [(a,Int)] -> [[a]]
permuteRepMAux [] = [[]]
permuteRepMAux xs =
do (ys,(a,n),zs) <- splitEverywhere xs
let m = pred n
liftM (a:)
(permuteRepMAux (ys ++ (m>0, (a, m)) ?: zs))
infixr 5 ?:
(?:) :: (Bool, a) -> [a] -> [a]
(True,a) ?: xs = a:xs
(False,_) ?: xs = xs
choose :: Int -> Int -> [[Bool]]
choose n k =
if k<0 || k>n
then []
else
if n==0
then [[]]
else
map (False:) (choose (pred n) k) ++
map (True:) (choose (pred n) (pred k))
chooseMSL :: Int -> Int -> [[Bool]]
chooseMSL n0 k0 =
flip MS.evalStateT k0 $ fmap catMaybes $ sequence $
intersperse (MS.StateT $ \k -> [(Just False, k), (Just True, pred k)]) $
flip map [n0,n0-1..0] $ \n ->
MS.gets (\k -> 0<=k && k<=n) >>= guard >> return Nothing
_chooseMSL :: Int -> Int -> [[Bool]]
_chooseMSL n0 k0 =
flip MS.evalStateT k0 $ do
count <-
forM [n0,n0-1..1] $ \n ->
MS.StateT $ \k ->
guard (0<=k && k<=n) >> [(False, k), (True, pred k)]
MS.gets (0==) >>= guard
return count
{- |
Generate all choices of n elements out of the list x with repetitions.
\"variation\" seems to be used historically,
but I like it more than \"k-permutation\".
-}
variateRep :: Int -> [a] -> [[a]]
variateRep n x = nest n (\y -> concatMap (\z -> map (z:) y) x) [[]]
variateRepMSL :: Int -> [a] -> [[a]]
variateRepMSL = replicateM
{- |
Generate all choices of n elements out of the list x without repetitions.
It holds
@ variate (length xs) xs == permute xs @
-}
variate :: Int -> [a] -> [[a]]
variate 0 _ = [[]]
variate n x =
concatMap (\(y, ys) -> map (y:) (variate (n-1) ys))
(removeEach x)
variateMSL :: Int -> [a] -> [[a]]
variateMSL n xs =
flip MS.evalStateT xs $ replicateM n $
MS.StateT removeEach
{- |
Generate all choices of n elements out of the list x
respecting the order in x and without repetitions.
-}
tuples :: Int -> [a] -> [[a]]
tuples 0 _ = [[]]
tuples r xs =
concatMap (\(y:ys) -> map (y:) (tuples (r-1) ys))
(init (tails xs))
tuplesMSL :: Int -> [a] -> [[a]]
tuplesMSL n xs =
flip MS.evalStateT xs $ replicateM n $
MS.StateT $ mapMaybe viewL . tails
_tuplesMSL :: Int -> [a] -> [[a]]
_tuplesMSL n xs =
flip MS.evalStateT xs $
replicateM n $ do
yl <- MS.get
(y:ys) <- MT.lift $ tails yl
MS.put ys
return y
tuplesRec :: Int -> [a] -> [[a]]
tuplesRec k xt =
if k<0
then []
else
case xt of
[] -> guard (k==0) >> [[]]
x:xs ->
tuplesRec k xs ++
map (x:) (tuplesRec (pred k) xs)
partitions :: [a] -> [([a],[a])]
partitions =
foldr
(\x -> concatMap (\(lxs,rxs) -> [(x:lxs,rxs), (lxs,x:rxs)]))
[([],[])]
{- |
Number of possibilities arising in rectification of a predicate
in deductive database theory.
Stefan Brass, \"Logische Programmierung und deduktive Datenbanken\", 2007,
page 7-60
This is isomorphic to the partition of @n@-element sets
into @k@ non-empty subsets.
> *Combinatorics> map (length . uncurry rectifications) $ do x<-[0..10]; y<-[0..x]; return (x,[1..y::Int])
> [1,0,1,0,1,1,0,1,3,1,0,1,7,6,1,0,1,15,25,10,1,0,1,31,90,65,15,1,0,1,63,301,350,140,21,1,0,1,127,966,1701,1050,266,28,1,0,1,255,3025,7770,6951,2646,462,36,1,0,1,511,9330,34105,42525,22827,5880,750,45,1]
-}
rectifications :: Int -> [a] -> [[a]]
rectifications =
let recourse _ 0 xt =
if null xt
then [[]]
else []
recourse ys n xt =
let n1 = pred n
in liftM2 (:) ys (recourse ys n1 xt) ++
case xt of
[] -> []
(x:xs) -> map (x:) (recourse (ys++[x]) n1 xs)
in recourse []
{- |
Their number is @k^n@.
-}
{-
setPartitionsEmpty :: Int -> [a] -> [[[a]]]
setPartitionsEmpty k =
let recourse [] = [replicate k []]
recourse (x:xs) =
map (\(ys0,y,ys1) -> ys0 ++ [x:y] ++ ys1) $
concatMap splitEverywhere (recourse xs)
{-
do xs1 <- recourse xs
(ys0,y,ys1) <- splitEverywhere xs1
return (ys0 ++ [x:y] ++ ys1)
-}
in recourse
-}
setPartitions :: Int -> [a] -> [[[a]]]
setPartitions 0 xs =
if null xs
then [[]]
else [ ]
setPartitions _ [] = []
setPartitions 1 xs = [[xs]] -- unnecessary for correctness, but useful for efficiency
setPartitions k (x:xs) =
do (rest, choosen) <- partitions xs
part <- setPartitions (pred k) rest
return ((x:choosen) : part)
{-* Compute the number of certain compositions from a number of elements. -}
{- |
@chooseFromIndex n k i == choose n k !! i@
-}
chooseFromIndex :: Integral a => a -> a -> a -> [Bool]
chooseFromIndex n 0 _ = genericReplicate n False
chooseFromIndex n k i =
let n1 = pred n
p = binomial n1 k
b = i>=p
in b :
if b
then chooseFromIndex n1 (pred k) (i-p)
else chooseFromIndex n1 k i
chooseFromIndexList :: Integral a => a -> a -> a -> [Bool]
chooseFromIndexList n k0 i0 =
-- (\((0,0), xs) -> xs) $
snd $
mapAccumL
(\(k,i) bins ->
let p = genericIndex (bins++[0]) k
b = i>=p
in (if b
then (pred k, i-p)
else (k, i),
b))
(k0,i0) $
reverse $
genericTake n binomials
chooseFromIndexMaybe :: Int -> Int -> Int -> Maybe [Bool]
chooseFromIndexMaybe n k i =
toMaybe
(0 <= i && i < binomial n k)
(chooseFromIndex n k i)
-- error ("chooseFromIndex: out of range " ++ show (n, k, i))
chooseToIndex :: Integral a => [Bool] -> (a, a, a)
chooseToIndex =
foldl
(\(n,k0,i0) (bins,b) ->
let (k1,i1) = if b then (succ k0, i0 + genericIndex (bins++[0]) k1) else (k0,i0)
in (succ n, k1, i1))
(0,0,0) .
zip binomials .
reverse
{-* Generate complete lists of combinatorial numbers. -}
factorial :: Integral a => a -> a
factorial n = product [1..n]
{-| Pascal's triangle containing the binomial coefficients. -}
binomial :: Integral a => a -> a -> a
binomial n k =
let bino n' k' =
if k'<0
then 0
else genericIndex (binomialSeq n') k'
in if n<2*k
then bino n (n-k)
else bino n k
binomialSeq :: Integral a => a -> [a]
binomialSeq n =
{- this does not work because the corresponding numbers are not always divisible
product (zipWith div [n', pred n' ..] [1..k'])
-}
scanl (\acc (num,den) -> div (acc*num) den) 1
(zip [n, pred n ..] [1..n])
binomialGen :: (Integral a, Fractional b) => b -> a -> b
binomialGen n k = genericIndex (binomialSeqGen n) k
binomialSeqGen :: (Fractional b) => b -> [b]
binomialSeqGen n =
scanl (\acc (num,den) -> acc*num / den) 1
(zip (iterate (subtract 1) n) (iterate (1+) 1))
multinomial :: Integral a => [a] -> a
multinomial =
product . mapAdjacent binomial . scanr1 (+)
{-* Generate complete lists of factorial numbers. -}
factorials :: Num a => [a]
factorials = scanl (*) 1 (iterate (+1) 1)
{-|
Pascal's triangle containing the binomial coefficients.
Only efficient if a prefix of all rows is required.
It is not efficient for picking particular rows
or even particular elements.
-}
binomials :: Num a => [[a]]
binomials =
let conv11 x = zipWith (+) ([0]++x) (x++[0])
in iterate conv11 [1]
{- |
@catalanNumber n@ computes the number of binary trees with @n@ nodes.
-}
catalanNumber :: Integer -> Integer
catalanNumber n =
let (c,r) = divMod (binomial (2*n) n) (n+1)
in if r==0
then c
else error "catalanNumber: Integer implementation broken"
{- |
Compute the sequence of Catalan numbers by recurrence identity.
It is @catalanNumbers !! n == catalanNumber n@
-}
catalanNumbers :: Num a => [a]
catalanNumbers =
let xs = 1 : PowerSeries.mul xs xs
in xs
derangementNumber :: Integer -> Integer
derangementNumber n =
sum (scanl (*) ((-1) ^ mod n 2) [-n,1-n..(-1)])
{- |
Number of fix-point-free permutations with @n@ elements.
-}
derangementNumbers :: Num a => [a]
derangementNumbers =
-- OEIS-A166: a(n) = n·a(n-1)+(-1)^n
-- y(x) = 1/(1+x) + x · (t -> y(t)·t)'(x)
let xs = PowerSeries.add
(cycle [1,-1])
(0 : PowerSeries.differentiate (0 : xs))
in xs
derangementNumbersAlt :: Num a => [a]
derangementNumbersAlt =
-- OEIS-A166: a(n) = (n-1)·(a(n-1)+a(n-2))
-- y(x) = 1 + x^2 · (t -> y(t)·(1+t))'(x)
let xs =
1 : 0 :
PowerSeries.differentiate
(PowerSeries.add xs (0 : xs))
in xs
derangementNumbersInclExcl :: Num a => [a]
derangementNumbersInclExcl =
let xs = zipWith (-) factorials (map (scalarProduct xs . init) binomials)
in xs
-- generation of all possibilities and computation of their number should be in different modules
{- |
Number of partitions of an @n@ element set into @k@ non-empty subsets.
Known as Stirling numbers .
-}
setPartitionNumbers :: Num a => [[a]]
setPartitionNumbers =
-- s_{n+1,k} = s_{n,k-1} + k·s_{n,k}
iterate (\x -> 0 : PowerSeries.add x (PowerSeries.differentiate x)) [1]
{- |
@surjectiveMappingNumber n k@ computes the number of surjective mappings
from a @n@ element set to a @k@ element set.
-}
surjectiveMappingNumber :: Integer -> Integer -> Integer
surjectiveMappingNumber n k =
foldl subtract 0 $
zipWith (*)
(map (^n) [0..])
(binomialSeq k)
surjectiveMappingNumbers :: Num a => [[a]]
surjectiveMappingNumbers =
iterate
(\x -> 0 : PowerSeries.differentiate
(PowerSeries.add x (0 : x))) [1]
surjectiveMappingNumbersStirling :: Num a => [[a]]
surjectiveMappingNumbersStirling =
map (zipWith (*) factorials) setPartitionNumbers
{- |
Multiply two Fibonacci matrices, that is matrices of the form
> /F[n-1] F[n] \
> \F[n] F[n+1]/
-}
fiboMul ::
(Integer,Integer,Integer) ->
(Integer,Integer,Integer) ->
(Integer,Integer,Integer)
fiboMul (f0,f1,f2) (g0,g1,g2) =
let h0 = f0*g0 + f1*g1
h1 = f0*g1 + f1*g2
-- h1 = f1*g0 + f2*g1
h2 = f1*g1 + f2*g2
in (h0,h1,h2)
{-
Fast computation using matrix power of
> /0 1\
> \1 1/
Hard-coded fast power with integer exponent.
Better use a generic algorithm.
-}
fibonacciNumber :: Integer -> Integer
fibonacciNumber x =
let aux 0 = (1,0,1)
aux (-1) = (-1,1,0)
aux n =
let (m,r) = divMod n 2
f = aux m
f2 = fiboMul f f
in if r==0
then f2
else fiboMul (0,1,1) f2
(_,y,_) = aux x
in y
{- |
Number of possibilities to compose a 2 x n rectangle of n bricks.
> ||| |-- --|
> ||| |-- --|
-}
fibonacciNumbers :: [Integer]
fibonacciNumbers =
let xs = 0 : ys
ys = 1 : zipWith (+) xs ys
in xs
{- * Auxiliary functions -}
{- candidates for Useful -}
{- | Create a list of all possible rotations of the input list. -}
allCycles :: [a] -> [[a]]
allCycles x =
Match.take x (map (Match.take x) (iterate tail (cycle x)))