combinat 0.2.3.1 → 0.2.4

raw patch · 4 files changed

+479/−20 lines, 4 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Math/Combinat/Numbers.hs +2/−18
• Math/Combinat/Numbers/Series.hs +452/−0
• Math/Combinat/Permutations.hs +21/−0
• combinat.cabal +4/−2

@@ -6,6 +6,8 @@  module Math.Combinat.Numbers where +--------------------------------------------------------------------------------+ import Data.Array  --------------------------------------------------------------------------------@@ -124,24 +126,6 @@     f k = toRational (paritySign (n+k) * factorial k * stirling2nd n k)          / toRational (k+1) ------------------------------------------------------------------------------------ * Power series---- | Power series expansion of --- --- > @1 / ( (1-x^a_1) * (1-x^a_2) * ... * (1-x^a_n) )@------ Example:------ @(coinSeries [2,3,5]) !! k@ is the number of ways --- to pay @k@ dollars with coins of two, three and five dollars.------ TODO: better name?-coinSeries :: [Int] -> [Integer]-coinSeries []  = 1 : repeat 0-coinSeries (k:ks) = xs where-  xs = zipWith (+) (coinSeries ks) (replicate k 0 ++ xs) -   --------------------------------------------------------------------------------   

@@ -0,0 +1,452 @@++-- | Some basic power series expansions.+-- This module is not re-exported by "Math.Combinat".+--+-- Note: the \"@convolveWithXXX@\" functions are much faster than the equivalent+-- @(XXX \`convolve\`)@!+-- +-- TODO: better names for these functions.+--++{-# LANGUAGE CPP, GeneralizedNewtypeDeriving #-}+module Math.Combinat.Numbers.Series where++--------------------------------------------------------------------------------++import Data.List++#ifdef QUICKCHECK+import System.Random+import Test.QuickCheck+#endif++--------------------------------------------------------------------------------++-- | The series [1,0,0,0,0,...], which is the neutral element for the convolution.+{-# SPECIALIZE unitSeries :: [Integer] #-}+unitSeries :: Num a => [a]+unitSeries = 1 : repeat 0++-- | Convolution of series. The result is always an infinite list. Warning: This is slow!+convolve :: Num a => [a] -> [a] -> [a]+convolve xs1 ys1 = res where+  res = [ foldl' (+) 0 (zipWith (*) xs (reverse (take n ys)))+        | n<-[1..] +        ]+  xs = xs1 ++ repeat 0+  ys = ys1 ++ repeat 0++-- | Convolution of many series. Still slow!+convolveMany :: Num a => [[a]] -> [a]+convolveMany []  = 1 : repeat 0+convolveMany xss = foldl1 convolve xss++--------------------------------------------------------------------------------+-- * \"Coin\" series++-- | Power series expansion of +-- +-- > 1 / ( (1-x^k_1) * (1-x^k_2) * ... * (1-x^k_n) )+--+-- Example:+--+-- @(coinSeries [2,3,5])!!k@ is the number of ways +-- to pay @k@ dollars with coins of two, three and five dollars.+--+-- TODO: better name?+coinSeries :: [Int] -> [Integer]+coinSeries [] = 1 : repeat 0+coinSeries (k:ks) = xs where+  xs = zipWith (+) (coinSeries ks) (replicate k 0 ++ xs) ++-- | Generalization of the above to include coefficients: expansion of +--  +-- > 1 / ( (1-a_1*x^k_1) * (1-a_2*x^k_2) * ... * (1-a_n*x^k_n) ) +-- +coinSeries' :: Num a => [(a,Int)] -> [a]+coinSeries' [] = 1 : repeat 0+coinSeries' ((a,k):aks) = xs where+  xs = zipWith (+) (coinSeries' aks) (replicate k 0 ++ map (*a) xs) ++convolveWithCoinSeries :: [Int] -> [Integer] -> [Integer]+convolveWithCoinSeries ks series1 = worker ks where+  series = series1 ++ repeat 0+  worker [] = series+  worker (k:ks) = xs where+    xs = zipWith (+) (worker ks) (replicate k 0 ++ xs)++convolveWithCoinSeries' :: Num a => [(a,Int)] -> [a] -> [a]+convolveWithCoinSeries' ks series1 = worker ks where+  series = series1 ++ repeat 0+  worker [] = series+  worker ((a,k):aks) = xs where+    xs = zipWith (+) (worker aks) (replicate k 0 ++ map (*a) xs)++--------------------------------------------------------------------------------+-- * Reciprocals of products of polynomials++-- | Convolution of many 'pseries', that is, the expansion of the reciprocal+-- of a product of polynomials+productPSeries :: [[Int]] -> [Integer]+productPSeries = foldl (flip convolveWithPSeries) unitSeries++-- | The same, with coefficients.+productPSeries' :: Num a => [[(a,Int)]] -> [a]+productPSeries' = foldl (flip convolveWithPSeries') unitSeries++convolveWithProductPSeries :: [[Int]] -> [Integer] -> [Integer]+convolveWithProductPSeries kss ser = foldl (flip convolveWithPSeries) ser kss++-- | This is the most general function in this module; all the others+-- are special cases of this one.  +convolveWithProductPSeries' :: Num a => [[(a,Int)]] -> [a] -> [a] +convolveWithProductPSeries' akss ser = foldl (flip convolveWithPSeries') ser akss+  +--------------------------------------------------------------------------------+-- * Reciprocals of polynomials++-- Reciprocals of polynomials, without coefficients++#ifdef QUICKCHECK+-- | Expansion of @1 / (1-x^k)@. Included for completeness only; +-- it equals to @coinSeries [k]@, and for example+-- for @k=4@ it is simply+-- +-- > [1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0...]+--+pseries1 :: Int -> [Integer]+pseries1 k1 = convolveWithPSeries1 k1 unitSeries ++-- | The expansion of @1 / (1-x^k_1-x^k_2)@+pseries2 :: Int -> Int -> [Integer]+pseries2 k1 k2 = convolveWithPSeries2 k1 k2 unitSeries ++-- | The expansion of @1 / (1-x^k_1-x^k_2-x^k_3)@+pseries3 :: Int -> Int -> Int -> [Integer]+pseries3 k1 k2 k3 = convolveWithPSeries3 k1 k2 k3 unitSeries+#endif++-- | The power series expansion of +--+-- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)+--+pseries :: [Int] -> [Integer]+pseries ks = convolveWithPSeries ks unitSeries++#ifdef QUICKCHECK+-- | Convolve with (the expansion of) @1 / (1-x^k1)@+convolveWithPSeries1 :: Int -> [Integer] -> [Integer]+convolveWithPSeries1 k1 series1 = xs where+  series = series1 ++ repeat 0 +  xs = zipWith (+) series ( replicate k1 0 ++ xs )++-- | Convolve with (the expansion of) @1 / (1-x^k1-x^k2)@+convolveWithPSeries2 :: Int -> Int -> [Integer] -> [Integer]+convolveWithPSeries2 k1 k2 series1 = xs where+  series = series1 ++ repeat 0 +  xs = zipWith3 (\x y z -> x + y + z)+    series+    ( replicate k1 0 ++ xs )+    ( replicate k2 0 ++ xs )+    +-- | Convolve with (the expansion of) @1 / (1-x^k_1-x^k_2-x^k_3)@+convolveWithPSeries3 :: Int -> Int -> Int -> [Integer] -> [Integer]+convolveWithPSeries3 k1 k2 k3 series1 = xs where+  series = series1 ++ repeat 0 +  xs = zipWith4 (\x y z w -> x + y + z + w)+    series+    ( replicate k1 0 ++ xs )+    ( replicate k2 0 ++ xs )+    ( replicate k3 0 ++ xs )+#endif++-- | Convolve with (the expansion of) +--+-- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)+--+convolveWithPSeries :: [Int] -> [Integer] -> [Integer]+convolveWithPSeries ks series1 = ys where +  series = series1 ++ repeat 0 +  ys = worker ks ys +  worker [] _ = series +  worker (k:ks) ys = xs where+    xs = zipWith (+) (replicate k 0 ++ ys) (worker ks ys)++--------------------------------------------------------------------------------+--  Reciprocals of polynomials, with coefficients++#ifdef QUICKCHECK+-- | @1 / (1 - a*x^k)@. +-- For example, for @a=3@ and @k=2@ it is just+-- +-- > [1,0,3,0,9,0,27,0,81,0,243,0,729,0,2187,0,6561,0,19683,0...]+--+pseries1' :: Num a => (a,Int) -> [a]+pseries1' ak1 = convolveWithPSeries1' ak1 unitSeries++-- | @1 / (1 - a_1*x^k_1 - a_2*x^k_2)@+pseries2' :: Num a => (a,Int) -> (a,Int) -> [a]+pseries2' ak1 ak2 = convolveWithPSeries2' ak1 ak2 unitSeries++-- | @1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3)@+pseries3' :: Num a => (a,Int) -> (a,Int) -> (a,Int) -> [a]+pseries3' ak1 ak2 ak3 = convolveWithPSeries3' ak1 ak2 ak3 unitSeries+#endif++-- | The expansion of +--+-- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)+--+pseries' :: Num a => [(a,Int)] -> [a]+pseries' aks = convolveWithPSeries' aks unitSeries++#ifdef QUICKCHECK+-- | Convolve with @1 / (1 - a*x^k)@. +convolveWithPSeries1' :: Num a => (a,Int) -> [a] -> [a]+convolveWithPSeries1' (a1,k1) series1 = xs where+  series = series1 ++ repeat 0 +  xs = zipWith (+)+    series+    ( replicate k1 0 ++ map (*a1) xs )++-- | Convolve with @1 / (1 - a_1*x^k_1 - a_2*x^k_2)@+convolveWithPSeries2' :: Num a => (a,Int) -> (a,Int) -> [a] -> [a]+convolveWithPSeries2' (a1,k1) (a2,k2) series1 = xs where+  series = series1 ++ repeat 0 +  xs = zipWith3 (\x y z -> x + y + z)+    series+    ( replicate k1 0 ++ map (*a1) xs )+    ( replicate k2 0 ++ map (*a2) xs )+    +-- | Convolve with @1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3)@+convolveWithPSeries3' :: Num a => (a,Int) -> (a,Int) -> (a,Int) -> [a] -> [a]+convolveWithPSeries3' (a1,k1) (a2,k2) (a3,k3) series1 = xs where+  series = series1 ++ repeat 0 +  xs = zipWith4 (\x y z w -> x + y + z + w)+    series+    ( replicate k1 0 ++ map (*a1) xs )+    ( replicate k2 0 ++ map (*a2) xs )+    ( replicate k3 0 ++ map (*a3) xs )+#endif++-- | Convolve with (the expansion of) +--+-- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)+--+convolveWithPSeries' :: Num a => [(a,Int)] -> [a] -> [a]+convolveWithPSeries' aks series1 = ys where +  series = series1 ++ repeat 0 +  ys = worker aks ys +  worker [] _ = series+  worker ((a,k):aks) ys = xs where+    xs = zipWith (+) (replicate k 0 ++ map (*a) ys) (worker aks ys)++data Sign = Plus | Minus deriving (Eq,Show)++signValue :: Num a => Sign -> a+signValue Plus  =  1+signValue Minus = -1++signedPSeries :: [(Sign,Int)] -> [Integer] +signedPSeries aks = convolveWithSignedPSeries aks unitSeries++-- | Convolve with (the expansion of) +--+-- > 1 / (1 +- x^k_1 +- x^k_2 +- x^k_3 +- ... +- x^k_n)+--+-- Should be faster than using `convolveWithPSeries'`.+-- Note: 'Plus' corresponds to the coefficient @-1@ in `pseries'` (since+-- there is a minus sign in the definition there)!+convolveWithSignedPSeries :: [(Sign,Int)] -> [Integer] -> [Integer]+convolveWithSignedPSeries aks series1 = ys where +  series = series1 ++ repeat 0 +  ys = worker aks ys +  worker [] _ = series+  worker ((a,k):aks) ys = xs where+    xs = case a of+      Minus -> zipWith (+) one two +      Plus  -> zipWith (-) one two+    one = worker aks ys+    two = replicate k 0 ++ ys+     +--------------------------------------------------------------------------------++#ifdef QUICKCHECK++swap :: (a,b) -> (b,a)+swap (x,y) = (y,x)++-- compare the first 1000 elements of the infinite lists+(=!=) :: (Eq a, Num a) => [a] -> [a] -> Bool+(=!=) xs1 ys1 = (take m xs == take m ys) where +  m = 1000+  xs = xs1 ++ repeat 0+  ys = ys1 ++ repeat 0+infix 4 =!=++newtype Nat = Nat { fromNat :: Int } deriving (Eq,Ord,Show,Num,Random)+newtype Ser = Ser { fromSer :: [Integer] } deriving (Eq,Ord,Show)+newtype Exp  = Exp  { fromExp  ::  Int  } deriving (Eq,Ord,Show,Num,Random)+newtype Exps = Exps { fromExps :: [Int] } deriving (Eq,Ord,Show)+newtype CoeffExp  = CoeffExp  { fromCoeffExp  ::  (Integer,Int)  } deriving (Eq,Ord,Show)+newtype CoeffExps = CoeffExps { fromCoeffExps :: [(Integer,Int)] } deriving (Eq,Ord,Show)++minSerSize = 0    :: Int+maxSerSize = 1000 :: Int++minSerValue = -10000 :: Integer+maxSerValue =  10000 :: Integer++rndList :: (RandomGen g, Random a) => Int -> (a, a) -> g -> ([a], g)+rndList n minmax g = swap $ mapAccumL f g [1..n] where+  f g _ = swap $ randomR minmax g ++instance Arbitrary Nat where+  arbitrary = choose (Nat 0 , Nat 750)++instance Arbitrary Exp where+  arbitrary = choose (Exp 1 , Exp 32)++instance Arbitrary CoeffExp where+  arbitrary = do+    coeff <- choose (minSerValue, maxSerValue) :: Gen Integer+    exp <- arbitrary :: Gen Exp+    return $ CoeffExp (coeff,fromExp exp)+   +instance Random Ser where+  random g = (Ser list, g2) where+    (size,g1) = randomR (minSerSize,maxSerSize) g+    (list,g2) = rndList size (minSerValue,maxSerValue) g1+  randomR _ = random++instance Random Exps where+  random g = (Exps list, g2) where+    (size,g1) = randomR (0,10) g+    (list,g2) = rndList size (1,32) g1+  randomR _ = random++instance Random CoeffExps where+  random g = (CoeffExps (zip list2 list1), g3) where+    (size,g1) = randomR (0,10) g+    (list1,g2) = rndList size (1,32) g1+    (list2,g3) = rndList size (minSerValue,maxSerValue) g2+  randomR _ = random+  +instance Arbitrary Ser where+  arbitrary = choose undefined++instance Arbitrary Exps where+  arbitrary = choose undefined++instance Arbitrary CoeffExps where+  arbitrary = choose undefined+  +-- TODO: quickcheck test properties++checkAll :: IO ()+checkAll = do+  let f :: Testable a => a -> IO ()+      f = quickCheck+{- +  -- these are very slow, because random is slow+  putStrLn "leftIdentity"  ; f prop_leftIdentity+  putStrLn "rightIdentity" ; f prop_rightIdentity+  putStrLn "commutativity" ; f prop_commutativity+  putStrLn "associativity" ; f prop_associativity+-}+  putStrLn "convPSeries1 vs generic" ; f prop_conv1_vs_gen+  putStrLn "convPSeries2 vs generic" ; f prop_conv2_vs_gen+  putStrLn "convPSeries3 vs generic" ; f prop_conv3_vs_gen+  putStrLn "convPSeries1' vs generic" ; f prop_conv1_vs_gen'+  putStrLn "convPSeries2' vs generic" ; f prop_conv2_vs_gen'+  putStrLn "convPSeries3' vs generic" ; f prop_conv3_vs_gen'+  putStrLn "convolve_pseries"  ; f prop_convolve_pseries +  putStrLn "convolve_pseries'" ; f prop_convolve_pseries' +  putStrLn "coinSeries vs pseries"  ; f prop_coin_vs_pseries+  putStrLn "coinSeries vs pseries'" ; f prop_coin_vs_pseries'+     +prop_leftIdentity ser = ( xs =!= unitSeries `convolve` xs ) where +  xs = fromSer ser ++prop_rightIdentity ser = ( unitSeries `convolve` xs =!= xs ) where +  xs = fromSer ser ++prop_commutativity ser1 ser2 = ( xs `convolve` ys =!= ys `convolve` xs ) where +  xs = fromSer ser1+  ys = fromSer ser2++prop_associativity ser1 ser2 ser3 = ( one =!= two ) where+  one = (xs `convolve` ys) `convolve` zs+  two = xs `convolve` (ys `convolve` zs)+  xs = fromSer ser1+  ys = fromSer ser2+  zs = fromSer ser3+  +prop_conv1_vs_gen exp1 ser = ( one =!= two ) where+  one = convolveWithPSeries1 k1 xs +  two = convolveWithPSeries [k1] xs+  k1 = fromExp exp1+  xs = fromSer ser  ++prop_conv2_vs_gen exp1 exp2 ser = (one =!= two) where+  one = convolveWithPSeries2 k1 k2 xs +  two = convolveWithPSeries [k2,k1] xs+  k1 = fromExp exp1+  k2 = fromExp exp2+  xs = fromSer ser  ++prop_conv3_vs_gen exp1 exp2 exp3 ser = (one =!= two) where+  one = convolveWithPSeries3 k1 k2 k3 xs +  two = convolveWithPSeries [k2,k3,k1] xs+  k1 = fromExp exp1+  k2 = fromExp exp2+  k3 = fromExp exp3+  xs = fromSer ser  ++prop_conv1_vs_gen' exp1 ser = ( one =!= two ) where+  one = convolveWithPSeries1' ak1 xs +  two = convolveWithPSeries' [ak1] xs+  ak1 = fromCoeffExp exp1+  xs = fromSer ser  ++prop_conv2_vs_gen' exp1 exp2 ser = (one =!= two) where+  one = convolveWithPSeries2' ak1 ak2 xs +  two = convolveWithPSeries' [ak2,ak1] xs+  ak1 = fromCoeffExp exp1+  ak2 = fromCoeffExp exp2+  xs = fromSer ser  ++prop_conv3_vs_gen' exp1 exp2 exp3 ser = (one =!= two) where+  one = convolveWithPSeries3' ak1 ak2 ak3 xs +  two = convolveWithPSeries' [ak2,ak3,ak1] xs+  ak1 = fromCoeffExp exp1+  ak2 = fromCoeffExp exp2+  ak3 = fromCoeffExp exp3+  xs = fromSer ser  ++prop_convolve_pseries exps1 ser = (one =!= two) where+  one = convolveWithPSeries ks1 xs +  two = xs `convolve` pseries ks1 +  ks1 = fromExps exps1+  xs = fromSer ser  ++prop_convolve_pseries' cexps1 ser = (one =!= two) where+  one = convolveWithPSeries' aks1 xs +  two = xs `convolve` pseries' aks1 +  aks1 = fromCoeffExps cexps1+  xs = fromSer ser  ++prop_coin_vs_pseries exps1 = (one =!= two) where+  one = coinSeries ks1 +  two = convolveMany (map pseries1 ks1)+  ks1 = fromExps exps1++prop_coin_vs_pseries' cexps1 = (one =!= two) where+  one = coinSeries' aks1 +  two = convolveMany (map pseries1' aks1)+  aks1 = fromCoeffExps cexps1+    +#endif ++--------------------------------------------------------------------------------+

@@ -45,6 +45,11 @@   , permuteMultiset   , countPermuteMultiset   , fasc2B_algorithm_L++#ifdef QUICKCHECK+    -- * QuickCheck +  , checkAll+#endif QUICKCHECK   )    where @@ -440,6 +445,22 @@ instance Arbitrary DisjointCycles    where arbitrary = choose undefined instance Arbitrary SameSize          where arbitrary = choose undefined +-- | Runs all quickCheck tests+checkAll :: IO ()+checkAll = do+  let f :: Testable a => a -> IO ()+      f = quickCheck+  f prop_disjcyc1+  f prop_disjcyc2 +  f prop_randCyclic+  f prop_inverse+  f prop_mulPerm+  f prop_mulSign      +  f prop_invMul+  f prop_cyclSign+  f prop_permIsPerm+  f prop_isEven+           prop_disjcyc1 perm = ( perm == disjointCyclesToPermutation n (permutationToDisjointCycles perm) )   where n = permutationSize perm prop_disjcyc2 k dcyc = ( dcyc == permutationToDisjointCycles (disjointCyclesToPermutation n dcyc) )

@@ -1,5 +1,5 @@ Name:                combinat-Version:             0.2.3.1+Version:             0.2.4 Synopsis:            Generation of various combinatorial objects. Description:         A collection of functions to generate combinatorial                      objects like partitions, combinations, permutations,@@ -43,6 +43,7 @@    Exposed-Modules:     Math.Combinat,                         Math.Combinat.Numbers,+                       Math.Combinat.Numbers.Series,                        Math.Combinat.Sets,                        Math.Combinat.Tuples,                         Math.Combinat.Combinations,@@ -57,7 +58,8 @@      Other-Modules:       Math.Combinat.Helper -  Extensions:          MultiParamTypeClasses, ScopedTypeVariables, CPP+  Extensions:          CPP, MultiParamTypeClasses, ScopedTypeVariables, +                       GeneralizedNewtypeDeriving     Hs-Source-Dirs:      .