packages feed

combinatorial (empty) → 0.0

raw patch · 18 files changed

+1994/−0 lines, 18 filesdep +QuickCheckdep +arraydep +basesetup-changed

Dependencies added: QuickCheck, array, base, combinatorial, containers, transformers, utility-ht

Files

+ Changes.md view
@@ -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
+ LICENSE view
@@ -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.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ combinatorial.cabal view
@@ -0,0 +1,77 @@+Name:             combinatorial+Version:          0.0+License:          BSD3+License-File:     LICENSE+Author:           Henning Thielemann <haskell@henning-thielemann.de>+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>+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
+ src/Combinatorics.hs view
@@ -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.+<http://oeis.org/A048993>++> *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.++<http://oeis.org/A000166>+-}+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 <http://oeis.org/A048993>.+-}+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.++<http://oeis.org/A019538>+-}+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)))
+ src/Combinatorics/BellNumbers.hs view
@@ -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
+ src/Combinatorics/CardPairs.hs view
@@ -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
+ src/Combinatorics/Coin.hs view
@@ -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
+ src/Combinatorics/Mastermind.hs view
@@ -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)
+ src/Combinatorics/MaxNim.hs view
@@ -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..]
+ src/Combinatorics/PaperStripGame.hs view
@@ -0,0 +1,87 @@+{- |+Number of possible games as described in+<http://projecteuler.net/problem=306>.+-}+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
+ src/Combinatorics/Partitions.hs view
@@ -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
+ src/Combinatorics/Permutation/WithoutSomeFixpoints.hs view
@@ -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++{- | <http://oeis.org/A047920> -}+numbers :: (Num a) => [[a]]+numbers =+   tail $ scanl (\row fac -> scanl (-) fac row) [] $+   scanl (*) 1 $ iterate (1+) 1
+ src/Combinatorics/TreeDepth.hs view
@@ -0,0 +1,130 @@+module Combinatorics.TreeDepth where++{-+Date: Mon, 18 Apr 2005 18:00:22 +0200+From: Daniel Beer <daniel.beer@informatik.tu-chemnitz.de>+To: Hellseher <lemming@henning-thielemann.de>+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
+ src/Combinatorics/Utility.hs view
@@ -0,0 +1,4 @@+module Combinatorics.Utility where++scalarProduct :: Num a => [a] -> [a] -> a+scalarProduct x y = sum (zipWith (*) x y)
+ src/Polynomial.hs view
@@ -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
+ src/PowerSeries.hs view
@@ -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)
+ test/Test.hs view
@@ -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) :+      []