combinat 0.2.8.2 → 0.2.9.0

raw patch · 30 files changed

+3556/−766 lines, 30 filesdep +tastydep +tasty-hunitdep +tasty-quickcheckdep ~basedep ~containersPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: tasty, tasty-hunit, tasty-quickcheck

Dependency ranges changed: base, containers

API changes (from Hackage documentation)

Files

• LICENSE +1/−1
• Math/Combinat/Helper.hs +10/−1
• Math/Combinat/Numbers.hs +6/−191
• Math/Combinat/Numbers/Integers.hs +113/−0
• Math/Combinat/Numbers/Primes.hs +1/−96
• Math/Combinat/Numbers/Sequences.hs +198/−0
• Math/Combinat/Numbers/Series.hs +76/−18
• Math/Combinat/Partitions/Integer.hs +94/−416
• Math/Combinat/Partitions/Integer/Compact.hs +819/−0
• Math/Combinat/Partitions/Integer/Count.hs +215/−0
• Math/Combinat/Partitions/Integer/IntList.hs +398/−0
• Math/Combinat/Partitions/Integer/Naive.hs +202/−0
• Math/Combinat/Partitions/Skew.hs +18/−0
• Math/Combinat/Partitions/Skew/Ribbon.hs +364/−0
• Math/Combinat/Permutations.hs +41/−0
• Math/Combinat/Sets/VennDiagrams.hs +150/−0
• Math/Combinat/Sign.hs +27/−1
• Math/Combinat/Tableaux.hs +3/−2
• Math/Combinat/Tableaux/Skew.hs +2/−1
• cbits/c_compact_partition.c +24/−0
• combinat.cabal +33/−9
• test/TestSuite.hs +3/−1
• test/Tests/Common.hs +4/−0
• test/Tests/Partitions/Compact.hs +390/−0
• test/Tests/Partitions/Integer.hs +29/−1
• test/Tests/Partitions/Ribbon.hs +86/−0
• test/Tests/Partitions/Skew.hs +31/−2
• test/Tests/Permutations.hs +28/−0
• test/Tests/Series.hs +184/−21
• test/Tests/SkewTableaux.hs +6/−5

@@ -1,4 +1,4 @@-Copyright (c) 2008-2016, Balazs Komuves+Copyright (c) 2008-2018, Balazs Komuves All rights reserved.  Redistribution and use in source and binary forms, with or without

@@ -52,6 +52,12 @@ sum' :: Num a => [a] -> a sum' = foldl' (+) 0 +interleave :: [a] -> [a] -> [a]+interleave (x:xs) (y:ys) = x : y : interleave xs ys+interleave [x]    []     = x : []+interleave []     []     = []+interleave _      _      = error "interleave: shouldn't happen"+ -------------------------------------------------------------------------------- -- * equality and ordering  @@ -63,8 +69,11 @@ reverseOrdering GT = LT reverseOrdering EQ = EQ +reverseComparing :: Ord b => (a -> b) -> a -> a -> Ordering+reverseComparing f x y = compare (f y) (f x)+ reverseCompare :: Ord a => a -> a -> Ordering-reverseCompare x y = reverseOrdering $ compare x y+reverseCompare x y = compare y x   -- reverseOrdering $ compare x y  reverseSort :: Ord a => [a] -> [a] reverseSort = sortBy reverseCompare

@@ -1,194 +1,9 @@ --- | A few important number sequences. ---  --- See the \"On-Line Encyclopedia of Integer Sequences\",--- <https://oeis.org> .--module Math.Combinat.Numbers where------------------------------------------------------------------------------------import Data.Array--import Math.Combinat.Helper ( sum' )-import Math.Combinat.Sign-------------------------------------------------------------------------------------- | A000142.-factorial :: Integral a => a -> Integer-factorial n-  | n <  0    = error "factorial: input should be nonnegative"-  | n == 0    = 1-  | otherwise = product [1..fromIntegral n]---- | A006882.-doubleFactorial :: Integral a => a -> Integer-doubleFactorial n-  | n <  0    = error "doubleFactorial: input should be nonnegative"-  | n == 0    = 1-  | odd n     = product [1,3..fromIntegral n]-  | otherwise = product [2,4..fromIntegral n]---- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below.-binomial :: Integral a => a -> a -> Integer-binomial n k -  | k > n = 0-  | k < 0 = 0-  | k > (n `div` 2) = binomial n (n-k)-  | otherwise = (product [n'-k'+1 .. n']) `div` (product [1..k'])-  where -    k' = fromIntegral k-    n' = fromIntegral n---- | The extension of the binomial function to negative inputs. This should satisfy the following properties:------ > for n,k >=0 : signedBinomial n k == binomial n k--- > for any n,k : signedBinomial n k == signedBinomial n (n-k) --- > for k >= 0  : signedBinomial (-n) k == (-1)^k * signedBinomial (n+k-1) k------ Note: This is compatible with Mathematica's @Binomial@ function.----signedBinomial :: Int -> Int -> Integer-signedBinomial n k-  | n >= 0     = binomial n k-  | k >= 0     = negateIfOdd    k  $ binomial (k-n-1)   k  -  | otherwise  = negateIfOdd (n+k) $ binomial (-k-1) (-n-1)--{--test_signed_0 = [ signedBinomial ( n) k == signedBinomial ( n) ( n-k)                | n<-[-30..40] , k<-[-30..40] ]-test_signed_1 = [ signedBinomial (-n) k == signedBinomial (-n) (-n-k)                | n<-[-30..40] , k<-[-30..40] ]-test_signed_2 = [ signedBinomial (-n) k == negateIfOdd k $ signedBinomial (n+k-1) k  | n<-[-30..40] , k<-[0..30] ]--}---- | A given row of the Pascal triangle; equivalent to a sequence of binomial --- numbers, but much more efficient. You can also left-fold over it.------ > pascalRow n == [ binomial n k | k<-[0..n] ]-pascalRow :: Integral a => a -> [Integer]-pascalRow n' = worker 0 1 where-  n = fromIntegral n'-  worker j x-    | j>n   = [] -    | True  = let j'=j+1 in x : worker j' (div (x*(n-j)) j') --multinomial :: Integral a => [a] -> Integer-multinomial xs = div-  (factorial (sum xs))-  (product [ factorial x | x<-xs ])  -  ------------------------------------------------------------------------------------ * Catalan numbers---- | Catalan numbers. OEIS:A000108.-catalan :: Integral a => a -> Integer-catalan n -  | n < 0     = 0-  | otherwise = binomial (n+n) n `div` fromIntegral (n+1)---- | Catalan's triangle. OEIS:A009766.--- Note:------ > catalanTriangle n n == catalan n--- > catalanTriangle n k == countStandardYoungTableaux (toPartition [n,k])----catalanTriangle :: Integral a => a -> a -> Integer-catalanTriangle n k-  | k > n     = 0-  | k < 0     = 0-  | otherwise = (binomial (n+k) n * fromIntegral (n-k+1)) `div` fromIntegral (n+1)------------------------------------------------------------------------------------- * Stirling numbers---- | Rows of (signed) Stirling numbers of the first kind. OEIS:A008275.--- Coefficients of the polinomial @(x-1)*(x-2)*...*(x-n+1)@.--- This function uses the recursion formula.-signedStirling1stArray :: Integral a => a -> Array Int Integer-signedStirling1stArray n-  | n <  1    = error "stirling1stArray: n should be at least 1"-  | n == 1    = listArray (1,1 ) [1]-  | otherwise = listArray (1,n') [ lkp (k-1) - fromIntegral (n-1) * lkp k | k<-[1..n'] ] +module Math.Combinat.Numbers +  ( module Math.Combinat.Numbers.Sequences+  , module Math.Combinat.Numbers.Integers+  )    where-    prev = signedStirling1stArray (n-1)-    n' = fromIntegral n :: Int-    lkp j | j <  1    = 0-          | j >= n'   = 0-          | otherwise = prev ! j -        --- | (Signed) Stirling numbers of the first kind. OEIS:A008275.--- This function uses 'signedStirling1stArray', so it shouldn't be used--- to compute /many/ Stirling numbers.------ Argument order: @signedStirling1st n k@----signedStirling1st :: Integral a => a -> a -> Integer-signedStirling1st n k -  | k==0 && n==0 = 1-  | k < 1        = 0-  | k > n        = 0-  | otherwise    = signedStirling1stArray n ! (fromIntegral k) --- | (Unsigned) Stirling numbers of the first kind. See 'signedStirling1st'.-unsignedStirling1st :: Integral a => a -> a -> Integer-unsignedStirling1st n k = abs (signedStirling1st n k)---- | Stirling numbers of the second kind. OEIS:A008277.--- This function uses an explicit formula.--- --- Argument order: @stirling2nd n k@----stirling2nd :: Integral a => a -> a -> Integer-stirling2nd n k -  | k==0 && n==0 = 1-  | k < 1        = 0-  | k > n        = 0-  | otherwise = sum xs `div` factorial k where-      xs = [ negateIfOdd (k-i) $ binomial k i * (fromIntegral i)^n | i<-[0..k] ]------------------------------------------------------------------------------------- * Bernoulli numbers---- | Bernoulli numbers. @bernoulli 1 == -1%2@ and @bernoulli k == 0@ for--- k>2 and /odd/. This function uses the formula involving Stirling numbers--- of the second kind. Numerators: A027641, denominators: A027642.-bernoulli :: Integral a => a -> Rational-bernoulli n -  | n <  0    = error "bernoulli: n should be nonnegative"-  | n == 0    = 1-  | n == 1    = -1/2-  | otherwise = sum [ f k | k<-[1..n] ] -  where-    f k = toRational (negateIfOdd (n+k) $ factorial k * stirling2nd n k) -        / toRational (k+1)------------------------------------------------------------------------------------- * Bell numbers---- | Bell numbers (Sloane's A000110) from B(0) up to B(n). B(0)=B(1)=1, B(2)=2, etc. ------ The Bell numbers count the number of /set partitions/ of a set of size @n@--- --- See <http://en.wikipedia.org/wiki/Bell_number>----bellNumbersArray :: Integral a => a -> Array Int Integer-bellNumbersArray nn = arr where-  arr = array (0::Int,n) kvs -  n = fromIntegral nn :: Int-  kvs = (0,1) : [ (k, f k) | k<-[1..n] ] -  f n = sum' [ binomial (n-1) k * arr ! k | k<-[0..n-1] ]---- | The n-th Bell number B(n), using the Stirling numbers of the second kind.--- This may be slower than using 'bellNumbersArray'.-bellNumber :: Integral a => a -> Integer-bellNumber nn-  | n <  0     = error "bellNumber: expecting a nonnegative index"-  | n == 0     = 1-  | otherwise  = sum' [ stirling2nd n k | k<-[1..n] ] -  where-    n = fromIntegral nn :: Int------------------------------------------------------------------------------------- +import Math.Combinat.Numbers.Sequences+import Math.Combinat.Numbers.Integers

@@ -0,0 +1,113 @@++-- | Operations on integers++module Math.Combinat.Numbers.Integers +  ( -- * Integer logarithm+    integerLog2+  , ceilingLog2+    -- * Integer square root+  , isSquare+  , integerSquareRoot+  , ceilingSquareRoot+  , integerSquareRoot' +  , integerSquareRootNewton'+  )+  where++--------------------------------------------------------------------------------++-- import Math.Combinat.Numbers++import Data.List ( group , sort )+import Data.Bits++import System.Random++--------------------------------------------------------------------------------+-- Integer logarithm++-- | Largest integer @k@ such that @2^k@ is smaller or equal to @n@+integerLog2 :: Integer -> Integer+integerLog2 n = go n where+  go 0 = -1+  go k = 1 + go (shiftR k 1)++-- | Smallest integer @k@ such that @2^k@ is larger or equal to @n@+ceilingLog2 :: Integer -> Integer+ceilingLog2 0 = 0+ceilingLog2 n = 1 + go (n-1) where+  go 0 = -1+  go k = 1 + go (shiftR k 1)+  +--------------------------------------------------------------------------------+-- Integer square root++isSquare :: Integer -> Bool+isSquare n = +  if (fromIntegral $ mod n 32) `elem` rs +    then snd (integerSquareRoot' n) == 0+    else False+  where+    rs = [0,1,4,9,16,17,25] :: [Int]+    +-- | Integer square root (largest integer whose square is smaller or equal to the input)+-- using Newton's method, with a faster (for large numbers) inital guess based on bit shifts.+integerSquareRoot :: Integer -> Integer+integerSquareRoot = fst . integerSquareRoot'++-- | Smallest integer whose square is larger or equal to the input+ceilingSquareRoot :: Integer -> Integer+ceilingSquareRoot n = (if r>0 then u+1 else u) where (u,r) = integerSquareRoot' n ++-- | We also return the excess residue; that is+--+-- > (a,r) = integerSquareRoot' n+-- +-- means that+--+-- > a*a + r = n+-- > a*a <= n < (a+1)*(a+1)+integerSquareRoot' :: Integer -> (Integer,Integer)+integerSquareRoot' n+  | n<0 = error "integerSquareRoot: negative input"+  | n<2 = (n,0)+  | otherwise = go firstGuess +  where+    k = integerLog2 n+    firstGuess = 2^(div (k+2) 2) -- !! note that (div (k+1) 2) is NOT enough !!+    go a = +      if m < a+        then go a' +        else (a, r + a*(m-a))+      where+        (m,r) = divMod n a+        a' = div (m + a) 2++-- | Newton's method without an initial guess. For very small numbers (<10^10) it+-- is somewhat faster than the above version.+integerSquareRootNewton' :: Integer -> (Integer,Integer)+integerSquareRootNewton' n+  | n<0 = error "integerSquareRootNewton: negative input"+  | n<2 = (n,0)+  | otherwise = go (div n 2) +  where+    go a = +      if m < a+        then go a' +        else (a, r + a*(m-a))+      where+        (m,r) = divMod n a+        a' = div (m + a) 2++{-+-- brute force test of integer square root+isqrt_test n1 n2 = +  [ k +  | k<-[n1..n2] +  , let (a,r) = integerSquareRoot' k+  , (a*a+r/=k) || (a*a>k) || (a+1)*(a+1)<=k +  ]+-}++--------------------------------------------------------------------------------+

@@ -9,15 +9,6 @@     -- * Prime factorization   , groupIntegerFactors   , integerFactorsTrialDivision-    -- * Integer logarithm-  , integerLog2-  , ceilingLog2-    -- * Integer square root-  , isSquare-  , integerSquareRoot-  , ceilingSquareRoot-  , integerSquareRoot' -  , integerSquareRootNewton'     -- * Modulo @m@ arithmetic   , powerMod     -- * Prime testing@@ -29,7 +20,7 @@  -------------------------------------------------------------------------------- --- import Math.Combinat.Numbers+import Math.Combinat.Numbers.Integers  import Data.List ( group , sort ) import Data.Bits@@ -119,92 +110,6 @@ -- brute force testing of factors ifactorsTest :: (Integer -> [Integer]) -> Integer -> Bool ifactorsTest alg n = and [ product (alg k) == k | k<-[1..n] ]   --}------------------------------------------------------------------------------------- Integer logarithm---- | Largest integer @k@ such that @2^k@ is smaller or equal to @n@-integerLog2 :: Integer -> Integer-integerLog2 n = go n where-  go 0 = -1-  go k = 1 + go (shiftR k 1)---- | Smallest integer @k@ such that @2^k@ is larger or equal to @n@-ceilingLog2 :: Integer -> Integer-ceilingLog2 0 = 0-ceilingLog2 n = 1 + go (n-1) where-  go 0 = -1-  go k = 1 + go (shiftR k 1)-  ------------------------------------------------------------------------------------ Integer square root--isSquare :: Integer -> Bool-isSquare n = -  if (fromIntegral $ mod n 32) `elem` rs -    then snd (integerSquareRoot' n) == 0-    else False-  where-    rs = [0,1,4,9,16,17,25] :: [Int]-    --- | Integer square root (largest integer whose square is smaller or equal to the input)--- using Newton's method, with a faster (for large numbers) inital guess based on bit shifts.-integerSquareRoot :: Integer -> Integer-integerSquareRoot = fst . integerSquareRoot'---- | Smallest integer whose square is larger or equal to the input-ceilingSquareRoot :: Integer -> Integer-ceilingSquareRoot n = (if r>0 then u+1 else u) where (u,r) = integerSquareRoot' n ---- | We also return the excess residue; that is------ > (a,r) = integerSquareRoot' n--- --- means that------ > a*a + r = n--- > a*a <= n < (a+1)*(a+1)-integerSquareRoot' :: Integer -> (Integer,Integer)-integerSquareRoot' n-  | n<0 = error "integerSquareRoot: negative input"-  | n<2 = (n,0)-  | otherwise = go firstGuess -  where-    k = integerLog2 n-    firstGuess = 2^(div (k+2) 2) -- !! note that (div (k+1) 2) is NOT enough !!-    go a = -      if m < a-        then go a' -        else (a, r + a*(m-a))-      where-        (m,r) = divMod n a-        a' = div (m + a) 2---- | Newton's method without an initial guess. For very small numbers (<10^10) it--- is somewhat faster than the above version.-integerSquareRootNewton' :: Integer -> (Integer,Integer)-integerSquareRootNewton' n-  | n<0 = error "integerSquareRootNewton: negative input"-  | n<2 = (n,0)-  | otherwise = go (div n 2) -  where-    go a = -      if m < a-        then go a' -        else (a, r + a*(m-a))-      where-        (m,r) = divMod n a-        a' = div (m + a) 2--{---- brute force test of integer square root-isqrt_test n1 n2 = -  [ k -  | k<-[n1..n2] -  , let (a,r) = integerSquareRoot' k-  , (a*a+r/=k) || (a*a>k) || (a+1)*(a+1)<=k -  ] -}  --------------------------------------------------------------------------------

@@ -0,0 +1,198 @@++-- | Some important number sequences. +--  +-- See the \"On-Line Encyclopedia of Integer Sequences\",+-- <https://oeis.org> .++module Math.Combinat.Numbers.Sequences where++--------------------------------------------------------------------------------++import Data.Array++import Math.Combinat.Helper ( sum' )+import Math.Combinat.Sign++--------------------------------------------------------------------------------+-- * Factorial++-- | A000142.+factorial :: Integral a => a -> Integer+factorial n+  | n <  0    = error "factorial: input should be nonnegative"+  | n == 0    = 1+  | otherwise = product [1..fromIntegral n]++-- | A006882.+doubleFactorial :: Integral a => a -> Integer+doubleFactorial n+  | n <  0    = error "doubleFactorial: input should be nonnegative"+  | n == 0    = 1+  | odd n     = product [1,3..fromIntegral n]+  | otherwise = product [2,4..fromIntegral n]++--------------------------------------------------------------------------------+-- * Binomial and multinomial++-- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below.+binomial :: Integral a => a -> a -> Integer+binomial n k +  | k > n = 0+  | k < 0 = 0+  | k > (n `div` 2) = binomial n (n-k)+  | otherwise = (product [n'-k'+1 .. n']) `div` (product [1..k'])+  where +    k' = fromIntegral k+    n' = fromIntegral n++-- | The extension of the binomial function to negative inputs. This should satisfy the following properties:+--+-- > for n,k >=0 : signedBinomial n k == binomial n k+-- > for any n,k : signedBinomial n k == signedBinomial n (n-k) +-- > for k >= 0  : signedBinomial (-n) k == (-1)^k * signedBinomial (n+k-1) k+--+-- Note: This is compatible with Mathematica's @Binomial@ function.+--+signedBinomial :: Int -> Int -> Integer+signedBinomial n k+  | n >= 0     = binomial n k+  | k >= 0     = negateIfOdd    k  $ binomial (k-n-1)   k  +  | otherwise  = negateIfOdd (n+k) $ binomial (-k-1) (-n-1)++{-+test_signed_0 = [ signedBinomial ( n) k == signedBinomial ( n) ( n-k)                | n<-[-30..40] , k<-[-30..40] ]+test_signed_1 = [ signedBinomial (-n) k == signedBinomial (-n) (-n-k)                | n<-[-30..40] , k<-[-30..40] ]+test_signed_2 = [ signedBinomial (-n) k == negateIfOdd k $ signedBinomial (n+k-1) k  | n<-[-30..40] , k<-[0..30] ]+-}++-- | A given row of the Pascal triangle; equivalent to a sequence of binomial +-- numbers, but much more efficient. You can also left-fold over it.+--+-- > pascalRow n == [ binomial n k | k<-[0..n] ]+pascalRow :: Integral a => a -> [Integer]+pascalRow n' = worker 0 1 where+  n = fromIntegral n'+  worker j x+    | j>n   = [] +    | True  = let j'=j+1 in x : worker j' (div (x*(n-j)) j') ++multinomial :: Integral a => [a] -> Integer+multinomial xs = div+  (factorial (sum xs))+  (product [ factorial x | x<-xs ])  +  +--------------------------------------------------------------------------------+-- * Catalan numbers++-- | Catalan numbers. OEIS:A000108.+catalan :: Integral a => a -> Integer+catalan n +  | n < 0     = 0+  | otherwise = binomial (n+n) n `div` fromIntegral (n+1)++-- | Catalan's triangle. OEIS:A009766.+-- Note:+--+-- > catalanTriangle n n == catalan n+-- > catalanTriangle n k == countStandardYoungTableaux (toPartition [n,k])+--+catalanTriangle :: Integral a => a -> a -> Integer+catalanTriangle n k+  | k > n     = 0+  | k < 0     = 0+  | otherwise = (binomial (n+k) n * fromIntegral (n-k+1)) `div` fromIntegral (n+1)++--------------------------------------------------------------------------------+-- * Stirling numbers++-- | Rows of (signed) Stirling numbers of the first kind. OEIS:A008275.+-- Coefficients of the polinomial @(x-1)*(x-2)*...*(x-n+1)@.+-- This function uses the recursion formula.+signedStirling1stArray :: Integral a => a -> Array Int Integer+signedStirling1stArray n+  | n <  1    = error "stirling1stArray: n should be at least 1"+  | n == 1    = listArray (1,1 ) [1]+  | otherwise = listArray (1,n') [ lkp (k-1) - fromIntegral (n-1) * lkp k | k<-[1..n'] ] +  where+    prev = signedStirling1stArray (n-1)+    n' = fromIntegral n :: Int+    lkp j | j <  1    = 0+          | j >= n'   = 0+          | otherwise = prev ! j +        +-- | (Signed) Stirling numbers of the first kind. OEIS:A008275.+-- This function uses 'signedStirling1stArray', so it shouldn't be used+-- to compute /many/ Stirling numbers.+--+-- Argument order: @signedStirling1st n k@+--+signedStirling1st :: Integral a => a -> a -> Integer+signedStirling1st n k +  | k==0 && n==0 = 1+  | k < 1        = 0+  | k > n        = 0+  | otherwise    = signedStirling1stArray n ! (fromIntegral k)++-- | (Unsigned) Stirling numbers of the first kind. See 'signedStirling1st'.+unsignedStirling1st :: Integral a => a -> a -> Integer+unsignedStirling1st n k = abs (signedStirling1st n k)++-- | Stirling numbers of the second kind. OEIS:A008277.+-- This function uses an explicit formula.+-- +-- Argument order: @stirling2nd n k@+--+stirling2nd :: Integral a => a -> a -> Integer+stirling2nd n k +  | k==0 && n==0 = 1+  | k < 1        = 0+  | k > n        = 0+  | otherwise = sum xs `div` factorial k where+      xs = [ negateIfOdd (k-i) $ binomial k i * (fromIntegral i)^n | i<-[0..k] ]++--------------------------------------------------------------------------------+-- * Bernoulli numbers++-- | Bernoulli numbers. @bernoulli 1 == -1%2@ and @bernoulli k == 0@ for+-- k>2 and /odd/. This function uses the formula involving Stirling numbers+-- of the second kind. Numerators: A027641, denominators: A027642.+bernoulli :: Integral a => a -> Rational+bernoulli n +  | n <  0    = error "bernoulli: n should be nonnegative"+  | n == 0    = 1+  | n == 1    = -1/2+  | otherwise = sum [ f k | k<-[1..n] ] +  where+    f k = toRational (negateIfOdd (n+k) $ factorial k * stirling2nd n k) +        / toRational (k+1)++--------------------------------------------------------------------------------+-- * Bell numbers++-- | Bell numbers (Sloane's A000110) from B(0) up to B(n). B(0)=B(1)=1, B(2)=2, etc. +--+-- The Bell numbers count the number of /set partitions/ of a set of size @n@+-- +-- See <http://en.wikipedia.org/wiki/Bell_number>+--+bellNumbersArray :: Integral a => a -> Array Int Integer+bellNumbersArray nn = arr where+  arr = array (0::Int,n) kvs +  n = fromIntegral nn :: Int+  kvs = (0,1) : [ (k, f k) | k<-[1..n] ] +  f n = sum' [ binomial (n-1) k * arr ! k | k<-[0..n-1] ]++-- | The n-th Bell number B(n), using the Stirling numbers of the second kind.+-- This may be slower than using 'bellNumbersArray'.+bellNumber :: Integral a => a -> Integer+bellNumber nn+  | n <  0     = error "bellNumber: expecting a nonnegative index"+  | n == 0     = 1+  | otherwise  = sum' [ stirling2nd n k | k<-[1..n] ] +  where+    n = fromIntegral nn :: Int++--------------------------------------------------------------------------------+++ 

@@ -8,7 +8,7 @@ -- TODO: better names for these functions. -- -{-# LANGUAGE CPP, GeneralizedNewtypeDeriving #-}+{-# LANGUAGE CPP, BangPatterns, GeneralizedNewtypeDeriving #-} module Math.Combinat.Numbers.Series where  --------------------------------------------------------------------------------@@ -63,9 +63,17 @@ scaleSeries :: Num a => a -> [a] -> [a] scaleSeries s = map (*s) +-- | A different implementation, taken from:+--+-- M. Douglas McIlroy: Power Series, Power Serious  mulSeries :: Num a => [a] -> [a] -> [a]-mulSeries = convolve+mulSeries xs ys = go (xs ++ repeat 0) (ys ++ repeat 0) where+  go (f:fs) ggs@(g:gs) = f*g : (scaleSeries f gs) `addSeries` go fs ggs +-- | Multiplication of power series. This implementation is a synonym for 'convolve'+mulSeriesNaive :: Num a => [a] -> [a] -> [a]+mulSeriesNaive = convolve+ productOfSeries :: Num a => [[a]] -> [a] productOfSeries = convolveMany @@ -90,6 +98,14 @@ -------------------------------------------------------------------------------- -- * Reciprocals of general power series +-- | Division of series.+--+-- Taken from: M. Douglas McIlroy: Power Series, Power Serious +divSeries :: (Eq a, Fractional a) => [a] -> [a] -> [a]+divSeries xs ys = go (xs ++ repeat 0) (ys ++ repeat 0) where+  go (0:fs)     (0:gs) = go fs gs+  go (f:fs) ggs@(g:gs) = let q = f/g in q : go (fs `subSeries` scaleSeries q gs) ggs+ -- | Given a power series, we iteratively compute its multiplicative inverse reciprocalSeries :: (Eq a, Fractional a) => [a] -> [a] reciprocalSeries series = case series of@@ -119,9 +135,13 @@ -- | @g \`composeSeries\` f@ is the power series expansion of @g(f(x))@. -- This is a synonym for @flip substitute@. ----- We require that the constant term of @f@ is zero.+-- This implementation is taken from+--+-- M. Douglas McIlroy: Power Series, Power Serious  composeSeries :: (Eq a, Num a) => [a] -> [a] -> [a]-composeSeries g f = substitute f g+composeSeries xs ys = go (xs ++ repeat 0) (ys ++ repeat 0) where+  go (f:fs) (0:gs) = f : mulSeries gs (go fs (0:gs))+  go (f:fs) (_:gs) = error "PowerSeries/composeSeries: we expect the the constant term of the inner series to be zero"  -- | @substitute f g@ is the power series corresponding to @g(f(x))@.  -- Equivalently, this is the composition of univariate functions (in the \"wrong\" order).@@ -129,10 +149,18 @@ -- Note: for this to be meaningful in general (not depending on convergence properties), -- we need that the constant term of @f@ is zero. substitute :: (Eq a, Num a) => [a] -> [a] -> [a]-substitute as_ bs_ = +substitute f g = composeSeries g f++-- | Naive implementation of 'composeSeries' (via 'substituteNaive')+composeSeriesNaive :: (Eq a, Num a) => [a] -> [a] -> [a]+composeSeriesNaive g f = substituteNaive f g++-- | Naive implementation of 'substitute'+substituteNaive :: (Eq a, Num a) => [a] -> [a] -> [a]+substituteNaive as_ bs_ =    case head as of     0 -> [ f n | n<-[0..] ]-    _ -> error "PowerSeries/substitute: we expect the the constant term of the inner series to be zero"+    _ -> error "PowerSeries/substituteNaive: we expect the the constant term of the inner series to be zero"   where     as = as_ ++ repeat 0     bs = bs_ ++ repeat 0@@ -148,6 +176,19 @@ -------------------------------------------------------------------------------- -- * Lagrange inversions +-- | We expect the input series to match @(0:a1:_)@. with a1 nonzero The following is true for the result (at least with exact arithmetic):+--+-- > substitute f (lagrangeInversion f) == (0 : 1 : repeat 0)+-- > substitute (lagrangeInversion f) f == (0 : 1 : repeat 0)+--+-- This implementation is taken from:+--+-- M. Douglas McIlroy: Power Series, Power Serious +lagrangeInversion :: (Eq a, Fractional a) => [a] -> [a]+lagrangeInversion xs = go (xs ++ repeat 0) where+  go (0:fs) = rs where rs = 0 : divSeries unitSeries (composeSeries fs rs)+  go (_:fs) = error "lagrangeInversion: the series should start with (0 + a1*x + a2*x^2 + ...) where a1 is non-zero"+ -- | Coefficients of the Lagrange inversion lagrangeCoeff :: Partition -> Integer lagrangeCoeff p = div numer denom where@@ -162,11 +203,11 @@ -- > substitute f (integralLagrangeInversion f) == (0 : 1 : repeat 0) -- > substitute (integralLagrangeInversion f) f == (0 : 1 : repeat 0) ---integralLagrangeInversion :: (Eq a, Num a) => [a] -> [a]-integralLagrangeInversion series_ = +integralLagrangeInversionNaive :: (Eq a, Num a) => [a] -> [a]+integralLagrangeInversionNaive series_ =    case series of     (0:1:rest) -> 0 : 1 : [ f n | n<-[1..] ]-    _ -> error "integralLagrangeInversion: the series should start with (0 + x + a2*x^2 + ...)"+    _ -> error "integralLagrangeInversionNaive: the series should start with (0 + x + a2*x^2 + ...)"   where     series = series_ ++ repeat 0     as  = tail series @@ -175,20 +216,16 @@               | p <- partitions n               ]  --- | We expect the input series to match @(0:a1:_)@. with a1 nonzero The following is true for the result (at least with exact arithmetic):------ > substitute f (lagrangeInversion f) == (0 : 1 : repeat 0)--- > substitute (lagrangeInversion f) f == (0 : 1 : repeat 0)----lagrangeInversion :: (Eq a, Fractional a) => [a] -> [a]-lagrangeInversion series_ = +-- | Naive implementation of 'lagrangeInversion'+lagrangeInversionNaive :: (Eq a, Fractional a) => [a] -> [a]+lagrangeInversionNaive series_ =    case series of     (0:a1:rest) -> if a1 ==0        then err        else 0 : (1/a1) : [ f n / a1^(n+1) | n<-[1..] ]     _ -> err   where-    err    = error "lagrangeInversion: the series should start with (0 + a1*x + a2*x^2 + ...) where a1 is non-zero"+    err    = error "lagrangeInversionNaive: the series should start with (0 + a1*x + a2*x^2 + ...) where a1 is non-zero"     series = series_ ++ repeat 0     a1  = series !! 1     as  = map (/a1) (tail series)@@ -196,8 +233,22 @@     f n = sum [ fromInteger (lagrangeCoeff p) * product [ (a i)^j | (i,j) <- toExponentialForm p ]               | p <- partitions n               ] -  ++ --------------------------------------------------------------------------------+-- * Differentiation and integration++differentiateSeries :: Num a => [a] -> [a]+differentiateSeries (y:ys) = go (1::Int) ys where+  go !n (x:xs) = fromIntegral n * x : go (n+1) xs+  go _  []     = []++integrateSeries :: Fractional a => [a] -> [a]+integrateSeries ys = 0 : go (1::Int) ys where+  go !n (x:xs) = x / (fromIntegral n) : go (n+1) xs+  go _  []     = []++-------------------------------------------------------------------------------- -- * Power series expansions of elementary functions  -- | Power series expansion of @exp(x)@@@ -214,6 +265,13 @@ sinSeries :: Fractional a => [a] sinSeries = go 1 1 where   go i e = 0 : e : go (i+2) (-e / ((i+1)*(i+2)))++-- | Alternative implementation using differential equations.+--+-- Taken from: M. Douglas McIlroy: Power Series, Power Serious+cosSeries2, sinSeries2 :: Fractional a => [a]+cosSeries2 = unitSeries `subSeries` integrateSeries sinSeries2+sinSeries2 =                        integrateSeries cosSeries2  -- | Power series expansion of @cosh(x)@ coshSeries :: Fractional a => [a]

@@ -18,7 +18,54 @@ --   {-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables #-}-module Math.Combinat.Partitions.Integer where+module Math.Combinat.Partitions.Integer +  ( -- module Math.Combinat.Partitions.Integer.Count+    module Math.Combinat.Partitions.Integer.Naive+    -- * Types and basic stuff+  , Partition+    -- * Conversion to\/from lists+  , fromPartition +  , mkPartition +  , toPartition +  , toPartitionUnsafe +  , isPartition +    -- * Union and sum+  , unionOfPartitions+  , sumOfPartitions+    -- * Generating partitions+  , partitions +  , partitions'+  , allPartitions +  , allPartitionsGrouped +  , allPartitions'  +  , allPartitionsGrouped'  +    -- * Counting partitions+  , countPartitions+  , countPartitions'+  , countAllPartitions+  , countAllPartitions'+  , countPartitionsWithKParts +    -- * Random partitions+  , randomPartition+  , randomPartitions+    -- * Dominating \/ dominated partitions+  , dominatedPartitions +  , dominatingPartitions +    -- * Partitions with given number of parts+  , partitionsWithKParts+    -- * Partitions with only odd\/distinct parts+  , partitionsWithOddParts +  , partitionsWithDistinctParts+    -- * Sub- and super-partitions of a given partition+  , subPartitions +  , allSubPartitions +  , superPartitions +    -- * ASCII Ferrers diagrams+  , PartitionConvention(..)+  , asciiFerrersDiagram +  , asciiFerrersDiagram'+  )+  where  -------------------------------------------------------------------------------- @@ -36,31 +83,29 @@ import Data.Array import System.Random ------------------------------------------------------------------------------------ * Type and basic stuff---- | A partition of an integer. The additional invariant enforced here is that partitions --- are monotone decreasing sequences of /positive/ integers. The @Ord@ instance is lexicographical.-newtype Partition = Partition [Int] deriving (Eq,Ord,Show,Read)--instance HasNumberOfParts Partition where-  numberOfParts (Partition p) = length p+import Math.Combinat.Partitions.Integer.Naive+import Math.Combinat.Partitions.Integer.IntList+import Math.Combinat.Partitions.Integer.Count  ---------------------------------------------------------------------------------+-- * Conversion to\/from lists++fromPartition :: Partition -> [Int]+fromPartition (Partition_ part) = part    -- | Sorts the input, and cuts the nonpositive elements. mkPartition :: [Int] -> Partition-mkPartition xs = Partition $ sortBy (reverseCompare) $ filter (>0) xs---- | Assumes that the input is decreasing.-toPartitionUnsafe :: [Int] -> Partition-toPartitionUnsafe = Partition+mkPartition xs = toPartitionUnsafe $ sortBy (reverseCompare) $ filter (>0) xs  -- | Checks whether the input is an integer partition. See the note at 'isPartition'! toPartition :: [Int] -> Partition toPartition xs = if isPartition xs   then toPartitionUnsafe xs   else error "toPartition: not a partition"++-- | Assumes that the input is decreasing.+toPartitionUnsafe :: [Int] -> Partition+toPartitionUnsafe = Partition_    -- | This returns @True@ if the input is non-increasing sequence of  -- /positive/ integers (possibly empty); @False@ otherwise.@@ -70,211 +115,43 @@ isPartition [x] = x > 0 isPartition (x:xs@(y:_)) = (x >= y) && isPartition xs -isEmptyPartition :: Partition -> Bool-isEmptyPartition (Partition p) = null p--emptyPartition :: Partition-emptyPartition = Partition []--instance CanBeEmpty Partition where-  empty   = emptyPartition-  isEmpty = isEmptyPartition--fromPartition :: Partition -> [Int]-fromPartition (Partition part) = part---- | The first element of the sequence.-partitionHeight :: Partition -> Int-partitionHeight (Partition part) = case part of-  (p:_) -> p-  []    -> 0-  --- | The length of the sequence (that is, the number of parts).-partitionWidth :: Partition -> Int-partitionWidth (Partition part) = length part--instance HasHeight Partition where-  height = partitionHeight- -instance HasWidth Partition where-  width = partitionWidth--heightWidth :: Partition -> (Int,Int)-heightWidth part = (height part, width part)---- | The weight of the partition ---   (that is, the sum of the corresponding sequence).-partitionWeight :: Partition -> Int-partitionWeight (Partition part) = sum' part--instance HasWeight Partition where -  weight = partitionWeight---- | The dual (or conjugate) partition.-dualPartition :: Partition -> Partition-dualPartition (Partition part) = Partition (_dualPartition part)--instance HasDuality Partition where -  dual = dualPartition--data Pair = Pair !Int !Int--_dualPartition :: [Int] -> [Int]-_dualPartition [] = []-_dualPartition xs = go 0 (diffSequence xs) [] where-  go !i (d:ds) acc = go (i+1) ds (d:acc)-  go n  []     acc = finish n acc -  finish !j (k:ks) = replicate k j ++ finish (j-1) ks-  finish _  []     = []--{---- more variations:--_dualPartition_b :: [Int] -> [Int]-_dualPartition_b [] = []-_dualPartition_b xs = go 1 (diffSequence xs) [] where-  go !i (d:ds) acc = go (i+1) ds ((d,i):acc)-  go _  []     acc = concatMap (\(d,i) -> replicate d i) acc--_dualPartition_c :: [Int] -> [Int]-_dualPartition_c [] = []-_dualPartition_c xs = reverse $ concat $ zipWith f [1..] (diffSequence xs) where-  f _ 0 = []-  f k d = replicate d k--}---- | A simpler, but bit slower (about twice?) implementation of dual partition-_dualPartitionNaive :: [Int] -> [Int]-_dualPartitionNaive [] = []-_dualPartitionNaive xs@(k:_) = [ length $ filter (>=i) xs | i <- [1..k] ]---- | From a sequence @[a1,a2,..,an]@ computes the sequence of differences--- @[a1-a2,a2-a3,...,an-0]@-diffSequence :: [Int] -> [Int]-diffSequence = go where-  go (x:ys@(y:_)) = (x-y) : go ys -  go [x] = [x]-  go []  = []+--------------------------------------------------------------------------------+-- * Union and sum --- | Example:+-- | This is simply the union of parts. For example  ----- > elements (toPartition [5,4,1]) ==--- >   [ (1,1), (1,2), (1,3), (1,4), (1,5)--- >   , (2,1), (2,2), (2,3), (2,4)--- >   , (3,1)--- >   ]+-- > Partition [4,2,1] `unionOfPartitions` Partition [4,3,1] == Partition [4,4,3,2,1,1] ---elements :: Partition -> [(Int,Int)]-elements (Partition part) = _elements part--_elements :: [Int] -> [(Int,Int)]-_elements shape = [ (i,j) | (i,l) <- zip [1..] shape, j<-[1..l] ] -------------------------------------------------------------------------------------- * Exponential form+-- Note: This is the dual of pointwise sum, 'sumOfPartitions'+--+unionOfPartitions :: Partition -> Partition -> Partition +unionOfPartitions (Partition_ xs) (Partition_ ys) = mkPartition (xs ++ ys) --- | We convert a partition to exponential form.--- @(i,e)@ mean @(i^e)@; for example @[(1,4),(2,3)]@ corresponds to @(1^4)(2^3) = [2,2,2,1,1,1,1]@. Another example:+-- | Pointwise sum of the parts. For example: ----- > toExponentialForm (Partition [5,5,3,2,2,2,2,1,1]) == [(1,2),(2,4),(3,1),(5,2)]+-- > Partition [3,2,1,1] `sumOfPartitions` Partition [4,3,1] == Partition [7,5,2,1] ---toExponentialForm :: Partition -> [(Int,Int)]-toExponentialForm = _toExponentialForm . fromPartition--_toExponentialForm :: [Int] -> [(Int,Int)]-_toExponentialForm = reverse . map (\xs -> (head xs,length xs)) . group--fromExponentialFrom :: [(Int,Int)] -> Partition-fromExponentialFrom = Partition . sortBy reverseCompare . go where-  go ((j,e):rest) = replicate e j ++ go rest-  go []           = []   -------------------------------------------------------------------------------------- * Automorphisms ---- | Computes the number of \"automorphisms\" of a given integer partition.-countAutomorphisms :: Partition -> Integer  -countAutomorphisms = _countAutomorphisms . fromPartition--_countAutomorphisms :: [Int] -> Integer-_countAutomorphisms = multinomial . map length . group+-- Note: This is the dual of 'unionOfPartitions'+--+sumOfPartitions :: Partition -> Partition -> Partition +sumOfPartitions (Partition_ xs) (Partition_ ys) = Partition_ (longZipWith 0 0 (+) xs ys) ----------------------------------------------------------------------------------+-------------------------------------------------------------------------------- -- * Generating partitions  -- | Partitions of @d@. partitions :: Int -> [Partition]-partitions = map Partition . _partitions---- | Partitions of @d@, as lists-_partitions :: Int -> [[Int]]-_partitions d = go d d where-  go _  0  = [[]]-  go !h !n = [ a:as | a<-[1..min n h], as <- go a (n-a) ]---- | Number of partitions of @n@-countPartitions :: Int -> Integer-countPartitions n = partitionCountList !! n---- | This uses 'countPartitions'', and thus is slow-countPartitionsNaive :: Int -> Integer-countPartitionsNaive d = countPartitions' (d,d) d-------------------------------------------------------------------------------------- | Infinite list of number of partitions of @0,1,2,...@------ This uses the infinite product formula the generating function of partitions, recursively--- expanding it; it is quite fast.------ > partitionCountList == map countPartitions [0..]----partitionCountList :: [Integer]-partitionCountList = final where--  final = go 1 (1:repeat 0) --  go !k (x:xs) = x : go (k+1) ys where-    ys = zipWith (+) xs (take k final ++ ys)-    -- explanation:-    --   xs == drop k $ f (k-1)-    --   ys == drop k $ f (k  )  --{---Full explanation of 'partitionCountList':-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~--let f k = productPSeries $ map (:[]) [1..k]--f 0 = [1,0,0,0,0,0,0,0...]-f 1 = [1,1,1,1,1,1,1,1...]-f 2 = [1,1,2,2,3,3,4,4...]-f 3 = [1,1,2,3,4,5,7,8...]--observe: --* take (k+1) (f k) == take (k+1) partitionCountList-* f (k+1) == zipWith (+) (f k) (replicate (k+1) 0 ++ f (k+1))--now apply (drop (k+1)) to the second one : --* drop (k+1) (f (k+1)) == zipWith (+) (drop (k+1) $ f k) (f (k+1))-* f (k+1) = take (k+1) final ++ drop (k+1) (f (k+1))+partitions = map toPartitionUnsafe . _partitions --}+-- | Partitions of d, fitting into a given rectangle. The order is again lexicographic.+partitions'  +  :: (Int,Int)     -- ^ (height,width)+  -> Int           -- ^ d+  -> [Partition]+partitions' hw d = map toPartitionUnsafe $ _partitions' hw d          -------------------------------------------------------------------------------- --- | Naive infinite list of number of partitions of @0,1,2,...@------ > partitionCountListNaive == map countPartitionsNaive [0..]------ This is much slower than the power series expansion above.----partitionCountListNaive :: [Integer]-partitionCountListNaive = map countPartitionsNaive [0..]- -- | All integer partitions up to a given degree (that is, all integer partitions whose sum is less or equal to @d@) allPartitions :: Int -> [Partition] allPartitions d = concat [ partitions i | i <- [0..d] ]@@ -296,48 +173,14 @@   -> [[Partition]] allPartitionsGrouped' (h,w) = [ partitions' (h,w) i | i <- [0..d] ] where d = h*w --- | # = \\binom { h+w } { h }-countAllPartitions' :: (Int,Int) -> Integer-countAllPartitions' (h,w) = -  binomial (h+w) (min h w)-  --sum [ countPartitions' (h,w) i | i <- [0..d] ] where d = h*w -countAllPartitions :: Int -> Integer-countAllPartitions d = sum' [ countPartitions i | i <- [0..d] ]---- | Integer partitions of @d@, fitting into a given rectangle, as lists.-_partitions' -  :: (Int,Int)     -- ^ (height,width)-  -> Int           -- ^ d-  -> [[Int]]        -_partitions' _ 0 = [[]] -_partitions' ( 0 , _) d = if d==0 then [[]] else []-_partitions' ( _ , 0) d = if d==0 then [[]] else []-_partitions' (!h ,!w) d = -  [ i:xs | i <- [1..min d h] , xs <- _partitions' (i,w-1) (d-i) ]---- | Partitions of d, fitting into a given rectangle. The order is again lexicographic.-partitions'  -  :: (Int,Int)     -- ^ (height,width)-  -> Int           -- ^ d-  -> [Partition]-partitions' hw d = map toPartitionUnsafe $ _partitions' hw d        --countPartitions' :: (Int,Int) -> Int -> Integer-countPartitions' _ 0 = 1-countPartitions' (0,_) d = if d==0 then 1 else 0-countPartitions' (_,0) d = if d==0 then 1 else 0-countPartitions' (h,w) d = sum-  [ countPartitions' (i,w-1) (d-i) | i <- [1..min d h] ] -- --------------------------------------------------------------------------------- -- * Random partitions  -- | Uniformly random partition of the given weight.  -- -- NOTE: This algorithm is effective for small @n@-s (say @n@ up to a few hundred \/ one thousand it should work nicely),--- and the first time it is executed may be slower (as it needs to build the table 'partitionCountList' first)+-- and the first time it is executed may be slower (as it needs to build the table of partitions counts first) -- -- Algorithm of Nijenhuis and Wilf (1975); see --@@ -359,8 +202,7 @@   -> g -> ([Partition], g) randomPartitions howmany n = runRand $ replicateM howmany (worker n []) where -  table = listArray (0,n) $ take (n+1) partitionCountList :: Array Int Integer-  cnt k = table ! k+  cnt = countPartitions     finish :: [(Int,Int)] -> Partition   finish = mkPartition . concatMap f where f (j,d) = replicate j d@@ -387,21 +229,8 @@     let jd@(!j,!d) = find_jd m capm     worker (m - j*d) (jd:acc) -------------------------------------------------------------------------------------- * Dominance order ---- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions--- (this is partial ordering on the set of partitions of @n@).------ See <http://en.wikipedia.org/wiki/Dominance_order>----dominates :: Partition -> Partition -> Bool-dominates (Partition qs) (Partition ps) -  = and $ zipWith (>=) (sums (qs ++ repeat 0)) (sums ps)-  where-    sums = scanl (+) 0-+--------------------------------------------------------------------------------+-- * Dominating \/ dominated partitions  -- | Lists all partitions of the same weight as @lambda@ and also dominated by @lambda@ -- (that is, all partial sums are less or equal):@@ -409,21 +238,7 @@ -- > dominatedPartitions lam == [ mu | mu <- partitions (weight lam), lam `dominates` mu ] --  dominatedPartitions :: Partition -> [Partition]    -dominatedPartitions (Partition lambda) = map Partition (_dominatedPartitions lambda)--_dominatedPartitions :: [Int] -> [[Int]]-_dominatedPartitions []     = [[]]-_dominatedPartitions lambda = go (head lambda) w dsums 0 where--  n = length lambda-  w = sum    lambda-  dsums = scanl1 (+) (lambda ++ repeat 0)--  go _   0 _       _  = [[]]-  go !h !w (!d:ds) !e  -    | w >  0  = [ (a:as) | a <- [1..min h (d-e)] , as <- go a (w-a) ds (e+a) ] -    | w == 0  = [[]]-    | w <  0  = error "_dominatedPartitions: fatal error; shouldn't happen"+dominatedPartitions (Partition_ lambda) = map Partition_ (_dominatedPartitions lambda)  -- | Lists all partitions of the sime weight as @mu@ and also dominating @mu@ -- (that is, all partial sums are greater or equal):@@ -431,21 +246,7 @@ -- > dominatingPartitions mu == [ lam | lam <- partitions (weight mu), lam `dominates` mu ] --  dominatingPartitions :: Partition -> [Partition]    -dominatingPartitions (Partition mu) = map Partition (_dominatingPartitions mu)--_dominatingPartitions :: [Int] -> [[Int]]-_dominatingPartitions []     = [[]]-_dominatingPartitions mu     = go w w dsums 0 where--  n = length mu-  w = sum    mu-  dsums = scanl1 (+) (mu ++ repeat 0)--  go _   0 _       _  = [[]]-  go !h !w (!d:ds) !e  -    | w >  0  = [ (a:as) | a <- [max 0 (d-e)..min h w] , as <- go a (w-a) ds (e+a) ] -    | w == 0  = [[]]-    | w <  0  = error "_dominatingPartitions: fatal error; shouldn't happen"+dominatingPartitions (Partition_ mu) = map Partition_ (_dominatingPartitions mu)  -------------------------------------------------------------------------------- -- * Partitions with given number of parts@@ -460,7 +261,7 @@   :: Int    -- ^ @k@ = number of parts   -> Int    -- ^ @n@ = the integer we partition   -> [Partition]-partitionsWithKParts k n = map Partition $ go n k n where+partitionsWithKParts k n = map Partition_ $ go n k n where {-   h = max height   k = number of parts@@ -472,23 +273,12 @@     | k == 1     = if h>=n && n>=1 then [[n]] else []     | otherwise  = [ a:p | a <- [1..(min h (n-k+1))] , p <- go a (k-1) (n-a) ] -countPartitionsWithKParts -  :: Int    -- ^ @k@ = number of parts-  -> Int    -- ^ @n@ = the integer we partition-  -> Integer-countPartitionsWithKParts k n = go n k n where-  go !h !k !n -    | k <  0     = 0-    | k == 0     = if h>=0 && n==0 then 1 else 0-    | k == 1     = if h>=n && n>=1 then 1 else 0-    | otherwise  = sum' [ go a (k-1) (n-a) | a<-[1..(min h (n-k+1))] ]- -------------------------------------------------------------------------------- -- * Partitions with only odd\/distinct parts  -- | Partitions of @n@ with only odd parts partitionsWithOddParts :: Int -> [Partition]-partitionsWithOddParts d = map Partition (go d d) where+partitionsWithOddParts d = map Partition_ (go d d) where   go _  0  = [[]]   go !h !n = [ a:as | a<-[1,3..min n h], as <- go a (n-a) ] @@ -510,143 +300,31 @@ -- > length (partitionsWithDistinctParts d) == length (partitionsWithOddParts d) -- partitionsWithDistinctParts :: Int -> [Partition]-partitionsWithDistinctParts d = map Partition (go d d) where+partitionsWithDistinctParts d = map Partition_ (go d d) where   go _  0  = [[]]   go !h !n = [ a:as | a<-[1..min n h], as <- go (a-1) (n-a) ]  -------------------------------------------------------------------------------- -- * Sub- and super-partitions of a given partition --- | Returns @True@ of the first partition is a subpartition (that is, fit inside) of the second.--- This includes equality-isSubPartitionOf :: Partition -> Partition -> Bool-isSubPartitionOf (Partition ps) (Partition qs) = and $ zipWith (<=) ps (qs ++ repeat 0)---- | This is provided for convenience\/completeness only, as:------ > isSuperPartitionOf q p == isSubPartitionOf p q----isSuperPartitionOf :: Partition -> Partition -> Bool-isSuperPartitionOf (Partition qs) (Partition ps) = and $ zipWith (<=) ps (qs ++ repeat 0)-- -- | Sub-partitions of a given partition with the given weight: -- -- > sort (subPartitions d q) == sort [ p | p <- partitions d, isSubPartitionOf p q ] -- subPartitions :: Int -> Partition -> [Partition]-subPartitions d (Partition ps) = map Partition (_subPartitions d ps)--_subPartitions :: Int -> [Int] -> [[Int]]-_subPartitions d big-  | null big       = if d==0 then [[]] else []-  | d > sum' big   = []-  | d < 0          = []-  | otherwise      = go d (head big) big-  where-    go :: Int -> Int -> [Int] -> [[Int]]-    go !k !h []      = if k==0 then [[]] else []-    go !k !h (b:bs) -      | k<0 || h<0   = []-      | k==0         = [[]]-      | h==0         = []-      | otherwise    = [ this:rest | this <- [1..min h b] , rest <- go (k-this) this bs ]------------------------------------------+subPartitions d (Partition_ ps) = map Partition_ (_subPartitions d ps)  -- | All sub-partitions of a given partition allSubPartitions :: Partition -> [Partition]-allSubPartitions (Partition ps) = map Partition (_allSubPartitions ps)--_allSubPartitions :: [Int] -> [[Int]]-_allSubPartitions big -  | null big   = [[]]-  | otherwise  = go (head big) big-  where-    go _  [] = [[]]-    go !h (b:bs) -      | h==0         = []-      | otherwise    = [] : [ this:rest | this <- [1..min h b] , rest <- go this bs ]------------------------------------------+allSubPartitions (Partition_ ps) = map Partition_ (_allSubPartitions ps)  -- | Super-partitions of a given partition with the given weight: -- -- > sort (superPartitions d p) == sort [ q | q <- partitions d, isSubPartitionOf p q ] -- superPartitions :: Int -> Partition -> [Partition]-superPartitions d (Partition ps) = map Partition (_superPartitions d ps)--_superPartitions :: Int -> [Int] -> [[Int]]-_superPartitions dd small-  | dd < w0     = []-  | null small  = _partitions dd-  | otherwise   = go dd w1 dd (small ++ repeat 0)-  where-    w0 = sum' small-    w1 = w0 - head small-    -- d = remaining weight of the outer partition we are constructing-    -- w = remaining weight of the inner partition (we need to reserve at least this amount)-    -- h = max height (decreasing)-    go !d !w !h (!a:as@(b:_)) -      | d < 0     = []-      | d == 0    = if a == 0 then [[]] else []-      | otherwise = [ this:rest | this <- [max 1 a .. min h (d-w)] , rest <- go (d-this) (w-b) this as ]+superPartitions d (Partition_ ps) = map toPartitionUnsafe (_superPartitions d ps)     ------------------------------------------------------------------------------------ * The Pieri rule---- | The Pieri rule computes @s[lambda]*h[n]@ as a sum of @s[mu]@-s (each with coefficient 1).------ See for example <http://en.wikipedia.org/wiki/Pieri's_formula>----pieriRule :: Partition -> Int -> [Partition] -pieriRule (Partition lambda) n = map Partition (_pieriRule lambda n) where--  -- | We assume here that @lambda@ is a partition (non-increasing sequence of /positive/ integers)! -  _pieriRule :: [Int] -> Int -> [[Int]] -  _pieriRule lambda n-    | n == 0     = [lambda]-    | n <  0     = [] -    | otherwise  = go n diffs dsums (lambda++[0]) -    where-      diffs = n : diffSequence lambda                 -- maximum we can add to a given row-      dsums = reverse $ scanl1 (+) (reverse diffs)    -- partial sums of remaining total we can add-      go !k (d:ds) (p:ps@(q:_)) (l:ls) -        | k > p     = []-        | otherwise = [ h:tl | a <- [ max 0 (k-q) .. min d k ] , let h = l+a , tl <- go (k-a) ds ps ls ]-      go !k [d]    _      [l]    = if k <= d -                                     then if l+k>0 then [[l+k]] else [[]]-                                     else []-      go !k []     _      _      = if k==0 then [[]] else []---- | The dual Pieri rule computes @s[lambda]*e[n]@ as a sum of @s[mu]@-s (each with coefficient 1)-dualPieriRule :: Partition -> Int -> [Partition] -dualPieriRule lam n = map dualPartition $ pieriRule (dualPartition lam) n---{- --- moved to "Math.Combinat.Tableaux.GelfandTsetlin"---- | Computes the Schur expansion of @h[n1]*h[n2]*h[n3]*...*h[nk]@ via iterating the Pieri rule-iteratedPieriRule :: Num coeff => [Int] -> Map Partition coeff-iteratedPieriRule = iteratedPieriRule' (Partition [])---- | Iterating the Pieri rule, we can compute the Schur expansion of--- @h[lambda]*h[n1]*h[n2]*h[n3]*...*h[nk]@-iteratedPieriRule' :: Num coeff => Partition -> [Int] -> Map Partition coeff-iteratedPieriRule' plambda ns = iteratedPieriRule'' (plambda,1) ns--{-# SPECIALIZE iteratedPieriRule'' :: (Partition,Int    ) -> [Int] -> Map Partition Int     #-}-{-# SPECIALIZE iteratedPieriRule'' :: (Partition,Integer) -> [Int] -> Map Partition Integer #-}-iteratedPieriRule'' :: Num coeff => (Partition,coeff) -> [Int] -> Map Partition coeff-iteratedPieriRule'' (plambda,coeff0) ns = worker (Map.singleton plambda coeff0) ns where-  worker old []     = old-  worker old (n:ns) = worker new ns where-    stuff = [ (coeff, pieriRule lam n) | (lam,coeff) <- Map.toList old ] -    new   = foldl' f Map.empty stuff -    f t0 (c,ps) = foldl' (\t p -> Map.insertWith (+) p c t) t0 ps  --}  -------------------------------------------------------------------------------- -- * ASCII Ferrers diagrams

@@ -0,0 +1,819 @@++{- | Compact representation of integer partitions.++Partitions are conceptually nonincreasing sequences of /positive/ integers.++When the partition fits into a 15x15 rectangle, we encode the parts as nibbles in a single 64-bit word.+The most significant nibble is the first element, and the least significant nibble is used to encode+the length. This way equality and comparison of 64-bit words is the same as the corresponding operations+for partitions (lexicographic ordering).++This will make working with small partitions much more memory efficient (very helpful when+building tables indexed by partitions, for example!) and hopefully quite a bit faster, too.++When they do not fit into a 15x15 rectangle, but fit into 255x7, 255x15, 255x23 or 255x31, respectively,+then we extend the above to use the bytes of 1, 2, 3 or 4 64-bit words.++In the general case, we encode the partition as a list of 64-bit words, each encoding 4 16-bit parts.++Partitions with elements bigger than 65535 are not supported.++Note: This is an internal module, you are not supposed to import it directly.+-}++{-# LANGUAGE BangPatterns, PatternSynonyms, ViewPatterns, ForeignFunctionInterface #-}+module Math.Combinat.Partitions.Integer.Compact where++--------------------------------------------------------------------------------++import Data.Bits+import Data.Word+import Data.Ord+import Data.List ( intercalate , group , sort , sortBy , foldl' , scanl' ) ++import Math.Combinat.Compositions ( compositions' )+++--------------------------------------------------------------------------------+-- * The compact partition data type++data Partition+  = Nibble   {-# UNPACK #-} !Word64+  | Medium1  {-# UNPACK #-} !Word64+  | Medium2  {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64+  | Medium3  {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64+  | Medium4  {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64+  | WordList {-# UNPACK #-} !Int ![Word64]+  deriving (Eq,Show)++--------------------------------------------------------------------------------+  +-- | for debugging+partitionPrefixChar :: Partition -> Char+partitionPrefixChar p = case p of+  Nibble   {} -> 'N'+  Medium1  {} -> '1'+  Medium2  {} -> '2' +  Medium3  {} -> '3' +  Medium4  {} -> '4' +  WordList {} -> 'L'++{- +instance Show Partition where+  show compact = partitionPrefixChar compact +               : '<' : intercalate "," (map show $ toList compact) ++ ">"+-}++instance Ord Partition where+  compare = cmp+               +--------------------------------------------------------------------------------+-- * Pattern synonyms ++-- | Pattern sysnonyms allows us to use existing code with minimal modifications+pattern Nil :: Partition+pattern Nil <- (isEmpty -> True) where+        Nil =  empty++pattern Cons :: Int -> Partition -> Partition+pattern Cons x xs <- (uncons -> Just (x,xs)) where+        Cons x xs = cons x xs++-- | Simulated newtype constructor +pattern Partition_ :: [Int] -> Partition+pattern Partition_ xs <- (toList -> xs) where+        Partition_ xs = fromDescList xs++pattern Head :: Int -> Partition +pattern Head h <- (height -> h)++pattern Tail :: Partition -> Partition+pattern Tail xs <- (partitionTail -> xs)++pattern Length :: Int -> Partition +pattern Length n <- (width -> n)        ++--------------------------------------------------------------------------------+-- * Lexicographic comparison++-- | The lexicographic ordering+cmp :: Partition -> Partition -> Ordering+cmp (Nibble  a)           (Nibble  b)           = compare  a             b+cmp (Medium1 a1)          (Medium1 b1)          = compare  a1            b1+cmp (Medium2 a1 a2)       (Medium2 b1 b2)       = compare (a1,a2)       (b1,b2)+cmp (Medium3 a1 a2 a3)    (Medium3 b1 b2 b3)    = compare (a1,a2,a3)    (b1,b2,b3)+cmp (Medium4 a1 a2 a3 a4) (Medium4 b1 b2 b3 b4) = compare (a1,a2,a3,a4) (b1,b2,b3,b4)+cmp (WordList _ as)       (WordList _ bs)       = compare  as            bs+cmp p                     q                     = compare (toList p)    (toList q)+  +--------------------------------------------------------------------------------+-- * Basic (de)constructrion++empty :: Partition+empty = Nibble 0++isEmpty :: Partition -> Bool+isEmpty compact = case compact of+  Nibble x -> (x == 0)+  _        -> False++--------------------------------------------------------------------------------++singleton :: Int -> Partition+singleton x+  | x == 0      = Nibble 0+  | x <= 15     = Nibble     $ shiftL (i2w x) 60 + 1+  | x <= 255    = Medium1    $ shiftL (i2w x) 56 + 1+  | x <= 65535  = WordList 1 [ shiftL (i2w x) 48 ]+  | otherwise   = error "singleton: partitions with elements bigger than 65535 are not supported"++--------------------------------------------------------------------------------++uncons :: Partition -> Maybe (Int,Partition)+uncons compact = case compact of+  Nibble  0           -> Nothing+  Nibble  w           -> Just ( w2i (shiftR w  60) , Nibble $ shiftL (w .&. 0x0ffffffffffffff0) 4 + ((w .&. 15) - 1) )+  Medium1 w1          -> Just ( w2i (shiftR w1 56) , partitionTail compact )+  Medium2 w1 w2       -> Just ( w2i (shiftR w1 56) , partitionTail compact )+  Medium3 w1 w2 w3    -> Just ( w2i (shiftR w1 56) , partitionTail compact )+  Medium4 w1 w2 w3 w4 -> Just ( w2i (shiftR w1 56) , partitionTail compact )+  WordList n (w:rest) -> Just ( w2i (shiftR w  48) , partitionTail compact )++--------------------------------------------------------------------------------++-- | @partitionTail p == snd (uncons p)@+partitionTail :: Partition -> Partition+partitionTail compact = case compact of++  Nibble  0 -> Nibble 0+  Nibble  w -> Nibble $ shiftL (w .&. 0x0ffffffffffffff0) 4 + ((w .&. 15) - 1) ++  Medium1 w1 ->+    let !y = (shiftR w1 48) .&. 255     -- next element+        !n = w1 .&. 15+    in  if y <= 15 +          then makeNibble (w2i $ n-1) $ safeTail $ toList compact+          else Medium1    $ shiftL (w1 .&. 0x00ffffffffffff00) 8 + (n-1) ++  Medium2 w1 w2 ->      +    let !y = (shiftR w1 48) .&. 255+        !n = w2 .&. 255+    in  if y <= 15 +          then makeNibble (w2i $ n-1) $ safeTail $ toList compact+          else if n <= 8+            then Medium1 $ shiftL (w1 .&. 0x00ffffffffffffff) 8 + shiftL (shiftR w2 56) 8 + (n-1) ++            else Medium2 ( shiftL  w1 8 + shiftR w2 56 ) +                         ( shiftL (w2 .&. 0x00ffffffffffff00) 8 + (n-1) )++  Medium3 w1 w2 w3 ->   +    let !y = (shiftR w1 48) .&. 255+        !n = w3 .&. 255+    in  if y <= 15 && n <= 16+          then makeNibble (w2i $ n-1) $ safeTail $ toList compact+          else if n <= 16+            then Medium2 ( shiftL  w1 8 + shiftR w2 56 ) +                         ( shiftL  w2 8 + shiftR w3 56 + shiftL (shiftR w3 56) 8 + (n-1) )+ +            else Medium3 ( shiftL  w1 8 + shiftR w2 56 ) +                         ( shiftL  w2 8 + shiftR w3 56 ) +                         ( shiftL (w3 .&. 0x00ffffffffffff00) 8 + (n-1) )+ +  _ -> +    let n = width compact+    in  fromDescList' (n-1) $ safeTail $ toList compact ++--------------------------------------------------------------------------------++-- | We assume that @x >= partitionHeight p@!+cons :: Int -> Partition -> Partition+cons !x !compact = case compact of++  Nibble 0                -> singleton x+  +  Nibble word+    | x <= 15  && n < 15  -> Nibble $ shiftR word 4 + shiftL xw 60 + (n+1)+    | x <= 255            -> makeMedium   (w2i $ n+1) (x : toList compact)+    | otherwise           -> makeWordList (w2i $ n+1) (x : toList compact)+    where  +      n  = word .&. 15+      xw = i2w x+      +  Medium1 w1+    | x <= 255 && n < 7   -> Medium1 (shiftR w1 8 + shiftL xw 56 + (n+1))+    | x <= 255            -> Medium2 (shiftR w1 8 + shiftL xw 56        ) 8+    | otherwise           -> makeWordList (w2i $ n+1) (x : toList compact)+    where  +      n  = w1 .&. 255+      xw = i2w x++  Medium2 w1 w2+    | x <= 255 && n < 15  -> Medium2 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56 + (n+1))+    | x <= 255            -> Medium3 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56        ) 16+    | otherwise           -> makeWordList (w2i $ n+1) (x : toList compact)+    where  +      n  = w2 .&. 255+      xw = i2w x++  Medium3 w1 w2 w3+    | x <= 255 && n < 23  -> Medium3 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56) (shiftR w3 8 + shiftL (w2 .&. 255) 56 + (n+1))+    | x <= 255            -> Medium4 (shiftR w1 8 + shiftL xw 56) (shiftR w2 8 + shiftL (w1 .&. 255) 56) (shiftR w3 8 + shiftL (w2 .&. 255) 56        ) 24+    | otherwise           -> makeWordList (w2i $ n+1) (x : toList compact)+    where  +      n  = w3 .&. 255+      xw = i2w x++  Medium4 w1 w2 w3 w4+    | x <= 255 && n < 31  -> Medium4 (shiftR w1 8 + shiftL  xw          56) +                                     (shiftR w2 8 + shiftL (w1 .&. 255) 56) +                                     (shiftR w3 8 + shiftL (w2 .&. 255) 56) +                                     (shiftR w4 8 + shiftL (w3 .&. 255) 56 + (n+1))+    | otherwise           -> makeWordList (w2i $ n+1) (x : toList compact)+    where  +      n = w4 .&. 255+      xw = i2w x+      +  _ -> +    let n = width compact+    in  fromDescList' (n+1) (x : toList compact)++--------------------------------------------------------------------------------++-- | We assume that the element is not bigger than the last element!+snoc :: Partition -> Int -> Partition+snoc !compact  0 = compact+snoc !compact !x = case compact of++  Nibble 0 -> singleton x++  Nibble word+    | n < 15    -> Nibble $ (word + 1) .|. shiftL (i2w x) ((15-n)*4)+    | otherwise -> makeMedium (n+1) (toList compact ++ [x])+    where  +      n = w2i (word .&. 15)+      +  Medium1 w1+    | n < 7     -> Medium1 $ (w1 + 1) .|. shiftL (i2w x) ((7-n)*8)+    | otherwise -> Medium2 ((w1 .&. 0xffffffffffffff00) + i2w x) 8+    where  +      n = w2i (w1 .&. 255)++  Medium2 w1 w2+    | n < 15    -> Medium2 w1 $ (w2 + 1) .|. shiftL (i2w x) ((15-n)*8)+    | otherwise -> Medium3 w1 ((w2 .&. 0xffffffffffffff00) + i2w x) 16+    where  +      n = w2i (w2 .&. 255)++  Medium3 w1 w2 w3+    | n < 23    -> Medium3 w1 w2 $ (w3 + 1) .|. shiftL (i2w x) ((23-n)*8)+    | otherwise -> Medium4 w1 w2 ((w3 .&. 0xffffffffffffff00) + i2w x) 24+    where  +      n = w2i (w3 .&. 255)++  Medium4 w1 w2 w3 w4+    | n < 31    -> Medium4 w1 w2 w3 $ (w4 + 1) .|. shiftL (i2w x) ((31-n)*8)+    | otherwise -> makeWordList (n + 1) (toList compact ++ [x])+    where  +      n = w2i (w4 .&. 255)+  +  WordList n list -> WordList (n+1) (go list) where+    go :: [Word64] -> [Word64]+    go (w:[]) = case mod n 4 of+                  0 -> w : shiftL (i2w x) 48 : []+                  1 -> w + shiftL (i2w x) 32 : []+                  2 -> w + shiftL (i2w x) 16 : []+                  3 -> w +        (i2w x)    : []+    go (w:ws) = w : go ws+    go []     = shiftL (i2w x) 48 : []+ +{-   +  _ -> +    let n = width compact+    in  makeWordList (n+1) (toList compact ++ [x])+-}++--------------------------------------------------------------------------------+-- * exponential form++toExponentialForm :: Partition -> [(Int,Int)]+toExponentialForm = map (\xs -> (head xs,length xs)) . group . toAscList++fromExponentialForm :: [(Int,Int)] -> Partition+fromExponentialForm = fromDescList . concatMap f . sortBy g where+  f (!i,!e) = replicate e i+  g (!i, _) (!j,_) = compare j i++--------------------------------------------------------------------------------+-- * Width and height of the bounding rectangle++-- | Width, or the number of parts+width :: Partition -> Int+width compact = case compact of+  Nibble        word -> w2i (word .&.  15)+  Medium1       word -> w2i (word .&. 255)+  Medium2 _     word -> w2i (word .&. 255)+  Medium3 _ _   word -> w2i (word .&. 255)+  Medium4 _ _ _ word -> w2i (word .&. 255)+  WordList n    _    -> n++-- | Height, or the first (that is, the largest) element+height :: Partition -> Int+height compact = case compact of+  Nibble  word        -> w2i (shiftR word 60)+  Medium1 word        -> w2i (shiftR word 56)+  Medium2 word _      -> w2i (shiftR word 56)+  Medium3 word _ _    -> w2i (shiftR word 56)+  Medium4 word _ _ _  -> w2i (shiftR word 56)+  WordList _ (word:_) -> w2i (shiftR word 48)++-- | Width and height +widthHeight :: Partition -> (Int,Int)+widthHeight compact = case compact of +  Nibble  word            -> ( w2i (word  .&.  15) , w2i (shiftR word  60) )+  Medium1 word            -> ( w2i (word  .&. 255) , w2i (shiftR word  56) )+  Medium2 word1     word2 -> ( w2i (word2 .&. 255) , w2i (shiftR word1 56) )+  Medium3 word1 _   word3 -> ( w2i (word3 .&. 255) , w2i (shiftR word1 56) )+  Medium4 word1 _ _ word4 -> ( w2i (word4 .&. 255) , w2i (shiftR word1 56) )+  WordList n (word:_)     -> ( n                   , w2i (shiftR word  48) )++--------------------------------------------------------------------------------+-- * Differential sequence++-- | From a non-increasing sequence @[a1,a2,..,an]@ this computes the sequence of differences+-- @[a1-a2,a2-a3,...,an-0]@+diffSequence :: Partition -> [Int]+diffSequence compact = case compact of++  Nibble 0 -> []++  Nibble w -> +    let !nw = (w .&. 15) +        !w' = w - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w (60 - i*4) - shiftR w' (56 - i*4)) .&. 15 | i<-[0..n-1] ]++  Medium1 w -> +    let !nw = (w .&. 255) +        !w' = w - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w (56 - i*8) - shiftR w' (48 - i*8)) .&. 255 | i<-[0..n-1] ]++  Medium2 w1 w2 -> +    let !nw = (w2 .&. 255) +        !w2' = w2 - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1  (48 - i*8)) .&. 255 | i<-[0..6]   ] ++ +        ( w2i $ (       w1            - shiftR w2   56       ) .&. 255               ) : +        [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2' (48 - i*8)) .&. 255 | i<-[0..n-9] ] ++  Medium3 w1 w2 w3 -> +    let !nw = (w3 .&. 255) +        !w3' = w3 - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1  (48 - i*8)) .&. 255 | i<-[0..6]    ] ++ +        ( w2i $ (       w1            - shiftR w2   56       ) .&. 255                ) : +        [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2  (48 - i*8)) .&. 255 | i<-[0..6]    ] +++        ( w2i $ (       w2            - shiftR w3   56       ) .&. 255                ) : +        [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3' (48 - i*8)) .&. 255 | i<-[0..n-17] ] ++  Medium4 w1 w2 w3 w4 -> +    let !nw = (w4 .&. 255) +        !w4' = w4 - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1  (48 - i*8)) .&. 255 | i<-[0..6]    ] ++ +        ( w2i $ (       w1            - shiftR w2   56       ) .&. 255                ) : +        [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2  (48 - i*8)) .&. 255 | i<-[0..6]    ] +++        ( w2i $ (       w2            - shiftR w3   56       ) .&. 255                ) : +        [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3  (48 - i*8)) .&. 255 | i<-[0..6]    ] +++        ( w2i $ (       w3            - shiftR w4   56       ) .&. 255                ) : +        [ w2i $ (shiftR w4 (56 - i*8) - shiftR w4' (48 - i*8)) .&. 255 | i<-[0..n-25] ] ++  WordList {} -> go (toList compact) where+    go (x:ys@(y:_)) = (x-y) : go ys +    go [x] = [x]+    go []  = []++----------------------------------------++-- | From a non-increasing sequence @[a1,a2,..,an]@ this computes the reversed sequence of differences+-- @[ a[n]-0 , a[n-1]-a[n] , ... , a[2]-a[3] , a[1]-a[2] ] @+reverseDiffSequence :: Partition -> [Int]+reverseDiffSequence compact = case compact of++  Nibble 0 -> []++  Nibble w -> +    let !nw = (w .&. 15) +        !w' = w - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w (60 - i*4) - shiftR w' (56 - i*4)) .&. 15 | i<-toZero (n-1) ]++  Medium1 w -> +    let !nw = (w .&. 255) +        !w' = w - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w (56 - i*8) - shiftR w' (48 - i*8)) .&. 255 | i<-toZero (n-1) ]++  Medium2 w1 w2 -> +    let !nw = (w2 .&. 255) +        !w2' = w2 - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2' (48 - i*8)) .&. 255 | i<-toZero (n-9) ] +++        ( w2i $ (       w1            - shiftR w2   56       ) .&. 255                   ) : +        [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1  (48 - i*8)) .&. 255 | i<-toZero 6     ]  +        +  Medium3 w1 w2 w3 -> +    let !nw = (w3 .&. 255) +        !w3' = w3 - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3' (48 - i*8)) .&. 255 | i<-toZero (n-17) ] +++        ( w2i $ (       w2            - shiftR w3   56       ) .&. 255                    ) : +        [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2  (48 - i*8)) .&. 255 | i<-toZero 6      ] +++        ( w2i $ (       w1            - shiftR w2   56       ) .&. 255                    ) : +        [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1  (48 - i*8)) .&. 255 | i<-toZero 6      ] ++  Medium4 w1 w2 w3 w4 -> +    let !nw = (w4 .&. 255) +        !w4' = w4 - nw+        !n  = w2i nw+    in  [ w2i $ (shiftR w4 (56 - i*8) - shiftR w4' (48 - i*8)) .&. 255 | i<-toZero (n-25) ] +++        ( w2i $ (       w3            - shiftR w4   56       ) .&. 255                    ) : +        [ w2i $ (shiftR w3 (56 - i*8) - shiftR w3  (48 - i*8)) .&. 255 | i<-toZero 6      ] +++        ( w2i $ (       w2            - shiftR w3   56       ) .&. 255                    ) : +        [ w2i $ (shiftR w2 (56 - i*8) - shiftR w2  (48 - i*8)) .&. 255 | i<-toZero 6      ] +++        ( w2i $ (       w1            - shiftR w2   56       ) .&. 255                    ) : +        [ w2i $ (shiftR w1 (56 - i*8) - shiftR w1  (48 - i*8)) .&. 255 | i<-toZero 6      ] ++  WordList {} -> (h : go asclist) where+    asclist@(h:_) = toAscList compact+    go (x:ys@(y:_)) = (y-x) : go ys +    go [_] = []+    go []  = []++--------------------------------------------------------------------------------+-- *  Dual partition++foreign import ccall unsafe "c_dual_nibble" c_dual_nibble :: Word64 -> Word64++dualPartition :: Partition -> Partition+dualPartition compact = case compact of++  Nibble 0 -> Nibble 0+  Nibble w -> Nibble (c_dual_nibble w)  +  _        -> if (w <= 255 && h <= 31)+                then makeMedium   h dualList+                else makeWordList h dualList+  where+    (w,h) = widthHeight compact+    dualList = concat+      [ replicate d j+      | (j,d) <- zip (toOne w) (reverseDiffSequence compact)+      ]++--------------------------------------------------------------------------------+-- * Conversion to list++toList :: Partition -> [Int]+toList = toDescList++-- | returns a descending (non-increasing) list+toDescList :: Partition -> [Int]+toDescList compact = case compact of++  Nibble 0 -> []++  Nibble word -> +    let !n = w2i (word .&. 15) +    in  [ w2i (shiftR word  (60 - i*4) .&. 15 ) | i<-[0..n-1] ]++  Medium1 word1 -> +    let !n = w2i (word1 .&. 255)+    in  [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..n-1] ]++  Medium2 word1 word2 -> +    let !n = w2i (word2 .&. 255) +    in  [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..7]   ] +++        [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-[0..n-9] ] ++  Medium3 word1 word2 word3 -> +    let !n = w2i (word3 .&. 255) +    in  [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..7]    ] +++        [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-[0..7]    ] +++        [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-[0..n-17] ] ++  Medium4 word1 word2 word3 word4 -> +    let !n = w2i (word4 .&. 255) +    in  [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-[0..7]    ] +++        [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-[0..7]    ] +++        [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-[0..7]    ] +++        [ w2i (shiftR word4 (56 - i*8) .&. 255) | i<-[0..n-25] ] ++  WordList _ list -> go list where+    go :: [Word64] -> [Int]+    go !wlist = case wlist of+      (!w):(!ws) -> case ws of +        (_:_)      -> w2i (shiftR w 48          ) :+                      w2i (shiftR w 32 .&. 65535) :+                      w2i (shiftR w 16 .&. 65535) :+                      w2i (       w    .&. 65535) : go ws+        []         -> takeWhile (/=0) (fromWord w)+      []         -> []++    fromWord :: Word64 -> [Int]+    fromWord !word = +      [ w2i (shiftR word 48          )+      , w2i (shiftR word 32 .&. 65535)+      , w2i (shiftR word 16 .&. 65535)+      , w2i (       word    .&. 65535)+      ]++----------------------------------------++-- | Returns a reversed (ascending; non-decreasing) list+toAscList :: Partition -> [Int]+toAscList compact = case compact of++  Nibble 0 -> []++  Nibble word -> +    let !n = w2i (word .&. 15) +    in  [ w2i (shiftR word  (60 - i*4) .&. 15 ) | i<-toZero (n-1) ]++  Medium1 word1 -> +    let !n = w2i (word1 .&. 255)+    in  [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero (n-1) ]++  Medium2 word1 word2 -> +    let !n = w2i (word2 .&. 255) +    in  [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-toZero (n-9) ] +++        [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero 7     ] ++  Medium3 word1 word2 word3 -> +    let !n = w2i (word3 .&. 255) +    in  [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-toZero (n-17) ] +++        [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-toZero 7      ] +++        [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero 7      ]+        +  Medium4 word1 word2 word3 word4 -> +    let !n = w2i (word4 .&. 255) +    in  [ w2i (shiftR word4 (56 - i*8) .&. 255) | i<-toZero (n-25) ] +++        [ w2i (shiftR word3 (56 - i*8) .&. 255) | i<-toZero 7      ] +++        [ w2i (shiftR word2 (56 - i*8) .&. 255) | i<-toZero 7      ] +++        [ w2i (shiftR word1 (56 - i*8) .&. 255) | i<-toZero 7      ]++  WordList _ list -> dropWhile (==0) $ go (reverse list) where+    go :: [Word64] -> [Int]+    go !wlist = case wlist of+      (!w):ws -> w2i (       w    .&. 65535) : +                 w2i (shiftR w 16 .&. 65535) :+                 w2i (shiftR w 32 .&. 65535) :+                 w2i (shiftR w 48          ) : go ws+      [] -> []++{-+    go :: [Word64] -> [Int]+    go (w:[]) = fromWord w+    go (w:ws) = fromWord w ++ go ws+    go []     = []+    fromWord :: Word64 -> [Int]+    fromWord word = [ w2i (shiftR word (48 - i*16) .&. 65535) | i<-toZero 3 ] +-}++--------------------------------------------------------------------------------+-- * Conversion from list++fromDescList :: [Int] -> Partition+fromDescList list = fromDescList' (length list) list++-- | We assume that the input is a non-increasing list of /positive/ integers!+fromDescList' +  :: Int          -- ^ length+  -> [Int]        -- ^ the list+  -> Partition+fromDescList' !n !list =+  case list of+    []                           -> empty+    (h:_) | h <= 0               -> empty+          | h <= 15 && n <= 15   -> makeNibble   n list+          | h >  65535           -> error "partitions with elements bigger than 65535 are not supported"+          | h >  255 || n > 31   -> makeWordList n list+          | otherwise            -> makeMedium   n list++makeNibble :: Int -> [Int] -> Partition+makeNibble !n list = Nibble $ go (i2w n) 60 list where+  go !acc !k (x:xs) = go (acc + shiftL (i2w x) k) (k-4) xs+  go !acc _  []     = acc+{-   +makeNibble :: Int -> [Int] -> Partition+makeNibble !n list = Nibble +  $ sum' [ shiftL (i2w x) (60 - 4*i) | (i,x) <- zip [0..] list ] +  + i2w n+-}++makeMedium :: Int -> [Int] -> Partition+makeMedium !n list +  | n <= 7   = makeMedium1 n list+  | n <= 15  = makeMedium2 n list+  | n <= 23  = makeMedium3 n list+  | n <= 31  = makeMedium4 n list+  | otherwise = error "makeMedium: input list too big (should be smaller than 32)"++makeMedium1 :: Int -> [Int] -> Partition+makeMedium1 !n list = Medium1 +  $ sum' [ shiftL (fromIntegral x) (56 - 8*i) | (i,x) <- zip [0..] list ] +  + fromIntegral n++makeMedium2 :: Int -> [Int] -> Partition+makeMedium2 !n list = Medium2 word1 word2 where+  (list1,list2) = splitAt 8 list+  word1 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list1 ] +  word2 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list2 ] +        + fromIntegral n++makeMedium3 :: Int -> [Int] -> Partition+makeMedium3 !n list = Medium3 word1 word2 word3 where+  (list1,tmp  ) = splitAt 8 list+  (list2,list3) = splitAt 8 tmp+  word1 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list1 ] +  word2 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list2 ] +  word3 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list3 ] +        + i2w n++makeMedium4 :: Int -> [Int] -> Partition+makeMedium4 !n list = Medium4 word1 word2 word3 word4 where+  (list1,tmp1 ) = splitAt 8 list+  (list2,tmp2 ) = splitAt 8 tmp1+  (list3,list4) = splitAt 8 tmp2+  word1 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list1 ] +  word2 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list2 ] +  word3 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list3 ] +  word4 = sum' [ shiftL (i2w x) (56 - 8*i) | (i,x) <- zip [0..] list4 ] +        + i2w n+    +makeWordList :: Int -> [Int] -> Partition+makeWordList !n list = WordList n (go list) where   +  go :: [Int] -> [Word64]+  go !xs = case xs of+    (x:y:z:w:rest) -> makeWord x y z w : go rest+    (x:y:z:  []  ) -> makeWord x y z 0 : []+    (x:y:    []  ) -> makeWord x y 0 0 : []+    (x:      []  ) -> makeWord x 0 0 0 : []+    []             -> []+  makeWord !x !y !z !w = shiftL (i2w x) 48  +                       + shiftL (i2w y) 32  +                       + shiftL (i2w z) 16 +                       +        (i2w w)+{-+  go [] = []+  go xs = case splitAt 4 xs of+    (this,rest) -> case rest of+      [] -> makeWord (take 4 $ this ++ repeat 0) : []+      _  -> makeWord this : go rest+  makeWord [x,y,z,w] = shiftL (i2w x) 48  +                     + shiftL (i2w y) 32  +                     + shiftL (i2w z) 16 +                     +        (i2w w)+-}++--------------------------------------------------------------------------------+-- * Partial orderings++isSubPartitionOf :: Partition -> Partition -> Bool+isSubPartitionOf p q = case (p,q) of++  (Nibble 0 , _       ) -> True+  +  (Nibble u , Nibble v) -> let !n = w2i (u .&. 15) +                           in  and [    (shiftR u (60 - i*4) .&. 15)+                                     <= (shiftR v (60 - i*4) .&. 15) +                                   | i<-[0..n-1] +                                   ]++  _                     -> and $ zipWith (<=) (toList p) (toList q ++ repeat 0)++dominates :: Partition -> Partition -> Bool+dominates q p = case (q,p) of++  (_        , Nibble 0 ) -> True++  (Nibble v , Nibble u ) -> go 60 0 0 where+                              n = u .&. 15                              +                              klimit = w2i (4*(15-n))+                              go !k !b !a = if k <= klimit +                                then True+                                else let !b' = b + (shiftR v k .&. 15)+                                         !a' = a + (shiftR u k .&. 15)+                                     in  if b' < a' +                                           then False +                                           else go (k-4) b' a'++  _                      -> and $ zipWith (>=) (sums $ toList q ++ repeat 0) (sums $ toList p) where+                              sums = tail . scanl' (+) 0++--------------------------------------------------------------------------------+-- * Pieri rule++-- | Expands to product @s[lambda]*h[1] = s[lambda]*e[1]@ as a sum of @s[mu]@-s. See <https://en.wikipedia.org/wiki/Pieri's_formula>+pieriRuleSingleBox :: Partition -> [Partition]+pieriRuleSingleBox !compact = case compact of++  Nibble 0 -> [ singleton 1 ]++  Nibble w | h < 15 -> +    [ Nibble  (w + shiftL 1 (60-4*i)) | (i,d)<-zip [0..n-1] diffs1 , d>0 ] ++ [ snoc compact 1 ]++  Medium1 w | h < 255 -> +    [ Medium1 (w + shiftL 1 (56-8*i)) | (i,d)<-zip [0..n-1] diffs1 , d>0 ] ++ [ snoc compact 1 ]++  Medium2 w1 w2 | h < 255 -> +    let (diffs1a,diffs1b) = splitAt 8 diffs1 +    in  [ Medium2    (w1 + shiftL 1 (56-8*i)) w2 | (i,d)<-zip [0..7  ] diffs1a , d>0 ] +++        [ Medium2 w1 (w2 + shiftL 1 (56-8*i))    | (i,d)<-zip [0..n-9] diffs1b , d>0 ] +++        [ snoc compact 1 ]++  Medium3 w1 w2 w3 | h < 255 -> +    let (diffs1a,tmp    ) = splitAt 8 diffs1 +        (diffs1b,diffs1c) = splitAt 8 tmp+    in  [ Medium3       (w1 + shiftL 1 (56-8*i)) w2 w3 | (i,d)<-zip [0..7   ] diffs1a , d>0 ] +++        [ Medium3    w1 (w2 + shiftL 1 (56-8*i)) w3    | (i,d)<-zip [0..7   ] diffs1b , d>0 ] +++        [ Medium3 w1 w2 (w3 + shiftL 1 (56-8*i))       | (i,d)<-zip [0..n-17] diffs1c , d>0 ] +++        [ snoc compact 1 ]+    +  _ -> genericSingleBox++  where+    (n,h)  =     widthHeight  compact+    list   =     toDescList   compact+    diffs1 = 1 : diffSequence compact++    genericSingleBox :: [Partition]+    genericSingleBox = map (fromDescList' n) (go list diffs1) ++ [ fromDescList' (n+1) (list ++ [1]) ] where+      go :: [Int] -> [Int] -> [[Int]]+      go (a:as) (d:ds) = if d > 0 then ((a+1):as) : map (a:) (go as ds) +                                  else              map (a:) (go as ds)+      go []     _      = []++-- | Expands to product @s[lambda]*h[k]@ as a sum of @s[mu]@-s. See <https://en.wikipedia.org/wiki/Pieri's_formula>+pieriRule :: Partition -> Int -> [Partition]+pieriRule !compact !k +  | k <  0                  = []+  | k == 0                  = [ compact ]+  | k == 1                  = pieriRuleSingleBox compact+  | h == 0                  = [ singleton k ]+  | h + k <= 15  && n < 15  = case compact of { Nibble w -> +                              [ Nibble (w + encode c)  | c <- comps ] }+  | otherwise               = [ fromDescList' (n+b) xs | c <- comps , let (b,xs) = add c ] ++  where+    (n,h)  = widthHeight compact+    list   = toDescList compact+    bounds = k : {- map (min k) -} (diffSequence compact) +    comps = compositions' bounds k++    add clist = go list clist where+      go (!p:ps) (!c:cs) = let (b,rest) = go ps cs in (b, (p+c):rest)+      go []      [c]     = if c>0 then (1,[c]) else (0,[])+      go _       _       = error "Compact/pieriRule/add: shouldn't happen"++    encode :: [Int] -> Word64+    encode = go 60 where+      go !k [c]    = if c==0 then 0 else shiftL (i2w c) k + 1+      go !k (c:cs) = shiftL (i2w c) k + go (k-4) cs+      go !k []     = error "Compact/pieriRule/encode: shouldn't happen"++--------------------------------------------------------------------------------+-- * local (internally used) utility functions++{-# INLINE i2w #-}+i2w :: Int -> Word64+i2w = fromIntegral++{-# INLINE w2i #-}+w2i :: Word64 -> Int+w2i = fromIntegral++{-# INLINE sum' #-}+sum' :: [Word64] -> Word64+sum' = foldl' (+) 0++{-# INLINE safeTail #-}+safeTail :: [Int] -> [Int]+safeTail xs = case xs of { [] -> [] ; _ -> tail xs }++{-# INLINE toZero #-}+toZero :: Int -> [Int]+toZero !n+  | n >  0  = n : toZero (n-1) +  | n == 0  = [0]+  | n <  0  = []++{-# INLINE toOne #-}+toOne :: Int -> [Int]+toOne !n+  | n >  1  = n : toOne (n-1) +  | n == 1  = [1]+  | n <  1  = []++--------------------------------------------------------------------------------++

@@ -0,0 +1,215 @@++-- | Counting partitions of integers.++{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables #-}+module Math.Combinat.Partitions.Integer.Count where++--------------------------------------------------------------------------------++import Data.List+import Control.Monad ( liftM , replicateM )++-- import Data.Map (Map)+-- import qualified Data.Map as Map++import Math.Combinat.Numbers ( factorial , binomial , multinomial )+import Math.Combinat.Numbers.Integers -- Primes+import Math.Combinat.Helper++import Data.Array+import System.Random++--------------------------------------------------------------------------------+-- * Infinite tables of integers++-- | A data structure which is essentially an infinite list of @Integer@-s,+-- but fast lookup (for reasonable small inputs)+newtype TableOfIntegers = TableOfIntegers [Array Int Integer]++lookupInteger :: TableOfIntegers -> Int -> Integer+lookupInteger (TableOfIntegers table) !n +  | n >= 0  = (table !! k) ! r+  | n <  0  = 0+  where+    (k,r) = divMod n 1024++makeTableOfIntegers+  :: ((Int -> Integer) -> (Int -> Integer))+  -> TableOfIntegers+makeTableOfIntegers user = table where+  calc  = user lkp+  lkp   = lookupInteger table+  table = TableOfIntegers+    [ listArray (0,1023) (map calc [a..b]) +    | k<-[0..] +    , let a = 1024*k +    , let b = 1024*(k+1) - 1 +    ]++--------------------------------------------------------------------------------+-- * Counting partitions++-- | Number of partitions of @n@ (looking up a table built using Euler's algorithm)+countPartitions :: Int -> Integer+countPartitions = lookupInteger partitionCountTable ++-- | This uses the power series expansion of the infinite product. It is slower than the above.+countPartitionsInfiniteProduct :: Int -> Integer+countPartitionsInfiniteProduct k = partitionCountListInfiniteProduct !! k++-- | This uses 'countPartitions'', and is (very) slow+countPartitionsNaive :: Int -> Integer+countPartitionsNaive d = countPartitions' (d,d) d++--------------------------------------------------------------------------------++-- | This uses Euler's algorithm to compute p(n)+--+-- See eg.:+-- NEIL CALKIN, JIMENA DAVIS, KEVIN JAMES, ELIZABETH PEREZ, AND CHARLES SWANNACK+-- COMPUTING THE INTEGER PARTITION FUNCTION+-- <http://www.math.clemson.edu/~kevja/PAPERS/ComputingPartitions-MathComp.pdf>+--+partitionCountTable :: TableOfIntegers+partitionCountTable = table where++  table = makeTableOfIntegers fun++  fun lkp !n +    | n >  1 = foldl' (+) 0 +             [ (if even k then negate else id) +                 ( lkp (n - div (k*(3*k+1)) 2)+                 + lkp (n - div (k*(3*k-1)) 2)+                 )+             | k <- [1..limit n]+             ]+    | n <  0 = 0+    | n == 0 = 1+    | n == 1 = 1++  limit :: Int -> Int+  limit !n = fromInteger $ ceilingSquareRoot (1 + div (nn+nn+1) 3) where+    nn = fromIntegral n :: Integer++-- | An infinite list containing all @p(n)@, starting from @p(0)@.+partitionCountList :: [Integer]+partitionCountList = map countPartitions [0..]++--------------------------------------------------------------------------------++-- | Infinite list of number of partitions of @0,1,2,...@+--+-- This uses the infinite product formula the generating function of partitions, +-- recursively expanding it; it is reasonably fast for small numbers.+--+-- > partitionCountListInfiniteProduct == map countPartitions [0..]+--+partitionCountListInfiniteProduct :: [Integer]+partitionCountListInfiniteProduct = final where++  final = go 1 (1:repeat 0) ++  go !k (x:xs) = x : go (k+1) ys where+    ys = zipWith (+) xs (take k final ++ ys)+    -- explanation:+    --   xs == drop k $ f (k-1)+    --   ys == drop k $ f (k  )  ++{-++Full explanation of 'partitionCountList':+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~++let f k = productPSeries $ map (:[]) [1..k]++f 0 = [1,0,0,0,0,0,0,0...]+f 1 = [1,1,1,1,1,1,1,1...]+f 2 = [1,1,2,2,3,3,4,4...]+f 3 = [1,1,2,3,4,5,7,8...]++observe: ++* take (k+1) (f k) == take (k+1) partitionCountList+* f (k+1) == zipWith (+) (f k) (replicate (k+1) 0 ++ f (k+1))++now apply (drop (k+1)) to the second one : ++* drop (k+1) (f (k+1)) == zipWith (+) (drop (k+1) $ f k) (f (k+1))+* f (k+1) = take (k+1) final ++ drop (k+1) (f (k+1))++-}++--------------------------------------------------------------------------------++-- | Naive infinite list of number of partitions of @0,1,2,...@+--+-- > partitionCountListNaive == map countPartitionsNaive [0..]+--+-- This is very slow.+--+partitionCountListNaive :: [Integer]+partitionCountListNaive = map countPartitionsNaive [0..]++--------------------------------------------------------------------------------+-- * Counting all partitions++countAllPartitions :: Int -> Integer+countAllPartitions d = sum' [ countPartitions i | i <- [0..d] ]++-- | Count all partitions fitting into a rectangle.+-- # = \\binom { h+w } { h }+countAllPartitions' :: (Int,Int) -> Integer+countAllPartitions' (h,w) = +  binomial (h+w) (min h w)+  --sum [ countPartitions' (h,w) i | i <- [0..d] ] where d = h*w++--------------------------------------------------------------------------------+-- * Counting fitting into a rectangle++-- | Number of of d, fitting into a given rectangle. Naive recursive algorithm.+countPartitions' :: (Int,Int) -> Int -> Integer+countPartitions' _ 0 = 1+countPartitions' (0,_) d = if d==0 then 1 else 0+countPartitions' (_,0) d = if d==0 then 1 else 0+countPartitions' (h,w) d = sum+  [ countPartitions' (i,w-1) (d-i) | i <- [1..min d h] ] ++--------------------------------------------------------------------------------+-- * Partitions with given number of parts++-- | Count partitions of @n@ into @k@ parts.+--+-- Naive recursive algorithm.+--+countPartitionsWithKParts +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = the integer we partition+  -> Integer+countPartitionsWithKParts k n = go n k n where+  go !h !k !n +    | k <  0     = 0+    | k == 0     = if h>=0 && n==0 then 1 else 0+    | k == 1     = if h>=n && n>=1 then 1 else 0+    | otherwise  = sum' [ go a (k-1) (n-a) | a<-[1..(min h (n-k+1))] ]++--------------------------------------------------------------------------------+-- Partitions with only odd\/distinct parts++{-+-- | Partitions of @n@ with only odd parts+partitionsWithOddParts :: Int -> [Partition]+partitionsWithOddParts d = map Partition (go d d) where+  go _  0  = [[]]+  go !h !n = [ a:as | a<-[1,3..min n h], as <- go a (n-a) ]+-}++{-+-- > length (partitionsWithDistinctParts d) == length (partitionsWithOddParts d)+--+partitionsWithDistinctParts :: Int -> [Partition]+partitionsWithDistinctParts d = map Partition (go d d) where+  go _  0  = [[]]+  go !h !n = [ a:as | a<-[1..min n h], as <- go (a-1) (n-a) ]+-}++--------------------------------------------------------------------------------

@@ -0,0 +1,398 @@++-- | Partition functions working on lists of integers.+-- +-- It's not recommended to use this module directly.++{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables #-}+module Math.Combinat.Partitions.Integer.IntList where++--------------------------------------------------------------------------------++import Data.List+import Control.Monad ( liftM , replicateM )++import Math.Combinat.Numbers ( factorial , binomial , multinomial )+import Math.Combinat.Helper++import Data.Array+import System.Random++import Math.Combinat.Partitions.Integer.Count ( countPartitions )++--------------------------------------------------------------------------------+-- * Type and basic stuff++-- | Sorts the input, and cuts the nonpositive elements.+_mkPartition :: [Int] -> [Int]+_mkPartition xs = sortBy (reverseCompare) $ filter (>0) xs+ +-- | This returns @True@ if the input is non-increasing sequence of +-- /positive/ integers (possibly empty); @False@ otherwise.+--+_isPartition :: [Int] -> Bool+_isPartition []  = True+_isPartition [x] = x > 0+_isPartition (x:xs@(y:_)) = (x >= y) && _isPartition xs+++_dualPartition :: [Int] -> [Int]+_dualPartition [] = []+_dualPartition xs = go 0 (_diffSequence xs) [] where+  go !i (d:ds) acc = go (i+1) ds (d:acc)+  go n  []     acc = finish n acc +  finish !j (k:ks) = replicate k j ++ finish (j-1) ks+  finish _  []     = []++--------------------------------------------------------------------------------++{-+-- more variations:++_dualPartition_b :: [Int] -> [Int]+_dualPartition_b [] = []+_dualPartition_b xs = go 1 (diffSequence xs) [] where+  go !i (d:ds) acc = go (i+1) ds ((d,i):acc)+  go _  []     acc = concatMap (\(d,i) -> replicate d i) acc++_dualPartition_c :: [Int] -> [Int]+_dualPartition_c [] = []+_dualPartition_c xs = reverse $ concat $ zipWith f [1..] (diffSequence xs) where+  f _ 0 = []+  f k d = replicate d k+-}++-- | A simpler, but bit slower (about twice?) implementation of dual partition+_dualPartitionNaive :: [Int] -> [Int]+_dualPartitionNaive [] = []+_dualPartitionNaive xs@(k:_) = [ length $ filter (>=i) xs | i <- [1..k] ]++-- | From a sequence @[a1,a2,..,an]@ computes the sequence of differences+-- @[a1-a2,a2-a3,...,an-0]@+_diffSequence :: [Int] -> [Int]+_diffSequence = go where+  go (x:ys@(y:_)) = (x-y) : go ys +  go [x] = [x]+  go []  = []++-- | Example:+--+-- > _elements [5,4,1] ==+-- >   [ (1,1), (1,2), (1,3), (1,4), (1,5)+-- >   , (2,1), (2,2), (2,3), (2,4)+-- >   , (3,1)+-- >   ]+--++_elements :: [Int] -> [(Int,Int)]+_elements shape = [ (i,j) | (i,l) <- zip [1..] shape, j<-[1..l] ] ++---------------------------------------------------------------------------------+-- * Exponential form++-- | We convert a partition to exponential form.+-- @(i,e)@ mean @(i^e)@; for example @[(1,4),(2,3)]@ corresponds to @(1^4)(2^3) = [2,2,2,1,1,1,1]@. Another example:+--+-- > toExponentialForm (Partition [5,5,3,2,2,2,2,1,1]) == [(1,2),(2,4),(3,1),(5,2)]+--+_toExponentialForm :: [Int] -> [(Int,Int)]+_toExponentialForm = reverse . map (\xs -> (head xs,length xs)) . group++_fromExponentialForm :: [(Int,Int)] -> [Int]+_fromExponentialForm = sortBy reverseCompare . go where+  go ((j,e):rest) = replicate e j ++ go rest+  go []           = []   ++---------------------------------------------------------------------------------+-- * Generating partitions++-- | Partitions of @d@, as lists+_partitions :: Int -> [[Int]]+_partitions d = go d d where+  go _  0  = [[]]+  go !h !n = [ a:as | a<-[1..min n h], as <- go a (n-a) ]++-- | All integer partitions up to a given degree (that is, all integer partitions whose sum is less or equal to @d@)+_allPartitions :: Int -> [[Int]]+_allPartitions d = concat [ _partitions i | i <- [0..d] ]++-- | All integer partitions up to a given degree (that is, all integer partitions whose sum is less or equal to @d@),+-- grouped by weight+_allPartitionsGrouped :: Int -> [[[Int]]]+_allPartitionsGrouped d = [ _partitions i | i <- [0..d] ]++---------------------------------------------------------------------------------++-- | Integer partitions of @d@, fitting into a given rectangle, as lists.+_partitions' +  :: (Int,Int)     -- ^ (height,width)+  -> Int           -- ^ d+  -> [[Int]]        +_partitions' _ 0 = [[]] +_partitions' ( 0 , _) d = if d==0 then [[]] else []+_partitions' ( _ , 0) d = if d==0 then [[]] else []+_partitions' (!h ,!w) d = +  [ i:xs | i <- [1..min d h] , xs <- _partitions' (i,w-1) (d-i) ]++---------------------------------------------------------------------------------+-- * Random partitions++-- | Uniformly random partition of the given weight. +--+-- NOTE: This algorithm is effective for small @n@-s (say @n@ up to a few hundred \/ one thousand it should work nicely),+-- and the first time it is executed may be slower (as it needs to build the table 'partitionCountList' first)+--+-- Algorithm of Nijenhuis and Wilf (1975); see+--+-- * Knuth Vol 4A, pre-fascicle 3B, exercise 47;+--+-- * Nijenhuis and Wilf: Combinatorial Algorithms for Computers and Calculators, chapter 10+--+_randomPartition :: RandomGen g => Int -> g -> ([Int], g)+_randomPartition n g = (p, g') where+  ([p], g') = _randomPartitions 1 n g++-- | Generates several uniformly random partitions of @n@ at the same time.+-- Should be a little bit faster then generating them individually.+--+_randomPartitions +  :: forall g. RandomGen g +  => Int   -- ^ number of partitions to generate+  -> Int   -- ^ the weight of the partitions+  -> g -> ([[Int]], g)+_randomPartitions howmany n = runRand $ replicateM howmany (worker n []) where++  cnt = countPartitions+ +  finish :: [(Int,Int)] -> [Int]+  finish = _mkPartition . concatMap f where f (j,d) = replicate j d++  fi :: Int -> Integer +  fi = fromIntegral++  find_jd :: Int -> Integer -> (Int,Int)+  find_jd m capm = go 0 [ (j,d) | j<-[1..n], d<-[1..div m j] ] where+    go :: Integer -> [(Int,Int)] -> (Int,Int)+    go !s []   = (1,1)       -- ??+    go !s [jd] = jd          -- ??+    go !s (jd@(j,d):rest) = +      if s' > capm +        then jd +        else go s' rest+      where+        s' = s + fi d * cnt (m - j*d)++  worker :: Int -> [(Int,Int)] -> Rand g [Int]+  worker  0 acc = return $ finish acc+  worker !m acc = do+    capm <- randChoose (0, (fi m) * cnt m - 1)+    let jd@(!j,!d) = find_jd m capm+    worker (m - j*d) (jd:acc)+++---------------------------------------------------------------------------------+-- * Dominance order ++-- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions+-- (this is partial ordering on the set of partitions of @n@).+--+-- See <http://en.wikipedia.org/wiki/Dominance_order>+--+_dominates :: [Int] -> [Int] -> Bool+_dominates qs ps+  = and $ zipWith (>=) (sums (qs ++ repeat 0)) (sums ps)+  where+    sums = scanl (+) 0++-- | Lists all partitions of the same weight as @lambda@ and also dominated by @lambda@+-- (that is, all partial sums are less or equal):+--+-- > dominatedPartitions lam == [ mu | mu <- partitions (weight lam), lam `dominates` mu ]+-- +_dominatedPartitions :: [Int] -> [[Int]]+_dominatedPartitions []     = [[]]+_dominatedPartitions lambda = go (head lambda) w dsums 0 where++  n = length lambda+  w = sum    lambda+  dsums = scanl1 (+) (lambda ++ repeat 0)++  go _   0 _       _  = [[]]+  go !h !w (!d:ds) !e  +    | w >  0  = [ (a:as) | a <- [1..min h (d-e)] , as <- go a (w-a) ds (e+a) ] +    | w == 0  = [[]]+    | w <  0  = error "_dominatedPartitions: fatal error; shouldn't happen"++-- | Lists all partitions of the sime weight as @mu@ and also dominating @mu@+-- (that is, all partial sums are greater or equal):+--+-- > dominatingPartitions mu == [ lam | lam <- partitions (weight mu), lam `dominates` mu ]+-- +_dominatingPartitions :: [Int] -> [[Int]]+_dominatingPartitions []     = [[]]+_dominatingPartitions mu     = go w w dsums 0 where++  n = length mu+  w = sum    mu+  dsums = scanl1 (+) (mu ++ repeat 0)++  go _   0 _       _  = [[]]+  go !h !w (!d:ds) !e  +    | w >  0  = [ (a:as) | a <- [max 0 (d-e)..min h w] , as <- go a (w-a) ds (e+a) ] +    | w == 0  = [[]]+    | w <  0  = error "_dominatingPartitions: fatal error; shouldn't happen"++--------------------------------------------------------------------------------+-- * Partitions with given number of parts++-- | Lists partitions of @n@ into @k@ parts.+--+-- > sort (partitionsWithKParts k n) == sort [ p | p <- partitions n , numberOfParts p == k ]+--+-- Naive recursive algorithm.+--+_partitionsWithKParts +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = the integer we partition+  -> [[Int]]+_partitionsWithKParts k n = go n k n where+{-+  h = max height+  k = number of parts+  n = integer+-}+  go !h !k !n +    | k <  0     = []+    | k == 0     = if h>=0 && n==0 then [[] ] else []+    | k == 1     = if h>=n && n>=1 then [[n]] else []+    | otherwise  = [ a:p | a <- [1..(min h (n-k+1))] , p <- go a (k-1) (n-a) ]++--------------------------------------------------------------------------------+-- * Partitions with only odd\/distinct parts++-- | Partitions of @n@ with only odd parts+_partitionsWithOddParts :: Int -> [[Int]]+_partitionsWithOddParts d = (go d d) where+  go _  0  = [[]]+  go !h !n = [ a:as | a<-[1,3..min n h], as <- go a (n-a) ]++{-+-- | Partitions of @n@ with only even parts+--+-- Note: this is not very interesting, it's just @(map.map) (2*) $ _partitions (div n 2)@+--+_partitionsWithEvenParts :: Int -> [[Int]]+_partitionsWithEvenParts d = (go d d) where+  go _  0  = [[]]+  go !h !n = [ a:as | a<-[2,4..min n h], as <- go a (n-a) ]+-}++-- | Partitions of @n@ with distinct parts.+-- +-- Note:+--+-- > length (partitionsWithDistinctParts d) == length (partitionsWithOddParts d)+--+_partitionsWithDistinctParts :: Int -> [[Int]]+_partitionsWithDistinctParts d = (go d d) where+  go _  0  = [[]]+  go !h !n = [ a:as | a<-[1..min n h], as <- go (a-1) (n-a) ]++--------------------------------------------------------------------------------+-- * Sub- and super-partitions of a given partition++-- | Returns @True@ of the first partition is a subpartition (that is, fit inside) of the second.+-- This includes equality+_isSubPartitionOf :: [Int] -> [Int] -> Bool+_isSubPartitionOf ps qs = and $ zipWith (<=) ps (qs ++ repeat 0)++-- | This is provided for convenience\/completeness only, as:+--+-- > isSuperPartitionOf q p == isSubPartitionOf p q+--+_isSuperPartitionOf :: [Int] -> [Int] -> Bool+_isSuperPartitionOf qs ps = and $ zipWith (<=) ps (qs ++ repeat 0)+++-- | Sub-partitions of a given partition with the given weight:+--+-- > sort (subPartitions d q) == sort [ p | p <- partitions d, isSubPartitionOf p q ]+--+_subPartitions :: Int -> [Int] -> [[Int]]+_subPartitions d big+  | null big       = if d==0 then [[]] else []+  | d > sum' big   = []+  | d < 0          = []+  | otherwise      = go d (head big) big+  where+    go :: Int -> Int -> [Int] -> [[Int]]+    go !k !h []      = if k==0 then [[]] else []+    go !k !h (b:bs) +      | k<0 || h<0   = []+      | k==0         = [[]]+      | h==0         = []+      | otherwise    = [ this:rest | this <- [1..min h b] , rest <- go (k-this) this bs ]++----------------------------------------++-- | All sub-partitions of a given partition+_allSubPartitions :: [Int] -> [[Int]]+_allSubPartitions big +  | null big   = [[]]+  | otherwise  = go (head big) big+  where+    go _  [] = [[]]+    go !h (b:bs) +      | h==0         = []+      | otherwise    = [] : [ this:rest | this <- [1..min h b] , rest <- go this bs ]++----------------------------------------++-- | Super-partitions of a given partition with the given weight:+--+-- > sort (superPartitions d p) == sort [ q | q <- partitions d, isSubPartitionOf p q ]+--+_superPartitions :: Int -> [Int] -> [[Int]]+_superPartitions dd small+  | dd < w0     = []+  | null small  = _partitions dd+  | otherwise   = go dd w1 dd (small ++ repeat 0)+  where+    w0 = sum' small+    w1 = w0 - head small+    -- d = remaining weight of the outer partition we are constructing+    -- w = remaining weight of the inner partition (we need to reserve at least this amount)+    -- h = max height (decreasing)+    go !d !w !h (!a:as@(b:_)) +      | d < 0     = []+      | d == 0    = if a == 0 then [[]] else []+      | otherwise = [ this:rest | this <- [max 1 a .. min h (d-w)] , rest <- go (d-this) (w-b) this as ]+    +--------------------------------------------------------------------------------+-- * The Pieri rule++-- | The Pieri rule computes @s[lambda]*h[n]@ as a sum of @s[mu]@-s (each with coefficient 1).+--+-- See for example <http://en.wikipedia.org/wiki/Pieri's_formula>+--+-- | We assume here that @lambda@ is a partition (non-increasing sequence of /positive/ integers)! +_pieriRule :: [Int] -> Int -> [[Int]] +_pieriRule lambda n+    | n == 0     = [lambda]+    | n <  0     = [] +    | otherwise  = go n diffs dsums (lambda++[0]) +    where+      diffs = n : _diffSequence lambda                 -- maximum we can add to a given row+      dsums = reverse $ scanl1 (+) (reverse diffs)    -- partial sums of remaining total we can add+      go !k (d:ds) (p:ps@(q:_)) (l:ls) +        | k > p     = []+        | otherwise = [ h:tl | a <- [ max 0 (k-q) .. min d k ] , let h = l+a , tl <- go (k-a) ds ps ls ]+      go !k [d]    _      [l]    = if k <= d +                                     then if l+k>0 then [[l+k]] else [[]]+                                     else []+      go !k []     _      _      = if k==0 then [[]] else []++-- | The dual Pieri rule computes @s[lambda]*e[n]@ as a sum of @s[mu]@-s (each with coefficient 1)+_dualPieriRule :: [Int] -> Int -> [[Int]] +_dualPieriRule lam n = map _dualPartition $ _pieriRule (_dualPartition lam) n++--------------------------------------------------------------------------------

@@ -0,0 +1,202 @@++-- | Naive implementation of partitions of integers, encoded as list of @Int@-s.+--+-- Integer partitions are nonincreasing sequences of positive integers.+--+-- This is an internal module, you are not supposed to import it directly.+--+ ++{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, PatternSynonyms, ViewPatterns #-}+module Math.Combinat.Partitions.Integer.Naive where++--------------------------------------------------------------------------------++import Data.List +import Control.Monad ( liftM , replicateM )++-- import Data.Map (Map)+-- import qualified Data.Map as Map++import Math.Combinat.Classes+import Math.Combinat.ASCII as ASCII+import Math.Combinat.Numbers (factorial,binomial,multinomial)+import Math.Combinat.Helper++import Data.Array+import System.Random++import Math.Combinat.Partitions.Integer.IntList+import Math.Combinat.Partitions.Integer.Count ( countPartitions )++--------------------------------------------------------------------------------+-- * Type and basic stuff++-- | A partition of an integer. The additional invariant enforced here is that partitions +-- are monotone decreasing sequences of /positive/ integers. The @Ord@ instance is lexicographical.+newtype Partition = Partition [Int] deriving (Eq,Ord,Show,Read)++instance HasNumberOfParts Partition where+  numberOfParts (Partition p) = length p++---------------------------------------------------------------------------------++isEmptyPartition :: Partition -> Bool+isEmptyPartition (Partition p) = null p++emptyPartition :: Partition+emptyPartition = Partition []++instance CanBeEmpty Partition where+  empty   = emptyPartition+  isEmpty = isEmptyPartition++-- | The first element of the sequence.+partitionHeight :: Partition -> Int+partitionHeight (Partition part) = case part of+  (p:_) -> p+  []    -> 0+  +-- | The length of the sequence (that is, the number of parts).+partitionWidth :: Partition -> Int+partitionWidth (Partition part) = length part++instance HasHeight Partition where+  height = partitionHeight+ +instance HasWidth Partition where+  width = partitionWidth++heightWidth :: Partition -> (Int,Int)+heightWidth part = (height part, width part)++-- | The weight of the partition +--   (that is, the sum of the corresponding sequence).+partitionWeight :: Partition -> Int+partitionWeight (Partition part) = sum' part++instance HasWeight Partition where +  weight = partitionWeight++-- | The dual (or conjugate) partition.+dualPartition :: Partition -> Partition+dualPartition (Partition part) = Partition (_dualPartition part)++instance HasDuality Partition where +  dual = dualPartition++-- | Example:+--+-- > elements (toPartition [5,4,1]) ==+-- >   [ (1,1), (1,2), (1,3), (1,4), (1,5)+-- >   , (2,1), (2,2), (2,3), (2,4)+-- >   , (3,1)+-- >   ]+--+elements :: Partition -> [(Int,Int)]+elements (Partition part) = _elements part++--------------------------------------------------------------------------------+-- * Pattern synonyms ++-- | Pattern sysnonyms allows us to use existing code with minimal modifications+pattern Nil :: Partition+pattern Nil <- (isEmpty -> True) where+        Nil =  empty++pattern Cons :: Int -> Partition -> Partition+pattern Cons x xs  <- (unconsPartition -> Just (x,xs)) where+        Cons x (Partition xs) = Partition (x:xs)++-- | Simulated newtype constructor +pattern Partition_ :: [Int] -> Partition+pattern Partition_ xs = Partition xs++pattern Head :: Int -> Partition +pattern Head h <- (head . toDescList -> h)++pattern Tail :: Partition -> Partition+pattern Tail xs <- (Partition . tail . toDescList -> xs)++pattern Length :: Int -> Partition +pattern Length n <- (partitionWidth -> n)        + +---------------------------------------------------------------------------------+-- * Exponential form++-- | We convert a partition to exponential form.+-- @(i,e)@ mean @(i^e)@; for example @[(1,4),(2,3)]@ corresponds to @(1^4)(2^3) = [2,2,2,1,1,1,1]@. Another example:+--+-- > toExponentialForm (Partition [5,5,3,2,2,2,2,1,1]) == [(1,2),(2,4),(3,1),(5,2)]+--+toExponentialForm :: Partition -> [(Int,Int)]+toExponentialForm = _toExponentialForm . toDescList++fromExponentialForm :: [(Int,Int)] -> Partition+fromExponentialForm = Partition . _fromExponentialForm where++--------------------------------------------------------------------------------+-- * List-like operations++-- | From a sequence @[a1,a2,..,an]@ computes the sequence of differences+-- @[a1-a2,a2-a3,...,an-0]@+diffSequence :: Partition -> [Int]+diffSequence = go . toDescList where+  go (x:ys@(y:_)) = (x-y) : go ys +  go [x] = [x]+  go []  = []++unconsPartition :: Partition -> Maybe (Int,Partition)+unconsPartition (Partition xs) = case xs of+  (y:ys) -> Just (y, Partition ys)+  []     -> Nothing++toDescList :: Partition -> [Int]+toDescList (Partition xs) = xs++---------------------------------------------------------------------------------+-- * Dominance order ++-- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions+-- (this is partial ordering on the set of partitions of @n@).+--+-- See <http://en.wikipedia.org/wiki/Dominance_order>+--+dominates :: Partition -> Partition -> Bool+dominates (Partition qs) (Partition ps) +  = and $ zipWith (>=) (sums (qs ++ repeat 0)) (sums ps)+  where+    sums = scanl (+) 0++--------------------------------------------------------------------------------+-- * Containment partial ordering++-- | Returns @True@ of the first partition is a subpartition (that is, fit inside) of the second.+-- This includes equality+isSubPartitionOf :: Partition -> Partition -> Bool+isSubPartitionOf (Partition ps) (Partition qs) = and $ zipWith (<=) ps (qs ++ repeat 0)++-- | This is provided for convenience\/completeness only, as:+--+-- > isSuperPartitionOf q p == isSubPartitionOf p q+--+isSuperPartitionOf :: Partition -> Partition -> Bool+isSuperPartitionOf (Partition qs) (Partition ps) = and $ zipWith (<=) ps (qs ++ repeat 0)+    +--------------------------------------------------------------------------------+-- * The Pieri rule++-- | The Pieri rule computes @s[lambda]*h[n]@ as a sum of @s[mu]@-s (each with coefficient 1).+--+-- See for example <http://en.wikipedia.org/wiki/Pieri's_formula>+--+pieriRule :: Partition -> Int -> [Partition] +pieriRule (Partition lambda) n = map Partition (_pieriRule lambda n) where++-- | The dual Pieri rule computes @s[lambda]*e[n]@ as a sum of @s[mu]@-s (each with coefficient 1)+dualPieriRule :: Partition -> Int -> [Partition] +dualPieriRule lam n = map dualPartition $ pieriRule (dualPartition lam) n++--------------------------------------------------------------------------------++

@@ -82,6 +82,13 @@ instance HasDuality SkewPartition where   dual = dualSkewPartition +-- | See "partitionElements"+skewPartitionElements :: SkewPartition -> [(Int, Int)]+skewPartitionElements (SkewPartition abs) = concat+  [ [ (i,j) | j <- [a+1 .. a+b] ]+  | (i,(a,b)) <- zip [1..] abs+  ]+ -------------------------------------------------------------------------------- -- * Listing skew partitions @@ -109,6 +116,17 @@   where     outerWeight = innerWeight + skewWeight      innerWeight = weight inner ++--------------------------------------------------------------------------------+-- connected components++{-+connectedComponents :: SkewPartition -> [((Int,Int),SkewPartition)]+connectedComponents = error "connectedComponents: not implemented yet"++isConnectedSkewPartition :: SkewPartition -> Bool+isConnectedSkewPartition skewp = length (connectedComponents skewp) == 1+-}  -------------------------------------------------------------------------------- -- * ASCII

@@ -0,0 +1,364 @@++-- | Ribbons (also called border strips, skew hooks, skew rim hooks, etc...).+--+-- Ribbons are skew partitions that are 1) connected, 2) do not contain+-- 2x2 blocks. Intuitively, they are 1-box wide continuous strips on+-- the boundary.+--+-- An alternative definition that they are skew partitions whose projection+-- to the diagonal line is a continuous segment of width 1.++{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}+module Math.Combinat.Partitions.Skew.Ribbon where++--------------------------------------------------------------------------------++import Data.Array+import Data.List+import Data.Maybe++import qualified Data.Map as Map++import Math.Combinat.Sets+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Integer.IntList ( _diffSequence )+import Math.Combinat.Partitions.Skew+import Math.Combinat.Tableaux+import Math.Combinat.Tableaux.LittlewoodRichardson+import Math.Combinat.Tableaux.GelfandTsetlin+import Math.Combinat.Helper++--------------------------------------------------------------------------------+-- * Corners (TODO: move to Partitions - but we also want to refactor that)++-- | The coordinates of the outer corners +outerCorners :: Partition -> [(Int,Int)]+outerCorners = outerCornerBoxes++-- | The coordinates of the inner corners, including the two on the two coordinate+-- axes. For the partition @[5,4,1]@ the result should be @[(0,5),(1,4),(2,1),(3,0)]@+extendedInnerCorners:: Partition -> [(Int,Int)]+extendedInnerCorners (Partition_ ps) = (0, head ps') : catMaybes mbCorners where+  ps' = ps ++ [0]+  mbCorners = zipWith3 f [1..] (tail ps') (_diffSequence ps') +  f !y !x !k = if k > 0 then Just (y,x) else Nothing++-- | Sequence of all the (extended) corners+extendedCornerSequence :: Partition -> [(Int,Int)]+extendedCornerSequence (Partition_ ps) = {- if null ps then [(0,0)] else -} interleave inner outer where+  inner = (0, head ps') : [ (y,x) | (y,x,k) <- zip3 [1..] (tail ps') diff , k>0 ]+  outer =                 [ (y,x) | (y,x,k) <- zip3 [1..] ps'        diff , k>0 ]+  diff = _diffSequence ps'+  ps' = ps ++ [0]++-- | The inner corner /boxes/ of the partition. Coordinates are counted from 1+-- (cf.the 'elements' function), and the first coordinate is the row, the second+-- the column (in English notation).+--+-- For the partition @[5,4,1]@ the result should be @[(1,4),(2,1)]@+--+-- > innerCornerBoxes lambda == (tail $ init $ extendedInnerCorners lambda)+--+innerCornerBoxes :: Partition -> [(Int,Int)]+innerCornerBoxes (Partition_ ps) = +  case ps of+    []  -> []+    _   -> catMaybes mbCorners +  where+    mbCorners = zipWith3 f [1..] (tail ps) (_diffSequence ps) +    f !y !x !k = if k > 0 then Just (y,x) else Nothing++-- | The outer corner /boxes/ of the partition. Coordinates are counted from 1+-- (cf.the 'elements' function), and the first coordinate is the row, the second+-- the column (in English notation).+--+-- For the partition @[5,4,1]@ the result should be @[(1,5),(2,4),(3,1)]@+outerCornerBoxes :: Partition -> [(Int,Int)]+outerCornerBoxes (Partition_ ps) = catMaybes mbCorners where+  mbCorners = zipWith3 f [1..] ps (_diffSequence ps) +  f !y !x !k = if k > 0 then Just (y,x) else Nothing++-- | The outer and inner corner boxes interleaved, so together they form +-- the turning points of the full border strip+cornerBoxSequence :: Partition -> [(Int,Int)]+cornerBoxSequence (Partition_ ps) = if null ps then [] else interleave outer inner where+  inner = [ (y,x) | (y,x,k) <- zip3 [1..] tailps diff , k>0 ]+  outer = [ (y,x) | (y,x,k) <- zip3 [1..] ps     diff , k>0 ]+  diff = _diffSequence ps+  tailps = case ps of { [] -> [] ; _-> tail ps }++--------------------------------------------------------------------------------++-- | Naive (and very slow) implementation of @innerCornerBoxes@, for testing purposes+innerCornerBoxesNaive :: Partition -> [(Int,Int)]+innerCornerBoxesNaive part = filter f boxes where+  boxes = elements part+  f (y,x) =       elem (y+1,x  ) boxes+          &&      elem (y  ,x+1) boxes+          && not (elem (y+1,x+1) boxes)++-- | Naive (and very slow) implementation of @outerCornerBoxes@, for testing purposes+outerCornerBoxesNaive :: Partition -> [(Int,Int)]+outerCornerBoxesNaive part = filter f boxes where+  boxes = elements part+  f (y,x) =  not (elem (y+1,x  ) boxes)+          && not (elem (y  ,x+1) boxes)+          && not (elem (y+1,x+1) boxes)++--------------------------------------------------------------------------------+-- * Ribbon++-- | A skew partition is a a ribbon (or border strip) if and only if projected+-- to the diagonals the result is an interval.+isRibbon :: SkewPartition -> Bool+isRibbon skewp = go Nothing proj where+  proj = Map.toList +       $ Map.fromListWith (+) [ (x-y , 1) | (y,x) <- skewPartitionElements skewp ]+  go Nothing   []            = False+  go (Just _)  []            = True+  go Nothing   ((a,h):rest)  = (h == 1) &&               go (Just a) rest  +  go (Just b)  ((a,h):rest)  = (h == 1) && (a == b+1) && go (Just a) rest++{-+-- | Naive (and slow) reference implementation of "isRibbon"+isRibbonNaive :: SkewPartition -> Bool+isRibbonNaive skewp = isConnectedSkewPartition skewp && no2x2 where+  boxes = skewPartitionElements skewp+  no2x2 = and +    [ not ( elem (y+1,x  ) boxes &&             +            elem (y  ,x+1) boxes &&  +            elem (y+1,x+1) boxes )        -- no 2x2 blocks +    | (y,x) <- boxes +    ]+-}++toRibbon :: SkewPartition -> Maybe Ribbon+toRibbon skew = +  if not (isRibbon skew)+    then Nothing+    else Just ribbon +  where+    ribbon =  Ribbon+      { rbShape  = skew+      , rbLength = skewPartitionWeight skew+      , rbHeight = height+      , rbWidth  = width+      }+    elems  = skewPartitionElements skew+    height = (length $ group $ sort $ map fst elems) - 1    -- TODO: optimize these+    width  = (length $ group $ sort $ map snd elems) - 1++-- | Border strips (or ribbons) are defined to be skew partitions which are +-- connected and do not contain 2x2 blocks.+-- +-- The /length/ of a border strip is the number of boxes it contains,+-- and its /height/ is defined to be one less than the number of rows+-- (in English notation) it occupies. The /width/ is defined symmetrically to +-- be one less than the number of columns it occupies.+--+data Ribbon = Ribbon+  { rbShape  :: SkewPartition+  , rbLength :: Int+  , rbHeight :: Int+  , rbWidth  :: Int+  }+  deriving (Eq,Ord,Show)++--------------------------------------------------------------------------------+-- * Inner border strips++-- | Ribbons (or border strips) are defined to be skew partitions which are +-- connected and do not contain 2x2 blocks. This function returns the+-- border strips whose outer partition is the given one.+innerRibbons :: Partition -> [Ribbon]+innerRibbons part@(Partition ps) = if null ps then [] else strips where++  strips  = [ mkStrip i j +            | i<-[1..n] , _canStartStrip (annArr!i)+            , j<-[i..n] , _canEndStrip   (annArr!j)+            ]++  n       = length annList+  annList = annotatedInnerBorderStrip part+  annArr  = listArray (1,n) annList++  mkStrip !i1 !i2 = Ribbon shape len height width where+    ps'   = ps ++ [0]+    shape = SkewPartition [ (p-k,k) | (i,p,q) <- zip3 [1..] ps (tail ps') , let k = indent i p q ] +    indent !i !p !q +      | i <  y1    = 0+      | i >  y2    = 0+      | i == y2    = p - x2 + 1     -- the order is important here !!!+      | otherwise  = p - q  + 1     -- because of the case y1 == y2 == i++    len    = i2 - i1 + 1+    height = y2 - y1+    width  = x1 - x2+    BorderBox _ _ y1 x1 = annArr ! i1+    BorderBox _ _ y2 x2 = annArr ! i2++-- | Inner border strips (or ribbons) of the given length+innerRibbonsOfLength :: Partition -> Int -> [Ribbon]+innerRibbonsOfLength part@(Partition ps) givenLength = if null ps then [] else strips where++  strips  = [ mkStrip i j +            | i<-[1..n] , _canStartStrip (annArr!i)+            , j<-[i..n] , _canEndStrip   (annArr!j)+            , j-i+1 == givenLength+            ]++  n       = length annList+  annList = annotatedInnerBorderStrip part+  annArr  = listArray (1,n) annList++  mkStrip !i1 !i2 = Ribbon shape givenLength height width where+    ps'   = ps ++ [0]+    shape = SkewPartition [ (p-k,k) | (i,p,q) <- zip3 [1..] ps (tail ps') , let k = indent i p q ] +    indent !i !p !q +      | i <  y1    = 0+      | i >  y2    = 0+      | i == y2    = p - x2 + 1     -- the order is important here !!!+      | otherwise  = p - q  + 1     -- because of the case y1 == y2 == i++    height = y2 - y1+    width  = x1 - x2+    BorderBox _ _ y1 x1 = annArr ! i1+    BorderBox _ _ y2 x2 = annArr ! i2+++--------------------------------------------------------------------------------+-- * Outer border strips++-- | Hooks of length @n@ (TODO: move to the partition module)+listHooks :: Int -> [Partition]+listHooks 0 = []+listHooks 1 = [ Partition [1] ]+listHooks n = [ Partition (k : replicate (n-k) 1) | k<-[1..n] ]++-- | Outer border strips (or ribbons) of the given length+outerRibbonsOfLength :: Partition -> Int -> [Ribbon]+outerRibbonsOfLength part@(Partition ps) givenLength = result where++  result = if null ps +    then [ Ribbon shape givenLength ht wd+         | p <- listHooks givenLength+         , let shape = mkSkewPartition (p,part)+         , let ht = partitionWidth  p - 1        -- pretty inconsistent names here :(((+         , let wd = partitionHeight p - 1+         ]+    else strips ++  strips  = [ mkStrip i j +            | i<-[1..n] , _canStartStrip (annArr!i)+            , j<-[i..n] , _canEndStrip   (annArr!j)+            , j-i+1 == givenLength+            ]+ +  ysize = partitionWidth  part+  xsize = partitionHeight part+ +  annList  =  [ BorderBox True False 1 x | x <- reverse [xsize+2 .. xsize+givenLength ] ]+           ++ annList0 +           ++ [ BorderBox False True y 1 | y <-         [ysize+2 .. ysize+givenLength ] ]+ +  n        = length annList+  annList0 = annotatedOuterBorderStrip part+  annArr   = listArray (1,n) annList++  mkStrip !i1 !i2 = Ribbon shape len height width where+    ps'   = (-666) : ps ++ replicate (givenLength) 0+    shape = SkewPartition [ (p,k) | (i,p,q) <- zip3 [1..max ysize y2] (tail ps') ps' , let k = indent i p q ] +    indent !i !p !q +      | i <  y1    = 0+      | i >  y2    = 0+      | i == y1    = x1 - p    -- the order is important here !!!+--      | i == y2    = x2 - p     +      | otherwise  = q - p  + 1   ++    len    = i2 - i1 + 1+    height = y2 - y1+    width  = x1 - x2+    BorderBox _ _ y1 x1 = annArr ! i1+    BorderBox _ _ y2 x2 = annArr ! i2++--------------------------------------------------------------------------------+-- * Naive implementations (for testing)++-- | Naive (and slow) implementation listing all inner border strips+innerRibbonsNaive :: Partition -> [Ribbon]+innerRibbonsNaive outer = list where+  list = [ Ribbon skew (len skew) (ht skew) (wt skew)+         | skew <- allSkewPartitionsWithOuterShape outer+         , isRibbon skew+         ]+  len skew = length (skewPartitionElements skew)+  ht  skew = (length $ group $ sort $ map fst $ skewPartitionElements skew) - 1+  wt  skew = (length $ group $ sort $ map snd $ skewPartitionElements skew) - 1+++-- | Naive (and slow) implementation listing all inner border strips of the given length+innerRibbonsOfLengthNaive :: Partition -> Int -> [Ribbon]+innerRibbonsOfLengthNaive outer givenLength = list where+  pweight = partitionWeight outer+  list = [ Ribbon skew (len skew) (ht skew) (wt skew)+         | skew <- skewPartitionsWithOuterShape outer givenLength+         , isRibbon skew+         ]+  len skew = length (skewPartitionElements skew)+  ht  skew = (length $ group $ sort $ map fst $ skewPartitionElements skew) - 1+  wt  skew = (length $ group $ sort $ map snd $ skewPartitionElements skew) - 1++-- | Naive (and slow) implementation listing all outer border strips of the given length+outerRibbonsOfLengthNaive :: Partition -> Int -> [Ribbon]+outerRibbonsOfLengthNaive inner givenLength = list where+  pweight = partitionWeight inner+  list = [ Ribbon skew (len skew) (ht skew) (wt skew)+         | skew <- skewPartitionsWithInnerShape inner givenLength+         , isRibbon skew+         ]+  len skew = length (skewPartitionElements skew)+  ht  skew = (length $ group $ sort $ map fst $ skewPartitionElements skew) - 1+  wt  skew = (length $ group $ sort $ map snd $ skewPartitionElements skew) - 1++--------------------------------------------------------------------------------+-- * Annotated borders++-- | A box on the border of a partition+data BorderBox = BorderBox+  { _canStartStrip :: !Bool+  , _canEndStrip   :: !Bool+  , _yCoord :: !Int+  , _xCoord :: !Int+  }+  deriving Show+ +-- | The boxes of the full inner border strip, annotated with whether a border strip +-- can start or end there.+annotatedInnerBorderStrip :: Partition -> [BorderBox]+annotatedInnerBorderStrip partition = if isEmptyPartition partition then [] else list where+  list    = goVert (head corners) (tail corners) +  corners = extendedCornerSequence partition  ++  goVert  (y1,x ) ((y2,_ ):rest) = [ BorderBox True (y==y2) y x | y<-[y1+1..y2] ] ++ goHoriz (y2,x) rest+  goVert  _       []             = [] ++  goHoriz (y ,x1) ((_, x2):rest) = case rest of+    [] -> [ BorderBox False True    y x | x<-[x1-1,x1-2..x2+1] ]+    _  -> [ BorderBox False (x/=x2) y x | x<-[x1-1,x1-2..x2  ] ] ++ goVert (y,x2) rest++-- | The boxes of the full outer border strip, annotated with whether a border strip +-- can start or end there.+annotatedOuterBorderStrip :: Partition -> [BorderBox]+annotatedOuterBorderStrip partition = if isEmptyPartition partition then [] else list where+  list    = goVert (head corners) (tail corners) +  corners = extendedCornerSequence partition  ++  goVert  (y1,x ) ((y2,_ ):rest) = [ BorderBox (y==y1) (y/=y2) (y+1) (x+1) | y<-[y1..y2] ] ++ goHoriz (y2,x) rest+  goVert  _       []             = [] ++  goHoriz (y ,x1) ((_, x2):rest) = case rest of+    [] -> [ BorderBox True (x==0) (y+1) (x+1) | x<-[x1-1,x1-2..x2  ] ]+    _  -> [ BorderBox True False  (y+1) (x+1) | x<-[x1-1,x1-2..x2+1] ] ++ goVert (y,x2) rest+++--------------------------------------------------------------------------------

@@ -28,6 +28,7 @@   , permutationToDisjointCycles   , disjointCyclesToPermutation   , numberOfCycles+  , concatPermutations     -- * Queries   , isIdentityPermutation   , isReversePermutation@@ -63,6 +64,9 @@   , permuteList   , permuteLeft , permuteRight   , permuteLeftList , permuteRightList+    -- * Sorting+  , sortingPermutationAsc +  , sortingPermutationDesc     -- * ASCII drawing   , asciiPermutation   , asciiDisjointCycles@@ -224,6 +228,17 @@ isIdentityPermutation (Permutation ar) = (elems ar == [1..n]) where   (1,n) = bounds ar +-- | Given a permutation of @n@ and a permutation of @m@, we return+-- a permutation of @n+m@ resulting by putting them next to each other.+-- This should satisfy+--+-- > permuteList p1 xs ++ permuteList p2 ys == permuteList (concatPermutations p1 p2) (xs++ys)+--+concatPermutations :: Permutation -> Permutation -> Permutation +concatPermutations perm1 perm2 = toPermutationUnsafe list where+  n    = permutationSize perm1+  list = fromPermutation perm1 ++ map (+n) (fromPermutation perm2)+ -------------------------------------------------------------------------------- -- * ASCII @@ -738,6 +753,32 @@ permuteLeftList perm xs = elems $ permuteLeft perm $ arr where   arr = listArray (1,n) xs :: Array Int a   n   = permutationSize perm++--------------------------------------------------------------------------------++-- | Given a list of things, we return a permutation which sorts them into+-- ascending order, that is:+--+-- > permuteList (sortingPermutationAsc xs) xs == sort xs+--+-- Note: if the things are not unique, then the sorting permutations is not+-- unique either; we just return one of them.+--+sortingPermutationAsc :: Ord a => [a] -> Permutation+sortingPermutationAsc xs = toPermutation (map fst sorted) where+  sorted = sortBy (comparing snd) $ zip [1..] xs++-- | Given a list of things, we return a permutation which sorts them into+-- descending order, that is:+--+-- > permuteList (sortingPermutationDesc xs) xs == reverse (sort xs)+--+-- Note: if the things are not unique, then the sorting permutations is not+-- unique either; we just return one of them.+--+sortingPermutationDesc :: Ord a => [a] -> Permutation+sortingPermutationDesc xs = toPermutation (map fst sorted) where+  sorted = sortBy (reverseComparing snd) $ zip [1..] xs  -------------------------------------------------------------------------------- -- * Permutations of distinct elements

@@ -0,0 +1,150 @@++-- | Venn diagrams. See <https://en.wikipedia.org/wiki/Venn_diagram>+--+-- TODO: write a more efficient implementation (for example an array of size @2^n@)+--++{-# LANGUAGE BangPatterns #-}+module Math.Combinat.Sets.VennDiagrams where++--------------------------------------------------------------------------------++import Data.List++import GHC.TypeLits+import Data.Proxy++import qualified Data.Map as Map+import Data.Map (Map)++import Math.Combinat.Compositions+import Math.Combinat.ASCII++--------------------------------------------------------------------------------++-- | Venn diagrams of @n@ sets. Each possible zone is annotated with a value+-- of type @a@. A typical use case is to annotate with the cardinality of the+-- given zone.+--+-- Internally this is representated by a map from @[Bool]@, where @True@ means element +-- of the set, @False@ means not.+--+-- TODO: write a more efficient implementation (for example an array of size 2^n)+newtype VennDiagram a = VennDiagram { vennTable :: Map [Bool] a } deriving (Eq,Ord,Show)++-- | How many sets are in the Venn diagram+vennDiagramNumberOfSets :: VennDiagram a -> Int+vennDiagramNumberOfSets (VennDiagram table) = length $ fst $ Map.findMin table++-- | How many zones are in the Venn diagram+--+-- > vennDiagramNumberOfZones v == 2 ^ (vennDiagramNumberOfSets v)+--+vennDiagramNumberOfZones :: VennDiagram a -> Int+vennDiagramNumberOfZones venn = 2 ^ (vennDiagramNumberOfSets venn)++-- | How many /nonempty/ zones are in the Venn diagram+vennDiagramNumberOfNonemptyZones :: VennDiagram Int -> Int+vennDiagramNumberOfNonemptyZones (VennDiagram table) = length $ filter (/=0) $ Map.elems table++unsafeMakeVennDiagram :: [([Bool],a)] -> VennDiagram a+unsafeMakeVennDiagram = VennDiagram . Map.fromList++-- | We call venn diagram trivial if all the intersection zones has zero cardinality+-- (that is, the original sets are all disjoint)+isTrivialVennDiagram :: VennDiagram Int -> Bool+isTrivialVennDiagram (VennDiagram table) = and [ c == 0 | (bs,c) <- Map.toList table , isIntersection bs ] where+  isIntersection bs = case filter id bs of+    []  -> False+    [_] -> False+    _   -> True++printVennDiagram :: Show a => VennDiagram a -> IO ()+printVennDiagram = putStrLn . prettyVennDiagram++prettyVennDiagram :: Show a => VennDiagram a -> String+prettyVennDiagram = unlines . asciiLines . asciiVennDiagram++asciiVennDiagram :: Show a => VennDiagram a -> ASCII+asciiVennDiagram (VennDiagram table) = asciiFromLines $ map f (Map.toList table) where+  f (bs,a) = "{" ++ extendTo (length bs) [ if b then z else ' ' | (b,z) <- zip bs abc ] ++ "} -> " ++ show a+  extendTo k str = str ++ replicate (k - length str) ' '+  abc = ['A'..'Z']++instance Show a => DrawASCII (VennDiagram a) where+  ascii = asciiVennDiagram++-- | Given a Venn diagram of cardinalities, we compute the cardinalities of the+-- original sets (note: this is slow!)+vennDiagramSetCardinalities :: VennDiagram Int -> [Int]+vennDiagramSetCardinalities (VennDiagram table) = go n list where+  list = Map.toList table+  n = length $ fst $ head list+  go :: Int -> [([Bool],Int)] -> [Int]+  go !0 _  = []+  go !k xs = this : go (k-1) (map xtail xs) where+    this = foldl' (+) 0 [ c | ((True:_) , c) <- xs ]+  xtail (bs,c) = (tail bs,c)++--------------------------------------------------------------------------------++-- | Given the cardinalities of some finite sets, we list all possible+-- Venn diagrams.+--+-- Note: we don't include the empty zone in the tables, because it's always empty.+--+-- Remark: if each sets is a singleton set, we get back set partitions:+--+-- > > [ length $ enumerateVennDiagrams $ replicate k 1 | k<-[1..8] ]+-- > [1,2,5,15,52,203,877,4140]+-- >+-- > > [ countSetPartitions k | k<-[1..8] ]+-- > [1,2,5,15,52,203,877,4140]+--+-- Maybe this could be called multiset-partitions?+--+-- Example:+--+-- > autoTabulate RowMajor (Right 6) $ map ascii $ enumerateVennDiagrams [2,3,3]+--+enumerateVennDiagrams :: [Int] -> [VennDiagram Int]+enumerateVennDiagrams dims = +  case dims of+    []     -> []+    [d]    -> venns1 d+    (d:ds) -> concatMap (worker (length ds) d) $ enumerateVennDiagrams ds+  where++    worker !n !d (VennDiagram table) = result where++      list   = Map.toList table+      falses = replicate n False++      comps k = compositions' (map snd list) k+      result = +        [ unsafeMakeVennDiagram $ +            [ (False:tfs    , m-c) | ((tfs,m),c) <- zip list comp ] +++            [ (True :tfs    ,   c) | ((tfs,m),c) <- zip list comp ] +++            [ (True :falses , d-k) ]+        | k <- [0..d]+        , comp <- comps k+        ]++    venns1 :: Int -> [VennDiagram Int]+    venns1 p = [ theVenn ] where +      theVenn = unsafeMakeVennDiagram [ ([True],p) ] ++--------------------------------------------------------------------------------++{-++-- | for testing only+venns2 :: Int -> Int -> [Venn Int]+venns2 p q = +  [ mkVenn [ ([t,f],p-k) , ([f,t],q-k) , ([t,t],k) ]+  | k <- [0..min p q] +  ]+  where+    t = True+    f = False+-}

@@ -1,12 +1,19 @@  -- | Signs -{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP, BangPatterns #-} module Math.Combinat.Sign where  --------------------------------------------------------------------------------  import Data.Monoid++-- Semigroup became a superclass of Monoid+#if MIN_VERSION_base(4,11,0)     +import Data.Foldable+import Data.Semigroup+#endif+ import System.Random  --------------------------------------------------------------------------------@@ -16,10 +23,29 @@   | Minus   deriving (Eq,Ord,Show,Read) +--------------------------------------------------------------------------------++-- Semigroup became a superclass of Monoid+#if MIN_VERSION_base(4,11,0)        ++instance Semigroup Sign where+  (<>)    = mulSign+  sconcat = foldl1 mulSign+ instance Monoid Sign where   mempty  = Plus+  mconcat = productOfSigns++#else++instance Monoid Sign where+  mempty  = Plus   mappend = mulSign   mconcat = productOfSigns++#endif++--------------------------------------------------------------------------------  instance Random Sign where   random        g = let (b,g') = random g in (if b    then Plus else Minus, g')

@@ -31,8 +31,9 @@ import Data.List  import Math.Combinat.Classes-import Math.Combinat.Numbers (factorial,binomial)-import Math.Combinat.Partitions+import Math.Combinat.Numbers ( factorial , binomial )+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Integer.IntList ( _dualPartition ) import Math.Combinat.ASCII import Math.Combinat.Helper 

@@ -16,6 +16,7 @@  import Math.Combinat.Classes import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Integer.IntList ( _diffSequence ) import Math.Combinat.Partitions.Skew import Math.Combinat.Tableaux import Math.Combinat.ASCII@@ -118,7 +119,7 @@    stuff = worker as bs ds (repeat 1)    (as,bs) = unzip abs-  ds = diffSequence as+  ds = _diffSequence as      -- | @worker inner outerMinusInner innerdiffs lowerbound   worker :: [Int] -> [Int] -> [Int] -> [Int] -> [[(Int,[Int])]]

@@ -0,0 +1,24 @@++#include <stdint.h>++// -----------------------------------------------------------------------------++uint64_t c_dual_nibble(uint64_t word)+{+  uint64_t n    = (word & 15);    // length+  uint64_t h    = (word >> 60);   // height+  uint64_t dual = h;              // length of dual = height of original+  uint64_t w    = word - n;       // zero out the low nibble++  uint64_t o = 60;+  for(uint64_t i=0; i<n; i++)+  { uint64_t k = ( (w >> (64-4*(n-i))) +                 - (w >> (60-4*(n-i))) ) & 15 ;   // diff+    for(uint64_t j=0;j<k;j++) { dual |= (n-i) << o ; o -= 4 ; } +  }++  return dual;+}++// -----------------------------------------------------------------------------+

@@ -1,5 +1,5 @@ Name:                combinat-Version:             0.2.8.2+Version:             0.2.9.0 Synopsis:            Generate and manipulate various combinatorial objects. Description:         A collection of functions to generate, count, manipulate                      and visualize all kinds of combinatorial objects like @@ -8,13 +8,13 @@ License:             BSD3 License-file:        LICENSE Author:              Balazs Komuves-Copyright:           (c) 2008-2016 Balazs Komuves+Copyright:           (c) 2008-2018 Balazs Komuves Maintainer:          bkomuves (plus) hackage (at) gmail (dot) com-Homepage:            http://code.haskell.org/~bkomuves/+Homepage:            http://moire.be/haskell/ Stability:           Experimental Category:            Math-Tested-With:         GHC == 7.10.3-Cabal-Version:       >= 1.18+Tested-With:         GHC == 8.0.2+Cabal-Version:       1.24 Build-Type:          Simple  extra-doc-files:     svg/*.svg @@ -22,7 +22,12 @@ extra-source-files:  svg/*.svg                      svg/src/gen_figures.hs                      - +source-repository head+  type:                darcs +  location:            http://moire.be/haskell/projects/combinat/++--------------------------------------------------------------------------------+ Library    Build-Depends:       base >= 4 && < 5, array >= 0.5, containers, random, transformers@@ -30,10 +35,13 @@   Exposed-Modules:     Math.Combinat                        Math.Combinat.Classes                        Math.Combinat.Numbers-                       Math.Combinat.Numbers.Series+                       Math.Combinat.Numbers.Integers+                       Math.Combinat.Numbers.Sequences                        Math.Combinat.Numbers.Primes+                       Math.Combinat.Numbers.Series                        Math.Combinat.Sign                        Math.Combinat.Sets+                       Math.Combinat.Sets.VennDiagrams                        Math.Combinat.Tuples                         Math.Combinat.Compositions                        Math.Combinat.Groups.Thompson.F@@ -42,7 +50,12 @@                        Math.Combinat.Groups.Braid.NF                        Math.Combinat.Partitions                        Math.Combinat.Partitions.Integer+                       Math.Combinat.Partitions.Integer.Count+                       Math.Combinat.Partitions.Integer.Compact+                       Math.Combinat.Partitions.Integer.Naive+                       Math.Combinat.Partitions.Integer.IntList                        Math.Combinat.Partitions.Skew+                       Math.Combinat.Partitions.Skew.Ribbon                        Math.Combinat.Partitions.Set                        Math.Combinat.Partitions.NonCrossing                        Math.Combinat.Partitions.Plane@@ -72,8 +85,12 @@    Hs-Source-Dirs:      . +  C-Sources:           cbits/c_compact_partition.c+   ghc-options:         -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports+     +--------------------------------------------------------------------------------  test-suite combinat-tests                       @@ -89,12 +106,19 @@                        Tests.SkewTableaux                        Tests.Thompson                        Tests.Partitions.Integer+                       Tests.Partitions.Compact                        Tests.Partitions.Skew+                       Tests.Partitions.Ribbon -  build-depends:       base >= 4 && < 5, array >= 0.5, containers, random, transformers,+  build-depends:       base >= 4 && < 5, array >= 0.5, containers >= 0.5, random, transformers,                        combinat,-                       QuickCheck >= 2, test-framework, test-framework-quickcheck2+                       test-framework, +                       test-framework-quickcheck2, QuickCheck >= 2,+                       tasty, tasty-quickcheck, tasty-hunit    Default-Language:    Haskell2010   Default-Extensions:  CPP, BangPatterns++--------------------------------------------------------------------------------+ 

@@ -9,6 +9,7 @@ import Tests.Permutations       ( testgroup_Permutations      ) import Tests.Partitions.Integer ( testgroup_IntegerPartitions ) import Tests.Partitions.Skew    ( testgroup_SkewPartitions    )+import Tests.Partitions.Ribbon  ( testgroup_Ribbon      ) import Tests.Braid              ( testgroup_Braid                                  , testgroup_Braid_NF          ) import Tests.Series             ( testgroup_PowerSeries       )@@ -24,9 +25,11 @@ tests :: [Test] tests =    [ testgroup_Permutations+  , testgroup_PowerSeries     , testGroup "Partitions"        [ testgroup_IntegerPartitions       , testgroup_SkewPartitions+      , testgroup_Ribbon       ]   , testgroup_SkewTableaux   , testgroup_ThompsonF@@ -35,7 +38,6 @@       [ testgroup_Braid        , testgroup_Braid_NF        ]-  , testgroup_PowerSeries     ]  --------------------------------------------------------------------------------

@@ -29,6 +29,10 @@ myMkGen fun = MkGen (\r _ -> let (x,_) = fun r in x)
 
 -- | Generates a random element.
+myMkGen' :: (a -> b) -> (forall g. RandomGen g => g -> (a,g)) -> Gen b
+myMkGen' h fun = MkGen (\r _ -> let (x,_) = fun r in h x)
+
+-- | Generates a random element.
 myMkSizedGen :: (forall g. RandomGen g => Int -> g -> (a,g)) -> Gen a
 myMkSizedGen fun = MkGen (\r siz -> let (x,_) = fun siz r in x)
 

@@ -0,0 +1,390 @@++module Tests.Partitions.Compact where ++--------------------------------------------------------------------------------++import Data.List hiding ( uncons )+import Data.Ord++import Test.Tasty+import Test.Tasty.HUnit      as U+import Test.Tasty.QuickCheck as Q ++import Math.Combinat.Partitions.Integer.Compact+import Math.Combinat.Partitions.Integer as P+-- import qualified Math.Combinat.Parititions.Integer.IntList as P ++--------------------------------------------------------------------------------++main = defaultMain tests++tests :: TestTree+tests = testGroup "Tests" [properties, unitTests]++--------------------------------------------------------------------------------++unitTests :: TestTree+unitTests = testGroup "Unit tests"+  [ testCase "toList . fromList == id /1" $ allTrue [ xs == toList (fromDescList xs) | xs <- _testPartitions ]+  , testCase "toList . fromList == id /2" $ allTrue [ xs == toList (fromDescList xs) | xs <- _allparts 18    ]+  , testCase "toAscList . fromList == reverse /1" $ allTrue [ reverse xs == toAscList (fromDescList xs) | xs <- _testPartitions ]+  , testCase "toAscList . fromList == reverse /2" $ allTrue [ reverse xs == toAscList (fromDescList xs) | xs <- _allparts 18    ]+  , testCase "fromList . toList == id"    $ allTrue [ p ==  fromDescList (toList p ) | p  <- testPartitions  ]+  , testCase "singleton"               $ allTrue [ toList (singleton n) == [n] | n <- [1..300] ]+  , testCase "singleton 0 is empty"    $ allTrue [ toList (singleton 0) == [] ]+  , testCase "uncons empty"            $ allTrue [ uncons empty == Nothing ]+  , testCase "uncons singleton"        $ allTrue [ uncons (singleton x) == Just (x,empty) | x <- [1..300] ]+  , testCase "cons/snoc 0 empty"       $ allTrue [ (cons 0 empty) == empty , (snoc empty 0) == empty ] +  , testCase "cons empty"              $ allTrue [ toList (cons n empty) == [n] | n <- [1..300] ]+  , testCase "snoc empty"              $ allTrue [ toList (snoc empty n) == [n] | n <- [1..300] ]+  , testCase "width/height of empty"   $ allTrue [ width empty == 0 , height empty == 0 ]+  , testCase "width of all "           $ allTrue [ length   xs == width  p | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "height of all"           $ allTrue [ safeHead xs == height p | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "(width,height)"          $ allTrue [ widthHeight p == (width p, height p) | p <- testPartitions ]+  , testCase "tail of all"             $ allTrue [ safeTail xs == toList (partitionTail p) | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "toList using uncons"     $ allTrue [ xs == toListViaUncons p                      | xs <- _testPartitionsSmall , let p = fromDescList xs ] +  , testCase "fromList using cons"     $ allTrue [ xs == toList (fromListViaCons (DescList xs)) | xs <- _testPartitionsSmall ]+  , testCase "fromList using snoc"     $ allTrue [ xs == toList (fromListViaSnoc (DescList xs)) | xs <- _testPartitionsSmall ]+  , testCase "reflexivity"             $ allTrue [ (fromDescList xs == p) | xs <- _testPartitions , let p = fromDescList xs]+  , testCase "snoc1"                   $ allTrue [ toList (snoc     p 1)  == xs ++ [1]           | xs <- _testPartitions , let p = fromDescList xs]+  , testCase "snocN/2..5"              $ allTrue [ toList (snocN n (p,1)) == xs ++ replicate n 1 | xs <- _testPartitions , let p = fromDescList xs , n <- [2..5] ]+  , testCase "compare/staircase"       $ allTrue [ compare xs1 xs2 == cmp p1 p2 | n1<-[0..100] , n2<-[0..100], let xs1 = _staircase n1 , let xs2 = _staircase n2 , let p1 = staircase n1 , let p2 = staircase n2 ]+  , testCase "compare/slope"           $ allTrue [ compare xs1 xs2 == cmp p1 p2 | n1<-[0..100] , n2<-[0..100], let xs1 = _slope n1 , let xs2 = _slope n2 , let p1 = slope n1 , let p2 = slope n2 ]+  , testCase "compare/steep"           $ allTrue [ compare xs1 xs2 == cmp p1 p2 | n1<-[0..100] , n2<-[0..100], let xs1 = _steep n1 , let xs2 = _steep n2 , let p1 = steep n1 , let p2 = steep n2 ]+  , testCase "compare/small"           $ allTrue [ compare xs1 xs2 == cmp p1 p2 | xs1 <- _allparts 12 , xs2 <- _allparts 12 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "compare/15"              $ allTrue [ compare xs1 xs2 == cmp p1 p2 | xs1 <- _parts    15 , xs2 <- _parts    15 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "compare/16"              $ allTrue [ compare xs1 xs2 == cmp p1 p2 | xs1 <- _parts    16 , xs2 <- _parts    16 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "compare/17"              $ allTrue [ compare xs1 xs2 == cmp p1 p2 | xs1 <- _parts    17 , xs2 <- _parts    17 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "ineq/small"              $ allTrue [ ineqTestPartition p1 p2 | xs1 <- _allparts 13 , xs2 <- _allparts 13 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "consN 15 / staircase"    $ allTrue [ (replicate k 15  ++ xs) == toList (consN k (15,p)) | n<-[0..15] , let xs = _staircase n , let p = fromDescList xs , k <- [1..80] ]+  , testCase "consN 15 / slope"        $ allTrue [ (replicate k 15  ++ xs) == toList (consN k (15,p)) | n<-[0..15] , let xs = _slope n , let p = fromDescList xs , k <- [1..80] ]+  , testCase "consN 15 / steep"        $ allTrue [ (replicate k 15  ++ xs) == toList (consN k (15,p)) | n<-[0..15] , let xs = _steep n , let p = fromDescList xs , k <- [1..40] ]+  , testCase "consN 16 / staircase"    $ allTrue [ (replicate k  16 ++ xs) == toList (consN k (16,p)) | n<-[0..16] , let xs = _staircase n , let p = fromDescList xs , k <- [1..80] ]+  , testCase "consN 16 / slope"        $ allTrue [ (replicate k  16 ++ xs) == toList (consN k (16,p)) | n<-[0..16] , let xs = _slope n , let p = fromDescList xs , k <- [1..80] ]+  , testCase "consN 16 / steep"        $ allTrue [ (replicate k  16 ++ xs) == toList (consN k (16,p)) | n<-[0..16] , let xs = _steep n , let p = fromDescList xs , k <- [1..40] ]+  , testCase "consN 256 / staircase"   $ allTrue [ (replicate k 256 ++ xs) == toList (consN k (256,p)) | n<-[0..40] , let xs = _staircase n , let p = fromDescList xs , k <- [1..40] ]+  , testCase "consN 256 / slope"       $ allTrue [ (replicate k 256 ++ xs) == toList (consN k (256,p)) | n<-[0..40] , let xs = _slope n , let p = fromDescList xs , k <- [1..35] ]+  , testCase "consN 256 / steep"       $ allTrue [ (replicate k 256 ++ xs) == toList (consN k (256,p)) | n<-[0..40] , let xs = _steep n , let p = fromDescList xs , k <- [1..35] ]+  , testCase "diffSequence"            $ allTrue [ diffSequence p == refDiffSeq xs | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "reverseDiffSequence"     $ allTrue [ reverseDiffSequence p == reverse (refDiffSeq xs) | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "dual . dual == id"       $ allTrue [ dualPartition (dualPartition p) == p | p <- testPartitions ]+  , testCase "dual == reference impl." $ allTrue [ toList (dualPartition p) == P._dualPartition xs | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "toExponentialForm"       $ allTrue [ toExponentialForm p == P._toExponentialForm xs | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "fromExponentialForm"     $ allTrue [ toList (fromExponentialForm ef) == xs | xs <- _testPartitions , let p = fromDescList xs , let ef = P._toExponentialForm xs ]+  , testCase "to / from expo. form"    $ allTrue [ toList (fromExponentialForm $ toExponentialForm p) == xs | xs <- _testPartitions , let p = fromDescList xs ]+  , testCase "isSubPartitionOf/small"  $ allTrue [ P._isSubPartitionOf xs1 xs2 == isSubPartitionOf p1 p2 | xs1 <- _allparts 12 , xs2 <- _allparts 12 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "isSubPartitionOf/15"     $ allTrue [ P._isSubPartitionOf xs1 xs2 == isSubPartitionOf p1 p2 | xs1 <- _parts    15 , xs2 <- _parts    15 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "isSubPartitionOf/16"     $ allTrue [ P._isSubPartitionOf xs1 xs2 == isSubPartitionOf p1 p2 | xs1 <- _parts    16 , xs2 <- _parts    16 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "isSubPartitionOf/17"     $ allTrue [ P._isSubPartitionOf xs1 xs2 == isSubPartitionOf p1 p2 | xs1 <- _parts    17 , xs2 <- _parts    17 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "dominates/small"         $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _allparts 12 , xs2 <- _allparts 12 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "dominates/15"            $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _parts    15 , xs2 <- _parts    15 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "dominates/16"            $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _parts    16 , xs2 <- _parts    16 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "dominates/17"            $ allTrue [ P._dominates xs1 xs2 == dominates p1 p2 | xs1 <- _parts    17 , xs2 <- _parts    17 , let p1 = fromDescList xs1 , let p2 = fromDescList xs2 ]+  , testCase "pieriRuleSingleBox"      $ allTrue [ pieriRuleSingleBox p =%%= map fromDescList (P._pieriRule xs 1) | xs <- _testPartitions      , let p = fromDescList xs ] +  , testCase "pieriRule"               $ allTrue [ pieriRule p k        =%%= map fromDescList (P._pieriRule xs k) | xs <- every10th _testPartitionsSmall , let p = fromDescList xs , k <- [1..2] ] +  ]++--------------------------------------------------------------------------------++properties :: TestTree+properties = localOption (QuickCheckTests 1000)  -- 200+           $ testGroup "Properties" +  [ prop "toList . fromList == id" $ \(DescList xs) -> (toList (fromDescList xs) == xs)+  , prop "toAscList . fromList == reverse" $ \(DescList xs) -> (toAscList (fromDescList xs) == reverse xs)+  , prop "snoc1/list"              $ \(DescList xs) -> toList (snoc (fromDescList xs) 1) == xs ++ [1]+  , prop "snocN/list"              $ \(DescList xs) (SmallN n) -> toList (snocN n (fromDescList xs,1)) == xs ++ replicate n 1+  , prop "fromList . toList == id" $ \p    -> (fromDescList (toList p ) == p )+  , prop "compare"                 $ \p q  -> cmp p q == compare (toList p) (toList q)+  , prop "uncons"                  $ \p    -> (unconsTest p) == unconsList (toList p)+  , prop "width"                   $ \p    -> width  p == length   (toList p)+  , prop "height"                  $ \p    -> height p == safeHead (toList p)+  , prop "(width,height)"          $ \p    -> let xs = toList p in widthHeight p == (length xs, safeHead xs)+  , prop "snoc1"                   $ \p    -> toList (snoc p 1) == toList p ++ [1]+  , prop "snocN"                   $ \p (SmallN n) -> toList (snocN n (p,1)) == toList p ++ replicate n 1+  , prop "cons/head"               $ \p            -> let a = max 1 (height p) in toList (cons a p) == a : toList p +  , prop "cons/head+1"             $ \p            -> let a = height p + 1     in toList (cons a p) == a : toList p +  , prop "cons/head+k"             $ \p (SmallN k) -> let a = height p + k     in toList (cons a p) == a : toList p +  , prop "consN/head"              $ \p (SmallN n) -> let a = max 1 (height p) in toList (consN n (a,p)) == replicate n a ++ toList p +  , prop "consN/head+1"            $ \p (SmallN n) -> let a = height p + 1     in toList (consN n (a,p)) == replicate n a ++ toList p +  , prop "consN/head+k"            $ \p (SmallN n) (SmallN k) -> let a = height p + k  in toList (consN n (a,p)) == replicate n a ++ toList p +  , prop "tailN"                   $ \p (SmallN n) -> toList (tailN n p) == drop n (toList p)+  , prop "isEmpty of tail length  "  $ \p ->                      isEmpty (tailN (width p    ) p)  +  , prop "isEmpty of tail length-1"  $ \p -> width p == 0 || not (isEmpty (tailN (width p - 1) p))+  , prop "toList using uncons"     $ \p  -> toList p == toListViaUncons p+  , prop "fromList using cons"     $ \dlist@(DescList xs) -> fromDescList xs == fromListViaCons dlist+  , prop "fromList using snoc"     $ \dlist@(DescList xs) -> fromDescList xs == fromListViaSnoc dlist+  , prop "cmp a b vs. cmp b a"     $ \p q -> cmp p q == reverseOrdering (cmp q p)+  , prop "ineq"                    $ \p q -> ineqTestPartition p q+  , prop "diffSequence"            $ \p -> diffSequence p == refDiffSeq (toList p)+  , prop "reverseDiffSequence"     $ \p -> reverseDiffSequence p == reverse (refDiffSeq (toList p))+  , prop "dual . dual == id"       $ \p -> dualPartition (dualPartition p) == p+  , prop "dual = reference impl."  $ \p -> toList (dualPartition p) == P._dualPartition (toList p)+  , prop "fromExpo . toExpo == id" $ \p -> fromExponentialForm (toExponentialForm p) == p+  , prop "isSubPartitionOf"        $ \p q -> isSubPartitionOf p q == P._isSubPartitionOf (toList p) (toList q)+  , prop "dominates"               $ \p q -> dominates p q == P._dominates (toList p) (toList q)+  , prop "pieriRuleSingleBox"      $ \p -> map toList (pieriRuleSingleBox p) =%%= P._pieriRule (toList p) 1 ++  , localOption (QuickCheckTests 100) +  $ prop "pieriRule /1"            $ \p (PieriK k) -> map toList (pieriRule p k) =%%= P._pieriRule (toList p) k +  , localOption (QuickCheckTests 100) +  $ prop "pieriRule /2"            $ \p (PieriK k) -> (pieriRule p k) =%%= map fromDescList (P._pieriRule (toList p) k) +  ]++ineqTestPartition :: Partition -> Partition -> Bool+ineqTestPartition = ineqTest++--------------------------------------------------------------------------------++ineqTest :: Ord a => a -> a -> Bool+ineqTest a b = case (a<b , a==b , a>b) of+  (True ,False,False) -> True+  (False,True ,False) -> True+  (False,False,True ) -> True+  _                   -> False++every10th :: [a] -> [a]+every10th = everyNth 10++everyNth :: Int -> [a] -> [a]+everyNth k = go where+  go []     = []+  go (x:xs) = x : drop (k-1) xs++infix 4 =%%=+(=%%=) :: Ord a => [a] -> [a] -> Bool+(=%%=) xs ys = (sort xs == sort ys)++allTrue :: [Bool] -> Assertion+allTrue bools = (and bools @=? True)++prop :: Testable a => TestName -> a -> TestTree+prop = Q.testProperty++(<#>) :: (a -> c) -> (b -> d) -> (a,b) -> (c,d)+(<#>) f g (x,y) = (f x, g y)++reverseOrdering :: Ordering -> Ordering+reverseOrdering ord = case ord of+  LT -> GT+  GT -> LT+  EQ -> EQ++--------------------------------------------------------------------------------++newtype SmallN = SmallN Int deriving (Eq,Ord,Show)++instance Arbitrary SmallN where+  arbitrary = SmallN <$> choose (1,80)++newtype PieriK = PieriK Int deriving (Eq,Ord,Show)++instance Arbitrary PieriK where+  arbitrary = PieriK <$> choose (1,3)++newtype DescList = DescList [Int] deriving (Eq,Ord,Show)++instance Arbitrary DescList where+  arbitrary = (DescList . reverse . sort . map getPositive) <$> arbitrary++instance Arbitrary Partition where+  arbitrary = do+    DescList xs <- arbitrary+    return $ fromDescList xs++--------------------------------------------------------------------------------++_allparts n = P._allPartitions n+_parts    n = P._partitions    n ++_staircase n   = [n,n-1..1]+_rectangle h n = replicate n h++_slope  n =         concatMap (\x->[x,x]) [n,n-1..1]+_slope1 n =         concatMap (\x->[x,x]) [n,n-1..1] ++ [1]+_slope2 n = (n+1) : concatMap (\x->[x,x]) [n,n-1..1] +_slope3 n = (n+1) : concatMap (\x->[x,x]) [n,n-1..1] ++ [1]++_steep  n = [n,n-2..1]++--------------------------------------------------------------------------------++allparts = map fromDescList . _allparts+parts    = map fromDescList . _parts++staircase n   = fromDescList $ _staircase n +rectangle h n = fromDescList $ _rectangle h n ++slope n = fromDescList $ _slope n+steep n = fromDescList $ _steep n++--------------------------------------------------------------------------------++testPartitionsSmall = map fromDescList _testPartitionsSmall+testPartitions      = map fromDescList _testPartitions++_testPartitionsSmall = concat+  [ _allparts 25+  , [ _rectangle 1   n | n <- [1..75] ]+  , [ _rectangle 15  n | n <- [1..75] ]+  , [ _rectangle 16  n | n <- [1..75] ]+  , [ _rectangle 255 n | n <- [1..75] ]+  , [ _rectangle 256 n | n <- [1..75] ]+  , [ _staircase n     | n <- [1..75] ]+  , [ _slope  n        | n <- [1..75] ]+  , [ _slope1 n        | n <- [1..75] ]+  , [ _slope2 n        | n <- [1..75] ]+  , [ _slope3 n        | n <- [1..75] ]+  , [ _steep  n        | n <- [1..75] ] +  , [ _staircase n     | n <- [200,255,256,257,258] ]+  ]++_testPartitions = _testPartitionsSmall ++ _testPartRandom++_testPartRandom = concat+  [ [ drop n seq       | seq <- randomSequences , n<-[0..79] ] +  , [ take n seq       | seq <- randomSequences , n<-[0..79] ] +  ]++--------------------------------------------------------------------------------+-- * reference++snocN :: Int -> (Partition,Int) -> Partition+snocN n (p,x) +  | n <= 0    = p+  | otherwise = snocN (n-1) (snoc p x, x)++consN :: Int -> (Int,Partition) -> Partition+consN n (x,p) +  | n <= 0    = p+  | otherwise = consN (n-1) (x, cons x p)++tailN :: Int -> Partition -> Partition+tailN n p +  | n <= 0    = p +  | otherwise = tailN (n-1) (partitionTail p)++unconsTest :: Partition -> Maybe (Int,[Int])+unconsTest xs = case uncons xs of+  Nothing    -> Nothing+  Just (x,p) -> Just (x, toList p)++unconsList :: [Int] -> Maybe (Int,[Int])+unconsList xs = case xs of+  []     -> Nothing+  (x:xs) -> Just (x,xs)++safeHead :: [Int] -> Int+safeHead xs = case xs of+  []     -> 0+  (x:xs) -> x++refDiffSeq :: [Int] -> [Int]+refDiffSeq = go where+  go (x:ys@(y:_)) = (x-y) : go ys +  go [x] = [x]+  go []  = []++--------------------------------------------------------------------------------++toListViaUncons :: Partition -> [Int]+toListViaUncons = go where+  go p = case uncons p of+    Nothing    -> []+    Just (x,p) -> x : go p+  +fromListViaCons :: DescList -> Partition+fromListViaCons (DescList list) = go list where+  go xs = case xs of+    []     -> empty+    (x:xs) -> cons x (go xs)++fromListViaSnoc :: DescList -> Partition+fromListViaSnoc (DescList list) = go (reverse list) where+  go xs = case xs of+    []     -> empty+    (x:xs) -> snoc (go xs) x++--------------------------------------------------------------------------------++{- +-- generated by:+import Data.List ; import Control.Monad ; import System.Random+genseq = (reverse . sort) <$> (replicateM 80 $ randomRIO (1::Int,255))+main = do+  seqs <- replicateM 64 genseq+  writeFile "rnd.txt" $ unlines $ map show seqs+-}++randomSequences =   +  [ [255,255,251,249,245,239,238,235,233,232,224,222,214,214,213,209,202,197,197,194,193,189,181,168,165,164,164,162,156,152,152,152,151,148,147,142,136,132,130,127,124,123,120,119,117,115,112,112,107,105,90,84,81,75,74,69,64,57,56,55,53,49,48,47,44,43,39,38,31,20,19,18,18,17,17,15,14,10,9,7]+  , [252,246,241,240,238,233,232,230,229,228,227,221,219,214,212,212,204,203,201,195,192,192,188,186,186,182,180,179,176,174,170,166,162,162,153,149,147,141,140,140,139,139,135,133,128,127,126,126,125,125,124,124,114,114,113,112,108,102,94,91,82,82,80,73,71,67,57,56,51,41,36,34,24,22,18,9,8,7,1,1]+  , [255,254,250,250,248,247,245,243,239,236,223,217,215,209,208,205,203,200,199,198,196,194,191,187,184,174,174,172,164,161,159,158,157,153,152,150,149,149,144,135,133,128,127,122,119,119,118,115,115,113,111,103,98,91,88,86,86,74,72,72,67,65,65,56,55,52,51,49,48,38,38,37,34,32,20,17,14,13,8,1]+  , [249,224,224,218,210,206,204,196,193,189,188,188,182,177,174,174,173,173,170,163,151,151,150,147,147,145,142,139,138,134,132,129,126,124,123,121,120,113,112,112,112,109,109,108,106,103,98,89,89,85,85,85,82,82,82,79,75,72,70,66,64,63,56,56,47,45,45,42,27,24,22,20,20,19,16,13,4,4,1,1]+  , [255,254,252,247,246,240,238,233,232,227,226,225,223,213,213,212,209,203,200,195,194,193,193,188,188,185,182,180,175,172,171,170,166,161,158,154,150,146,146,144,142,140,133,131,128,125,119,118,117,115,114,110,109,103,96,91,89,88,81,76,74,59,57,46,37,32,31,24,24,23,16,12,12,9,8,6,6,6,5,2]+  , [253,248,241,234,230,228,226,225,224,222,220,219,217,213,212,207,200,196,183,179,179,173,173,165,162,157,154,147,147,142,137,136,135,134,134,132,128,127,118,113,106,105,104,104,102,94,90,89,81,77,73,68,65,63,63,59,56,54,46,45,44,40,32,32,31,28,26,26,25,22,21,21,16,13,12,11,11,10,5,3]+  , [253,253,253,252,252,246,245,243,241,236,235,234,232,227,224,222,219,210,209,205,200,198,198,195,194,194,188,184,177,175,171,170,169,167,163,162,162,155,138,138,131,125,119,114,106,99,97,95,87,85,85,83,79,71,71,69,65,57,56,49,46,40,37,33,31,29,28,24,23,16,16,16,15,13,12,11,8,8,4,2]+  , [254,252,252,248,233,231,230,228,225,224,218,194,190,190,186,184,177,177,176,171,170,166,165,164,151,149,147,143,136,134,134,133,128,124,122,121,115,114,107,106,105,101,100,98,94,92,92,91,86,84,81,80,76,76,74,71,70,67,67,67,65,58,54,54,53,53,52,52,48,46,46,43,41,37,36,24,19,14,12,7]+  , [255,254,253,251,249,248,248,245,242,236,231,224,221,214,212,205,203,198,195,194,194,190,184,183,183,178,175,163,154,154,149,148,147,145,144,136,134,127,120,117,115,113,112,109,107,107,104,102,91,86,86,84,83,81,81,80,78,77,76,74,73,69,58,57,55,46,39,39,39,38,37,33,32,31,29,27,21,21,15,8]+  , [253,252,252,252,247,242,236,227,227,224,218,214,213,211,209,207,203,202,198,196,191,189,183,183,176,170,169,168,166,164,163,159,157,157,151,145,138,133,130,128,123,119,112,111,109,106,105,101,100,99,96,94,91,76,74,69,69,69,66,66,65,64,61,60,59,45,44,43,41,39,30,28,26,24,21,16,11,6,1,1]+  , [254,253,250,248,238,236,235,232,231,228,228,216,213,213,211,210,207,205,201,200,198,187,182,180,180,178,175,171,168,159,157,157,152,150,143,139,138,138,133,124,123,122,120,120,118,117,115,114,112,111,111,105,105,98,92,92,86,83,80,78,76,70,69,66,65,57,55,54,45,41,36,36,34,21,20,19,12,6,5,3]+  , [252,252,250,249,245,239,229,228,223,218,211,210,207,204,202,201,200,197,189,187,185,181,180,177,175,173,170,168,159,158,158,156,152,145,140,140,137,135,134,133,132,111,108,105,100,97,91,90,89,83,79,78,76,68,67,66,61,58,57,56,50,48,48,47,41,37,33,33,27,27,26,14,13,10,9,9,8,8,8,5]+  , [255,255,246,238,236,235,231,227,226,224,217,214,214,207,204,203,200,200,198,196,191,189,189,187,185,181,179,177,173,172,167,164,164,163,159,154,149,149,143,142,140,136,135,135,126,125,120,112,112,101,100,98,96,96,93,92,89,89,80,79,77,76,76,56,50,48,41,41,40,37,32,29,21,19,19,19,17,11,7,4]+  , [248,246,242,237,235,233,233,232,231,229,229,227,222,222,214,213,208,204,203,199,197,194,192,192,192,190,189,184,183,179,175,175,174,173,171,166,162,160,155,155,153,149,146,145,144,137,127,119,103,100,99,99,96,95,93,87,86,86,84,82,80,80,79,74,66,65,63,59,57,49,47,45,39,32,31,29,25,13,9,3]+  , [252,251,246,243,236,232,228,226,225,217,214,212,206,204,202,201,198,195,194,193,189,187,184,181,177,170,170,165,164,161,157,156,156,145,142,138,136,136,129,128,118,117,116,115,106,97,96,96,96,94,88,84,84,83,83,78,76,75,74,71,70,67,67,64,54,51,47,37,34,31,26,19,14,10,9,9,6,5,3,1]+  , [254,253,247,245,244,244,241,240,237,234,234,233,231,222,217,205,202,199,196,195,192,192,190,184,181,180,179,166,166,162,161,160,157,156,155,149,140,137,134,130,130,129,122,122,111,109,101,100,98,96,96,90,84,83,83,81,74,69,68,64,54,51,51,50,41,41,37,36,34,30,23,21,20,20,16,12,9,5,4,2]+  , [247,241,241,240,240,237,236,231,230,222,216,206,202,201,200,187,185,179,174,170,169,168,165,160,158,155,154,153,147,145,140,137,136,133,132,132,116,111,108,100,98,94,92,90,89,87,87,87,87,84,77,74,74,71,71,70,67,65,63,58,57,56,49,48,47,46,43,43,42,38,37,33,26,19,17,17,17,11,6,3]+  , [253,250,250,248,246,244,235,229,228,225,223,222,219,219,218,218,217,217,214,213,212,210,210,208,206,197,194,190,184,181,179,173,171,169,164,162,159,143,133,132,130,114,112,107,107,106,105,100,88,87,86,85,83,79,78,78,78,78,66,64,62,61,61,52,49,35,31,28,25,24,22,21,20,18,18,16,15,8,2,1]+  , [252,243,243,234,232,231,219,219,218,210,208,194,193,191,186,184,181,178,172,171,169,168,162,154,148,148,144,143,131,129,126,125,121,119,115,114,111,106,103,101,99,97,92,88,84,82,79,74,71,66,66,63,58,47,47,47,44,43,42,41,39,38,33,33,32,31,29,29,27,24,24,20,20,19,17,16,15,13,11,9]+  , [255,255,245,240,240,236,236,236,236,234,232,227,225,221,218,214,213,207,206,199,192,190,190,186,181,176,168,167,165,164,162,162,156,155,155,151,150,149,143,143,143,142,139,136,136,132,131,130,125,121,121,120,117,100,95,88,88,88,87,82,81,73,70,67,66,64,57,54,46,44,41,35,30,30,19,15,13,12,9,6]+  , [254,252,249,249,248,247,247,241,239,232,231,228,224,220,214,214,213,206,202,200,198,195,194,192,185,184,183,180,179,176,157,149,144,142,138,138,136,132,130,125,123,123,123,117,117,117,113,112,105,103,97,84,82,79,78,72,71,67,61,57,53,53,52,51,44,41,36,34,30,28,27,27,26,23,22,21,13,5,5,3]+  , [253,250,243,239,235,228,228,225,218,217,217,211,209,209,206,202,199,196,188,185,181,179,178,174,173,172,170,170,168,166,153,153,149,146,146,145,143,139,138,122,118,117,110,109,107,106,98,95,93,88,86,85,83,77,68,66,63,62,62,62,59,55,50,49,48,47,47,47,41,35,34,33,32,26,24,16,15,12,11,9]+  , [254,252,248,245,244,243,238,238,238,232,231,227,225,220,219,218,217,216,214,211,211,209,207,205,204,203,201,200,195,193,186,180,178,174,173,173,170,170,165,157,151,150,148,140,134,123,116,109,105,105,105,96,95,94,91,89,86,81,76,73,71,57,53,50,47,47,47,47,46,45,42,38,34,27,25,19,17,14,7,7]+  , [253,250,246,245,241,241,239,238,235,231,230,230,223,213,205,199,199,198,197,195,193,190,187,184,184,175,174,173,172,170,169,168,167,165,163,154,152,140,140,139,134,131,130,128,126,119,118,112,112,101,95,94,94,91,88,86,86,83,76,75,75,75,70,58,58,57,49,46,42,41,41,34,33,32,31,29,28,17,16,7]+  , [253,253,251,246,246,240,232,228,225,223,218,218,212,211,210,202,197,192,192,192,192,189,186,186,182,180,177,168,163,163,158,155,150,149,141,134,130,130,130,126,118,115,115,110,103,103,103,92,90,85,81,79,75,70,69,69,54,54,51,49,47,44,42,38,35,35,30,27,18,13,11,11,10,8,7,6,4,3,3,3]+  , [254,251,242,239,238,233,230,227,225,225,215,215,212,209,207,205,205,203,195,195,189,187,187,187,186,182,175,172,171,169,165,160,158,157,157,149,148,143,143,142,142,136,132,129,120,119,117,112,111,110,108,107,100,99,93,91,88,88,86,80,80,77,76,73,70,70,69,67,62,60,51,49,29,25,19,18,16,9,8,1]+  , [251,250,250,249,239,238,233,233,227,227,224,222,221,220,220,218,212,203,199,199,195,193,193,189,186,182,176,175,168,163,162,147,145,145,139,138,136,130,128,128,126,124,122,114,109,106,99,98,90,87,87,83,82,76,72,71,68,65,64,63,62,62,59,56,55,46,41,40,35,31,29,29,28,14,13,10,10,8,7,6]+  , [250,245,239,237,235,235,232,231,224,224,222,219,218,215,211,206,191,190,190,188,179,177,173,170,169,167,166,162,162,160,159,158,155,145,142,140,140,135,131,129,121,116,116,108,106,102,101,96,95,91,89,88,84,80,77,75,70,70,68,60,58,55,54,52,49,44,43,42,39,36,33,32,32,29,23,21,18,17,14,11]+  , [254,254,250,247,241,241,239,232,230,223,221,217,214,200,197,196,194,194,188,185,184,183,183,174,172,171,168,164,164,158,148,146,144,141,133,131,126,125,124,124,123,122,120,114,113,111,110,106,104,104,99,99,92,88,84,74,71,70,67,66,66,59,57,56,55,52,51,50,44,41,39,37,32,29,28,27,25,9,6,3]+  , [248,243,241,232,227,226,225,225,220,218,218,218,214,211,206,203,201,198,196,195,190,189,186,175,174,163,155,150,146,142,141,139,137,136,135,135,133,126,123,118,118,113,110,108,107,106,105,102,102,102,101,100,100,99,96,95,95,91,87,84,80,79,74,73,73,72,60,55,54,54,52,49,44,38,20,19,19,18,14,5]+  , [255,253,253,250,249,245,243,242,239,237,233,233,224,224,222,219,218,210,208,204,192,180,176,168,164,160,160,158,156,155,154,151,150,149,149,142,141,139,138,131,129,126,123,123,114,114,110,110,110,103,99,98,92,92,88,88,87,85,81,79,76,72,67,65,63,63,62,59,57,52,51,50,46,44,41,36,23,15,7,1]+  , [251,248,246,243,236,235,230,227,226,224,221,219,218,216,215,214,214,207,204,202,202,196,192,186,181,179,179,178,176,175,174,169,164,152,149,147,144,139,136,135,134,134,129,125,123,114,105,102,99,89,84,84,81,80,78,75,62,62,59,59,58,58,57,56,55,54,45,38,35,30,29,22,22,16,15,13,11,7,6,4]+  , [255,254,253,243,241,230,230,228,228,227,225,223,223,213,213,205,204,203,198,192,192,190,186,180,180,177,177,174,173,161,156,155,134,127,125,123,120,111,108,102,100,100,97,89,86,80,78,76,75,74,74,74,71,71,70,70,66,66,62,60,59,56,49,47,44,42,38,35,30,30,17,17,17,15,13,13,9,9,6,3]+  , [255,251,249,245,239,239,237,235,230,229,225,224,221,219,213,211,209,207,204,201,194,194,190,186,185,185,183,171,168,166,164,160,159,153,148,140,137,136,133,128,125,121,121,121,118,112,112,109,108,108,106,100,95,94,90,88,88,87,85,84,81,70,68,67,67,66,56,48,43,38,36,35,33,32,32,24,22,19,17,5]+  , [252,252,252,249,245,243,240,236,228,227,224,222,222,216,212,209,206,204,201,192,191,190,185,179,173,172,169,169,165,162,158,158,149,148,147,140,133,128,124,122,122,121,116,115,114,108,103,102,101,98,97,93,93,92,88,87,86,86,78,69,64,63,57,52,45,43,37,35,34,31,26,22,21,18,16,13,9,5,3,2]+  , [254,254,251,248,248,247,245,235,223,222,207,207,206,204,197,194,192,188,184,182,180,178,175,175,173,167,167,165,164,163,163,160,157,155,150,144,141,138,135,131,125,124,121,114,109,101,99,96,91,90,88,85,79,76,67,66,65,59,57,57,57,53,47,34,32,27,26,26,24,22,20,18,18,10,6,5,4,3,2,1]+  , [254,253,250,249,241,240,239,238,234,233,230,229,226,226,226,225,222,219,218,212,207,207,202,199,196,194,189,189,188,184,175,170,166,165,161,160,158,157,154,136,132,128,125,124,122,119,118,116,112,108,108,104,94,91,90,86,80,78,76,72,72,69,58,57,56,52,52,45,44,41,34,27,26,24,18,12,11,5,3,3]+  , [251,245,242,240,237,236,234,227,226,225,225,225,224,223,221,219,209,206,206,206,200,200,198,197,195,195,195,194,192,190,189,186,180,172,169,163,161,155,150,150,148,146,145,135,135,128,120,118,116,113,112,108,108,105,103,102,102,94,92,86,82,80,76,72,61,61,54,52,52,49,46,42,38,34,31,30,26,24,16,14]+  , [255,254,246,245,242,241,239,236,234,233,231,229,227,226,219,217,217,216,214,213,212,205,203,197,195,194,192,191,186,182,179,166,160,158,157,156,156,148,146,138,133,131,129,119,113,112,111,108,107,106,105,103,101,97,86,86,78,77,77,74,63,61,59,55,47,45,45,40,39,36,35,33,31,26,24,18,17,15,11,8]+  , [255,245,245,244,240,239,233,232,228,227,225,224,217,211,211,209,206,205,203,201,200,198,197,192,187,174,172,161,161,157,153,152,151,146,146,145,144,142,139,135,131,128,127,123,121,114,113,112,110,97,90,87,85,82,78,78,78,76,75,70,69,69,68,68,66,56,49,45,45,43,43,34,34,33,32,27,26,19,15,7]+  , [252,246,244,234,223,218,218,217,215,209,206,202,201,198,197,196,196,192,190,183,179,178,176,176,172,163,162,162,159,153,150,149,142,141,134,134,133,128,127,117,114,112,110,109,108,105,102,97,92,87,84,79,77,75,74,71,70,66,65,56,56,53,53,50,44,36,36,33,27,22,20,14,10,9,9,8,7,5,4,1]+  , [254,252,251,251,250,248,247,238,235,235,231,227,225,216,212,212,208,193,189,188,186,179,162,160,158,157,156,150,148,148,141,141,140,134,133,133,131,123,118,116,112,109,106,106,106,106,94,91,88,86,84,80,79,72,71,71,61,56,55,53,49,43,40,36,33,33,29,29,29,22,21,18,16,11,11,8,7,4,4,4]+  , [254,250,248,246,238,238,237,237,229,225,225,224,221,220,216,216,214,212,210,209,208,197,193,192,181,181,179,170,168,166,161,159,156,156,148,146,139,127,125,125,124,120,115,113,109,101,96,85,82,81,77,71,71,67,65,59,59,57,51,51,48,47,46,45,45,44,43,43,41,38,36,31,24,24,11,9,8,7,4,3]+  , [252,237,235,230,228,223,221,220,209,206,206,203,200,195,191,189,185,176,174,168,166,166,159,159,159,156,155,151,149,146,146,143,143,138,137,131,129,120,118,117,116,115,111,108,105,104,103,102,98,98,97,96,95,91,85,70,62,60,55,54,54,47,45,38,37,35,26,25,25,21,18,17,14,12,12,11,10,6,4,3]+  , [247,247,247,245,238,234,234,230,229,227,225,224,223,221,219,215,209,207,206,204,202,200,199,199,193,185,184,182,181,174,171,171,166,162,162,158,156,155,153,151,148,145,136,134,132,122,116,115,111,106,104,97,95,94,93,89,82,73,70,63,63,60,56,51,47,47,45,42,40,39,38,33,31,29,24,18,18,18,12,3]+  , [254,252,249,247,241,239,236,230,228,228,228,228,223,222,222,219,212,205,203,202,201,199,196,195,194,194,191,183,177,169,168,167,166,164,159,157,155,153,151,147,145,144,136,132,129,121,118,110,110,108,106,104,102,98,96,96,95,93,84,80,73,65,64,62,55,53,48,46,44,43,43,40,33,30,27,26,20,18,14,9]+  , [253,245,245,245,244,238,236,235,230,230,230,223,221,220,220,217,213,210,203,201,193,182,165,163,162,161,161,155,154,153,152,148,147,140,133,132,128,120,118,116,114,114,113,96,91,89,81,80,79,76,76,68,68,67,65,63,55,51,50,49,46,44,41,40,39,37,32,32,30,29,21,21,20,18,15,13,13,10,7,5]+  , [253,242,239,238,238,237,234,233,229,228,227,225,222,222,221,217,212,208,203,197,197,191,190,190,187,185,184,182,178,175,173,169,164,164,163,158,146,145,144,140,134,131,127,125,113,111,107,106,98,97,97,97,96,95,89,89,83,82,77,71,60,60,59,55,54,51,39,37,35,31,28,22,22,21,19,14,14,8,4,3]+  , [253,245,243,237,237,227,226,224,222,221,218,217,216,214,213,206,201,192,192,191,190,189,188,188,182,179,169,167,159,148,140,139,138,135,132,126,125,124,122,117,116,114,100,96,90,85,82,81,80,79,77,75,72,69,69,64,61,58,53,53,51,50,49,49,47,44,43,39,38,38,37,37,33,26,17,16,11,8,5,2]+  , [254,252,248,245,244,242,238,235,235,233,233,228,227,219,208,205,202,200,199,194,190,190,188,183,181,178,177,173,166,163,162,153,143,138,136,128,126,120,119,118,117,114,110,109,108,107,104,103,96,95,94,92,84,80,78,72,70,69,68,67,64,64,63,59,54,50,48,43,40,35,28,27,26,25,22,18,17,16,6,3]+  , [254,250,249,242,239,239,234,233,219,219,218,215,215,205,201,201,199,199,198,198,195,194,191,187,184,184,184,183,182,179,175,175,162,161,161,158,158,154,153,153,149,147,138,127,126,123,119,114,113,108,106,106,105,104,101,101,100,99,98,92,90,88,88,86,86,82,77,71,68,68,60,53,44,42,37,30,23,9,4,1]+  , [255,255,250,249,247,243,242,234,233,229,228,228,227,225,223,218,213,212,211,210,210,198,190,186,184,184,175,173,171,166,160,156,155,152,151,151,147,143,138,134,131,131,118,117,115,115,105,97,95,95,91,90,83,77,73,70,70,68,66,54,50,49,45,38,33,32,30,25,24,20,17,13,13,13,10,8,7,6,5,2]+  , [253,251,249,247,245,243,237,231,230,224,223,222,218,217,214,212,212,209,202,201,186,179,174,171,169,166,164,164,162,157,155,152,148,142,139,134,128,123,120,110,106,104,104,101,100,98,97,97,93,85,80,80,78,77,76,74,70,65,65,54,54,54,43,42,37,32,28,27,26,26,25,22,16,15,14,14,11,5,3,1]+  , [252,250,243,241,238,232,232,223,211,211,211,208,207,206,205,204,201,201,201,199,190,189,188,179,166,158,156,155,152,150,148,148,147,146,144,139,135,132,131,128,118,116,115,114,114,107,99,94,93,89,85,83,82,80,77,73,71,70,68,66,60,58,57,55,53,49,43,40,39,33,28,28,27,27,13,7,6,6,3,1]+  , [255,252,250,250,249,249,244,237,233,232,231,231,228,227,227,218,217,211,201,199,198,198,196,189,189,189,189,181,181,178,173,168,168,167,162,160,160,158,156,155,149,145,145,143,141,137,137,130,130,126,120,116,112,109,105,103,102,100,92,90,89,88,86,82,73,71,71,70,67,66,65,63,50,49,46,26,22,18,17,8]+  , [255,252,246,245,239,237,237,233,233,232,231,231,220,220,218,214,210,206,206,203,199,197,196,194,194,188,185,184,179,156,156,154,149,146,146,143,138,136,136,135,133,130,130,110,109,108,108,102,98,94,93,91,88,88,87,83,79,79,78,78,75,75,74,74,69,58,53,51,50,47,33,31,31,22,15,10,8,5,1,1]+  , [253,247,246,244,244,242,241,235,234,228,225,212,210,209,207,204,203,201,199,194,192,187,185,184,180,180,177,174,172,171,170,167,165,155,155,154,150,148,144,141,139,137,130,118,117,116,114,112,102,102,96,90,83,81,79,74,70,68,61,57,49,48,47,46,45,45,42,40,34,32,28,27,19,17,14,9,8,3,2,1]+  , [242,239,237,234,232,229,225,221,214,212,210,202,202,199,196,196,194,193,184,183,182,181,177,176,175,170,167,167,166,164,161,157,155,155,152,146,144,138,137,132,130,130,124,115,115,114,114,111,101,97,92,87,78,78,78,73,70,68,62,58,57,55,51,51,50,43,43,40,39,36,34,33,32,31,24,23,22,16,11,9]+  , [255,253,251,251,251,249,249,246,245,244,242,236,224,216,210,210,206,205,199,199,197,195,193,192,191,190,189,182,181,180,173,171,170,169,160,159,149,144,142,140,139,136,130,125,123,117,95,92,91,85,79,77,76,75,63,61,59,57,57,53,48,47,41,41,38,29,28,27,27,24,17,17,17,15,13,9,6,5,4,4]+  , [252,252,246,245,237,232,232,228,221,218,215,214,211,208,206,202,201,200,198,198,194,191,191,180,178,174,169,166,166,164,163,162,159,159,158,149,146,144,138,132,127,125,123,122,121,118,115,113,107,107,106,105,102,101,100,92,86,83,80,75,74,73,72,72,65,61,59,55,54,47,47,43,31,24,17,14,14,9,8,4]+  , [253,250,249,246,242,241,236,230,226,223,219,213,201,199,198,191,182,182,179,174,171,168,168,166,161,156,155,153,149,144,140,137,136,132,129,129,126,123,121,118,106,102,101,100,99,95,93,89,89,87,86,81,80,78,74,73,73,72,66,62,59,57,55,49,49,43,41,38,35,32,28,27,24,20,17,15,11,6,5,2]+  , [254,253,249,246,240,229,223,216,212,210,207,206,206,203,203,201,197,195,189,184,183,182,178,175,172,170,166,164,154,151,145,142,141,140,138,125,124,118,117,116,114,107,103,103,94,92,89,86,84,83,81,78,68,66,64,56,56,54,50,47,46,44,40,35,29,27,27,26,26,25,24,24,18,17,12,8,8,7,3,2]+  , [255,248,244,243,238,237,237,231,229,228,228,219,214,211,208,202,202,199,195,192,191,184,183,179,174,170,168,166,163,159,158,158,154,133,130,130,127,126,125,123,121,119,114,98,98,89,89,88,87,83,79,73,69,65,62,58,57,56,51,49,49,48,43,43,40,36,34,33,31,26,23,22,20,19,18,17,17,14,11,7]+  , [254,247,246,246,242,240,237,236,236,229,224,224,219,217,217,215,214,202,201,197,193,189,177,172,172,170,168,164,161,156,155,153,152,152,146,145,144,144,140,138,127,124,123,121,115,110,106,99,99,98,94,91,90,90,89,77,77,73,72,70,69,68,66,65,63,60,58,56,54,47,47,41,40,35,34,22,19,18,12,2]+    ----+  , [510,509,496,489,477,476,462,455,452,443,442,426,424,422,407,406,394,380,377,376,375,335,333,328,323,314,309,299,292,288,287,285,271,265,232,231,212,204,192,191,190,184,182,181,181,164,163,156,154,145,141,141,139,136,128,123,122,112,112,97,97,94,91,86,84,72,67,65,58,56,56,55,50,49,48,47,38,24,14,13]+  , [505,488,475,470,467,466,462,461,437,423,419,399,395,390,390,386,384,381,379,378,372,371,369,365,355,344,344,336,332,322,306,298,296,292,285,278,268,264,252,241,236,229,227,225,219,213,211,205,205,192,189,189,185,172,167,161,156,151,150,150,146,144,139,132,128,123,117,96,80,72,63,43,42,37,37,37,24,22,20,4]+  , [508,498,496,482,476,473,468,460,435,433,427,423,423,402,393,392,387,381,367,360,357,353,353,348,343,335,324,313,311,299,298,297,291,281,263,263,258,258,246,245,239,223,220,213,211,205,198,195,190,181,180,151,150,147,133,125,122,111,109,108,98,89,87,84,82,77,72,67,64,63,62,57,42,34,33,32,15,9,5,4]+  , [506,505,495,488,487,485,485,483,450,448,443,440,435,427,426,422,400,398,396,389,375,369,367,366,358,358,353,350,338,334,326,316,316,306,296,291,288,287,271,237,234,228,221,214,210,207,202,192,192,190,190,183,182,181,175,170,165,159,155,150,145,139,136,134,132,114,113,105,103,83,82,71,70,59,34,24,12,7,5,4]+  ]++--------------------------------------------------------------------------------

@@ -1,7 +1,7 @@  -- | Tests for integer partitions. -{-# LANGUAGE CPP, BangPatterns #-}+{-# LANGUAGE CPP, BangPatterns, DataKinds, KindSignatures, ScopedTypeVariables #-} module Tests.Partitions.Integer where  --------------------------------------------------------------------------------@@ -13,6 +13,7 @@ import Tests.Common  import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Integer.Count  import Data.List import Control.Monad@@ -24,7 +25,30 @@ import Math.Combinat.Numbers ( factorial , binomial , multinomial ) import Math.Combinat.Helper +import Data.Proxy+import GHC.TypeLits+ --------------------------------------------------------------------------------++-- | Partitions of size at most n+newtype Part (n :: Nat) = Part (Partition) deriving (Eq,Show)++-- | usage: fromPart @20+fromPart :: Part n -> Partition+fromPart (Part p) = p++fromPart20 :: Part 20 -> Partition+fromPart20 (Part p) = p++fromPart30 :: Part 30 -> Partition+fromPart30 (Part p) = p ++instance forall n. KnownNat n => Arbitrary (Part n) where+  arbitrary = do+    n <- choose (0, fromInteger (natVal (Proxy :: Proxy n)))+    myMkGen' Part (randomPartition n)++-------------------------------------------------------------------------------- -- * Types and instances  newtype PartitionWeight     = PartitionWeight     Int              deriving (Eq,Show)@@ -69,6 +93,7 @@   , testProperty "dominated partitions"            prop_dominated_list   , testProperty "dominating partitions"           prop_dominating_list   , testProperty "counting partitions"             prop_countParts+  , testProperty "union/sum duality"               prop_union_sum_duality   ]  --------------------------------------------------------------------------------@@ -102,6 +127,9 @@  prop_countParts :: Bool prop_countParts = (take 50 partitionCountList == take 50 partitionCountListNaive)++prop_union_sum_duality :: Partition -> Partition -> Bool+prop_union_sum_duality p q = dualPartition (sumOfPartitions p q) == unionOfPartitions (dualPartition p) (dualPartition q)  -------------------------------------------------------------------------------- 

@@ -0,0 +1,86 @@++-- | Tests for ribbons (border strip skew partitions).+--++{-# LANGUAGE CPP, BangPatterns #-}+module Tests.Partitions.Ribbon where++--------------------------------------------------------------------------------++import Test.Framework+import Test.Framework.Providers.QuickCheck2+import Test.QuickCheck++import Tests.Common+import Tests.Partitions.Integer ( Part(..) , fromPart20 , fromPart30 )     -- Arbitrary instances++import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Skew+import Math.Combinat.Partitions.Skew.Ribbon++import Data.List++import Math.Combinat.Classes++--------------------------------------------------------------------------------+-- * instances++data Inner = Inner Partition Int deriving (Eq,Show)++data Outer = Outer Partition Int deriving (Eq,Show)++instance Arbitrary Inner where+  arbitrary = do+    p <- arbitrary+    let (w,h) = (partitionWidth p, partitionHeight p)+        n = w+h-1+    d <- choose (1,n)+    return $ Inner p d++instance Arbitrary Outer where+  arbitrary = do+    pp <- arbitrary+    let p = fromPart30 pp    -- smaller partitions+    let (w,h) = (partitionWidth p, partitionHeight p)+        n = w+h+10   +    d <- choose (1,n)+    return $ Outer p d++--------------------------------------------------------------------------------+-- * test group++testgroup_Ribbon :: Test+testgroup_Ribbon = testGroup "Ribbons and Corners" +  [ testGroup "Ribbons"  +      [ testProperty "all inner ribbons vs. naive"        prop_inner_all+      , testProperty "inner ribbons of length vs. naive"  prop_inner_length+      , testProperty "outer ribbons of length vs. naive"  prop_outer_length+      ]+  , testGroup "Corners"+      [ testProperty "inner corner boxes vs. naive" prop_innerCornerBoxes+      , testProperty "outer corner boxes vs. naive" prop_outerCornerBoxes +      ]+  ]++--------------------------------------------------------------------------------+-- * ribbon properties++prop_inner_all :: Partition -> Bool+prop_inner_all p = sort (innerRibbons p) == sort (innerRibbonsNaive p)++prop_inner_length :: Inner -> Bool+prop_inner_length (Inner p n) = sort (innerRibbonsOfLength p n) == sort (innerRibbonsOfLengthNaive p n)++prop_outer_length :: Outer -> Bool+prop_outer_length (Outer p n) = sort (outerRibbonsOfLength p n) == sort (outerRibbonsOfLengthNaive p n)++--------------------------------------------------------------------------------+-- * corner properties++prop_innerCornerBoxes :: Partition -> Bool+prop_innerCornerBoxes p  =  (innerCornerBoxes p == innerCornerBoxesNaive p)++prop_outerCornerBoxes :: Partition -> Bool+prop_outerCornerBoxes p  =  (outerCornerBoxes p == outerCornerBoxesNaive p)++--------------------------------------------------------------------------------

@@ -2,7 +2,7 @@ -- | Tests for skew partitions. -- -{-# LANGUAGE CPP, BangPatterns #-}+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, DataKinds, KindSignatures #-} module Tests.Partitions.Skew where  --------------------------------------------------------------------------------@@ -12,7 +12,7 @@ import Test.QuickCheck  import Tests.Common-import Tests.Partitions.Integer ()     -- Arbitrary instances+import Tests.Partitions.Integer ( Part(..) , fromPart20 , fromPart30 )     -- Arbitrary instances  import Math.Combinat.Partitions.Integer import Math.Combinat.Partitions.Skew@@ -21,6 +21,9 @@  import Math.Combinat.Classes +import Data.Proxy+import GHC.TypeLits+ -------------------------------------------------------------------------------- -- * instances @@ -34,6 +37,32 @@     k <- choose (0,ln-1)     let q = qs !! k     return $ mkSkewPartition (p,q) ++--------------------------------------------------------------------------------++-- | Skew partitions of size at most n+newtype Skew (n :: Nat) = Skew (SkewPartition) deriving (Eq,Show)++-- | usage: fromSkew @20+fromSkew :: Skew n -> SkewPartition+fromSkew (Skew p) = p++fromSkew20 :: Skew 20 -> SkewPartition+fromSkew20 (Skew p) = p++fromSkew30 :: Skew 30 -> SkewPartition+fromSkew30 (Skew p) = p ++instance forall nn. KnownNat nn => Arbitrary (Skew nn) where+  arbitrary = do+    Part p <- arbitrary :: Gen (Part nn)+    let n = partitionWeight p+    d <- choose (0,n)+    let qs = subPartitions d p+        ln = length qs+    k <- choose (0,ln-1)+    let q = qs !! k+    return $ Skew $ mkSkewPartition (p,q)   -------------------------------------------------------------------------------- -- * test group

@@ -56,6 +56,8 @@  data SameSize = SameSize Permutation Permutation deriving Show +data PermWithList = PWL Permutation [Int] deriving (Show)+ instance Random Permutation where   random g = randomPermutation size g1 where     (size,g1) = randomR (minPermSize,maxPermSize) g@@ -81,6 +83,20 @@     (prm2,g3) = randomPermutation size g2   randomR _ = random +randomRList :: (RandomGen g, Random a) => Int -> (a, a) -> g -> ([a],g)+randomRList n ab g0 = go n g0 where+  go 0   g = ([],g)+  go !k !g = let (x ,g' ) = randomR ab g +                 (xs,g'') = go (k-1) g'+             in  (x:xs,g'')++instance Random PermWithList where+  random g = (PWL prm xs, g3) where+    (size,g1) = randomR (minPermSize,maxPermSize) g+    (prm ,g2) = randomPermutation size g1 +    (xs  ,g3) = randomRList size (-100,100) g2+  randomR _ = random+ instance Arbitrary Nat where   arbitrary = choose (Nat 0 , Nat 50)      @@ -88,6 +104,7 @@ instance Arbitrary CyclicPermutation where arbitrary = choose undefined instance Arbitrary DisjointCycles    where arbitrary = choose undefined instance Arbitrary SameSize          where arbitrary = choose undefined+instance Arbitrary PermWithList      where arbitrary = choose undefined  -------------------------------------------------------------------------------- -- * test group@@ -129,6 +146,9 @@   , testProperty "number of inversions is the same for the inverse permutation"  prop_ninversions_inverse   , testProperty "merge sort algorithm = naive inversion count"                  prop_merge_inversions +  , testProperty "sortingPermutationAsc"    prop_sortingPermAsc+  , testProperty "sortingPermutationDesc"   prop_sortingPermDesc+  , testProperty "concatPermutations"       prop_concatPerm   ]  --------------------------------------------------------------------------------@@ -214,6 +234,14 @@ prop_ninversions_inverse perm = numberOfInversions perm == numberOfInversions (inverse perm)  prop_merge_inversions perm = (numberOfInversionsMerge perm == numberOfInversionsNaive perm)++prop_sortingPermAsc :: [Int] -> Bool +prop_sortingPermAsc xs = permuteList (sortingPermutationAsc xs) xs == sort xs++prop_sortingPermDesc :: [Int] ->  Bool+prop_sortingPermDesc xs = permuteList (sortingPermutationDesc xs) xs == reverse (sort xs)++prop_concatPerm (PWL p1 xs) (PWL p2 ys) = permuteList p1 xs ++ permuteList p2 ys == permuteList (concatPermutations p1 p2) (xs++ys)  -------------------------------------------------------------------------------- 

@@ -2,7 +2,7 @@ -- | Tests for power series -- -{-# LANGUAGE CPP, GeneralizedNewtypeDeriving #-}+{-# LANGUAGE CPP, GeneralizedNewtypeDeriving, DataKinds, KindSignatures #-} module Tests.Series where  --------------------------------------------------------------------------------@@ -15,7 +15,11 @@ import System.Random  import Data.List+import Data.Ratio +import GHC.TypeLits+import Data.Proxy+ import Math.Combinat.Sign import Math.Combinat.Numbers import Math.Combinat.Partitions.Integer@@ -123,25 +127,89 @@ swap (x,y) = (y,x) -} --- compare the first 500 elements of the infinite lists+-- compare the first N elements of the infinite lists+(=..=) :: (Eq a, Num a) => Int -> [a] -> [a] -> Bool+(=..=) n xs1 ys1 = take n xs == take n ys where+  xs = xs1 ++ repeat 0+  ys = ys1 ++ repeat 0++infix 4 =..=++-- compare the first 100 elements of the infinite lists (=!=) :: (Eq a, Num a) => [a] -> [a] -> Bool (=!=) xs1 ys1 = (take m xs == take m ys) where -  m = 500+  m = 100   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)+-- compare the first 500 elements of the infinite lists+(=!!=) :: (Eq a, Num a) => [a] -> [a] -> Bool+(=!!=) xs1 ys1 = (take m xs == take m ys) where +  m = 500+  xs = xs1 ++ repeat 0+  ys = ys1 ++ repeat 0++infix 4 =!!=++newtype XNat   = XNat   { fromXNat   :: Int      } deriving (Eq,Ord,Show,Num,Random)++newtype Rat    = Rat    { fromRat    :: Rational } deriving (Eq,Ord,Show,Num,Fractional)+newtype NZRat  = NZRat  { fromNZRat  :: Rational } deriving (Eq,Ord,Show,Num,Fractional)++-- type parameter is for controlling the size (length), because some tests are too slow+newtype Ser  (n :: Nat) = Ser  { fromSer'  :: [Integer] } deriving (Eq,Ord,Show)+newtype SerR (n :: Nat) = SerR { fromSerR' :: [Rational] } 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+---------------------------------------- +serProxy :: f (n :: Nat) -> Proxy n+serProxy _ = Proxy++seriesSize :: KnownNat (n :: Nat) => f (n :: Nat) -> Int+seriesSize ser = fromInteger $ natVal (serProxy ser) where ++----------------------------------------++fromSer  = fromSer500+fromSerR = fromSerR500++fromSer25 :: Ser 25 -> [Integer]+fromSer25 = fromSer'++fromSer100 :: Ser 100 -> [Integer]+fromSer100 = fromSer'++fromSer500 :: Ser 500 -> [Integer]+fromSer500 = fromSer'++----------------------------------------++fromSerR25 :: SerR 25 -> [Rational]+fromSerR25 = fromSerR'++fromSerR50 :: SerR 50 -> [Rational]+fromSerR50 = fromSerR'++fromSerR100 :: SerR 100 -> [Rational]+fromSerR100 = fromSerR'++fromSerR500 :: SerR 500 -> [Rational]+fromSerR500 = fromSerR'++----------------------------------------++{-+minSerSize = 0   :: Int+maxSerSize = 500 :: Int+-}+ minSerValue = -10000 :: Int maxSerValue =  10000 :: Int @@ -149,9 +217,20 @@ 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 Random Rat where+  random g = (Rat (fromIntegral x % fromIntegral y), g'') where+    (x,g' ) = randomR (-100,100::Int) g+    (y,g'') = randomR (   1, 25::Int) g'        -- hackety hack hack+  randomR _ g = random g +instance Random NZRat where+  random g = let (Rat q , g') = random g+             in  if q /= 0 then (NZRat q, g') else random g'            +  randomR _ g = random g++instance Arbitrary XNat where+  arbitrary = choose (XNat 0 , XNat 750)+ instance Arbitrary Exp where   arbitrary = choose (Exp 1 , Exp 32) @@ -161,14 +240,24 @@     exp   <- arbitrary :: Gen Exp     return $ CoeffExp (fromIntegral coeff, fromExp exp)    -instance Random Ser where-  random g = (Ser $ map fi list, g2) where-    (size,g1) = randomR (minSerSize,maxSerSize) g+instance KnownNat (n :: Nat) => Random (Ser n) where+  random g = (series, g2) where+    maxSerSize = seriesSize series+    series     = Ser (map fi list) +    (size,g1) = randomR (0,maxSerSize) g     (list,g2) = rndList size (minSerValue,maxSerValue) g1     fi :: Int -> Integer     fi = fromIntegral    randomR _ = random +instance KnownNat (n :: Nat) => Random (SerR n) where+  random g = (series, g2) where+    maxSerSize = seriesSize series+    series    = SerR (map fromRat list) +    (size,g1) = randomR (0,maxSerSize) g+    (list,g2) = rndList size (fromIntegral minSerValue, fromIntegral maxSerValue) g1+  randomR _ = random+ instance Random Exps where   random g = (Exps list, g2) where     (size,g1) = randomR (0,10) g@@ -181,10 +270,19 @@     (list1,g2) = rndList size (1,32) g1     (list2,g3) = rndList size (minSerValue,maxSerValue) g2   randomR _ = random++instance Arbitrary Rat where+  arbitrary = choose undefined++instance Arbitrary NZRat where+  arbitrary = choose undefined   -instance Arbitrary Ser where+instance KnownNat n => Arbitrary (Ser n) where   arbitrary = choose undefined +instance KnownNat n => Arbitrary (SerR n) where+  arbitrary = choose undefined+ instance Arbitrary Exps where   arbitrary = choose undefined @@ -197,7 +295,26 @@ testgroup_PowerSeries :: Test testgroup_PowerSeries = testGroup "Power series"   [ -    testProperty "convPSeries1 vs generic"     prop_conv1_vs_gen+    testProperty "mulSeries  == mulSeriesNaive"   prop_mulSeries_vs_naive+  , testProperty "divSeries  == mulWithRecip"     prop_divSeries_vs_mult_with_recip+  , testProperty "recip xs   == 1 / xs"           prop_recipSeries_vs_one_over+  , testProperty "compose    == composeNaive"     prop_compose_vs_naive+  , testProperty "substitute == substituteNaive"  prop_substitute_vs_naive+  , testProperty "inversion  == inversionNaive"   prop_inversion_vs_naive++  , testProperty "lagrange inversion works /1"       prop_lagrange_inversion1+  , testProperty "lagrange inversion works /2"       prop_lagrange_inversion2+  , testProperty "naive lagrange inversion works /1"       prop_lagrange_inversion_naive1+  , testProperty "naive lagrange inversion works /2"       prop_lagrange_inversion_naive2+  , testProperty "integral naive lagrange inversion works /1"       prop_lagrange_inversion_int_naive1+  , testProperty "integral naive lagrange inversion works /2"       prop_lagrange_inversion_int_naive2++  , testProperty "diff . int == id"            prop_diff_integrate+  , testProperty "tail (int . diff) == tail"   prop_integrate_diff+  , testProperty "sin vs sin2"                 prop_sin_vs_sin2+  , testProperty "cos vs cos2"                 prop_cos_vs_cos2++  , testProperty "convPSeries1 vs generic"     prop_conv1_vs_gen   , testProperty "convPSeries2 vs generic"     prop_conv2_vs_gen   , testProperty "convPSeries3 vs generic"     prop_conv3_vs_gen   , testProperty "convPSeries1' vs generic"    prop_conv1_vs_gen'@@ -217,23 +334,69 @@  -------------------------------------------------------------------------------- -- * properties++prop_mulSeries_vs_naive ser1 ser2 = (mulSeries xs ys =!= mulSeriesNaive xs ys) where+  xs = fromSer ser1+  ys = fromSer ser2++prop_divSeries_vs_mult_with_recip (NZRat q) ser1 ser2 = (=..=) 60 (divSeries xs ys) (mulSeries xs (reciprocalSeries ys)) where+  xs =     fromSerR100 ser1+  ys = q : fromSerR100 ser2++prop_recipSeries_vs_one_over (NZRat q) ser = (reciprocalSeries xs =!= divSeries unitSeries xs) where+  xs = q : fromSerR100 ser++prop_compose_vs_naive ser1 ser2 = (=..=) 25 (composeSeries xs ys) (composeSeriesNaive xs ys) where+  xs =     fromSer25 ser1+  ys = 0 : fromSer25 ser2++prop_substitute_vs_naive ser1 ser2 = (=..=) 25 (substitute xs ys) (substituteNaive xs ys) where+  xs = 0 : fromSer25 ser1+  ys =     fromSer25 ser2++prop_inversion_vs_naive (NZRat q) ser = (=..=) 25 (lagrangeInversion xs) (lagrangeInversionNaive xs) where+  xs = 0 : q : fromSerR25 ser++prop_lagrange_inversion1 (NZRat q) ser = (=..=) 35 (substitute f (lagrangeInversion f)) (0 : 1 : repeat 0) where f = 0 : q : fromSerR50 ser+prop_lagrange_inversion2 (NZRat q) ser = (=..=) 35 (substitute (lagrangeInversion f) f) (0 : 1 : repeat 0) where f = 0 : q : fromSerR50 ser++prop_lagrange_inversion_naive1 (NZRat q) ser = (=..=) 20 (substituteNaive f (lagrangeInversionNaive f)) (0 : 1 : repeat 0) where f = 0 : q : fromSerR25 ser+prop_lagrange_inversion_naive2 (NZRat q) ser = (=..=) 20 (substituteNaive (lagrangeInversionNaive f) f) (0 : 1 : repeat 0) where f = 0 : q : fromSerR25 ser++prop_lagrange_inversion_int_naive1 ser = (=..=) 20 (substituteNaive f (integralLagrangeInversionNaive f)) (0 : 1 : repeat 0) where f = 0 : 1 : fromSer25 ser+prop_lagrange_inversion_int_naive2 ser = (=..=) 20 (substituteNaive (integralLagrangeInversionNaive f) f) (0 : 1 : repeat 0) where f = 0 : 1 : fromSer25 ser++--------------------------------------------------------------------------------++prop_diff_integrate ser = (xs =!= differentiateSeries (integrateSeries xs)) where+  xs = fromSerR ser++prop_integrate_diff ser = (0 : tail xs =!= integrateSeries (differentiateSeries xs)) where+  xs = fromSerR ser++prop_cos_vs_cos2 = (cosSeries =!= (cosSeries2 :: [Rational])) +prop_sin_vs_sin2 = (sinSeries =!= (sinSeries2 :: [Rational])) ++--------------------------------------------------------------------------------       prop_leftIdentity ser = ( xs =!= unitSeries `convolve` xs ) where -  xs = fromSer ser +  xs = fromSer100 ser   prop_rightIdentity ser = ( unitSeries `convolve` xs =!= xs ) where -  xs = fromSer ser +  xs = fromSer100 ser   prop_commutativity ser1 ser2 = ( xs `convolve` ys =!= ys `convolve` xs ) where -  xs = fromSer ser1-  ys = fromSer ser2+  xs = fromSer100 ser1+  ys = fromSer100 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+  xs = fromSer100 ser1+  ys = fromSer100 ser2+  zs = fromSer100 ser3++--------------------------------------------------------------------------------    prop_conv1_vs_gen exp1 ser = ( one =!= two ) where   one = convolveWithPSeries1 k1 xs 

@@ -1,7 +1,7 @@ 
 -- | Tests for skew tableaux
 
-{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances, TypeApplications, DataKinds #-}
 module Tests.SkewTableaux where
 
 --------------------------------------------------------------------------------
@@ -14,7 +14,7 @@ import Test.QuickCheck.Gen
 
 import Tests.Partitions.Integer ()
-import Tests.Partitions.Skew    ()      -- arbitrary instances
+import Tests.Partitions.Skew ( Skew(..) , fromSkew20 , fromSkew30 )     -- Arbitrary instances
 
 import Math.Combinat.Tableaux
 import Math.Combinat.Tableaux.Skew
@@ -52,7 +52,8 @@ 
 instance Arbitrary (SkewTableau Int) where
   arbitrary = do
-    shape <- arbitrary
+    pshape <- arbitrary
+    let shape = fromSkew20 pshape      -- skew partition of size at most 20
     let w = skewPartitionWeight shape
     content <- replicateM w $ choose (1,1000)
     return $ fillSkewPartitionWithRowWord shape content
@@ -90,8 +91,8 @@   tableau = fillSkewPartitionWithColumnWord shape [1..]
   shape'  = skewTableauShape tableau
 
-prop_semistandard :: SkewPartition -> Bool
-prop_semistandard shape = and 
+prop_semistandard :: Skew 20 -> Bool
+prop_semistandard (Skew shape) = and 
   [ isSemiStandardSkewTableau st 
   | n  <- [kk..nn] 
   , st <- take 500 (semiStandardSkewTableaux n shape)         -- we only take the first 500 because impossibly slow otherwise