diff --git a/Changes.md b/Changes.md new file mode 100644 --- /dev/null +++ b/Changes.md @@ -0,0 +1,7 @@ +# Change log for the `combinatorial` package + +## 0.0 + +* Tests: replaced `(==>)` and custom cardinal types by `QC.forAll`. + +* extracted from HTam package diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,27 @@ +Copyright (c) Henning Thielemann 2016 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +SUCH DAMAGE. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/combinatorial.cabal b/combinatorial.cabal new file mode 100644 --- /dev/null +++ b/combinatorial.cabal @@ -0,0 +1,77 @@ +Name: combinatorial +Version: 0.0 +License: BSD3 +License-File: LICENSE +Author: Henning Thielemann +Maintainer: Henning Thielemann +Homepage: http://hub.darcs.net/thielema/combinatorial/ +Category: Math, Statistics +Synopsis: Count, enumerate, rank and unrank combinatorial objects +Description: + Counting, enumerating, ranking and unranking of combinatorial objects. + Well-known and less well-known basic combinatoric problems and examples. + . + The functions are not implemented in obviously stupid ways, + but they are also not optimized to the maximum extent. + The package is plain Haskell 98. + . + See also: + . + * @exact-combinatorics@: + Efficient computations of large combinatoric numbers. + . + * @combinat@: + Library for a similar purpose + with a different structure and selection of problems. +Tested-With: GHC==7.4.2, GHC==7.8.4, GHC==8.0.1 +Cabal-Version: >=1.14 +Build-Type: Simple +Extra-Source-Files: + Changes.md + +Source-Repository this + Tag: 0.0 + Type: darcs + Location: http://hub.darcs.net/thielema/combinatorial/ + +Source-Repository head + Type: darcs + Location: http://hub.darcs.net/thielema/combinatorial/ + +Library + Build-Depends: + containers >=0.4.2 && <0.6, + array >=0.4 && <0.6, + transformers >=0.3 && <0.6, + utility-ht >=0.0.8 && <0.13, + base >=4.5 && <5 + + GHC-Options: -Wall -fwarn-missing-import-lists + Hs-Source-Dirs: src + Default-Language: Haskell98 + Exposed-Modules: + Combinatorics + Combinatorics.Mastermind + Combinatorics.PaperStripGame + Combinatorics.CardPairs + Combinatorics.MaxNim + Combinatorics.TreeDepth + Combinatorics.BellNumbers + Combinatorics.Coin + Combinatorics.Partitions + Combinatorics.Permutation.WithoutSomeFixpoints + Other-Modules: + Combinatorics.Utility + PowerSeries + Polynomial + +Test-Suite combinatorial-test + Type: exitcode-stdio-1.0 + Build-Depends: + combinatorial, + QuickCheck >=2.5 && <3.0, + utility-ht, + base + Main-Is: test/Test.hs + GHC-Options: -Wall + Default-Language: Haskell98 diff --git a/src/Combinatorics.hs b/src/Combinatorics.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics.hs @@ -0,0 +1,581 @@ +{- | +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))) diff --git a/src/Combinatorics/BellNumbers.hs b/src/Combinatorics/BellNumbers.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/BellNumbers.hs @@ -0,0 +1,19 @@ +module Combinatorics.BellNumbers where + +import Combinatorics (binomials, ) +import Combinatorics.Utility (scalarProduct, ) +import qualified PowerSeries + + +{- List of Bell numbers computed with the recursive formula given in + Wurzel 2004-06, page 136 -} +bellRec :: Num a => [a] +bellRec = + 1 : map (scalarProduct bellRec) binomials + +bellSeries :: (Floating a, Enum a) => Int -> a +bellSeries n = + scalarProduct + (map (^n) [0..]) + (take 30 PowerSeries.derivativeCoefficients) + / exp 1 diff --git a/src/Combinatorics/CardPairs.hs b/src/Combinatorics/CardPairs.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/CardPairs.hs @@ -0,0 +1,323 @@ +{- +Compute how often it happens +that a Queen and a King are adjacent in a randomly ordered card set. +-} +module Combinatorics.CardPairs where + +import qualified Combinatorics as Comb + +import Data.Array (Array, (!), array, ) +import Data.Ix (Ix, ) +import qualified Data.List.HT as ListHT + +import qualified Control.Monad.Trans.State as State +import Control.Monad (liftM, liftM2, liftM3, replicateM, ) + +import Data.Ratio ((%), ) + + +type CardSet a = [(a, Int)] + +data Card = Other | Queen | King + deriving (Eq, Ord, Enum, Show) + +charFromCard :: Card -> Char +charFromCard card = + case card of + Other -> ' ' + Queen -> 'q' + King -> 'k' + +removeEach :: State.StateT (CardSet a) [] a +removeEach = + State.StateT $ + map (\(pre,(x,n),post) -> + (x, pre ++ + let m = pred n + in (if m>0 then ((x,m):) else id) + post)) . + ListHT.splitEverywhere + +normalizeSet :: CardSet a -> CardSet a +normalizeSet = filter ((>0) . snd) + +allPossibilities :: CardSet a -> [[a]] +allPossibilities set = + State.evalStateT + (replicateM (sum (map snd set)) removeEach) + (normalizeSet set) + +allPossibilitiesSmall :: [[Card]] +allPossibilitiesSmall = + allPossibilities [(Other, 4), (Queen, 2), (King, 2)] + +allPossibilitiesMedium :: [[Card]] +allPossibilitiesMedium = + allPossibilities [(Other, 4), (Queen, 4), (King, 4)] + +allPossibilitiesSkat :: [[Card]] +allPossibilitiesSkat = + allPossibilities [(Other, 24), (Queen, 4), (King, 4)] + + +adjacentCouple :: [Card] -> Bool +adjacentCouple = + or . + ListHT.mapAdjacent + (\x y -> (x==Queen && y==King) || (x==King && y==Queen)) + +adjacentCouplesSmall :: [[Card]] +adjacentCouplesSmall = + filter adjacentCouple $ + allPossibilities [(Other, 4), (Queen, 2), (King, 2)] + +exampleOutput :: IO () +exampleOutput = + mapM_ (print . map charFromCard) allPossibilitiesSmall + + +{- | +Candidate for utility-ht: +-} +sample :: (a -> b) -> [a] -> [(a,b)] +sample f = map (\x -> (x, f x)) + + +data CardCount i = + CardCount {otherCount, queenCount, kingCount :: i} + deriving (Eq, Ord, Ix, Show) + + +possibilitiesCardsNaive :: + CardCount Int -> Integer +possibilitiesCardsNaive (CardCount no nq nk) = + fromIntegral $ length $ + filter adjacentCouple $ + allPossibilities [(Other,no), (Queen,nq), (King,nk)] + +possibilitiesCardsDynamic :: + CardCount Int -> Array (CardCount Int) Integer +possibilitiesCardsDynamic (CardCount mo mq mk) = + let border = + liftM3 CardCount [0,1] [0..mq] [0..mk] ++ + liftM3 CardCount [0..mo] [0,1] [0..mk] ++ + liftM3 CardCount [0..mo] [0..mq] [0,1] + p = + array (CardCount 0 0 0, CardCount mo mq mk) $ + sample possibilitiesCardsNaive border ++ + sample + (\(CardCount no nq nk) -> + -- " ******" + p!(CardCount (no-1) nq nk) + + -- "q *****" + p!(CardCount (no-1) (nq-1) nk) + + -- "k *****" + p!(CardCount (no-1) nq (nk-1)) + + -- The following case is not handled correctly, + -- because the second 'q' can be part of a "qk". + -- "qq*****" + p!(CardCount no (nq-2) nk) + + -- "kk*****" + p!(CardCount no nq (nk-2)) + + -- "kq*****" + -- "qk*****" + 2 * Comb.multinomial [fromIntegral no, fromIntegral nq-1, fromIntegral nk-1]) + (liftM3 CardCount [2..mo] [2..mq] [2..mk]) + in p + + +sumCard :: Num i => CardCount i -> i +sumCard (CardCount x y z) = x+y+z + +{- +Candidate for utility-ht: slice + +http://hackage.haskell.org/packages/archive/event-list/0.1/doc/html/Data-EventList-Relative-TimeBody.html#v:slice +could be rewritten for plain lists. +-} + +{- | +Count the number of card set orderings +with adjacent queen and king. +We return a triple where the elements count with respect to an additional condition: +(card set starts with an ordinary card ' ', + start with queen 'q', + start with king 'k') +-} +possibilitiesCardsBorderNaive :: + CardCount Int -> CardCount Integer +possibilitiesCardsBorderNaive (CardCount no nq nk) = + foldl (\n (card:_) -> + case card of + Other -> n{otherCount = 1 + otherCount n} + Queen -> n{queenCount = 1 + queenCount n} + King -> n{kingCount = 1 + kingCount n}) + (CardCount 0 0 0) $ + filter adjacentCouple $ + allPossibilities [(Other,no), (Queen,nq), (King,nk)] + +possibilitiesCardsBorderDynamic :: + CardCount Int -> Array (CardCount Int) (CardCount Integer) +possibilitiesCardsBorderDynamic (CardCount mo mq mk) = + let p = + array (CardCount 0 0 0, CardCount mo mq mk) $ + liftM (\ nq -> (CardCount 0 nq 0, CardCount 0 0 0)) [1..mq] ++ + liftM (\ nk -> (CardCount 0 0 nk, CardCount 0 0 0)) [1..mk] ++ + liftM2 (\ nq nk -> ((CardCount 0 nq nk), + let s = fromIntegral $ nq+nk-1 + in CardCount 0 + (Comb.binomial s (fromIntegral nk)) + (Comb.binomial s (fromIntegral nq)))) + [1..mq] [1..mk] ++ + -- (CardCount 0 0 0) is redundant in the list, + -- its number is not needed anyway + liftM2 (\ no nk -> (CardCount no 0 nk, CardCount 0 0 0)) [0..mo] [0..mk] ++ + liftM2 (\ no nq -> (CardCount no nq 0, CardCount 0 0 0)) [0..mo] [0..mq] ++ + sample + (\(CardCount no nq nk) -> + let allP = Comb.multinomial [fromIntegral no, fromIntegral nq-1, fromIntegral nk-1] + in CardCount + (-- " ******" + sumCard (p ! CardCount (no-1) nq nk)) + (-- "q *****" + otherCount (p ! CardCount no (nq-1) nk) + + -- "qq*****" + queenCount (p ! CardCount no (nq-1) nk) + + -- "qk*****" + allP) + (-- "k *****" + otherCount (p ! CardCount no nq (nk-1)) + + -- "kk*****" + kingCount (p ! CardCount no nq (nk-1)) + + -- "kq*****" + allP)) + (liftM3 CardCount [1..mo] [1..mq] [1..mk]) + in p + +possibilitiesCardsBorder2Dynamic :: + CardCount Int -> Array (CardCount Int) (CardCount Integer) +possibilitiesCardsBorder2Dynamic (CardCount mo mq mk) = + let p = + array (CardCount 0 0 0, CardCount mo mq mk) $ + flip sample (liftM3 CardCount [0..mo] [0..mq] [0..mk]) $ + \(CardCount no nq nk) -> + let allP = Comb.multinomial [fromIntegral no, fromIntegral nq-1, fromIntegral nk-1] + test0 n f g = + if n==0 + then 0 + else g $ p ! f (n-1) + in CardCount + (test0 no (\io -> CardCount io nq nk) $ + -- " ******" + sumCard) + (test0 nq (\iq -> CardCount no iq nk) $ \pc -> + -- "q *****" + otherCount pc + + -- "qq*****" + queenCount pc + + -- "qk*****" + allP) + (test0 nk (\ik -> CardCount no nq ik) $ \pc -> + -- "k *****" + otherCount pc + + -- "kk*****" + kingCount pc + + -- "kq*****" + allP) + in p + +{- +for \{o,q,k\} \subset \{1,2,\dots\} +O_{o,q,k} = O_{o-1,q,k} + Q_{o-1,q,k} + K_{o-1,q,k} +Q_{o,q,k} = O_{o,q-1,k} + Q_{o,q-1,k} + M(o,q-1,k-1) +K_{o,q,k} = O_{o,q,k-1} + K_{o,q,k-1} + M(o,q-1,k-1) + +O = (O+Q+K)->(1,0,0) +Q = (O+Q)->(0,1,0) + M->(0,1,1) +K = (O+K)->(0,0,1) + M->(0,1,1) + +O = (O+Q+K)·x +Q = (O+Q)·y + y·z/(1-x-y-z) +K = (O+K)·z + y·z/(1-x-y-z) + +Q·(1-y) = O·y + y·z/(1-x-y-z) +K·(1-z) = O·z + y·z/(1-x-y-z) + +O = (O + (O·y + y·z/(1-x-y-z))/(1-y) + (O·z + y·z/(1-x-y-z))/(1-z))·x +O·(1-x-y-z)·(1-x) + = ((O·y·(1-x-y-z) + y·z)/(1-y) + (O·z·(1-x-y-z) + y·z)/(1-z))·x +O·(1-x-y-z)·(1-x)·(1-y)·(1-z) + = ((O·(1-x-y-z) + z)·y·(1-z) + (O·(1-x-y-z) + y)·z·(1-y))·x +O·(1-x-y-z + (1+x)·y·z)·(1-x-y-z) = x·y·z·(2-y-z) + +O+Q+K = O/x + = y·z·(2-y-z) / (1-x-y-z + (1+x)·y·z) / (1-x-y-z) +-} + + +{- +Pascalsches Dreieck als Potenzreihe von 1/(1-x-y) +ausgerechnet mit Matrizen. + +/n_{0,2}\ /n_{0,1}\ +|n_{1,1}| = |n_{1,0}| +\n_{1,2}/ \n_{1,1}/ + +/n_{1,1}\ /n_{0,1}\ +|n_{2,0}| = |n_{1,0}| +\n_{2,1}/ \n_{1,1}/ +-} + +testCardsBorderDynamic :: + (CardCount Integer, CardCount Integer, CardCount Integer) +testCardsBorderDynamic = + (possibilitiesCardsBorderNaive (CardCount 2 3 5), + possibilitiesCardsBorderDynamic (CardCount 5 5 5) ! (CardCount 2 3 5), + possibilitiesCardsBorder2Dynamic (CardCount 5 5 5) ! (CardCount 2 3 5)) + + +numberOfAllPossibilities :: CardCount Int -> Integer +numberOfAllPossibilities (CardCount no nq nk) = + Comb.multinomial [fromIntegral no, fromIntegral nq, fromIntegral nk] + + +cardSetSizeSkat :: CardCount Int +cardSetSizeSkat = CardCount 24 4 4 + +numberOfPossibilitiesSkat :: Integer +numberOfPossibilitiesSkat = + sumCard $ possibilitiesCardsBorder2Dynamic cardSetSizeSkat ! cardSetSizeSkat + +probabilitySkat :: Double +probabilitySkat = + fromRational $ + numberOfPossibilitiesSkat % numberOfAllPossibilities cardSetSizeSkat + + +cardSetSizeRummy :: CardCount Int +cardSetSizeRummy = CardCount 44 4 4 + +numberOfPossibilitiesRummy :: Integer +numberOfPossibilitiesRummy = + sumCard $ possibilitiesCardsBorder2Dynamic cardSetSizeRummy ! cardSetSizeRummy + +probabilityRummy :: Double +probabilityRummy = + fromRational $ + numberOfPossibilitiesRummy % numberOfAllPossibilities cardSetSizeRummy + + +{- | +Allow both Jack and King adjacent to Queen. +-} +cardSetSizeRummyJK :: CardCount Int +cardSetSizeRummyJK = CardCount 40 4 8 + +numberOfPossibilitiesRummyJK :: Integer +numberOfPossibilitiesRummyJK = + sumCard $ possibilitiesCardsBorder2Dynamic cardSetSizeRummyJK ! cardSetSizeRummyJK + +probabilityRummyJK :: Double +probabilityRummyJK = + fromRational $ + numberOfPossibilitiesRummyJK % numberOfAllPossibilities cardSetSizeRummyJK diff --git a/src/Combinatorics/Coin.hs b/src/Combinatorics/Coin.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/Coin.hs @@ -0,0 +1,21 @@ +{- | +How many possibilities are there for representing an amount of n ct +by the Euro coins 1ct, 2ct, 5ct, 10ct, 20ct, 50ct, 100ct, 200ct? +-} +module Combinatorics.Coin where + +import qualified Data.List as List +import qualified PowerSeries as PS + + +values :: [Int] +values = 1 : 2 : 5 : 10 : 20 : 50 : 100 : 200 : [] + +representationNumbersSingle :: Int -> [Integer] +representationNumbersSingle n = + cycle (1 : List.replicate (n-1) 0) + +representationNumbers :: [Integer] +representationNumbers = + foldl PS.mul PS.one $ + map representationNumbersSingle values diff --git a/src/Combinatorics/Mastermind.hs b/src/Combinatorics/Mastermind.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/Mastermind.hs @@ -0,0 +1,83 @@ +module Combinatorics.Mastermind ( + Eval(..), + evaluate, + evaluateAll, + formatEvalHistogram, + numberDistinct, + ) where + +import qualified Combinatorics.Permutation.WithoutSomeFixpoints as PermWOFP +import Combinatorics (binomial) + +import Text.Printf (printf) + +import qualified Data.Map as Map; import Data.Map (Map) +import qualified Data.Foldable as Fold +import qualified Data.List.HT as ListHT +import Data.Tuple.HT (mapPair) + + +{- | +Cf. @board-games@ package. +-} +data Eval = Eval {black, white :: Int} + deriving (Eq, Ord, Show) + +{- | +Given the code and a guess, compute the evaluation. +-} +evaluate :: (Ord a) => [a] -> [a] -> Eval +evaluate code attempt = + uncurry Eval $ + mapPair + (length, + Fold.sum . uncurry (Map.intersectionWith min) . + mapPair (histogram,histogram) . unzip) $ + ListHT.partition (uncurry (==)) $ + zip code attempt + +{- +*Combinatorics.Mastermind> filter ((Eval 2 0 ==) . evaluate "aabbb") $ replicateM 5 ['a'..'c'] +["aaaaa","aaaac","aaaca","aaacc","aacaa","aacac","aacca","aaccc","acbcc","accbc","acccb","cabcc","cacbc","caccb","ccbbc","ccbcb","cccbb"] +-} + +evaluateAll :: (Ord a) => [[a]] -> [a] -> Map Eval Int +evaluateAll codes attempt = histogram $ map (evaluate attempt) codes + +formatEvalHistogram :: Map Eval Int -> String +formatEvalHistogram m = + let n = maximum $ map (\(Eval b w) -> b+w) $ Map.keys m + in unlines $ + zipWith + (\b -> + unwords . + map (\w -> printf "%6d" $ Map.findWithDefault 0 (Eval b w) m)) + [0..] (reverse $ tail $ ListHT.inits [0..n]) + + +histogram :: (Ord a) => [a] -> Map a Int +histogram = Map.fromListWith (+) . map (\a -> (a,1)) + + +{- | +@numberDistinct n k b w@ computes the number of matching codes, +given that all codes have distinct symbols. +@n@ is the alphabet size, @k@ the width of the code, +@b@ the number of black evaluation sticks and +@w@ the number of white evaluation sticks. +-} +numberDistinct :: Int -> Int -> Int -> Int -> Integer +numberDistinct n k b w = + binomial (toInteger k) (toInteger b) + * + numberDistinctWhite (n-b) (k-b) w + +{- | +@numberDistinctWhite n k w == numberDistinct n k 0 w@ +-} +numberDistinctWhite :: Int -> Int -> Int -> Integer +numberDistinctWhite n k w = + let ni = toInteger n + ki = toInteger k + wi = toInteger w + in binomial ki wi * PermWOFP.numbers !! k !! w * binomial (ni-ki) (ki-wi) diff --git a/src/Combinatorics/MaxNim.hs b/src/Combinatorics/MaxNim.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/MaxNim.hs @@ -0,0 +1,37 @@ +{- | +Simulation of a game with the following rules: + +Players A and B alternatingly take numbers from a set of 2*n numbers. +Player A can choose freely from the remaining numbers, +whereas player B always chooses the maximum remaining number. +How many possibly outcomes of the games exist? +The order in which the numbers are taken is not respected. + +E-Mail by Daniel Beer from 2011-10-24. +-} +module Combinatorics.MaxNim where + +import qualified Data.Set as Set + + +{- | +We only track the number taken by player A +because player B will automatically have the complement set. +-} +gameRound :: (Set.Set Int, Set.Set Int) -> [(Set.Set Int, Set.Set Int)] +gameRound (takenByA, remaining) = do + a <- Set.toList remaining + return (Set.insert a takenByA, Set.deleteMax $ Set.delete a remaining) + +possibilities :: Int -> Set.Set (Set.Set Int) +possibilities n = + Set.fromList $ map fst $ + foldl (>>=) [(Set.empty, Set.fromList [1 .. 2*n])] $ + replicate n gameRound + +{- +This turns out to be the sequence of Catalan numbers. +-} +numberOfPossibilities :: [Int] +numberOfPossibilities = + map (Set.size . possibilities) [0..] diff --git a/src/Combinatorics/PaperStripGame.hs b/src/Combinatorics/PaperStripGame.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/PaperStripGame.hs @@ -0,0 +1,87 @@ +{- | +Number of possible games as described in +. +-} +module Combinatorics.PaperStripGame where + +import qualified Combinatorics as Combi +import qualified PowerSeries as PS +import qualified Data.List.HT as ListHT +import qualified Data.Tree as Tree +import Data.Tree (Tree, ) +import Data.List (inits, tails, ) +import Control.Monad (guard, ) + + +{- +representation: +store the original position of every box +-} +cutEverywhere0 :: [Int] -> [[Int]] +cutEverywhere0 xs = do + (ys, z0:z1:zs) <- zip (inits xs) (tails xs) + guard $ succ z0 == z1 + return $ ys++zs + +{- +representation: +list the sizes of the parts + +cutEverywhere1 [10] ~ cutEverywhere [0..9] +cutEverywhere1 [2,5] ~ cutEverywhere [0,1,3,4,5,6,7] + or cutEverywhere [0,1,4,5,6,7,8] +-} +cutEverywhere1 :: [Int] -> [[Int]] +cutEverywhere1 zs = do + (xs,n,ys) <- ListHT.splitEverywhere zs + (a,b) <- cutPart n + return $ xs ++ filter (0/=) [a,b] ++ ys + +cutPart :: Int -> [(Int, Int)] +cutPart n = + zip [0..] $ takeWhile (>=0) $ iterate pred (n-2) + +treeOfGames :: Int -> Tree [Int] +treeOfGames n = + Tree.unfoldTree (\ns -> (ns, if null ns then [] else cutEverywhere1 ns)) [n] + +lengthOfGames :: Int -> [Int] +lengthOfGames = + let go n ls = + if all (<=1) ls + then [n] + else concatMap (go (succ n)) $ cutEverywhere1 ls + in go 0 . (:[]) + +{- +[1,1,1,2,3,6,12,26,60,144,366,960,2640,7464,21960,66240,206760,660240,2172240,7298640,... +-} +numbersOfGames :: [Int] +numbersOfGames = + map (length . lengthOfGames) [0..] + +{- +directions: + number of boxes -> + length of game v + +That is, the k-th column contains the histogram of (lengthOfGames n). + + | 0 1 2 3 4 5 6 7 8 9 10 +---------------------------------------------- +0 | 1 1 +1 | 1 2 1 +2 | 2 6 6 2 +3 | 6 24 36 24 6 +4 | 24 120 240 +5 | 120 + + +a_n_k = binomial (n+1) (k-2*n) * factorial k +-} + + +numbersOfGamesSeries :: [Integer] +numbersOfGamesSeries = + foldr (\(x0:x1:xs) ys -> x0 : x1 : PS.add xs ys) [] $ + zipWith PS.scale Combi.factorials $ tail Combi.binomials diff --git a/src/Combinatorics/Partitions.hs b/src/Combinatorics/Partitions.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/Partitions.hs @@ -0,0 +1,167 @@ +module Combinatorics.Partitions ( + pentagonalPowerSeries, + numPartitions, + partitionsInc, + partitionsDec, + allPartitionsInc, + + propInfProdLinearFactors, + propPentagonalPowerSeries, + propPentagonalsDifP, + propPentagonalsDifN, + propPartitions, + propNumPartitions, + ) where + +import qualified Data.List as List +import qualified PowerSeries as PS +import Data.Eq.HT (equating) + +{- + a(n) denotes the number in how many ways n can be presented as a sum of + positive integers: + a(n) n + 1 1 : 1 + 2 2 : 2, 1+1 + 3 3 : 3, 2+1, 1+1+1 + 5 4 : 4, 3+1, 2+2, 2+1+1, 1+1+1+1 + 7 5 : 5, 4+1, 3+2, 3+1+1, 2+2+1, 2+1+1+1, 1+1+1+1+1 + + Number of partitions: http://oeis.org/A000041 + Pentagonal numbers: http://oeis.org/A001318 +-} + +{- | +Pentagonal numbers are used to simplify the infinite product +\\prod_{i>0} (1-t^i) +It is known that the coefficients of the power series +are exclusively -1, 0 or 1. +The following is a very simple but inefficient implementation, +because of many multiplications with zero. +-} +prodLinearFactors :: Int -> PS.T Integer +prodLinearFactors n = + foldl PS.mul [1] $ take n $ map (1:) $ iterate (0:) [-1] + +infProdLinearFactors :: PS.T Integer +infProdLinearFactors = + zipWith (!!) + (scanl (\prod i -> delayedSub prod i prod) [1] [1..]) + [0..] + +propInfProdLinearFactors :: Int -> Bool +propInfProdLinearFactors n = + and $ + take (n+1) $ + zipWith (==) + infProdLinearFactors + (prodLinearFactors n) + + +pentagonalsP, pentagonalsN, + pentagonalsDifP, pentagonalsDifN :: [Int] + +pentagonalsP = map (\n -> div (n*(3*n-1)) 2) [0..] +pentagonalsN = map (\n -> div (n*(3*n+1)) 2) [0..] + +{- + (n+1)*(3*n+2) - n*(3*n-1) = 6*n+2 + (n+1)*(3*n+4) - n*(3*n+1) = 6*n+4 +-} +pentagonalsDifP = map (\n -> 3*n+1) [0..] +pentagonalsDifN = map (\n -> 3*n+2) [0..] + +propPentagonalsDifP :: Int -> Bool +propPentagonalsDifP n = + equating (take n) + pentagonalsDifP (zipWith (-) (tail pentagonalsP) pentagonalsP) + +propPentagonalsDifN :: Int -> Bool +propPentagonalsDifN n = + equating (take n) + pentagonalsDifN (zipWith (-) (tail pentagonalsN) pentagonalsN) + +{- + delay y by del and subtract it from x +-} +delayedSub :: [Integer] -> Int -> [Integer] -> [Integer] +delayedSub x del y = + let (a,b) = splitAt del x + in a ++ PS.sub b y + +{- + p00 p01 p02 p03 p04 p05 p06 p07 p08 p09 p10 p11 p12 p13 p14 p15 p16 p17 + - p00 p01 p02 p03 p04 p05 p06 p07 p08 p09 p10 p11 p12 p13 p14 p15 p16 + + p00 p01 p02 p03 p04 p05 p06 p07 p08 p09 p10 p11 p12 + - p00 p01 p02 p03 p04 p05 + ... + - p00 p01 p02 p03 p04 p05 p06 p07 p08 p09 p10 p11 p12 p13 p14 p15 + + p00 p01 p02 p03 p04 p05 p06 p07 p08 p09 p10 + - p00 p01 p02 + ... +-} +numPartitions :: [Integer] +numPartitions = + let accu = foldr (delayedSub numPartitions) (error "never evaluated") + ps = accu (tail pentagonalsDifP) + ns = accu (tail pentagonalsDifN) + in 1 : zipWith (+) ps (0:ns) + +{- | +This is a very efficient implementation of 'prodLinearFactors'. +-} +pentagonalPowerSeries :: [Integer] +pentagonalPowerSeries = + let make = concat . zipWith (\s n -> s : replicate (n-1) 0) (cycle [1,-1]) + in flip PS.sub [1] $ + PS.add + (make pentagonalsDifP) + (make pentagonalsDifN) + +propPentagonalPowerSeries :: Int -> Bool +propPentagonalPowerSeries n = + equating (take n) infProdLinearFactors pentagonalPowerSeries + + + +{- | Give all partitions of the natural number n + with summands which are at least k. + Not quite correct for k>n. -} +partitionsInc :: (Integral a) => a -> a -> [[a]] +partitionsInc k n = + concatMap (\y -> map (y:) (partitionsInc y (n-y))) [k .. div n 2] ++ [[n]] + +partitionsDec :: (Integral a) => a -> a -> [[a]] +partitionsDec 0 0 = [repeat 0] +partitionsDec _ 0 = [] +partitionsDec k n = + (if k>=n then [[n]] else []) ++ + concatMap (\y -> map (y:) (partitionsDec y (n-y))) + (takeWhile (>0) (iterate pred (min n k))) + +_partitionsInc :: (Integral a) => a -> a -> [[a]] +_partitionsInc k n = + if k>n + then [] + else concatMap (\y -> map (y:) (_partitionsInc y (n-y))) [k..(n-1)] + ++ [[n]] + +{- | it shall be k>0 && n>=0 ==> partitionsInc k n == allPartitionsInc !! k !! n + type Int is needed because of list node indexing -} +allPartitionsInc :: [[[[Int]]]] +allPartitionsInc = + let part :: Int -> Int -> [[Int]] + part k n = concatMap (\y -> map (y:) (xs !! y !! (n-y))) + [k .. div n 2] + ++ [[n]] + xs = repeat [[]] : map (\k -> map (part k) [0..]) [1..] + in xs + +propPartitions :: Int -> Int -> Bool +propPartitions k n = + partitionsInc k n == allPartitionsInc !! k !! n + +propNumPartitions :: Int -> Bool +propNumPartitions n = + equating (take n) + (map List.genericLength (allPartitionsInc !! 1)) numPartitions diff --git a/src/Combinatorics/Permutation/WithoutSomeFixpoints.hs b/src/Combinatorics/Permutation/WithoutSomeFixpoints.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/Permutation/WithoutSomeFixpoints.hs @@ -0,0 +1,19 @@ +module Combinatorics.Permutation.WithoutSomeFixpoints where + +import Combinatorics (permute) + +{- | +@enumerate n xs@ list all permutations of @xs@ +where the first @n@ elements do not keep there position +(i.e. are no fixpoints). + +Naive but comprehensible implementation. +-} +enumerate :: (Eq a) => Int -> [a] -> [[a]] +enumerate k xs = filter (and . zipWith (/=) xs . take k) $ permute xs + +{- | -} +numbers :: (Num a) => [[a]] +numbers = + tail $ scanl (\row fac -> scanl (-) fac row) [] $ + scanl (*) 1 $ iterate (1+) 1 diff --git a/src/Combinatorics/TreeDepth.hs b/src/Combinatorics/TreeDepth.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/TreeDepth.hs @@ -0,0 +1,130 @@ +module Combinatorics.TreeDepth where + +{- +Date: Mon, 18 Apr 2005 18:00:22 +0200 +From: Daniel Beer +To: Hellseher +Subject: Baum-Stochastik + + +Nimm folgenden Algorithmus, um einen zufälligen Baum mit n Knoten zu erzeugen: +Starte mit einem einzelnen Knoten (=Wurzel) +Schleife n-1 mal + wähle beliebigen Knoten v1 aus Graph + füge neuen Knoten v2 hinzu + füge Kante (v1,v2) hinzu + +So jetzt die Fragen: +a) Kann man den Erwartungswert für die Tiefe des Baums (also längster Pfad von Wurzel zu einem Blatt) +berechnen? +b) Kann man den Erwartungswert für die Anzahl der Blätter berechnen? +c) Erweiterung von (b). Kann man die zu erwartende Verteilung der Ausgangsgrade berechnen (so eine Art +Histogramm, das angibt wie oft welcher Ausgangsgrad erwartungsgemäß vorkommt)? + +Natürlich alles in Abhängigkeit von n versteht sich. +-} + +import qualified Polynomial as Poly +import qualified Data.Map as Map +import Data.Ratio ((%), ) + +{- Instead of handling probabilities + we make a complete case analysis and + talk only about the absolute frequencies. + That is we start with a one-node tree + then create a new two-node tree from it. + From (n-1)! n-node trees we create n! new (n+1)-node-trees. + + The expectation value of the depth of a node + is the n-th harmonic number. -} + +{-| @nodeDepth !! n !! k@ is the absolute frequency + of nodes with depth k in trees with n nodes. -} +nodeDepth :: [[Integer]] +nodeDepth = scanl (flip nodeDepthIt) [1] [1 ..] + +nodeDepthIt :: Integer -> [Integer] -> [Integer] +nodeDepthIt n = Poly.mul [n,1] + +{-| @treeDepth !! n !! m !! k@ is the absolute frequency + of nodes with depth k in trees with n nodes and depth m. + This can't work - the function carries not enough information + for recursive definition. +treeDepth :: [[[Integer]]] +treeDepth = iterate (\ls -> zipWith treeDepthIt ([[]]++ls) (ls++[[0]])) [[1]] + +treeDepthIt :: [Integer] -> [Integer] -> [Integer] +treeDepthIt nm0 nm1 = + foldl1 add [scale (if null nm0 then 0 else last nm0) (nm0 ++ [1]), + scale (sum (init nm1)) nm1, + 0 : init nm1] +-} + + +{-| + Trees are abstracted to lists of integers, + where each integer denotes the number of nodes + in the corresponding depth of the tree. + The number associated with each tree + is the frequency of this kind of tree + on random tree generation. +-} +type TreeFreq = Map.Map [Integer] Integer + +treeDepth :: [Rational] +treeDepth = + zipWith (%) + (map (sum . map (\(xs,c) -> fromIntegral (length xs) * c) . Map.toList) + treePrototypes) + (scanl (*) 1 [1 ..]) + +treeDepthSeq :: [[Integer]] +treeDepthSeq = + let count = map snd . Map.toList . Map.fromListWith (+) . + map (\(xs,c) -> (length xs, c)) . Map.toList + in map count treePrototypes + +treePrototypes :: [TreeFreq] +treePrototypes = + iterate treeDepthIt (Map.singleton [1] 1) + +extendTree :: [Integer] -> [[Integer]] +extendTree tree = + tail (snd (foldr + (\x (xs,ys) -> (x:xs, ((x+1):xs) : map (x:) ys)) ([],[]) tree)) ++ + [tree ++ [1]] + +treeDepthIt :: TreeFreq -> TreeFreq +treeDepthIt fm = + Map.fromListWith (+) + (concatMap (\(xs,c) -> zip (extendTree xs) (map (c*) xs)) + (Map.toList fm)) + + + +{-| @nodeDegree !! n !! k@ is the number of nodes + with outdegree k in a n-node tree. -} +nodeDegreeProb :: [[Rational]] +nodeDegreeProb = zipWith (\den -> map (%den)) (scanl1 (*) [1 ..]) nodeDegree + +nodeDegree :: [[Integer]] +nodeDegree = + scanl (flip (uncurry nodeDegreeIt)) [1] + (zip [0 ..] (scanl1 (*) [1 ..])) + +nodeDegreeIt :: Integer -> Integer -> [Integer] -> [Integer] +nodeDegreeIt n nFac = Poly.add [nFac] . Poly.mul [n,1] + +{-| expected value of node degree -} +nodeDegreeExpect :: [Rational] +nodeDegreeExpect = + zipWith (%) nodeDegreeExpectAux1 (scanl1 (*) [1 ..]) + +nodeDegreeExpectTrans :: Integer -> [Integer] -> [Integer] +nodeDegreeExpectTrans s x = + scanl (\acc (n,c) -> c + n*acc) s + (zip [1 ..] x) + +nodeDegreeExpectAux0, nodeDegreeExpectAux1 :: [Integer] +nodeDegreeExpectAux0 = nodeDegreeExpectTrans 1 (scanl1 (*) [1 ..]) +nodeDegreeExpectAux1 = nodeDegreeExpectTrans 0 nodeDegreeExpectAux0 diff --git a/src/Combinatorics/Utility.hs b/src/Combinatorics/Utility.hs new file mode 100644 --- /dev/null +++ b/src/Combinatorics/Utility.hs @@ -0,0 +1,4 @@ +module Combinatorics.Utility where + +scalarProduct :: Num a => [a] -> [a] -> a +scalarProduct x y = sum (zipWith (*) x y) diff --git a/src/Polynomial.hs b/src/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Polynomial.hs @@ -0,0 +1,48 @@ +module Polynomial ( + T, fromScalar, add, sub, neg, scale, mul, + differentiate, progression, + ) where + + +type T a = [a] + + +fromScalar :: a -> [a] +fromScalar = (:[]) + +-- | add two polynomials or series +add :: Num a => [a] -> [a] -> [a] +{- zipWith (+) would cut the resulting list + to the length of the shorter operand -} +add [] ys = ys +add xs [] = xs +add (x:xs) (y:ys) = x+y : add xs ys + +-- | subtract two polynomials or series +sub :: Num a => [a] -> [a] -> [a] +sub [] ys = map negate ys +sub xs [] = xs +sub (x:xs) (y:ys) = x-y : sub xs ys + +neg :: Num a => [a] -> [a] +neg = map negate + +-- | scale a polynomial or series by a factor +scale :: Num a => a -> [a] -> [a] +scale s = map (s*) + + +-- | multiply two polynomials or series +mul :: Num a => [a] -> [a] -> [a] +{- prevent from generation of many zeros + if the first operand is the empty list -} +mul [] = const [] +mul xs = foldr (\y zs -> add (scale y xs) (0:zs)) [] + + +progression :: Num a => [a] +progression = iterate (1+) 1 + + +differentiate :: (Num a) => [a] -> [a] +differentiate x = zipWith (*) (tail x) progression diff --git a/src/PowerSeries.hs b/src/PowerSeries.hs new file mode 100644 --- /dev/null +++ b/src/PowerSeries.hs @@ -0,0 +1,20 @@ +module PowerSeries ( + T, fromScalar, one, add, sub, neg, scale, mul, + derivativeCoefficients, differentiate, + ) where + +import Polynomial + (fromScalar, add, sub, neg, scale, mul, + differentiate, progression) + + +type T a = [a] + +one :: Num a => T a +one = fromScalar 1 + + +derivativeCoefficients :: Fractional a => T a +derivativeCoefficients = + scanl (/) 1 progression +-- map recip (scanl (*) 1 progression) diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -0,0 +1,341 @@ +module Main (main) where + +import qualified Combinatorics.Permutation.WithoutSomeFixpoints as PermWOFP +import qualified Combinatorics.Mastermind as Mastermind +import qualified Combinatorics.Partitions as Parts +import qualified Combinatorics.BellNumbers as Bell +import qualified Combinatorics as Comb + +import qualified Test.QuickCheck as QC +import Test.QuickCheck (Testable, quickCheck, ) + +import Control.Monad (liftM2, replicateM, ) +import Control.Applicative ((<$>), ) + +import qualified Data.List.Match as Match +import qualified Data.List.Key as Key +import qualified Data.List as List +import Data.Tuple.HT (uncurry3, ) +import Data.List.HT (allEqual, isAscending, ) +import Data.List (sort, nub, ) +import Data.Eq.HT (equating, ) + + + +permuteSum :: [Int] -> Bool +permuteSum xs = + sum (map sum (Comb.permute xs)) == + sum xs * Comb.factorial (length xs) + +permute :: Ord a => [a] -> Bool +permute xs = + allEqual $ + map (\p -> sort (p xs)) $ + Comb.permute : + Comb.permuteFast : + Comb.permuteShare : + [] + + +genPermuteRep :: QC.Gen [(Char, Int)] +genPermuteRep = do + xns <- QC.listOf $ liftM2 (,) QC.arbitrary $ QC.choose (0,10) + return $ Match.take (takeWhile (<=10) $ scanl1 (+) $ map snd xns) xns + +permuteRepM :: Eq a => [(a, Int)] -> Bool +permuteRepM xs = Comb.permuteRep xs == Comb.permuteRepM xs + +permuteRepNub :: Eq a => [(a, Int)] -> Bool +permuteRepNub xs' = + let xs = Key.nub fst xs' + perms = Comb.permuteRep xs + in perms == nub perms + +permuteRepMonotony :: Ord a => [(a, Int)] -> Bool +permuteRepMonotony = isAscending . Comb.permuteRep . Key.nub fst . sort + +permuteRepChoose :: Int -> Int -> Bool +permuteRepChoose n k = + Comb.choose n k == Comb.permuteRep [(False, n-k), (True, k)] + +chooseLength :: Int -> Int -> Bool +chooseLength n k = + all + (\x -> n == length x && k == length (filter id x)) + (Comb.choose n k) + + +genChooseIndex :: QC.Gen (Integer, Integer, Integer) +genChooseIndex = do + n <- QC.choose (0,25) + k <- QC.choose (0,n) + i <- QC.choose (0, Comb.binomial n k - 1) + return (n,k,i) + +chooseFromIndex :: Integer -> Integer -> Integer -> Bool +chooseFromIndex n k i = + Comb.chooseFromIndex n k i == Comb.chooseFromIndexList n k i + +chooseFromIndexSequence :: Int -> Int -> Bool +chooseFromIndexSequence n k = + map (Comb.chooseFromIndex n k) [0 .. Comb.binomial n k - 1] + == Comb.choose n k + +chooseToFromIndex :: Integer -> Integer -> Integer -> Bool +chooseToFromIndex n k i = + Comb.chooseToIndex (Comb.chooseFromIndex n k i) == (n, k, i) + +chooseFromToIndex :: [Bool] -> Bool +chooseFromToIndex bs = + uncurry3 Comb.chooseFromIndex + (Comb.chooseToIndex bs :: (Integer, Integer, Integer)) + == bs + + + +genVariate :: QC.Gen [Char] +genVariate = take 7 <$> QC.arbitrary + +variateRepMonad :: Eq a => Int -> [a] -> Bool +variateRepMonad n xs = + Comb.variateRep n xs == replicateM n xs + +variatePermute :: Eq a => [a] -> Bool +variatePermute xs = + Comb.variate (length xs) xs == Comb.permute xs + +variatePermuteClip :: Eq a => Int -> [a] -> Bool +variatePermuteClip n xs = + equating (take n) (Comb.variate (length xs) xs) (Comb.permute xs) + +_setPartitionsMonotony :: Ord a => Int -> [a] -> Bool +_setPartitionsMonotony k = + isAscending . Comb.setPartitions k . nub . sort + +rectificationsMonotony :: Ord a => Int -> [a] -> Bool +rectificationsMonotony k = + isAscending . Comb.rectifications k . nub . sort + + + +factorial :: [Char] -> Bool +factorial xs = + length (Comb.permute xs) == Comb.factorial (length xs) + + +binomial :: [Char] -> Int -> Bool +binomial xs k = + length (Comb.tuples k xs) == Comb.binomial (length xs) k + + +genBinomial :: QC.Gen (Integer, Integer) +genBinomial = do + n <- QC.choose (0,100) + k <- QC.choose (0,n) + return (n,k) + +binomialFactorial :: Integer -> Integer -> Bool +binomialFactorial n k = + let (q, r) = + divMod + (Comb.factorial n) + (Comb.factorial k * Comb.factorial (n-k)) + in r == 0 && Comb.binomial n k == q + + +binomialChoose :: Int -> Int -> Bool +binomialChoose n k = + length (Comb.choose n k) == Comb.binomial n k + +multinomialPermuteRep :: [(Char,Int)] -> Bool +multinomialPermuteRep xs = + length (Comb.permuteRep xs) == Comb.multinomial (map snd xs) + +multinomialCommutative :: [Integer] -> Bool +multinomialCommutative xs = + Comb.multinomial xs == Comb.multinomial (sort xs) + +setPartitionNumbers :: Int -> [Int] -> Bool +setPartitionNumbers k xs = + length (Comb.setPartitions k xs) == + (Comb.setPartitionNumbers !! length xs ++ repeat 0) !! k + +rectificationNumbers :: Int -> [Int] -> Bool +rectificationNumbers k xs = + length (Comb.rectifications k xs) == + (Comb.setPartitionNumbers !! k ++ repeat 0) !! length xs + + +surjectiveMappingNumber :: Int -> Bool +surjectiveMappingNumber = + equalFuncList2 Comb.surjectiveMappingNumber Comb.surjectiveMappingNumbers + +surjectiveMappingNumbers :: Int -> Bool +surjectiveMappingNumbers n = + allEqual $ map (take n) $ ( + Comb.surjectiveMappingNumbers : + Comb.surjectiveMappingNumbersStirling : + [] :: [[[Integer]]]) + + +equalFuncList :: (Integer -> Integer) -> [Integer] -> Int -> Bool +equalFuncList f xs n = + equating (take n) xs (map f $ iterate (1+) 0) + +factorials :: Int -> Bool +factorials = equalFuncList Comb.factorial Comb.factorials + +equalFuncList2 :: (Integer -> Integer -> Integer) -> [[Integer]] -> Int -> Bool +equalFuncList2 f xs n = + equating (take n) xs (zipWith (map . f) [0..] $ tail $ List.inits [0..]) + +binomials :: Int -> Bool +binomials = equalFuncList2 Comb.binomial Comb.binomials + +catalanNumbers :: Int -> Bool +catalanNumbers = equalFuncList Comb.catalanNumber Comb.catalanNumbers + +fibonacciNumbers :: Int -> Bool +fibonacciNumbers = equalFuncList Comb.fibonacciNumber Comb.fibonacciNumbers + +derangementNumber :: Int -> Bool +derangementNumber = equalFuncList Comb.derangementNumber Comb.derangementNumbers + +derangementNumbers :: Int -> Bool +derangementNumbers n = + allEqual $ map (take n) $ ( + Comb.derangementNumbers : + Comb.derangementNumbersAlt : + Comb.derangementNumbersInclExcl : + [] :: [[Integer]]) + + +bellSeries :: Int -> Bool +bellSeries = + equalFuncList + (\k -> round (Bell.bellSeries (fromInteger k) :: Double)) + (Bell.bellRec :: [Integer]) + + +genPermutationWOFP :: QC.Gen (Int, String) +genPermutationWOFP = do + xs <- take 6 . nub <$> QC.arbitrary + k <- QC.choose (0, length xs) + return (k,xs) + +permutationWOFP :: Int -> String -> Bool +permutationWOFP k xs = + PermWOFP.numbers !! length xs !! k == length (PermWOFP.enumerate k xs) + +permutationWOFPFactorial :: Int -> Bool +permutationWOFPFactorial k = + Comb.factorial (toInteger k) == PermWOFP.numbers !! k !! 0 + +permutationWOFPDerangement :: Int -> Bool +permutationWOFPDerangement k = + Comb.derangementNumber (toInteger k) == PermWOFP.numbers !! k !! k + + +genMastermindDistinct :: QC.Gen (Int, Int, Int, Int) +genMastermindDistinct = do + n <- QC.choose (0,12) + k <- QC.choose (0, min 5 n) + b <- QC.choose (0,k) + w <- QC.choose (0,k-b) + return (n,k,b,w) + +mastermindDistinct :: Int -> Int -> Int -> Int -> Bool +mastermindDistinct n k b w = + let alphabet = take n ['a'..] + code = take k alphabet + in Mastermind.numberDistinct n k b w == + (toInteger $ length $ + filter ((Mastermind.Eval b w ==) . Mastermind.evaluate code) $ + Comb.variate k alphabet) + + + +testUnit :: Testable prop => String -> prop -> IO () +testUnit label p = putStr (label++": ") >> quickCheck p + +main :: IO () +main = + sequence_ $ + testUnit "permutation sums" + (QC.forAll (take 6 <$> QC.arbitrary) permuteSum) : + testUnit "permutations" + (QC.forAll (take 6 <$> QC.arbitrary :: QC.Gen [Int]) permute) : + testUnit "permuteRepM" + (QC.forAll genPermuteRep permuteRepM) : + testUnit "permuteRepNub" + (QC.forAll genPermuteRep permuteRepNub) : + testUnit "permuteRepMonotony" + (QC.forAll genPermuteRep permuteRepMonotony) : + testUnit "permuteRepChoose" + (QC.forAll (QC.choose (0,10)) permuteRepChoose) : + testUnit "chooseLength" + (QC.forAll (QC.choose (0,10)) chooseLength) : + testUnit "chooseFromIndex" + (QC.forAll genChooseIndex $ uncurry3 chooseFromIndex) : + testUnit "chooseFromIndexSequence" + (QC.forAll (QC.choose (0,10)) chooseFromIndexSequence) : + testUnit "chooseToFromIndex" + (QC.forAll genChooseIndex $ uncurry3 chooseToFromIndex) : + testUnit "chooseFromToIndex" chooseFromToIndex : + testUnit "variation with repetitions with list monad" + (QC.forAll (QC.choose (0,6)) $ \n -> + QC.forAll genVariate $ variateRepMonad n) : + testUnit "variatePermute" (QC.forAll genVariate variatePermute) : + testUnit "permute expressed by variate" + (variatePermuteClip 1000 :: String -> Bool) : + testUnit "binomial vs. choose" + (QC.forAll (QC.choose (0,12)) binomialChoose) : + testUnit "multinomial vs. permutation with repetitions" + (QC.forAll genPermuteRep multinomialPermuteRep) : + testUnit "multinomial commutative" + (QC.forAll (QC.listOf $ QC.choose (0,300)) multinomialCommutative) : + testUnit "factorial vs. permute" + (QC.forAll (take 8 <$> QC.arbitrary) factorial) : + testUnit "binomial vs. tuples" + (QC.forAll (take 16 <$> QC.arbitrary) binomial) : + testUnit "binomial by factorial" + (QC.forAll genBinomial $ uncurry binomialFactorial) : + testUnit "factorial vs. factorials" (factorials 1000) : + testUnit "binomial vs. binomials" (binomials 100) : + testUnit "catalan numbers" (catalanNumbers 1000) : + testUnit "fibonacci numbers" (fibonacciNumbers 10000) : + testUnit "derangement number" (derangementNumber 1000) : + testUnit "derangement numbers" (derangementNumbers 1000) : + testUnit "set partition numbers" + (QC.forAll (QC.choose (0,10000)) $ \n -> + QC.forAll (take 7 <$> QC.arbitrary) $ setPartitionNumbers n) : + testUnit "rectification numbers" + (QC.forAll (QC.choose (0,7)) $ \n xs -> rectificationNumbers n xs) : + testUnit "rectification montony" + (QC.forAll (QC.choose (0,7)) $ \n xs -> + rectificationsMonotony n (xs::[Int])) : + testUnit "surjective mapping number" (surjectiveMappingNumber 20) : + testUnit "surjective mapping numbers" (surjectiveMappingNumbers 20) : + testUnit "bell series" (bellSeries 20) : + testUnit "permutation without some fixpoints" + (QC.forAll genPermutationWOFP $ uncurry permutationWOFP) : + testUnit "permutation without some fixpoints vs. factorial" + (QC.forAll (QC.choose (0,100)) permutationWOFPFactorial) : + testUnit "permutation without some fixpoints vs. derangement" + (QC.forAll (QC.choose (0,100)) permutationWOFPDerangement) : + testUnit "partitions infinite linear factors" + (QC.forAll (QC.choose (0,100)) Parts.propInfProdLinearFactors) : + testUnit "partitions pentagonal power series" + (Parts.propPentagonalPowerSeries 1000) : + testUnit "partitions positive pentagonal numbers" + (Parts.propPentagonalsDifP 10000) : + testUnit "partitions negative pentagonal numbers" + (Parts.propPentagonalsDifN 10000) : + testUnit "partitions" + (QC.forAll (QC.choose (1,10)) $ \k -> + QC.forAll (QC.choose (0,50)) $ \n -> Parts.propPartitions k n) : + testUnit "partitions count" (Parts.propNumPartitions 30) : + testUnit "mastermind with distinct symbols" + (QC.forAll genMastermindDistinct $ \(n,k,b,w) -> + mastermindDistinct n k b w) : + []