combinat 0.2.7.0 → 0.2.7.1
raw patch · 17 files changed
+1508/−298 lines, 17 filesdep ~containers
Dependency ranges changed: containers
Files
- LICENSE +1/−1
- Math/Combinat.hs +8/−5
- Math/Combinat/Helper.hs +14/−1
- Math/Combinat/Numbers.hs +4/−4
- Math/Combinat/Numbers/Series.hs +203/−3
- Math/Combinat/Partitions/Integer.hs +281/−46
- Math/Combinat/Partitions/Skew.hs +67/−0
- Math/Combinat/Permutations.hs +11/−5
- Math/Combinat/Sets.hs +29/−3
- Math/Combinat/Sign.hs +48/−0
- Math/Combinat/Tableaux.hs +6/−0
- Math/Combinat/Tableaux/GelfandTsetlin.hs +343/−0
- Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs +257/−0
- Math/Combinat/Tableaux/Kostka.hs +0/−222
- Math/Combinat/Tableaux/LittlewoodRichardson.hs +145/−0
- Math/Combinat/Tableaux/Skew.hs +78/−0
- combinat.cabal +13/−8
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2008-2014, Balazs Komuves+Copyright (c) 2008-2015, Balazs Komuves All rights reserved. Redistribution and use in source and binary forms, with or without
Math/Combinat.hs view
@@ -1,11 +1,11 @@ --- | A collection of functions to generate combinatorial--- objects like partitions, compositions, permutations,--- Young tableaux, various trees, etc etc.+-- | A collection of functions to generate and manipulate+-- combinatorial objects like partitions, compositions, +-- permutations, Young tableaux, various trees, etc etc. -- -- -- See also the @combinat-diagrams@ library for generating--- graphical representations of these structure using +-- graphical representations of (some of) these structure using -- the @diagrams@ library (<http://projects.haskell.org/diagrams>). -- --@@ -45,6 +45,7 @@ module Math.Combinat ( module Math.Combinat.Numbers+ , module Math.Combinat.Sign , module Math.Combinat.Sets , module Math.Combinat.Tuples , module Math.Combinat.Compositions@@ -54,9 +55,11 @@ , module Math.Combinat.Trees , module Math.Combinat.LatticePaths , module Math.Combinat.ASCII- ) where+ ) + where import Math.Combinat.Numbers+import Math.Combinat.Sign import Math.Combinat.Sets import Math.Combinat.Tuples import Math.Combinat.Compositions
Math/Combinat/Helper.hs view
@@ -187,5 +187,18 @@ return (s2, y:ys) --------------------------------------------------------------------------------- +-- * long zipwith ++longZipWith :: a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]+longZipWith a0 b0 f = go where+ go (x:xs) (y:ys) = f x y : go xs ys+ go [] ys = [ f a0 y | y<-ys ]+ go xs [] = [ f x b0 | x<-xs ]++{-+longZipWithZero :: (Num a, Num b) => (a -> b -> c) -> [a] -> [b] -> [c]+longZipWithZero = longZipWith 0 0 +-}++--------------------------------------------------------------------------------
Math/Combinat/Numbers.hs view
@@ -15,8 +15,8 @@ -------------------------------------------------------------------------------- -- | @(-1)^k@-paritySign :: Integral a => a -> Integer-paritySign k = if odd k then (-1) else 1+paritySignValue :: Integral a => a -> Integer+paritySignValue k = if odd k then (-1) else 1 -------------------------------------------------------------------------------- @@ -129,7 +129,7 @@ | k < 1 = 0 | k > n = 0 | otherwise = sum xs `div` factorial k where- xs = [ paritySign (k-i) * binomial k i * (fromIntegral i)^n | i<-[0..k] ]+ xs = [ paritySignValue (k-i) * binomial k i * (fromIntegral i)^n | i<-[0..k] ] -------------------------------------------------------------------------------- -- * Bernoulli numbers@@ -144,7 +144,7 @@ | n == 1 = -1/2 | otherwise = sum [ f k | k<-[1..n] ] where- f k = toRational (paritySign (n+k) * factorial k * stirling2nd n k) + f k = toRational (paritySignValue (n+k) * factorial k * stirling2nd n k) / toRational (k+1) --------------------------------------------------------------------------------
Math/Combinat/Numbers/Series.hs view
@@ -1,5 +1,5 @@ --- | Some basic power series expansions.+-- | Some basic univariate power series expansions. -- This module is not re-exported by "Math.Combinat". -- -- Note: the \"@convolveWithXXX@\" functions are much faster than the equivalent@@ -15,19 +15,66 @@ import Data.List +import Math.Combinat.Sign+import Math.Combinat.Numbers+import Math.Combinat.Partitions.Integer+import Math.Combinat.Helper+ #ifdef QUICKCHECK import System.Random import Test.QuickCheck #endif --------------------------------------------------------------------------------+-- * Trivial series -- | The series [1,0,0,0,0,...], which is the neutral element for the convolution. {-# SPECIALIZE unitSeries :: [Integer] #-} unitSeries :: Num a => [a] unitSeries = 1 : repeat 0 --- | Convolution of series. The result is always an infinite list. Warning: This is slow!+-- | Constant zero series+zeroSeries :: Num a => [a]+zeroSeries = repeat 0++-- | Power series representing a constant function+constSeries :: Num a => a -> [a]+constSeries x = x : repeat 0++-- | The power series representation of the identity function @x@+idSeries :: Num a => [a]+idSeries = 0 : 1 : repeat 0++-- | The power series representation of @x^n@+powerTerm :: Num a => Int -> [a]+powerTerm n = replicate n 0 ++ (1 : repeat 0)++--------------------------------------------------------------------------------+-- * Basic operations on power series++addSeries :: Num a => [a] -> [a] -> [a]+addSeries xs ys = longZipWith 0 0 (+) xs ys++subSeries :: Num a => [a] -> [a] -> [a]+subSeries xs ys = longZipWith 0 0 (-) xs ys++negateSeries :: Num a => [a] -> [a]+negateSeries = map negate++scaleSeries :: Num a => a -> [a] -> [a]+scaleSeries s = map (*s)++mulSeries :: Num a => [a] -> [a] -> [a]+mulSeries = convolve++productOfSeries :: Num a => [[a]] -> [a]+productOfSeries = convolveMany++--------------------------------------------------------------------------------+-- * Convolution (product)++-- | Convolution of series (that is, multiplication of power series). +-- The result is always an infinite list. Warning: This is slow! convolve :: Num a => [a] -> [a] -> [a] convolve xs1 ys1 = res where res = [ foldl' (+) 0 (zipWith (*) xs (reverse (take n ys)))@@ -36,12 +83,159 @@ xs = xs1 ++ repeat 0 ys = ys1 ++ repeat 0 --- | Convolution of many series. Still slow!+-- | Convolution (= product) of many series. Still slow! convolveMany :: Num a => [[a]] -> [a] convolveMany [] = 1 : repeat 0 convolveMany xss = foldl1 convolve xss --------------------------------------------------------------------------------+-- * Reciprocals of general power series++-- | Given a power series, we iteratively compute its multiplicative inverse+reciprocalSeries :: (Eq a, Fractional a) => [a] -> [a]+reciprocalSeries series = case series of+ [] -> error "reciprocalSeries: empty input series (const 0 function does not have an inverse)"+ (a:as) -> case a of+ 0 -> error "reciprocalSeries: input series has constant term 0"+ _ -> map (/a) $ integralReciprocalSeries $ map (/a) series++-- | Given a power series starting with @1@, we can compute its multiplicative inverse+-- without divisions.+--+{-# SPECIALIZE integralReciprocalSeries :: [Int] -> [Int] #-}+{-# SPECIALIZE integralReciprocalSeries :: [Integer] -> [Integer] #-}+integralReciprocalSeries :: (Eq a, Num a) => [a] -> [a]+integralReciprocalSeries series = case series of + [] -> error "integralReciprocalSeries: empty input series (const 0 function does not have an inverse)"+ (a:as) -> case a of+ 1 -> 1 : worker [1]+ _ -> error "integralReciprocalSeries: input series must start with 1"+ where+ worker bs = let b' = - sum (zipWith (*) (tail series) bs) + in b' : worker (b':bs)++--------------------------------------------------------------------------------+-- * Composition of formal power series++-- | @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.+composeSeries :: (Eq a, Num a) => [a] -> [a] -> [a]+composeSeries g f = substitute f g++-- | @substitute f g@ is the power series corresponding to @g(f(x))@. +-- Equivalently, this is the composition of univariate functions (in the \"wrong\" order).+--+-- 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_ = + 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"+ where+ as = as_ ++ repeat 0+ bs = bs_ ++ repeat 0+ a i = as !! i+ b j = bs !! j+ f n = sum+ [ b m * product [ (a i)^j | (i,j)<-es ] * fromInteger (multinomial (map snd es))+ | p <- partitions n + , let es = toExponentialForm p+ , let m = width p+ ]++--------------------------------------------------------------------------------+-- * Lagrange inversions++-- | Coefficients of the Lagrange inversion+lagrangeCoeff :: Partition -> Integer+lagrangeCoeff p = div numer denom where+ numer = (-1)^m * product (map fromIntegral [n+1..n+m])+ denom = fromIntegral (n+1) * product (map (factorial . snd) es)+ m = width p+ n = weight p+ es = toExponentialForm p++-- | We expect the input series to match @(0:1:_)@. The following is true for the result (at least with exact arithmetic):+--+-- > 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_ = + case series of+ (0:1:rest) -> 0 : 1 : [ f n | n<-[1..] ]+ _ -> error "integralLagrangeInversion: the series should start with (0 + x + a2*x^2 + ...)"+ where+ series = series_ ++ repeat 0+ as = tail series + a i = as !! i+ f n = sum [ fromInteger (lagrangeCoeff p) * product [ (a i)^j | (i,j) <- toExponentialForm p ]+ | 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_ = + 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"+ series = series_ ++ repeat 0+ a1 = series !! 1+ as = map (/a1) (tail series)+ a i = as !! i+ f n = sum [ fromInteger (lagrangeCoeff p) * product [ (a i)^j | (i,j) <- toExponentialForm p ]+ | p <- partitions n+ ] + +--------------------------------------------------------------------------------+-- * Power series expansions of elementary functions++-- | Power series expansion of @exp(x)@+expSeries :: Fractional a => [a]+expSeries = go 0 1 where+ go i e = e : go (i+1) (e / (i+1))++-- | Power series expansion of @cos(x)@+cosSeries :: Fractional a => [a]+cosSeries = go 0 1 where+ go i e = e : 0 : go (i+2) (-e / ((i+1)*(i+2)))++-- | Power series expansion of @sin(x)@+sinSeries :: Fractional a => [a]+sinSeries = go 1 1 where+ go i e = 0 : e : go (i+2) (-e / ((i+1)*(i+2)))++-- | Power series expansion of @cosh(x)@+coshSeries :: Fractional a => [a]+coshSeries = go 0 1 where+ go i e = e : 0 : go (i+2) (e / ((i+1)*(i+2)))++-- | Power series expansion of @sinh(x)@+sinhSeries :: Fractional a => [a]+sinhSeries = go 1 1 where+ go i e = 0 : e : go (i+2) (e / ((i+1)*(i+2)))++-- | Power series expansion of @log(1+x)@+log1Series :: Fractional a => [a]+log1Series = 0 : go 1 1 where+ go i e = (e/i) : go (i+1) (-e)++-- | Power series expansion of @(1-Sqrt[1-4x])/(2x)@ (the coefficients are the Catalan numbers)+dyckSeries :: Num a => [a]+dyckSeries = [ fromInteger (catalan i) | i<-[0..] ]++-------------------------------------------------------------------------------- -- * \"Coin\" series -- | Power series expansion of @@ -241,11 +435,13 @@ worker ((a,k):aks) ys = xs where xs = zipWith (+) (replicate k 0 ++ map (*a) ys) (worker aks ys) +{- data Sign = Plus | Minus deriving (Eq,Show) signValue :: Num a => Sign -> a signValue Plus = 1 signValue Minus = -1+-} signedPSeries :: [(Sign,Int)] -> [Integer] signedPSeries aks = convolveWithSignedPSeries aks unitSeries@@ -273,8 +469,12 @@ #ifdef QUICKCHECK +-- * Tests++{- swap :: (a,b) -> (b,a) swap (x,y) = (y,x)+-} -- compare the first 1000 elements of the infinite lists (=!=) :: (Eq a, Num a) => [a] -> [a] -> Bool
Math/Combinat/Partitions/Integer.hs view
@@ -17,13 +17,16 @@ -- <<svg/ferrers.svg>> -- -{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP, BangPatterns #-} module Math.Combinat.Partitions.Integer where -------------------------------------------------------------------------------- import Data.List +-- import Data.Map (Map)+-- import qualified Data.Map as Map+ import Math.Combinat.Helper import Math.Combinat.ASCII as ASCII import Math.Combinat.Numbers (factorial,binomial,multinomial)@@ -32,7 +35,7 @@ -- * 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.+-- are monotone decreasing sequences of /positive/ integers. The @Ord@ instance is lexicographical. newtype Partition = Partition [Int] deriving (Eq,Ord,Show,Read) ---------------------------------------------------------------------------------@@ -59,12 +62,12 @@ then toPartitionUnsafe xs else error "toPartition: not a partition" --- | Note: we only check that the sequence is ordered, but we /do not/ check for--- negative elements. This can be useful when working with symmetric functions.--- It may also change in the future...+-- | This returns @True@ if the input is non-increasing sequence of +-- /positive/ integers (possibly empty); @False@ otherwise.+-- isPartition :: [Int] -> Bool isPartition [] = True-isPartition [_] = True+isPartition [x] = x > 0 isPartition (x:xs@(y:_)) = (x >= y) && isPartition xs fromPartition :: Partition -> [Int]@@ -76,7 +79,7 @@ (p:_) -> p [] -> 0 --- | The length of the sequence.+-- | The length of the sequence (that is, the number of parts). width :: Partition -> Int width (Partition part) = length part @@ -92,11 +95,45 @@ dualPartition :: Partition -> Partition dualPartition (Partition part) = Partition (_dualPartition part) --- (we could be more efficient, but it hardly matters)+data Pair = Pair !Int !Int+ _dualPartition :: [Int] -> [Int] _dualPartition [] = []-_dualPartition xs@(k:_) = [ length $ filter (>=i) xs | i <- [1..k] ]+_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 (toPartition [5,4,1]) ==@@ -112,6 +149,25 @@ _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 :: 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.@@ -122,22 +178,51 @@ _countAutomorphisms = multinomial . map length . group ------------------------------------------------------------------------------------ * Dominance order +-- * Generating partitions --- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions--- (this is partial ordering on the set of partitions of @n@).------ See <http://en.wikipedia.org/wiki/Dominance_order>----dominates :: Partition -> Partition -> Bool-dominates (Partition qs) (Partition ps) - = and $ zipWith (>=) (sums (qs ++ repeat 0)) (sums ps)- where- sums = scanl (+) 0+-- | Partitions of @d@.+partitions :: Int -> [Partition]+partitions = map Partition . _partitions ------------------------------------------------------------------------------------- * 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) ] +countPartitions :: Int -> Integer+countPartitions d = countPartitions' (d,d) d++-- | 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] ]++-- | 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 -> [[Partition]]+allPartitionsGrouped d = [ partitions i | i <- [0..d] ]++-- | All integer partitions fitting into a given rectangle.+allPartitions' + :: (Int,Int) -- ^ (height,width)+ -> [Partition]+allPartitions' (h,w) = concat [ partitions' (h,w) i | i <- [0..d] ] where d = h*w++-- | All integer partitions fitting into a given rectangle, grouped by weight.+allPartitionsGrouped' + :: (Int,Int) -- ^ (height,width)+ -> [[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)@@ -163,36 +248,65 @@ countPartitions' (h,w) d = sum [ countPartitions' (i,w-1) (d-i) | i <- [1..min d h] ] --- | Partitions of @d@, as lists-_partitions :: Int -> [[Int]]-_partitions d = _partitions' (d,d) d+---------------------------------------------------------------------------------+-- * Dominance order --- | Partitions of @d@.-partitions :: Int -> [Partition]-partitions d = partitions' (d,d) d+-- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions+-- (this is partial ordering on the set of partitions of @n@).+--+-- See <http://en.wikipedia.org/wiki/Dominance_order>+--+dominates :: Partition -> Partition -> Bool+dominates (Partition qs) (Partition ps) + = and $ zipWith (>=) (sums (qs ++ repeat 0)) (sums ps)+ where+ sums = scanl (+) 0 -countPartitions :: Int -> Integer-countPartitions d = countPartitions' (d,d) d --- | All integer partitions fitting into a given rectangle.-allPartitions' - :: (Int,Int) -- ^ (height,width)- -> [[Partition]]-allPartitions' (h,w) = [ partitions' (h,w) i | i <- [0..d] ] where d = h*w+-- | 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 :: Partition -> [Partition] +dominatedPartitions (Partition lambda) = map Partition (_dominatedPartitions lambda) --- | 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 = [ partitions i | i <- [0..d] ]+_dominatedPartitions :: [Int] -> [[Int]]+_dominatedPartitions [] = [[]]+_dominatedPartitions lambda = go (head lambda) w dsums 0 where --- | # = \\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+ n = length lambda+ w = sum lambda+ dsums = scanl1 (+) (lambda ++ repeat 0) -countAllPartitions :: Int -> Integer-countAllPartitions d = sum' [ countPartitions i | i <- [0..d] ]+ 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 :: 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"+ -------------------------------------------------------------------------------- -- * Partitions with given number of parts @@ -229,7 +343,37 @@ | 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) ]++{-+-- | 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 -> [Partition]+partitionsWithEvenParts d = map Partition (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 -> [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) ]+ -------------------------------------------------------------------------------- -- * Sub-partitions of a given partition @@ -238,7 +382,6 @@ isSubPartitionOf :: Partition -> Partition -> Bool isSubPartitionOf (Partition ps) (Partition qs) = and $ zipWith (<=) ps (qs ++ repeat 0) ----------------------------------------- -- | Sub-partitions of a given partition with the given weight: --@@ -279,6 +422,58 @@ | otherwise = [] : [ this:rest | this <- [1..min h b] , rest <- go this bs ] --------------------------------------------------------------------------------+-- * The Pieri rule++-- | The Pieri rule computes @s[lambda]*h[n]@ as a sum of @s[mu]@-s (each with coefficient 1).+--+-- See for example <http://en.wikipedia.org/wiki/Pieri's_formula>+--+pieriRule :: Partition -> Int -> [Partition] +pieriRule (Partition lambda) n = map Partition (_pieriRule lambda n) where++ -- | We assume here that @lambda@ is a partition (non-increasing sequence of /positive/ integers)! + _pieriRule :: [Int] -> Int -> [[Int]] + _pieriRule lambda n+ | n == 0 = [lambda]+ | n < 0 = [] + | otherwise = go n diffs dsums (lambda++[0]) + where+ diffs = n : diffSequence lambda -- maximum we can add to a given row+ dsums = reverse $ scanl1 (+) (reverse diffs) -- partial sums of remaining total we can add+ go !k (d:ds) (p:ps@(q:_)) (l:ls) + | k > p = []+ | otherwise = [ h:tl | a <- [ max 0 (k-q) .. min d k ] , let h = l+a , tl <- go (k-a) ds ps ls ]+ go !k [d] _ [l] = if k <= d + then if l+k>0 then [[l+k]] else [[]]+ else []+ go !k [] _ _ = if k==0 then [[]] else []++-- | The dual Pieri rule computes @s[lambda]*e[n]@ as a sum of @s[mu]@-s (each with coefficient 1)+dualPieriRule :: Partition -> Int -> [Partition] +dualPieriRule lam n = map dualPartition $ pieriRule (dualPartition lam) n++{-+-- | 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 -- | Which orientation to draw the Ferrers diagrams.@@ -330,5 +525,45 @@ EnglishNotation -> fromPartition part EnglishNotationCCW -> reverse $ fromPartition $ dualPartition part FrenchNotation -> reverse $ fromPartition $ part++--------------------------------------------------------------------------------++{-++-- * Tests++-- we need some custom types for quickcheck generating not too large partitions...++newtype PartitionWeight = PartitionWeight Int+data PartitionWeightPair = PartitionWeightPair Int Int +data PartitionPair = PartitionPair Partition Int ++prop_partitions_in_bigbox :: PartitionWeight -> Bool+prop_partitions_in_bigbox (PartitionWeight n) = (partitions n == partitions' (n,n) n)++prop_kparts :: PartitionWeightPair -> Bool+prop_kparts (PartitionWeightPair n k) = (partitionsWithKParts k n == [ mu | mu <- partitions n, numberOfParts mu == k ])++prop_odd_partitions :: PartitionWeight -> Bool+prop_odd_partitions (PartitionWeight n) = + (partitionsWithOddParts n == [ mu | mu <- partitions n, and (map odd (fromPartition mu)) ])++prop_distinct_partitions :: PartitionWeight -> Bool+prop_distinct_partitions (PartitionWeight n) = + (partitionsWithDistinctParts n == [ mu | mu <- partitions n, let xs = fromPartition mu, xs == nub xs ])++prop_subparts :: PartitionPair -> Bool+prop_subparts (PartitionPair lam d) = (subPartitions d lam) == sort [ p | p <- partitions d, isSubPartitionOf p lam ]++prop_dual_dual :: Partition -> Bool+prop_dual_dual lam = (lam == dualPartition (dualPartition lam))++prop_dominated_list :: Partition -> Bool+prop_dominated_list lam = (dominatedPartitions lam == [ mu | mu <- partitions (weight lam), lam `dominates` mu ])++prop_dominating_list :: Partition -> Bool+prop_dominating_list mu = (dominatingPartitions mu == [ lam | lam <- partitions (weight mu ), lam `dominates` mu ])++-} --------------------------------------------------------------------------------
+ Math/Combinat/Partitions/Skew.hs view
@@ -0,0 +1,67 @@+ +-- | Skew partitions. +-- +-- Skew partitions are the difference of two integer partitions, denoted by @lambda/mu@. +-- + +{-# LANGUAGE BangPatterns #-} +module Math.Combinat.Partitions.Skew where + +-------------------------------------------------------------------------------- + +import Data.List + +import Math.Combinat.Partitions.Integer +import Math.Combinat.ASCII + +-------------------------------------------------------------------------------- + +-- | A skew partition @lambda/mu@ is represented by the list @[ (mu_i , lambda_i-mu_i) | i<-[1..n] ]@ +newtype SkewPartition = SkewPartition [(Int,Int)] deriving (Eq,Ord,Show) + +-- | @mkSkewPartition (lambda,mu)@ creates the skew partition @lambda/mu@ +mkSkewPartition :: (Partition,Partition) -> SkewPartition +mkSkewPartition ( lam@(Partition bs) , mu@(Partition as)) = if mu `isSubPartitionOf` lam + then SkewPartition $ zipWith (\b a -> (a,b-a)) bs (as ++ repeat 0) + else error "mkSkewPartition: mu should be a subpartition of lambda!" + +-- | This function \"cuts off\" the \"uninteresting parts\" of a skew partition +normalizeSkewPartition :: SkewPartition -> SkewPartition +normalizeSkewPartition (SkewPartition abs) = SkewPartition abs' where + (as,bs) = unzip abs + a0 = minimum as + k = length (takeWhile (==0) bs) + abs' = zip [ a-a0 | a <- drop k as ] (drop k bs) + +-- | Returns the outer and inner partition of a skew partition, respectively +fromSkewPartition :: SkewPartition -> (Partition,Partition) +fromSkewPartition (SkewPartition list) = (toPartition (zipWith (+) as bs) , toPartition (filter (>0) as)) where + (as,bs) = unzip list + +outerPartition :: SkewPartition -> Partition +outerPartition = fst . fromSkewPartition + +innerPartition :: SkewPartition -> Partition +innerPartition = snd . fromSkewPartition + +dualSkewPartition :: SkewPartition -> SkewPartition +dualSkewPartition = mkSkewPartition . f . fromSkewPartition where + f (lam,mu) = (dualPartition lam, dualPartition mu) + +asciiSkewFerrersDiagram :: SkewPartition -> ASCII +asciiSkewFerrersDiagram = asciiSkewFerrersDiagram' ('@','.') EnglishNotation + +asciiSkewFerrersDiagram' + :: (Char,Char) + -> PartitionConvention -- Orientation + -> SkewPartition + -> ASCII +asciiSkewFerrersDiagram' (outer,inner) orient (SkewPartition abs) = asciiFromLines stuff where + stuff = case orient of + EnglishNotation -> ls + EnglishNotationCCW -> reverse (transpose ls) + FrenchNotation -> reverse ls + ls = [ replicate a inner ++ replicate b outer | (a,b) <- abs ] + + +--------------------------------------------------------------------------------
Math/Combinat/Permutations.hs view
@@ -22,7 +22,7 @@ , disjointCyclesToPermutation , isEvenPermutation , isOddPermutation- , signOfPermutation+ , signOfPermutation -- , Sign(..) , isCyclicPermutation -- * Permutation groups , permute@@ -55,6 +55,8 @@ ) where +--------------------------------------------------------------------------------+ import Control.Monad import Control.Monad.ST @@ -69,6 +71,7 @@ import Data.Array.Unsafe import Math.Combinat.Helper+import Math.Combinat.Sign import Math.Combinat.Numbers (factorial,binomial) import System.Random@@ -212,11 +215,14 @@ isOddPermutation :: Permutation -> Bool isOddPermutation = not . isEvenPermutation --- | Plus 1 or minus 1.-signOfPermutation :: Num a => Permutation -> a+signOfPermutation :: Permutation -> Sign signOfPermutation perm = case isEvenPermutation perm of- True -> 1- False -> (-1)+ True -> Plus+ False -> Minus++-- | Plus 1 or minus 1.+signValueOfPermutation :: Num a => Permutation -> a+signValueOfPermutation = signValue . signOfPermutation isCyclicPermutation :: Permutation -> Bool isCyclicPermutation perm =
Math/Combinat/Sets.hs view
@@ -1,6 +1,7 @@ -- | Subsets. +{-# LANGUAGE BangPatterns, Rank2Types #-} module Math.Combinat.Sets ( -- * choices@@ -24,6 +25,11 @@ import Data.List ( sort , mapAccumL ) import System.Random +import Control.Monad+import Control.Monad.ST+import Data.Array.ST+import Data.Array.MArray+ -- import Data.Map (Map) -- import qualified Data.Map as Map @@ -122,15 +128,35 @@ randomChoice k n g0 = if k>n || k<0 then error "randomChoice: k out of range" - else (makeChoiceFromIndices as, g1) + else (makeChoiceFromIndicesKnuth n as, g1) where -- choose numbers from [1..n], [1..n-1], [1..n-2] etc (g1,as) = mapAccumL (\g m -> swap (randomR (1,m) g)) g0 [n,n-1..n-k+1] ++-------------------------------------------------------------------------------- -- | From a list of $k$ numbers, where the first is in the interval @[1..n]@, -- the second in @[1..n-1]@, the third in @[1..n-2]@, we create a size @k@ subset of @n@.-makeChoiceFromIndices :: [Int] -> [Int]-makeChoiceFromIndices = sort . go [] where+--+-- Knuth's method. The first argument is the number @n@.+--+makeChoiceFromIndicesKnuth :: Int -> [Int] -> [Int]+makeChoiceFromIndicesKnuth n xs = + runST $ do+ arr <- newArray_ (1,n) :: forall s. ST s (STUArray s Int Int)+ forM_ [1..n] $ \(!j) -> writeArray arr j j+ forM_ (zip [n,n-1..] xs) $ \(!j,!i) -> do+ a <- readArray arr j+ b <- readArray arr i+ writeArray arr j b+ writeArray arr i a+ sel <- forM (zip [n,n-1..] xs) $ \(!j,_) -> readArray arr j+ return (sort sel)++-- | From a list of $k$ numbers, where the first is in the interval @[1..n]@, +-- the second in @[1..n-1]@, the third in @[1..n-2]@, we create a size @k@ subset of @n@.+makeChoiceFromIndicesNaive :: [Int] -> [Int]+makeChoiceFromIndicesNaive = sort . go [] where go :: [Int] -> [Int] -> [Int] go acc (b:bs) = b' : go (insert b' acc) bs where b' = skip b acc
+ Math/Combinat/Sign.hs view
@@ -0,0 +1,48 @@+ +-- | Signs + +{-# LANGUAGE BangPatterns #-} +module Math.Combinat.Sign where + +-------------------------------------------------------------------------------- + +import Data.Monoid + +-------------------------------------------------------------------------------- + +data Sign + = Plus + | Minus + deriving (Eq,Ord,Show,Read) + +instance Monoid Sign where + mempty = Plus + mappend = mulSign + mconcat = productOfSigns + +signValue :: Num a => Sign -> a +signValue s = case s of + Plus -> 1 + Minus -> -1 + +paritySign :: Integral a => a -> Sign +paritySign x = if even x then Plus else Minus + +oppositeSign :: Sign -> Sign +oppositeSign s = case s of + Plus -> Minus + Minus -> Plus + +mulSign :: Sign -> Sign -> Sign +mulSign s1 s2 = case s1 of + Plus -> s2 + Minus -> oppositeSign s2 + +productOfSigns :: [Sign] -> Sign +productOfSigns = go Plus where + go !acc [] = acc + go !acc (x:xs) = case x of + Plus -> go acc xs + Minus -> go (oppositeSign acc) xs + +--------------------------------------------------------------------------------
Math/Combinat/Tableaux.hs view
@@ -29,11 +29,17 @@ import Math.Combinat.Helper import Math.Combinat.Numbers (factorial,binomial) import Math.Combinat.Partitions+import Math.Combinat.ASCII -------------------------------------------------------------------------------- -- * Basic stuff type Tableau a = [[a]]++asciiTableau :: Show a => Tableau a -> ASCII+asciiTableau t = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) + $ (map . map) asciiShow+ $ t _shape :: Tableau a -> [Int] _shape t = map length t
+ Math/Combinat/Tableaux/GelfandTsetlin.hs view
@@ -0,0 +1,343 @@+ +-- | Gelfand-Tsetlin patterns and Kostka numbers. +-- +-- Gelfand-Tsetlin patterns (or tableaux) are triangular arrays like +-- +-- > [ 3 ] +-- > [ 3 , 2 ] +-- > [ 3 , 1 , 0 ] +-- > [ 2 , 0 , 0 , 0 ] +-- +-- with both rows and columns non-increasing non-negative integers. +-- Note: these are in bijection with the semi-standard Young tableaux. +-- +-- If we add the further restriction that +-- the top diagonal reads @lambda@, +-- and the diagonal sums are partial sums of @mu@, where @lambda@ and @mu@ are two +-- partitions (in this case @lambda=[3,2]@ and @mu=[2,1,1,1]@), +-- then the number of the resulting patterns +-- or tableaux is the Kostka number @K(lambda,mu)@. +-- Actually @mu@ doesn't even need to the be non-increasing. +-- + +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} +module Math.Combinat.Tableaux.GelfandTsetlin where + +-------------------------------------------------------------------------------- + +import Data.List +import Data.Maybe +import Data.Monoid +import Data.Ord + +import Control.Monad +import Control.Monad.Trans.State + +import Data.Map (Map) +import qualified Data.Map as Map + +import Math.Combinat.Partitions.Integer +import Math.Combinat.Tableaux +import Math.Combinat.Helper +import Math.Combinat.ASCII + +-------------------------------------------------------------------------------- +-- * Kostka numbers + +-- | Kostka numbers (via counting Gelfand-Tsetlin patterns). See for example <http://en.wikipedia.org/wiki/Kostka_number> +-- +-- @K(lambda,mu)==0@ unless @lambda@ dominates @mu@: +-- +-- > [ mu | mu <- partitions (weight lam) , kostkaNumber lam mu > 0 ] == dominatedPartitions lam +-- +kostkaNumber :: Partition -> Partition -> Int +kostkaNumber = countKostkaGelfandTsetlinPatterns + +-- | Very naive (and slow) implementation of Kostka numbers, for reference. +kostkaNumberReferenceNaive :: Partition -> Partition -> Int +kostkaNumberReferenceNaive plambda pmu@(Partition mu) = length stuff where + stuff = [ 1 | t <- semiStandardYoungTableaux k plambda , cond t ] + k = length mu + cond t = [ (head xs, length xs) | xs <- group (sort $ concat t) ] == zip [1..] mu + +-------------------------------------------------------------------------------- + +-- | Lists all (positive) Kostka numbers @K(lambda,mu)@ with the given @lambda@: +-- +-- > kostkaNumbersWithGivenLambda lambda == Map.fromList [ (mu , kostkaNumber lambda mu) | mu <- dominatedPartitions lambda ] +-- +-- It's much faster than computing the individual Kostka numbers, but not as fast +-- as it could be. +-- +{-# SPECIALIZE kostkaNumbersWithGivenLambda :: Partition -> Map Partition Int #-} +{-# SPECIALIZE kostkaNumbersWithGivenLambda :: Partition -> Map Partition Integer #-} +kostkaNumbersWithGivenLambda :: forall coeff. Num coeff => Partition -> Map Partition coeff +kostkaNumbersWithGivenLambda plambda@(Partition lam) = evalState (worker lam) Map.empty where + + worker :: [Int] -> State (Map Partition (Map Partition coeff)) (Map Partition coeff) + worker unlam = case unlam of + [] -> return $ Map.singleton (Partition []) 1 + _ -> do + cache <- get + case Map.lookup (Partition unlam) cache of + Just sol -> return sol + Nothing -> do + let s = foldl' (+) 0 unlam + subsols <- forM (prevLambdas0 unlam) $ \p -> do + sub <- worker p + let t = s - foldl' (+) 0 p + f (Partition xs , c) = case xs of + (y:_) -> if t >= y then Just (Partition (t:xs) , c) else Nothing + [] -> if t > 0 then Just (Partition [t] , c) else Nothing + if t > 0 + then return $ Map.fromList $ mapMaybe f $ Map.toList sub + else return $ Map.empty + + let sol = Map.unionsWith (+) subsols + put $! (Map.insert (Partition unlam) sol cache) + return sol + + -- needs decreasing sequence + prevLambdas0 :: [Int] -> [[Int]] + prevLambdas0 (l:ls) = go l ls where + go b [a] = [ [x] | x <- [a..b] ] ++ [ [x,y] | x <- [a..b] , y<-[1..a] ] + go b (a:as) = [ x:xs | x <- [a..b] , xs <- go a as ] + go b [] = [] : [ [j] | j <- [1..b] ] + prevLambdas0 [] = [] + +-- | Lists all (positive) Kostka numbers @K(lambda,mu)@ with the given @mu@: +-- +-- > kostkaNumbersWithGivenMu mu == Map.fromList [ (lambda , kostkaNumber lambda mu) | lambda <- dominatingPartitions mu ] +-- +-- This function uses the iterated Pieri rule, and is relatively fast. +-- +kostkaNumbersWithGivenMu :: Partition -> Map Partition Int +kostkaNumbersWithGivenMu (Partition mu) = iteratedPieriRule (reverse mu) + +-------------------------------------------------------------------------------- +-- * Gelfand-Tsetlin patterns + +-- | A Gelfand-Tstetlin tableau +type GT = [[Int]] + + +asciiGT :: GT -> ASCII +asciiGT gt = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) + $ (map . map) asciiShow + $ gt + + +kostkaGelfandTsetlinPatterns :: Partition -> Partition -> [GT] +kostkaGelfandTsetlinPatterns lambda (Partition mu) = kostkaGelfandTsetlinPatterns' lambda mu + +-- | Generates all Kostka-Gelfand-Tsetlin tableau, that is, triangular arrays like +-- +-- > [ 3 ] +-- > [ 3 , 2 ] +-- > [ 3 , 1 , 0 ] +-- > [ 2 , 0 , 0 , 0 ] +-- +-- with both rows and column non-increasing such that +-- the top diagonal read lambda (in this case @lambda=[3,2]@) and the diagonal sums +-- are partial sums of mu (in this case @mu=[2,1,1,1]@) +-- +-- The number of such GT tableaux is the Kostka +-- number K(lambda,mu). +-- +kostkaGelfandTsetlinPatterns' :: Partition -> [Int] -> [GT] +kostkaGelfandTsetlinPatterns' plam@(Partition lambda0) mu0 + | minimum mu0 < 0 = [] + | wlam == 0 = if wmu == 0 then [ [] ] else [] + | wmu == wlam && plam `dominates` pmu = list + | otherwise = [] + where + + pmu = mkPartition mu0 + + nlam = length lambda0 + nmu = length mu0 + + n = max nlam nmu + + lambda = lambda0 ++ replicate (n - nlam) 0 + mu = mu0 ++ replicate (n - nmu ) 0 + + revlam = reverse lambda + + wmu = sum' mu + wlam = sum' lambda + + list = worker + revlam + (scanl1 (+) mu) + (replicate (n-1) 0) + (replicate (n ) 0) + [] + + worker + :: [Int] -- lambda_i in reverse order + -> [Int] -- partial sums of mu + -> [Int] -- sums of the tails of previous rows + -> [Int] -- last row + -> [[Int]] -- the lower part of GT tableau we accumulated so far (this is not needed if we only want to count) + -> [GT] + + worker (rl:rls) (smu:smus) (a:acc) (lastx0:lastrowt) table = stuff + where + x0 = smu - a + stuff = concat + [ worker rls smus (zipWith (+) acc (tail row)) (init row) (row:table) + | row <- boundedNonIncrSeqs' x0 (map (max rl) (max lastx0 x0 : lastrowt)) lambda + ] + worker [rl] _ _ _ table = [ [rl]:table ] + worker [] _ _ _ _ = [ [] ] + + boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]] + boundedNonIncrSeqs' = go where + go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ] + go _ [] _ = [[]] + go _ _ [] = [[]] + +-------------------------------------------------------------------------------- + +-- | This returns the corresponding Kostka number: +-- +-- > countKostkaGelfandTsetlinPatterns lambda mu == length (kostkaGelfandTsetlinPatterns lambda mu) +-- +countKostkaGelfandTsetlinPatterns :: Partition -> Partition -> Int +countKostkaGelfandTsetlinPatterns plam@(Partition lambda0) pmu@(Partition mu0) + | wlam == 0 = if wmu == 0 then 1 else 0 + | wmu == wlam && plam `dominates` pmu = cnt + | otherwise = 0 + where + + nlam = length lambda0 + nmu = length mu0 + + n = max nlam nmu + + lambda = lambda0 ++ replicate (n - nlam) 0 + mu = mu0 ++ replicate (n - nmu ) 0 + + revlam = reverse lambda + + wmu = sum' mu + wlam = sum' lambda + + cnt = worker + revlam + (scanl1 (+) mu) + (replicate (n-1) 0) + (replicate (n ) 0) + + worker + :: [Int] -- lambda_i in reverse order + -> [Int] -- partial sums of mu + -> [Int] -- sums of the tails of previous rows + -> [Int] -- last row + -> Int + + worker (rl:rls) (smu:smus) (a:acc) (lastx0:lastrowt) = stuff + where + x0 = smu - a + stuff = sum' + [ worker rls smus (zipWith (+) acc (tail row)) (init row) + | row <- boundedNonIncrSeqs' x0 (map (max rl) (max lastx0 x0 : lastrowt)) lambda + ] + worker [rl] _ _ _ = 1 + worker [] _ _ _ = 1 + + boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]] + boundedNonIncrSeqs' = go where + go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ] + go _ [] _ = [[]] + go _ _ [] = [[]] + +-------------------------------------------------------------------------------- + +{- + +-- | All non-increasing sentences between a lower and an upper bound +boundedNonIncrSeqs :: [Int] -> [Int] -> [[Int]] +boundedNonIncrSeqs as bs = case bs of + (h0:_) -> boundedNonIncrSeqs' h0 as bs + [] -> [[]] + +-- | All non-increasing sentences between a lower and an upper bound, and also less-or-equal than the given number +boundedNonIncrSeqs' :: Int -> [Int] -> [Int] -> [[Int]] +boundedNonIncrSeqs' = go where + go h0 (a:as) (b:bs) = [ h:hs | h <- [(max 0 a)..(min h0 b)] , hs <- go h as bs ] + go _ [] _ = [[]] + go _ _ [] = [[]] + +-- | All non-decreasing sentences between a lower and an upper bound +boundedNonDecrSeqs :: [Int] -> [Int] -> [[Int]] +boundedNonDecrSeqs = boundedNonDecrSeqs' 0 + +-- | All non-decreasing sentences between a lower and an upper bound, and also greator-or-equal then the given number +boundedNonDecrSeqs' :: Int -> [Int] -> [Int] -> [[Int]] +boundedNonDecrSeqs' h0 = go (max 0 h0) where + go h0 (a:as) (b:bs) = [ h:hs | h <- [(max h0 a)..b] , hs <- go h as bs ] + go _ [] _ = [[]] + go _ _ [] = [[]] + +-} + +-------------------------------------------------------------------------------- +-- * The iterated Pieri rule + +-- | Computes the Schur expansion of @h[n1]*h[n2]*h[n3]*...*h[nk]@ via iterating the Pieri rule. +-- Note: the coefficients are actually the Kostka numbers; the following is true: +-- +-- > Map.toList (iteratedPieriRule (fromPartition mu)) == [ (lam, kostkaNumber lam mu) | lam <- dominatingPartitions mu ] +-- +-- This should be faster than individually computing all these Kostka numbers. +-- +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 + +-------------------------------------------------------------------------------- + +-- | Computes the Schur expansion of @e[n1]*e[n2]*e[n3]*...*e[nk]@ via iterating the Pieri rule. +-- Note: the coefficients are actually the Kostka numbers; the following is true: +-- +-- > Map.toList (iteratedDualPieriRule (fromPartition mu)) == +-- > [ (dualPartition lam, kostkaNumber lam mu) | lam <- dominatingPartitions mu ] +-- +-- This should be faster than individually computing all these Kostka numbers. +-- It is a tiny bit slower than 'iteratedPieriRule'. +-- +iteratedDualPieriRule :: Num coeff => [Int] -> Map Partition coeff +iteratedDualPieriRule = iteratedDualPieriRule' (Partition []) + +-- | Iterating the Pieri rule, we can compute the Schur expansion of +-- @e[lambda]*e[n1]*e[n2]*e[n3]*...*e[nk]@ +iteratedDualPieriRule' :: Num coeff => Partition -> [Int] -> Map Partition coeff +iteratedDualPieriRule' plambda ns = iteratedDualPieriRule'' (plambda,1) ns + +{-# SPECIALIZE iteratedDualPieriRule'' :: (Partition,Int ) -> [Int] -> Map Partition Int #-} +{-# SPECIALIZE iteratedDualPieriRule'' :: (Partition,Integer) -> [Int] -> Map Partition Integer #-} +iteratedDualPieriRule'' :: Num coeff => (Partition,coeff) -> [Int] -> Map Partition coeff +iteratedDualPieriRule'' (plambda,coeff0) ns = worker (Map.singleton plambda coeff0) ns where + worker old [] = old + worker old (n:ns) = worker new ns where + stuff = [ (coeff, dualPieriRule 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 + +--------------------------------------------------------------------------------
+ Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs view
@@ -0,0 +1,257 @@++-- TODO: better name?++-- | This module contains a function to generate (equivalence classes of) +-- triangular tableaux of size /k/, strictly increasing to the right and +-- to the bottom. For example+-- +-- > 1 +-- > 2 4 +-- > 3 5 8 +-- > 6 7 9 10 +--+-- is such a tableau of size 4.+-- The numbers filling a tableau always consist of an interval @[1..c]@;+-- @c@ is called the /content/ of the tableaux. There is a unique tableau+-- of minimal content @2k-1@:+--+-- > 1 +-- > 2 3 +-- > 3 4 5 +-- > 4 5 6 7 +-- +-- Let us call the tableaux with maximal content (that is, @m = binomial (k+1) 2@)+-- /standard/. The number of such standard tableaux are+--+-- > 1, 1, 2, 12, 286, 33592, 23178480, ...+--+-- OEIS:A003121, \"Strict sense ballot numbers\", +-- <https://oeis.org/A003121>.+--+-- See +-- R. M. Thrall, A combinatorial problem, Michigan Math. J. 1, (1952), 81-88.+-- +-- The number of tableaux with content @c=m-d@ are+-- +-- > d= | 0 1 2 3 ...+-- > -----+----------------------------------------------+-- > k=2 | 1+-- > k=3 | 2 1+-- > k=4 | 12 18 8 1+-- > k=5 | 286 858 1001 572 165 22 1+-- > k=6 | 33592 167960 361114 436696 326196 155584 47320 8892 962 52 1 +--+-- We call these \"GT simplex tableaux\" (in the lack of a better name), since+-- they are in bijection with the simplicial cones in a canonical simplicial +-- decompositions of the Gelfand-Tsetlin cones (the content corresponds+-- to the dimension), which encode the combinatorics of Kostka numbers.+--++module Math.Combinat.Tableaux.GelfandTsetlin.Cone+ ( + -- * Types+ Tableau+ , Tri(..)+ , TriangularArray+ , fromTriangularArray+ , triangularArrayUnsafe+ -- * ASCII+ , asciiTriangularArray+ , asciiTableau+ -- * Content+ , gtSimplexContent+ , _gtSimplexContent+ , invertGTSimplexTableau+ , _invertGTSimplexTableau+ -- * Enumeration+ , gtSimplexTableaux+ , _gtSimplexTableaux+ , countGTSimplexTableaux+ ) + where++--------------------------------------------------------------------------------++import Data.Ix+import Data.Ord+import Data.List++import Control.Monad+import Control.Monad.ST+import Data.Array.IArray+import Data.Array.Unboxed+import Data.Array.ST++import Math.Combinat.Tableaux (Tableau)+import Math.Combinat.Helper+import Math.Combinat.ASCII++--------------------------------------------------------------------------------++-- | Triangular arrays+type TriangularArray a = Array Tri a++-- | Set of @(i,j)@ pairs with @i>=j>=1@.+newtype Tri = Tri { unTri :: (Int,Int) } deriving (Eq,Ord,Show)++binom2 :: Int -> Int+binom2 n = (n*(n-1)) `div` 2++index' :: Tri -> Int+index' (Tri (i,j)) = binom2 i + j - 1++-- it should be (1+8*m), +-- the 2 is a hack to be safe with the floating point stuff+deIndex' :: Int -> Tri +deIndex' m = Tri ( i+1 , m - binom2 (i+1) + 1 ) where+ i = ( (floor.sqrt.(fromIntegral::Int->Double)) (2+8*m) - 1 ) `div` 2 ++instance Ix Tri where+ index (a,b) x = index' x - index' a + inRange (a,b) x = (u<=j && j<=v) where+ u = index' a + v = index' b+ j = index' x+ range (a,b) = map deIndex' [ index' a .. index' b ] + rangeSize (a,b) = index' b - index' a + 1 ++{-# SPECIALIZE triangularArrayUnsafe :: Tableau Int -> TriangularArray Int #-}+triangularArrayUnsafe :: Tableau a -> TriangularArray a+triangularArrayUnsafe tableau = listArray (Tri (1,1),Tri (k,k)) (concat tableau) + where k = length tableau++{-# SPECIALIZE fromTriangularArray :: TriangularArray Int -> Tableau Int #-}+fromTriangularArray :: TriangularArray a -> Tableau a+fromTriangularArray arr = (map.map) snd $ groupBy (equating f) $ assocs arr+ where f = fst . unTri . fst++--------------------------------------------------------------------------------++asciiTriangularArray :: Show a => TriangularArray a -> ASCII+asciiTriangularArray = asciiTableau . fromTriangularArray++asciiTableau :: Show a => Tableau a -> ASCII+asciiTableau xxs = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) + $ (map . map) asciiShow+ $ xxs++--------------------------------------------------------------------------------++-- "fractional fillings"+data Hole = Hole Int Int deriving (Eq,Ord,Show)++type ReverseTableau = [[Int ]] +type ReverseHoleTableau = [[Hole]] ++toHole :: Int -> Hole+toHole k = Hole k 0++nextHole :: Hole -> Hole+nextHole (Hole k l) = Hole k (l+1)++{-# SPECIALIZE reverseTableau :: [[Int]] -> [[Int]] #-}+reverseTableau :: [[a]] -> [[a]]+reverseTableau = reverse . map reverse++--------------------------------------------------------------------------------++gtSimplexContent :: TriangularArray Int -> Int+gtSimplexContent arr = max (arr ! (fst (bounds arr))) (arr ! (snd (bounds arr))) -- we also handle inverted tableau++_gtSimplexContent :: Tableau Int -> Int+_gtSimplexContent t = max (head $ head t) (last $ last t) -- we also handle inverted tableau+ +normalize :: ReverseHoleTableau -> TriangularArray Int +normalize = snd . normalize'++-- returns ( content , tableau )+normalize' :: ReverseHoleTableau -> ( Int , TriangularArray Int ) +normalize' holes = ( c , array (Tri (1,1), Tri (k,k)) xys ) where+ k = length holes+ c = length sorted+ xys = concat $ zipWith hs [1..] sorted+ hs a xs = map (h a) xs+ h a (ij,_) = (Tri ij , a) + sorted = groupSortBy snd (concat withPos)+ withPos = zipWith f [1..] (reverseTableau holes) + f i xs = zipWith (g i) [1..] xs + g i j hole = ((i,j),hole) ++--------------------------------------------------------------------------------++startHole :: [Hole] -> [Int] -> Hole +startHole (t:ts) (p:ps) = max t (toHole p)+startHole (t:ts) [] = t+startHole [] (p:ps) = toHole p+startHole [] [] = error "startHole"++-- c is the "content" of the small tableau+enumHoles :: Int -> Hole -> [Hole]+enumHoles c start@(Hole k l) + = nextHole start + : [ Hole i 0 | i <- [k+1..c] ] ++ [ Hole i 1 | i <- [k+1..c] ]++helper :: Int -> [Int] -> [Hole] -> [[Hole]]+helper c [] this = [[]] +helper c prev@(p:ps) this = + [ t:rest | t <- enumHoles c (startHole this prev), rest <- helper c ps (t:this) ]++newLines' :: Int -> [Int] -> [[Hole]]+newLines' c lastReversed = helper c last [] + where+ top = head lastReversed+ last = reverse (top : lastReversed)++newLines :: [Int] -> [[Hole]]+newLines lastReversed = newLines' (head lastReversed) lastReversed++-- | Generates all tableaux of size @k@. Effective for @k<=6@.+gtSimplexTableaux :: Int -> [TriangularArray Int]+gtSimplexTableaux 0 = [ triangularArrayUnsafe [] ]+gtSimplexTableaux 1 = [ triangularArrayUnsafe [[1]] ]+gtSimplexTableaux k = map normalize $ concatMap f smalls where+ smalls :: [ [[Int]] ]+ smalls = map (reverseTableau . fromTriangularArray) $ gtSimplexTableaux (k-1)+ f :: [[Int]] -> [ [[Hole]] ]+ f small = map (:smallhole) $ map reverse $ newLines (head small) where+ smallhole = map (map toHole) small++_gtSimplexTableaux :: Int -> [Tableau Int]+_gtSimplexTableaux k = map fromTriangularArray $ gtSimplexTableaux k++--------------------------------------------------------------------------------++-- | Note: This is slow (it actually generates all the tableaux)+countGTSimplexTableaux :: Int -> [Int]+countGTSimplexTableaux = elems . sizes'++sizes' :: Int -> UArray Int Int+sizes' k = + runSTUArray $ do+ let (a,b) = ( 2*k-1 , binom2 (k+1) )+ ar <- newArray (a,b) 0 :: ST s (STUArray s Int Int) + mapM_ (worker ar) $ gtSimplexTableaux k + return ar+ where+ worker :: STUArray s Int Int -> TriangularArray Int -> ST s ()+ worker ar t = do+ let c = gtSimplexContent t + n <- readArray ar c + writeArray ar c (n+1)+ +--------------------------------------------------------------------------------++-- | We can flip the numbers in the tableau so that the interval @[1..c]@ becomes+-- @[c..1]@. This way we a get a maybe more familiar form, when each row and each+-- column is strictly /decreasing/ (to the right and to the bottom).+invertGTSimplexTableau :: TriangularArray Int -> TriangularArray Int +invertGTSimplexTableau t = amap f t where+ c = gtSimplexContent t + f x = c+1-x ++_invertGTSimplexTableau :: [[Int]] -> [[Int]]+_invertGTSimplexTableau t = (map . map) f t where+ c = _gtSimplexContent t + f x = c+1-x++--------------------------------------------------------------------------------
− Math/Combinat/Tableaux/Kostka.hs
@@ -1,222 +0,0 @@---- TODO: better name?---- | This module contains a function to generate (equivalence classes of) --- triangular tableaux of size /k/, strictly increasing to the right and --- to the bottom. For example--- --- > 1 --- > 2 4 --- > 3 5 8 --- > 6 7 9 10 ------ is such a tableau of size 4.--- The numbers filling a tableau always consist of an interval @[1..c]@;--- @c@ is called the /content/ of the tableaux. There is a unique tableau--- of minimal content @2k-1@:------ > 1 --- > 2 3 --- > 3 4 5 --- > 4 5 6 7 --- --- Let us call the tableaux with maximal content (that is, @m = binomial (k+1) 2@)--- /standard/. The number of standard tableaux are------ > 1, 1, 2, 12, 286, 33592, 23178480, ...------ OEIS:A003121, \"Strict sense ballot numbers\", --- <https://oeis.org/A003121>.------ See --- R. M. Thrall, A combinatorial problem, Michigan Math. J. 1, (1952), 81-88.--- --- The number of tableaux with content @c=m-d@ are--- --- > d= | 0 1 2 3 ...--- > -----+------------------------------------------------- > k=2 | 1--- > k=3 | 2 1--- > k=4 | 12 18 8 1--- > k=5 | 286 858 1001 572 165 22 1--- > k=6 | 33592 167960 361114 436696 326196 155584 47320 8892 962 52 1 ------ We call these \"Kostka tableaux\" (in the lack of a better name), since--- they are in bijection with the simplicial cones in a canonical simplicial --- decompositions of the Gelfand-Tsetlin cones (the content corresponds--- to the dimension), which encode the combinatorics of Kostka numbers.-----module Math.Combinat.Tableaux.Kostka - ( - Tableau- , Tri(..)- , TriangularArray- , fromTriangularArray- , triangularArrayUnsafe- , kostkaTableaux- , _kostkaTableaux- , countKostkaTableaux- , kostkaContent- , _kostkaContent- ) - where------------------------------------------------------------------------------------import Data.Ix-import Data.Ord-import Data.List--import Control.Monad-import Control.Monad.ST-import Data.Array.IArray-import Data.Array.Unboxed-import Data.Array.ST--import Math.Combinat.Tableaux (Tableau)-import Math.Combinat.Helper-------------------------------------------------------------------------------------- | Triangular arrays-type TriangularArray a = Array Tri a---- | Set of @(i,j)@ pairs with @i>=j>=1@.-newtype Tri = Tri { unTri :: (Int,Int) } deriving (Eq,Ord,Show)--binom2 :: Int -> Int-binom2 n = (n*(n-1)) `div` 2--index' :: Tri -> Int-index' (Tri (i,j)) = binom2 i + j - 1---- it should be (1+8*m), --- the 2 is a hack to be safe with the floating point stuff-deIndex' :: Int -> Tri -deIndex' m = Tri ( i+1 , m - binom2 (i+1) + 1 ) where- i = ( (floor.sqrt.(fromIntegral::Int->Double)) (2+8*m) - 1 ) `div` 2 --instance Ix Tri where- index (a,b) x = index' x - index' a - inRange (a,b) x = (u<=j && j<=v) where- u = index' a - v = index' b- j = index' x- range (a,b) = map deIndex' [ index' a .. index' b ] - rangeSize (a,b) = index' b - index' a + 1 --{-# SPECIALIZE triangularArrayUnsafe :: Tableau Int -> TriangularArray Int #-}-triangularArrayUnsafe :: Tableau a -> TriangularArray a-triangularArrayUnsafe tableau = listArray (Tri (1,1),Tri (k,k)) (concat tableau) - where k = length tableau--{-# SPECIALIZE fromTriangularArray :: TriangularArray Int -> Tableau Int #-}-fromTriangularArray :: TriangularArray a -> Tableau a-fromTriangularArray arr = (map.map) snd $ groupBy (equating f) $ assocs arr- where f = fst . unTri . fst- ------------------------------------------------------------------------------------- "fractional fillings"-data Hole = Hole Int Int deriving (Eq,Ord,Show)--type ReverseTableau = [[Int ]] -type ReverseHoleTableau = [[Hole]] --toHole :: Int -> Hole-toHole k = Hole k 0--nextHole :: Hole -> Hole-nextHole (Hole k l) = Hole k (l+1)--{-# SPECIALIZE reverseTableau :: [[Int]] -> [[Int]] #-}-reverseTableau :: [[a]] -> [[a]]-reverseTableau = reverse . map reverse------------------------------------------------------------------------------------kostkaContent :: TriangularArray Int -> Int-kostkaContent arr = arr ! (snd (bounds arr))--_kostkaContent :: Tableau Int -> Int-_kostkaContent = last . last- -normalize :: ReverseHoleTableau -> TriangularArray Int -normalize = snd . normalize'---- returns ( content , tableau )-normalize' :: ReverseHoleTableau -> ( Int , TriangularArray Int ) -normalize' holes = ( c , array (Tri (1,1), Tri (k,k)) xys ) where- k = length holes- c = length sorted- xys = concat $ zipWith hs [1..] sorted- hs a xs = map (h a) xs- h a (ij,_) = (Tri ij , a) - sorted = groupSortBy snd (concat withPos)- withPos = zipWith f [1..] (reverseTableau holes) - f i xs = zipWith (g i) [1..] xs - g i j hole = ((i,j),hole) ------------------------------------------------------------------------------------startHole :: [Hole] -> [Int] -> Hole -startHole (t:ts) (p:ps) = max t (toHole p)-startHole (t:ts) [] = t-startHole [] (p:ps) = toHole p-startHole [] [] = error "startHole"---- c is the "content" of the small tableau-enumHoles :: Int -> Hole -> [Hole]-enumHoles c start@(Hole k l) - = nextHole start - : [ Hole i 0 | i <- [k+1..c] ] ++ [ Hole i 1 | i <- [k+1..c] ]--helper :: Int -> [Int] -> [Hole] -> [[Hole]]-helper c [] this = [[]] -helper c prev@(p:ps) this = - [ t:rest | t <- enumHoles c (startHole this prev), rest <- helper c ps (t:this) ]--newLines' :: Int -> [Int] -> [[Hole]]-newLines' c lastReversed = helper c last [] - where- top = head lastReversed- last = reverse (top : lastReversed)--newLines :: [Int] -> [[Hole]]-newLines lastReversed = newLines' (head lastReversed) lastReversed---- | Generates all tableaux of size @k@. Effective for @k<=6@.-kostkaTableaux :: Int -> [TriangularArray Int]-kostkaTableaux 0 = [ triangularArrayUnsafe [] ]-kostkaTableaux 1 = [ triangularArrayUnsafe [[1]] ]-kostkaTableaux k = map normalize $ concatMap f smalls where- smalls :: [ [[Int]] ]- smalls = map (reverseTableau . fromTriangularArray) $ kostkaTableaux (k-1)- f :: [[Int]] -> [ [[Hole]] ]- f small = map (:smallhole) $ map reverse $ newLines (head small) where- smallhole = map (map toHole) small--_kostkaTableaux :: Int -> [Tableau Int]-_kostkaTableaux k = map fromTriangularArray $ kostkaTableaux k------------------------------------------------------------------------------------countKostkaTableaux :: Int -> [Int]-countKostkaTableaux = elems . sizes'--sizes' :: Int -> UArray Int Int-sizes' k = - runSTUArray $ do- let (a,b) = ( 2*k-1 , binom2 (k+1) )- ar <- newArray (a,b) 0 :: ST s (STUArray s Int Int) - mapM_ (worker ar) $ kostkaTableaux k - return ar- where- worker :: STUArray s Int Int -> TriangularArray Int -> ST s ()- worker ar t = do- let c = kostkaContent t - n <- readArray ar c - writeArray ar c (n+1)- ---------------------------------------------------------------------------------
+ Math/Combinat/Tableaux/LittlewoodRichardson.hs view
@@ -0,0 +1,145 @@++-- | The Littlewood-Richardson rule++module Math.Combinat.Tableaux.LittlewoodRichardson + ( lrRule , _lrRule+ ) + where++--------------------------------------------------------------------------------++import Data.List++import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Skew++import Data.Map.Strict (Map)+import qualified Data.Map.Strict as Map++--------------------------------------------------------------------------------++-- | @lrRule@ computes the expansion of a skew Schur function +-- @s[lambda/mu]@ via the Littlewood-Richardson rule.+--+-- Adapted from John Stembridge's Maple code: +-- <http://www.math.lsa.umich.edu/~jrs/software/SFexamples/LR_rule>+--+lrRule :: SkewPartition -> Map Partition Int+lrRule skew = _lrRule lam mu where+ (lam,mu) = fromSkewPartition skew++-- | @_lrRule lambda mu@ computes the expansion of the skew+-- Schur function @s[lambda/mu]@ via the Littlewood-Richardson rule.+---+{-# SPECIALIZE _lrRule :: Partition -> Partition -> Map Partition Int #-}+{-# SPECIALIZE _lrRule :: Partition -> Partition -> Map Partition Integer #-}+_lrRule :: Num coeff => Partition -> Partition -> Map Partition coeff+_lrRule plam@(Partition lam) pmu@(Partition mu0) = + if not (plam `dominates` pmu) + then Map.empty+ else foldl' f Map.empty [ nu | (nu,_) <- fillings n diagram ]+ where+ f old nu = Map.insertWith (+) (Partition nu) 1 old+ diagram = [ (i,j) | (i,a,b) <- reverse (zip3 [1..] lam mu) , j <- [b+1..a] ] + mu = mu0 ++ repeat 0+ n = sum $ zipWith (-) lam mu -- n == length diagram++{-+LR_rule:=proc(lambda) local l,mu,alpha,beta,i,j,dgrm;+ if not `LR_rule/fit`(lambda,args[2]) then RETURN(0) fi;+ l:=nops(lambda); mu:=[op(args[2]),0$l];+ dgrm:=[seq(seq([i,-j],j=-lambda[i]..-1-mu[i]),i=1..l)];+ if nargs>2 then alpha:=args[3];+ if nargs>3 then beta:=args[4] else beta:=[] fi;+ if not `LR_rule/fit`(alpha,beta) then RETURN(0) fi;+ l:=convert([op(lambda),op(beta)],`+`);+ if l<>convert([op(alpha),op(mu)],`+`) then RETURN(0) fi;+ nops(LR_fillings(dgrm,[alpha,beta]))+ else+ convert([seq(s[op(i[1])],i=LR_fillings(dgrm))],`+`)+ fi+end;+-}++--------------------------------------------------------------------------------++-- | A filling is a pair consisting a shape @nu@ and a lattice permutation @lp@+type Filling = ( [Int] , [Int] )++-- | A diagram is a set of boxes in a skew shape (in the right order)+type Diagram = [ (Int,Int) ]++-- | Note: we use reverse ordering in Diagrams compared the Stembridge's code.+-- Also, for performance reasons, we need the length of the diagram+fillings :: Int -> Diagram -> [Filling]+fillings _ [] = [ ([],[]) ]+fillings n diagram@((x,y):rest) = concatMap (nextLetter lower upper) (fillings (n-1) rest) where+ upper = case findIndex (==(x ,y+1)) diagram of { Just j -> n-j ; Nothing -> 0 }+ lower = case findIndex (==(x-1,y )) diagram of { Just j -> n-j ; Nothing -> 0 }++{-+LR_fillings:=proc(dgrm) local n,x,upper,lower;+ if dgrm=[] then+ if nargs=1 then x:=[] else x:=args[2][2] fi;+ RETURN([[x,[]]])+ fi;+ n:=nops(dgrm); x:=dgrm[n];+ if not member([x[1],x[2]+1],dgrm,'upper') then upper:=0 fi;+ if not member([x[1]-1,x[2]],dgrm,'lower') then lower:=0 fi;+ if nargs=1 then+ map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)]),lower,upper)+ else+ map(`LR/nextletter`,LR_fillings([op(1..n-1,dgrm)],args[2]),+ lower,upper,args[2][1])+ fi;+end:+-}++--------------------------------------------------------------------------------++nextLetter :: Int -> Int -> Filling -> [Filling]+nextLetter lower upper (nu,lpart) = stuff where+ stuff = [ ( incr i shape , lpart++[i] ) | i<-nlist ] + shape = nu ++ [0] + lb = if lower>0+ then lpart !! (lower-1)+ else 0+ ub = if upper>0 + then min (length shape) (lpart !! (upper-1)) + else length shape+ nlist = filter (>0) $ map f [lb+1..ub] + f j = if j==1 || shape!!(j-2) > shape!!(j-1) then j else 0++ -- increments the i-th element (starting from 1)+ incr :: Int -> [Int] -> [Int]+ incr i (x:xs) = case i of+ 0 -> finish (x:xs)+ 1 -> (x+1) : finish xs+ _ -> x : incr (i-1) xs+ incr _ [] = []++ -- removes tailing zeros+ finish :: [Int] -> [Int]+ finish (x:xs) = if x>0 then x : finish xs else [] + finish [] = [] ++{-+`LR/nextletter`:=proc(T) local shape,lp,lb,ub,i,nl;+ shape:=[op(T[1]),0]; lp:=T[2]; ub:=nops(shape);+ if nargs>3 then ub:=min(ub,nops(args[4])) fi;+ if args[2]=0 then lb:=0 else lb:=lp[args[2]] fi;+ if args[3]>0 then ub:=min(lp[args[3]],ub) fi;+ if nargs<4 then+ nl:=map(proc(x,y) if x=1 or y[x-1]>y[x] then x fi end,[$lb+1..ub],shape)+ else+ nl:=map(proc(x,y,z) if y[x]<z[x] and (x=1 or y[x-1]>y[x]) then x fi end,+ [$lb+1..ub],shape,args[4])+ fi;+ nl:=[seq([subsop(i=shape[i]+1,shape),[op(lp),i]],i=nl)];+ op(subs(0=NULL,nl))+end:+-}++--------------------------------------------------------------------------------++
+ Math/Combinat/Tableaux/Skew.hs view
@@ -0,0 +1,78 @@+ +-- | Skew tableaux are skew partitions filled with numbers. + +{-# LANGUAGE BangPatterns #-} + +module Math.Combinat.Tableaux.Skew where + +-------------------------------------------------------------------------------- + +import Data.List + +import Math.Combinat.Partitions.Integer +import Math.Combinat.Partitions.Skew +import Math.Combinat.Tableaux +import Math.Combinat.ASCII + +-------------------------------------------------------------------------------- + +-- | A skew tableau is represented by a list of offsets and entries +newtype SkewTableau a = SkewTableau [(Int,[a])] deriving (Eq,Ord,Show) + +instance Functor SkewTableau where + fmap f (SkewTableau t) = SkewTableau [ (a, map f xs) | (a,xs) <- t ] + +skewShape :: SkewTableau a -> SkewPartition +skewShape (SkewTableau list) = SkewPartition [ (o,length xs) | (o,xs) <- list ] + +-- | Semi-standard skew tableaux filled with numbers @[1..n]@ +semiStandardSkewTableaux :: Int -> SkewPartition -> [SkewTableau Int] +semiStandardSkewTableaux n (SkewPartition abs) = map SkewTableau stuff where + + stuff = worker as bs ds (repeat 1) + (as,bs) = unzip abs + ds = diffSequence as + + -- | @worker inner outerMinusInner innerdiffs lowerbound + worker :: [Int] -> [Int] -> [Int] -> [Int] -> [[(Int,[Int])]] + worker (a:as) (b:bs) (d:ds) lb = [ (a,this):rest + | this <- row b 1 lb + , let lb' = (replicate d 1 ++ map (+1) this) + , rest <- worker as bs ds lb' ] + worker [] _ _ _ = [ [] ] + + -- @row length minimum lowerbound@ + row 0 _ _ = [[]] + row _ _ [] = [] + row !k !m (!a:as) = [ x:xs | x <- [(max a m)..n] , xs <- row (k-1) x as ] + +{- +-- | 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 [] = [] +-} + +-------------------------------------------------------------------------------- + +asciiSkewTableau :: Show a => SkewTableau a -> ASCII +asciiSkewTableau = asciiSkewTableau' "." EnglishNotation + +asciiSkewTableau' + :: Show a + => String -- ^ string representing the elements of the inner (unfilled) partition + -> PartitionConvention -- Orientation + -> SkewTableau a + -> ASCII +asciiSkewTableau' innerstr orient (SkewTableau axs) = tabulate (HRight,VTop) (HSepSpaces 1, VSepEmpty) stuff where + stuff = case orient of + EnglishNotation -> es + EnglishNotationCCW -> reverse (transpose es) + FrenchNotation -> reverse es + inner = asciiFromString innerstr + es = [ replicate a inner ++ map asciiShow xs | (a,xs) <- axs ] + +--------------------------------------------------------------------------------
combinat.cabal view
@@ -1,14 +1,14 @@ Name: combinat-Version: 0.2.7.0-Synopsis: Generation of various combinatorial objects.-Description: A collection of functions to generate (and if there is - a formula, count) combinatorial objects like partitions, - compositions, permutations, Young tableaux, various trees, - etc.+Version: 0.2.7.1+Synopsis: Generate and manipulate various combinatorial objects.+Description: A collection of functions to generate, count and manipulate+ all kinds of combinatorial objects like partitions, + compositions, permutations, Young tableaux, binary trees, + and so on. License: BSD3 License-file: LICENSE Author: Balazs Komuves-Copyright: (c) 2008-2014 Balazs Komuves+Copyright: (c) 2008-2015 Balazs Komuves Maintainer: bkomuves (plus) hackage (at) gmail (dot) com Homepage: http://code.haskell.org/~bkomuves/ Stability: Experimental@@ -45,11 +45,13 @@ Math.Combinat.Numbers Math.Combinat.Numbers.Series Math.Combinat.Numbers.Primes+ Math.Combinat.Sign Math.Combinat.Sets Math.Combinat.Tuples Math.Combinat.Compositions Math.Combinat.Partitions Math.Combinat.Partitions.Integer+ Math.Combinat.Partitions.Skew Math.Combinat.Partitions.Set Math.Combinat.Partitions.NonCrossing Math.Combinat.Partitions.Plane@@ -57,7 +59,10 @@ Math.Combinat.Partitions.Vector Math.Combinat.Permutations Math.Combinat.Tableaux- Math.Combinat.Tableaux.Kostka+ Math.Combinat.Tableaux.Skew+ Math.Combinat.Tableaux.GelfandTsetlin+ Math.Combinat.Tableaux.GelfandTsetlin.Cone+ Math.Combinat.Tableaux.LittlewoodRichardson Math.Combinat.Trees Math.Combinat.Trees.Binary Math.Combinat.Trees.Nary