diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -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
diff --git a/Math/Combinat.hs b/Math/Combinat.hs
--- a/Math/Combinat.hs
+++ b/Math/Combinat.hs
@@ -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
diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs
--- a/Math/Combinat/Helper.hs
+++ b/Math/Combinat/Helper.hs
@@ -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 
+-}
+
+--------------------------------------------------------------------------------
   
diff --git a/Math/Combinat/Numbers.hs b/Math/Combinat/Numbers.hs
--- a/Math/Combinat/Numbers.hs
+++ b/Math/Combinat/Numbers.hs
@@ -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)
 
 --------------------------------------------------------------------------------
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
@@ -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
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
@@ -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 ])
+
+-}
 
 --------------------------------------------------------------------------------
diff --git a/Math/Combinat/Partitions/Skew.hs b/Math/Combinat/Partitions/Skew.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Partitions/Skew.hs
@@ -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 ]
+    
+  
+--------------------------------------------------------------------------------
diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs
--- a/Math/Combinat/Permutations.hs
+++ b/Math/Combinat/Permutations.hs
@@ -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 = 
diff --git a/Math/Combinat/Sets.hs b/Math/Combinat/Sets.hs
--- a/Math/Combinat/Sets.hs
+++ b/Math/Combinat/Sets.hs
@@ -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
diff --git a/Math/Combinat/Sign.hs b/Math/Combinat/Sign.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Sign.hs
@@ -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
+
+--------------------------------------------------------------------------------
diff --git a/Math/Combinat/Tableaux.hs b/Math/Combinat/Tableaux.hs
--- a/Math/Combinat/Tableaux.hs
+++ b/Math/Combinat/Tableaux.hs
@@ -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 
diff --git a/Math/Combinat/Tableaux/GelfandTsetlin.hs b/Math/Combinat/Tableaux/GelfandTsetlin.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Tableaux/GelfandTsetlin.hs
@@ -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  
+
+--------------------------------------------------------------------------------
diff --git a/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs b/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Tableaux/GelfandTsetlin/Cone.hs
@@ -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
+
+--------------------------------------------------------------------------------
diff --git a/Math/Combinat/Tableaux/Kostka.hs b/Math/Combinat/Tableaux/Kostka.hs
deleted file mode 100644
--- a/Math/Combinat/Tableaux/Kostka.hs
+++ /dev/null
@@ -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)
-     
---------------------------------------------------------------------------------
diff --git a/Math/Combinat/Tableaux/LittlewoodRichardson.hs b/Math/Combinat/Tableaux/LittlewoodRichardson.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Tableaux/LittlewoodRichardson.hs
@@ -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:
+-}
+
+--------------------------------------------------------------------------------
+
+
diff --git a/Math/Combinat/Tableaux/Skew.hs b/Math/Combinat/Tableaux/Skew.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Tableaux/Skew.hs
@@ -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 ]
+
+--------------------------------------------------------------------------------
diff --git a/combinat.cabal b/combinat.cabal
--- a/combinat.cabal
+++ b/combinat.cabal
@@ -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