diff --git a/Math/Combinat.hs b/Math/Combinat.hs
--- a/Math/Combinat.hs
+++ b/Math/Combinat.hs
@@ -26,17 +26,17 @@
 --  * \"count\" prefix: counting functions.
 
 module Math.Combinat 
-  ( module Math.Combinat.Sets
+  ( module Math.Combinat.Numbers
+  , module Math.Combinat.Sets
   , module Math.Combinat.Tuples
   , module Math.Combinat.Combinations
   , module Math.Combinat.Partitions
   , module Math.Combinat.Permutations
   , module Math.Combinat.Tableaux
   , module Math.Combinat.Trees
-  , binomial
-  , factorial
   ) where
 
+import Math.Combinat.Numbers
 import Math.Combinat.Sets
 import Math.Combinat.Tuples
 import Math.Combinat.Combinations
@@ -44,5 +44,3 @@
 import Math.Combinat.Permutations
 import Math.Combinat.Tableaux
 import Math.Combinat.Trees
-
-import Math.Combinat.Helper ( binomial , factorial )
diff --git a/Math/Combinat/Combinations.hs b/Math/Combinat/Combinations.hs
--- a/Math/Combinat/Combinations.hs
+++ b/Math/Combinat/Combinations.hs
@@ -3,7 +3,7 @@
 
 module Math.Combinat.Combinations where
 
-import Math.Combinat.Helper
+import Math.Combinat.Numbers (factorial,binomial)
 
 -------------------------------------------------------
 
diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs
--- a/Math/Combinat/Helper.hs
+++ b/Math/Combinat/Helper.hs
@@ -2,8 +2,15 @@
 module Math.Combinat.Helper where
 
 import Control.Monad
+
+import Data.List
+import Data.Ord
+import qualified Data.Set as Set
+
 import Debug.Trace
 
+--------------------------------------------------------------------------------
+
 debug :: Show a => a -> b -> b
 debug x y = trace ("-- " ++ show x ++ "\n") y
 
@@ -12,18 +19,10 @@
 swap :: (a,b) -> (b,a)
 swap (x,y) = (y,x)
 
--- helps testing the random rutines 
-count :: Eq a => a -> [a] -> Int
-count x xs = length $ filter (==x) xs
-
-fromJust :: Maybe a -> a
-fromJust (Just x) = x
-fromJust Nothing = error "fromJust: Nothing"
+--------------------------------------------------------------------------------
 
--- iterated function application
-nest :: Int -> (a -> a) -> a -> a
-nest 0 _ x = x
-nest n f x = nest (n-1) f (f x)
+equating :: Eq b => (a -> b) -> a -> a -> Bool
+equating f x y = (f x == f y)
 
 reverseOrdering :: Ordering -> Ordering
 reverseOrdering LT = GT
@@ -33,20 +32,30 @@
 reverseCompare :: Ord a => a -> a -> Ordering
 reverseCompare x y = reverseOrdering $ compare x y
 
-factorial :: Int -> Integer
-factorial 0 = 1
-factorial n = product [1..fromIntegral n]
+groupSortBy :: (Eq b, Ord b) => (a -> b) -> [a] -> [[a]]
+groupSortBy f = groupBy (equating f) . sortBy (comparing f) 
 
-binomial :: Int -> Int -> Integer
-binomial n k 
-  | k > n = 0
-  | k < 0 = 0
-  | k > (n `div` 2) = binomial n (n-k)
-  | otherwise = (product [n'-k'+1 .. n']) `div` (product [1..k'])
-  where 
-    k' = fromIntegral k
-    n' = fromIntegral n
+nubOrd :: Ord a => [a] -> [a]
+nubOrd = worker Set.empty where
+  worker _ [] = []
+  worker s (x:xs) 
+    | Set.member x s = worker s xs
+    | otherwise      = x : worker (Set.insert x s) xs
     
+--------------------------------------------------------------------------------
+
+-- helps testing the random rutines 
+count :: Eq a => a -> [a] -> Int
+count x xs = length $ filter (==x) xs
+
+--------------------------------------------------------------------------------
+
+fromJust :: Maybe a -> a
+fromJust (Just x) = x
+fromJust Nothing = error "fromJust: Nothing"
+
+--------------------------------------------------------------------------------
+
 intToBool :: Int -> Bool
 intToBool 0 = False
 intToBool 1 = True
@@ -56,6 +65,13 @@
 boolToInt False = 0
 boolToInt True  = 1
 
+--------------------------------------------------------------------------------
+    
+-- iterated function application
+nest :: Int -> (a -> a) -> a -> a
+nest 0 _ x = x
+nest n f x = nest (n-1) f (f x)
+
 unfold1 :: (a -> Maybe a) -> a -> [a]
 unfold1 f x = case f x of 
   Nothing -> [x] 
@@ -87,4 +103,6 @@
   (s2,ys) <- mapAccumM f s1 xs
   return (s2, y:ys)
 
+--------------------------------------------------------------------------------
+    
   
diff --git a/Math/Combinat/Numbers.hs b/Math/Combinat/Numbers.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Numbers.hs
@@ -0,0 +1,126 @@
+
+-- | A few important number sequences. 
+--  
+-- See the "On-Line Encyclopedia of Integer Sequences",
+-- <http://www.research.att.com/~njas/sequences/> .
+
+module Math.Combinat.Numbers where
+
+import Data.Array
+
+--------------------------------------------------------------------------------
+
+-- | @(-1)^k@
+paritySign :: Integral a => a -> Integer
+paritySign k = if odd k then (-1) else 1
+
+--------------------------------------------------------------------------------
+
+-- | A000142.
+factorial :: Integral a => a -> Integer
+factorial n
+  | n <  0    = error "factorial: input should be nonnegative"
+  | n == 0    = 1
+  | otherwise = product [1..fromIntegral n]
+
+-- | A006882.
+doubleFactorial :: Integral a => a -> Integer
+doubleFactorial n
+  | n <  0    = error "doubleFactorial: input should be nonnegative"
+  | n == 0    = 1
+  | odd n     = product [1,3..fromIntegral n]
+  | otherwise = product [2,4..fromIntegral n]
+
+-- | A007318.
+binomial :: Integral a => a -> a -> Integer
+binomial n k 
+  | k > n = 0
+  | k < 0 = 0
+  | k > (n `div` 2) = binomial n (n-k)
+  | otherwise = (product [n'-k'+1 .. n']) `div` (product [1..k'])
+  where 
+    k' = fromIntegral k
+    n' = fromIntegral n
+
+--------------------------------------------------------------------------------
+-- * Catalan numbers
+
+-- | Catalan numbers. OEIS:A000108.
+catalan :: Integral a => a -> Integer
+catalan n 
+  | n < 0     = 0
+  | otherwise = binomial (n+n) n `div` fromIntegral (n+1)
+
+-- | Catalan's triangle. OEIS:A009766.
+-- Note:
+--
+-- > catalanTriangle n n == catalan n
+-- > catalanTriangle n k == countStandardYoungTableaux (toPartition [n,k])
+--
+catalanTriangle :: Integral a => a -> a -> Integer
+catalanTriangle n k
+  | k > n     = 0
+  | k < 0     = 0
+  | otherwise = binomial (n+k) n * fromIntegral (n-k+1) `div` fromIntegral (n+1)
+
+--------------------------------------------------------------------------------
+-- * Stirling numbers
+
+-- | Rows of (signed) Stirling numbers of the first kind. OEIS:A008275.
+-- Coefficients of the polinomial @(x-1)*(x-2)*...*(x-n+1)@.
+-- This function uses the recursion formula.
+signedStirling1stArray :: Integral a => a -> Array Int Integer
+signedStirling1stArray n
+  | n <  1    = error "stirling1stArray: n should be at least 1"
+  | n == 1    = listArray (1,1 ) [1]
+  | otherwise = listArray (1,n') [ lkp (k-1) - fromIntegral (n-1) * lkp k | k<-[1..n'] ] 
+  where
+    prev = signedStirling1stArray (n-1)
+    n' = fromIntegral n :: Int
+    lkp j | j <  1    = 0
+          | j >= n'   = 0
+          | otherwise = prev ! j 
+        
+-- | (Signed) Stirling numbers of the first kind. OEIS:A008275.
+-- This function uses "signedStirling1stArray", so it shouldn't be used
+-- to compute /many/ Stirling numbers.
+signedStirling1st :: Integral a => a -> a -> Integer
+signedStirling1st n k 
+  | k < 1     = 0
+  | k > n     = 0
+  | otherwise = signedStirling1stArray n ! (fromIntegral k)
+
+-- | (Unsigned) Stirling numbers of the first kind. OEIS:A008275.
+-- This function uses "signedStirling1stArray", so it shouldn't be used
+-- to compute /many/ Stirling numbers.
+unsignedStirling1st :: Integral a => a -> a -> Integer
+unsignedStirling1st n k = abs (signedStirling1st n k)
+
+-- | Stirling numbers of the second kind. OEIS:A008277.
+-- This function uses an explicit formula.
+stirling2nd :: Integral a => a -> a -> Integer
+stirling2nd n k 
+  | 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] ]
+
+--------------------------------------------------------------------------------
+-- * Bernoulli numbers
+
+-- | Bernoulli numbers. @bernoulli 1 == -1%2@ and @bernoulli k == 0@ for
+-- k>2 and /odd/. This function uses the formula involving Stirling numbers
+-- of the second kind. Numerators: A027641, denominators: A027642.
+bernoulli :: Integral a => a -> Rational
+bernoulli n 
+  | n <  0    = error "bernoulli: n should be nonnegative"
+  | n == 0    = 1
+  | n == 1    = -1/2
+  | otherwise = sum [ f k | k<-[1..n] ] 
+  where
+    f k = toRational (paritySign (n+k) * factorial k * stirling2nd n k) 
+        / toRational (k+1)
+
+--------------------------------------------------------------------------------
+
+ 
diff --git a/Math/Combinat/Partitions.hs b/Math/Combinat/Partitions.hs
--- a/Math/Combinat/Partitions.hs
+++ b/Math/Combinat/Partitions.hs
@@ -15,6 +15,8 @@
   , weight
   , _dualPartition
   , dualPartition
+  , _elements
+  , elements
     -- * Generation
   , _partitions' 
   , partitions'  
@@ -31,6 +33,7 @@
 
 import Data.List
 import Math.Combinat.Helper
+import Math.Combinat.Numbers (factorial,binomial)
 
 -------------------------------------------------------
 
@@ -86,6 +89,19 @@
 _dualPartition :: [Int] -> [Int]
 _dualPartition [] = []
 _dualPartition xs@(k:_) = [ length $ filter (>=i) xs | i <- [1..k] ]
+
+-- Example:
+--
+-- > elements (toPartition [5,2,1]) ==
+-- > [ (1,1), (1,2), (1,3), (1,4), (1,5)
+-- > , (2,1), (2,2), (2,3), (2,4)
+-- > , (3,1)
+-- > ]
+elements :: Partition -> [(Int,Int)]
+elements (Partition part) = _elements part
+
+_elements :: [Int] -> [(Int,Int)]
+_elements shape = [ (i,j) | (i,l) <- zip [1..] shape, j<-[1..l] ] 
   
 -------------------------------------------------------
 
diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs
--- a/Math/Combinat/Permutations.hs
+++ b/Math/Combinat/Permutations.hs
@@ -51,11 +51,17 @@
 import Control.Monad
 import Control.Monad.ST
 
+#if BASE_VERSION < 4
+import Data.List 
+#else
 import Data.List hiding (permutations)
+#endif
+
 import Data.Array
 import Data.Array.ST
 
 import Math.Combinat.Helper
+import Math.Combinat.Numbers (factorial,binomial)
 
 import System.Random
 
diff --git a/Math/Combinat/Sets.hs b/Math/Combinat/Sets.hs
--- a/Math/Combinat/Sets.hs
+++ b/Math/Combinat/Sets.hs
@@ -2,22 +2,49 @@
 -- | Subsets. 
 
 module Math.Combinat.Sets 
-  ( kSublists
+  ( 
+    choose
+  , combine
+  , tuplesFromList
+  
+  , kSublists
   , sublists
   , countKSublists
   , countSublists
   ) 
   where
 
-import Math.Combinat.Helper
+import Math.Combinat.Numbers (factorial,binomial)
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 
+-- | All possible ways to choose @k@ elements from a list, without
+-- repetitions. \"Antisymmetric power\" for lists. Synonym for "kSublists".
+choose :: Int -> [a] -> [[a]]
+choose 0 _  = [[]]
+choose k [] = []
+choose k (x:xs) = map (x:) (choose (k-1) xs) ++ choose k xs  
+
+-- | All possible ways to choose @k@ elements from a list, /with repetitions/. 
+-- \"Symmetric power\" for lists. See also "Math.Combinat.Combinations".
+-- TODO: better name?
+combine :: Int -> [a] -> [[a]]
+combine 0 _  = [[]]
+combine k [] = []
+combine k xxs@(x:xs) = map (x:) (combine (k-1) xxs) ++ combine k xs  
+
+-- | \"Tensor power\" for lists.
+-- See also "Math.Combinat.Tuples".
+-- TODO: better name?
+tuplesFromList :: Int -> [a] -> [[a]]
+tuplesFromList 0 _  = [[]]
+tuplesFromList k xs = [ (y:ys) | y <- xs, ys <- tuplesFromList (k-1) xs ]
+ 
+--------------------------------------------------------------------------------
+
 -- | Sublists of a list having given number of elements.
 kSublists :: Int -> [a] -> [[a]]
-kSublists 0 _  = [[]]
-kSublists k [] = []
-kSublists k (x:xs) = map (x:) (kSublists (k-1) xs) ++ kSublists k xs  
+kSublists = choose
 
 -- | @# = \binom { n } { k }@.
 countKSublists :: Int -> Int -> Integer
@@ -32,4 +59,4 @@
 countSublists :: Int -> Integer
 countSublists n = 2 ^ n
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
diff --git a/Math/Combinat/Tableaux.hs b/Math/Combinat/Tableaux.hs
--- a/Math/Combinat/Tableaux.hs
+++ b/Math/Combinat/Tableaux.hs
@@ -25,9 +25,10 @@
 import Data.List
 
 import Math.Combinat.Helper
+import Math.Combinat.Numbers (factorial,binomial)
 import Math.Combinat.Partitions
 
--------------------------------------------------------
+--------------------------------------------------------------------------------
 -- * Basic stuff
 
 type Tableau a = [[a]]
@@ -41,13 +42,30 @@
 dualTableau :: Tableau a -> Tableau a
 dualTableau = transpose
 
-hooks :: Partition -> Tableau Int
+content :: Tableau a -> [a]
+content = concat
+
+-- | An element @(i,j)@ of the resulting tableau (which has shape of the
+-- given partition) means that the vertical part of the hook has length @i@,
+-- and the horizontal part @j@. The /hook length/ is thus @i+j-1@. 
+--
+-- Example:
+--
+-- > > mapM_ print $ hooks $ toPartition [5,4,1]
+-- > [(3,5),(2,4),(2,3),(2,2),(1,1)]
+-- > [(2,4),(1,3),(1,2),(1,1)]
+-- > [(1,1)]
+--
+hooks :: Partition -> Tableau (Int,Int)
 hooks part = zipWith f p [1..] where 
   p = fromPartition part
   q = _dualPartition p
-  f l i = zipWith (\x y -> x+y-i) q [l,l-1..1] 
+  f l i = zipWith (\x y -> (x-i+1,y)) q [l,l-1..1] 
 
--------------------------------------------------------
+hookLengths :: Partition -> Tableau Int
+hookLengths part = (map . map) (\(i,j) -> i+j-1) (hooks part) 
+
+--------------------------------------------------------------------------------
 -- * Row and column words
 
 rowWord :: Tableau a -> [a]
@@ -68,7 +86,7 @@
 columnWordToTableau :: Ord a => [a] -> Tableau a
 columnWordToTableau = transpose . rowWordToTableau
     
--------------------------------------------------------
+--------------------------------------------------------------------------------
 -- * Standard Young tableaux
 
 -- | Standard Young tableaux of a given shape.
@@ -109,10 +127,33 @@
   
 -- | hook-length formula
 countStandardYoungTableaux :: Partition -> Integer
-countStandardYoungTableaux part = {- debug (hooks part) $ -}
+countStandardYoungTableaux part = {- debug (hookLengths part) $ -}
   factorial n `div` h where
-    h = product $ map fromIntegral $ concat $ hooks part 
+    h = product $ map fromIntegral $ concat $ hookLengths part 
     n = weight part
         
--------------------------------------------------------
+--------------------------------------------------------------------------------
+-- * Semistandard Young tableaux
+   
+-- | Semistandard Young tableaux of given shape, \"naive\" algorithm    
+semiStandardYoungTableaux :: Int -> Partition -> [Tableau Int]
+semiStandardYoungTableaux n part = worker (repeat 0) shape where
+  shape = fromPartition part
+  worker _ [] = [[]] 
+  worker prevRow (s:ss) 
+    = [ (r:rs) | r <- row n s 1 prevRow, rs <- worker (map (+1) r) ss ]
+
+  -- weekly increasing lists of length @len@, pointwise at least @xs@, 
+  -- maximum value @n@, minimum value @prev@.
+  row :: Int -> Int -> Int -> [Int] -> [[Int]]
+  row _ 0   _    _      = [[]]
+  row n len prev (x:xs) = [ (a:as) | a <- [max x prev..n] , as <- row n (len-1) a xs ]
+
+-- | Stanley's hook formula (cf. Fulton page 55)
+countSemiStandardYoungTableaux :: Int -> Partition -> Integer
+countSemiStandardYoungTableaux n shape = k `div` h where
+  h = product $ map fromIntegral $ concat $ hookLengths shape 
+  k = product [ fromIntegral (n+j-i) | (i,j) <- elements shape ]
+  
+--------------------------------------------------------------------------------
     
diff --git a/Math/Combinat/Tableaux/Kostka.hs b/Math/Combinat/Tableaux/Kostka.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Tableaux/Kostka.hs
@@ -0,0 +1,222 @@
+
+-- 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\", 
+-- <http://www.research.att.com/~njas/sequences/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/Trees.hs b/Math/Combinat/Trees.hs
--- a/Math/Combinat/Trees.hs
+++ b/Math/Combinat/Trees.hs
@@ -50,6 +50,7 @@
 import System.Random
 
 import Math.Combinat.Helper
+import Math.Combinat.Numbers (factorial,binomial)
 
 -------------------------------------------------------
 -- * Types
diff --git a/combinat.cabal b/combinat.cabal
--- a/combinat.cabal
+++ b/combinat.cabal
@@ -1,5 +1,5 @@
 Name:                combinat
-Version:             0.2.1
+Version:             0.2.2
 Synopsis:            Generation of various combinatorial objects.
 Description:         A collection of functions to generate combinatorial
                      objects like partitions, combinations, permutations,
@@ -16,28 +16,40 @@
 Cabal-Version:       >= 1.2
 Build-Type:          Simple
 
+Flag withQuickCheck
+  Description: Compile with the QuickCheck tests. 
+  default: False
+
 Flag splitBase
   Description: Choose the new smaller, split-up base package.
 
-Flag withQuickCheck
-  Description: Compile with the QuickCheck tests. 
+Flag base4
+  Description: Base v4
   
 Library
   if flag(splitBase)
-    Build-Depends:       base >= 3, array, containers, random
+    if flag(base4)
+      Build-Depends:       base >= 4 && < 5, array, containers, random
+      cpp-options:         -DBASE_VERSION=4
+    else 
+      Build-Depends:       base >= 3 && < 4, array, containers, random
+      cpp-options:         -DBASE_VERSION=3
     if flag(withQuickCheck)
       Build-Depends:       QuickCheck
   else
-    Build-Depends:       base <  3
+    Build-Depends:       base < 3
+    cpp-options:         -DBASE_VERSION=2
 
 
   Exposed-Modules:     Math.Combinat, 
+                       Math.Combinat.Numbers,
                        Math.Combinat.Sets,
                        Math.Combinat.Tuples, 
                        Math.Combinat.Combinations,
                        Math.Combinat.Partitions,
                        Math.Combinat.Permutations,
                        Math.Combinat.Tableaux,
+                       Math.Combinat.Tableaux.Kostka,
                        Math.Combinat.Trees
   
   Other-Modules:       Math.Combinat.Helper