diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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 diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs --- a/Math/Combinat/Helper.hs +++ b/Math/Combinat/Helper.hs @@ -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 diff --git a/Math/Combinat/Numbers.hs b/Math/Combinat/Numbers.hs --- a/Math/Combinat/Numbers.hs +++ b/Math/Combinat/Numbers.hs @@ -1,194 +1,9 @@ --- | A few important number sequences. --- --- See the \"On-Line Encyclopedia of Integer Sequences\", --- . - -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 --- -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 diff --git a/Math/Combinat/Numbers/Integers.hs b/Math/Combinat/Numbers/Integers.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Numbers/Integers.hs @@ -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 + ] +-} + +-------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Numbers/Primes.hs b/Math/Combinat/Numbers/Primes.hs --- a/Math/Combinat/Numbers/Primes.hs +++ b/Math/Combinat/Numbers/Primes.hs @@ -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 - ] -} -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Numbers/Sequences.hs b/Math/Combinat/Numbers/Sequences.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Numbers/Sequences.hs @@ -0,0 +1,198 @@ + +-- | Some important number sequences. +-- +-- See the \"On-Line Encyclopedia of Integer Sequences\", +-- . + +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 +-- +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 + +-------------------------------------------------------------------------------- + + + diff --git a/Math/Combinat/Numbers/Series.hs b/Math/Combinat/Numbers/Series.hs --- a/Math/Combinat/Numbers/Series.hs +++ b/Math/Combinat/Numbers/Series.hs @@ -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] diff --git a/Math/Combinat/Partitions/Integer.hs b/Math/Combinat/Partitions/Integer.hs --- a/Math/Combinat/Partitions/Integer.hs +++ b/Math/Combinat/Partitions/Integer.hs @@ -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 --- -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 --- -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 diff --git a/Math/Combinat/Partitions/Integer/Compact.hs b/Math/Combinat/Partitions/Integer/Compact.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Partitions/Integer/Compact.hs @@ -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 +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 +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 = [] + +-------------------------------------------------------------------------------- + + diff --git a/Math/Combinat/Partitions/Integer/Count.hs b/Math/Combinat/Partitions/Integer/Count.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Partitions/Integer/Count.hs @@ -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 +-- +-- +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) ] +-} + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/Integer/IntList.hs b/Math/Combinat/Partitions/Integer/IntList.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Partitions/Integer/IntList.hs @@ -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 +-- +_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 +-- +-- | 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 + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Partitions/Integer/Naive.hs b/Math/Combinat/Partitions/Integer/Naive.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Partitions/Integer/Naive.hs @@ -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 +-- +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 +-- +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 + +-------------------------------------------------------------------------------- + + diff --git a/Math/Combinat/Partitions/Skew.hs b/Math/Combinat/Partitions/Skew.hs --- a/Math/Combinat/Partitions/Skew.hs +++ b/Math/Combinat/Partitions/Skew.hs @@ -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 diff --git a/Math/Combinat/Partitions/Skew/Ribbon.hs b/Math/Combinat/Partitions/Skew/Ribbon.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Partitions/Skew/Ribbon.hs @@ -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 + + +-------------------------------------------------------------------------------- diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs --- a/Math/Combinat/Permutations.hs +++ b/Math/Combinat/Permutations.hs @@ -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 diff --git a/Math/Combinat/Sets/VennDiagrams.hs b/Math/Combinat/Sets/VennDiagrams.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Sets/VennDiagrams.hs @@ -0,0 +1,150 @@ + +-- | Venn diagrams. See +-- +-- 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 +-} diff --git a/Math/Combinat/Sign.hs b/Math/Combinat/Sign.hs --- a/Math/Combinat/Sign.hs +++ b/Math/Combinat/Sign.hs @@ -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') diff --git a/Math/Combinat/Tableaux.hs b/Math/Combinat/Tableaux.hs --- a/Math/Combinat/Tableaux.hs +++ b/Math/Combinat/Tableaux.hs @@ -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 diff --git a/Math/Combinat/Tableaux/Skew.hs b/Math/Combinat/Tableaux/Skew.hs --- a/Math/Combinat/Tableaux/Skew.hs +++ b/Math/Combinat/Tableaux/Skew.hs @@ -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])]] diff --git a/cbits/c_compact_partition.c b/cbits/c_compact_partition.c new file mode 100644 --- /dev/null +++ b/cbits/c_compact_partition.c @@ -0,0 +1,24 @@ + +#include + +// ----------------------------------------------------------------------------- + +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> (64-4*(n-i))) + - (w >> (60-4*(n-i))) ) & 15 ; // diff + for(uint64_t j=0;j= 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 + +-------------------------------------------------------------------------------- + diff --git a/test/TestSuite.hs b/test/TestSuite.hs --- a/test/TestSuite.hs +++ b/test/TestSuite.hs @@ -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 ] -------------------------------------------------------------------------------- diff --git a/test/Tests/Common.hs b/test/Tests/Common.hs --- a/test/Tests/Common.hs +++ b/test/Tests/Common.hs @@ -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) diff --git a/test/Tests/Partitions/Compact.hs b/test/Tests/Partitions/Compact.hs new file mode 100644 --- /dev/null +++ b/test/Tests/Partitions/Compact.hs @@ -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 (ab) 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] + ] + +-------------------------------------------------------------------------------- diff --git a/test/Tests/Partitions/Integer.hs b/test/Tests/Partitions/Integer.hs --- a/test/Tests/Partitions/Integer.hs +++ b/test/Tests/Partitions/Integer.hs @@ -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) -------------------------------------------------------------------------------- diff --git a/test/Tests/Partitions/Ribbon.hs b/test/Tests/Partitions/Ribbon.hs new file mode 100644 --- /dev/null +++ b/test/Tests/Partitions/Ribbon.hs @@ -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) + +-------------------------------------------------------------------------------- diff --git a/test/Tests/Partitions/Skew.hs b/test/Tests/Partitions/Skew.hs --- a/test/Tests/Partitions/Skew.hs +++ b/test/Tests/Partitions/Skew.hs @@ -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 diff --git a/test/Tests/Permutations.hs b/test/Tests/Permutations.hs --- a/test/Tests/Permutations.hs +++ b/test/Tests/Permutations.hs @@ -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) -------------------------------------------------------------------------------- diff --git a/test/Tests/Series.hs b/test/Tests/Series.hs --- a/test/Tests/Series.hs +++ b/test/Tests/Series.hs @@ -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 diff --git a/test/Tests/SkewTableaux.hs b/test/Tests/SkewTableaux.hs --- a/test/Tests/SkewTableaux.hs +++ b/test/Tests/SkewTableaux.hs @@ -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