diff --git a/Math/Combinat/Groups/Braid.hs b/Math/Combinat/Groups/Braid.hs
--- a/Math/Combinat/Groups/Braid.hs
+++ b/Math/Combinat/Groups/Braid.hs
@@ -17,7 +17,8 @@
       CPP, BangPatterns, 
       ScopedTypeVariables, ExistentialQuantification,
       DataKinds, KindSignatures, Rank2Types,
-      TypeOperators, TypeFamilies #-}
+      TypeOperators, TypeFamilies,
+      StandaloneDeriving #-}
 
 module Math.Combinat.Groups.Braid where
 
@@ -44,24 +45,22 @@
 import Math.Combinat.ASCII
 import Math.Combinat.Sign
 import Math.Combinat.Helper
+import Math.Combinat.TypeLevel
+import Math.Combinat.Numbers.Series
 
 import Math.Combinat.Permutations ( Permutation(..) )
 import qualified Math.Combinat.Permutations as P
 
-#ifdef QUICKCHECK
-import Test.QuickCheck
-#endif
-
 --------------------------------------------------------------------------------
 -- * Artin generators
 
 -- | A standard Artin generator of a braid: @Sigma i@ represents twisting 
--- the neighbour strands @i@ and @(i+1)@, such that @#i@ goes /under/ @#(i+1)@
+-- the neighbour strands @i@ and @(i+1)@, such that strand @i@ goes /under/ strand @(i+1)@.
 --
 -- Note: The strands are numbered @1..n@.
 data BrGen
-  = Sigma    !Int
-  | SigmaInv !Int
+  = Sigma    !Int         -- ^ @i@ goes under @(i+1)@
+  | SigmaInv !Int         -- ^ @i@ goes above @(i+1)@
   deriving (Eq,Ord,Show)
  
 -- | The strand (more precisely, the first of the two strands) the generator twistes
@@ -118,6 +117,24 @@
 withSomeBraid :: SomeBraid -> (forall n. KnownNat n => Braid n -> a) -> a
 withSomeBraid sbraid f = case sbraid of SomeBraid braid -> f braid
 
+mkBraid :: (forall n. KnownNat n => Braid n -> a) -> Int -> [BrGen] -> a
+mkBraid f n w = y where
+  sb = someBraid n (Braid w)
+  y  = withSomeBraid sb f
+
+withBraid 
+  :: Int
+  -> (forall (n :: Nat). KnownNat n => Braid n)
+  -> (forall (n :: Nat). KnownNat n => Braid n -> a) 
+  -> a
+withBraid n polyBraid f = 
+  case snat of    
+    SomeNat pxy -> f (asProxyTypeOf1 polyBraid pxy)
+  where
+    snat = case someNatVal (fromIntegral n :: Integer) of
+      Just sn -> sn
+      Nothing -> error "withBraid: input is not a natural number"
+
 --------------------------------------------------------------------------------
 
 braidWord :: Braid n -> [BrGen]
@@ -201,7 +218,7 @@
 theGarsideBraid :: KnownNat n => Braid n
 theGarsideBraid = halfTwist 
 
--- | The inner automorphism defined by @X -> Delta^-1 X Delta@, 
+-- | The inner automorphism defined by @tau(X) = Delta^-1 X Delta@, 
 -- where @Delta@ is the positive half-twist.
 -- 
 -- This sends each generator @sigma_j@ to @sigma_(n-j)@.
@@ -212,6 +229,13 @@
   f (Sigma    i) = Sigma    (n-i)
   f (SigmaInv i) = SigmaInv (n-i)
 
+
+-- | The involution @tau@ on permutations (permutation braids)
+--
+tauPerm :: Permutation -> Permutation
+tauPerm (Permutation arr) = Permutation $ listArray (1,n) [ (n+1) - arr!(n-i) | i<-[0..n-1] ] where
+  (1,n) = bounds arr
+
 --------------------------------------------------------------------------------
 -- * Group operations
 
@@ -437,7 +461,40 @@
         (sgn,k) = brGenSignIdx g
         s = signValue sgn
 
+--------------------------------------------------------------------------------
+-- * Growth 
 
+-- | Bronfman's recursive formula for the reciprocial of the growth function 
+-- of /positive/ braids. It was already known (by Deligne) that these generating functions 
+-- are reciprocials of polynomials; Bronfman [1] gave a recursive formula for them.
+--
+-- > let count n l = length $ nub $ [ braidNormalForm w | w <- allPositiveBraidWords n l ]
+-- > let convertPoly (1:cs) = zip (map negate cs) [1..]
+-- > pseries' (convertPoly $ bronfmanH n) == expandBronfmanH n == [ count n l | l <- [0..] ] 
+--
+-- * [1] Aaron Bronfman: Growth functions of a class of monoids. Preprint, 2001
+--
+bronfmanH :: Int -> [Int]
+bronfmanH n = bronfmanHsList !! n
+
+-- | An infinite list containing the Bronfman polynomials:
+--
+-- > bronfmanH n = bronfmanHsList !! n
+--
+bronfmanHsList :: [[Int]]
+bronfmanHsList = list where
+  list = map go [0..]
+  go 0 = [1]
+  go n = sumSeries [ sgn i $ replicate (choose2 i) 0 ++ list !! (n-i) | i<-[1..n] ]
+  sgn i = if odd i then id else map negate
+  choose2 k = div (k*(k-1)) 2
+
+-- | Expands the reciprocial of @H(n)@ into an infinite power series,
+-- giving the growth function of the positive braids on @n@ strands.
+expandBronfmanH :: Int -> [Int]
+expandBronfmanH n = pseries' (convertPoly $ bronfmanH n) where
+  convertPoly (1:cs) = zip (map negate cs) [1..]
+   
 --------------------------------------------------------------------------------
 -- * ASCII diagram
 
@@ -534,29 +591,49 @@
 -}
 
 --------------------------------------------------------------------------------
+-- * List of all words
+
+-- | All positive braid words of the given length
+allPositiveBraidWords :: KnownNat n => Int -> [Braid n]
+allPositiveBraidWords l = braids where
+  n = numberOfStrands (head braids)
+  braids = map Braid $ _allPositiveBraidWords n l 
+
+-- | All braid words of the given length
+allBraidWords :: KnownNat n => Int -> [Braid n]
+allBraidWords l = braids where
+  n = numberOfStrands (head braids)
+  braids = map Braid $ _allBraidWords n l 
+
+-- | Untyped version of 'allPositiveBraidWords'
+_allPositiveBraidWords :: Int -> Int -> [[BrGen]]
+_allPositiveBraidWords n = go where
+  go 0 = [[]]
+  go k = [ Sigma i : rest | i<-[1..n-1] , rest <- go (k-1) ]
+
+-- | Untyped version of 'allBraidWords'
+_allBraidWords :: Int -> Int -> [[BrGen]]
+_allBraidWords n = go where
+  go 0 = [[]]
+  go k = [ gen : rest | gen <- gens , rest <- go (k-1) ]
+  gens = concat [ [ Sigma i , SigmaInv i ] | i<-[1..n-1] ]
+
+--------------------------------------------------------------------------------
 -- * Random braids  
 
 -- | Random braid word of the given length
 randomBraidWord :: (RandomGen g, KnownNat n) => Int -> g -> (Braid n, g)
-randomBraidWord len gen = (braid, gen') where
-  braid = Braid (map sig bjs)
-  n     = numberOfStrands braid
-  (gen',bjs) = mapAccumL worker gen [1..len]
-
-  worker !g _ = (g'',(b,j)) where
-    (j, g' ) = randomR (1,n-1) g
-    (b, g'') = random          g'
-
-  sig :: (Bool,Int) -> BrGen
-  sig (True ,j) = Sigma    j
-  sig (False,j) = SigmaInv j
+randomBraidWord len g = (braid, g') where
+  braid  = Braid w
+  n      = numberOfStrands braid
+  (w,g') = _randomBraidWord n len g
 
 -- | Random /positive/ braid word of the given length
 randomPositiveBraidWord :: (RandomGen g, KnownNat n) => Int -> g -> (Braid n, g)
-randomPositiveBraidWord len gen = (braid, gen') where
-  braid = Braid (map Sigma js)
-  n     = numberOfStrands braid
-  (gen',js) = mapAccumL (\(!g) _ -> swap (randomR (1,n-1) g)) gen [1..len]
+randomPositiveBraidWord len g = (braid, g') where
+  braid  = Braid w
+  n      = numberOfStrands braid
+  (w,g') = _randomPositiveBraidWord n len g
 
 --------------------------------------------------------------------------------
 
@@ -621,31 +698,47 @@
 
 --------------------------------------------------------------------------------
 
-#ifdef QUICKCHECK
-
--- | A permutation braid made convenient to use (type-level hackery)
-data PermBraid = forall n. KnownNat n => PermBraid Permutation (Braid n)
-
-mkPermBraid :: Permutation -> PermBraid
-mkPermBraid perm = 
-  case snat of    
-    SomeNat pxy -> PermBraid perm (asProxyTypeOf1 (permutationBraid perm) pxy)
-  where
-    n = P.permutationSize perm
-    Just snat = someNatVal (fromIntegral n :: Integer)
-
-prop_permBraid_perm :: PermBraid -> Bool
-prop_permBraid_perm (PermBraid perm braid) = (braidPermutation braid == perm)
-
-prop_permBraid_valid :: PermBraid -> Bool
-prop_permBraid_valid (PermBraid perm braid) = isPermutationBraid braid
+-- | This version of 'randomBraidWord' may be convenient to avoid the type level stuff
+withRandomBraidWord 
+  :: RandomGen g 
+  => (forall n. KnownNat n => Braid n -> a) 
+  -> Int                -- ^ number of strands
+  -> Int                -- ^ length of the random word
+  -> g -> (a, g)
+withRandomBraidWord f n len = runRand $ do
+  withSelectedM f (rand $ randomBraidWord len) n
 
-prop_braidPerm_comp :: KnownNat n => Braid n -> Braid n -> Bool
-prop_braidPerm_comp b1 b2 = (p == q) where
-  p = braidPermutation (compose b1 b2) 
-  q = braidPermutation b1 `P.multiply` braidPermutation b2
+-- | This version of 'randomPositiveBraidWord' may be convenient to avoid the type level stuff
+withRandomPositiveBraidWord 
+  :: RandomGen g 
+  => (forall n. KnownNat n => Braid n -> a) 
+  -> Int                -- ^ number of strands
+  -> Int                -- ^ length of the random word
+  -> g -> (a, g)
+withRandomPositiveBraidWord f n len = runRand $ do
+  withSelectedM f (rand $ randomPositiveBraidWord len) n
 
+-- | Untyped version of 'randomBraidWord'
+_randomBraidWord 
+  :: (RandomGen g) 
+  => Int                -- ^ number of strands
+  -> Int                -- ^ length of the random word
+  -> g -> ([BrGen], g)
+_randomBraidWord n len = runRand $ replicateM len $ do
+  k <- randChoose (1,n-1)
+  s <- randRoll
+  return $ case s of
+    Plus  -> Sigma k
+    Minus -> SigmaInv k
 
-#endif
+-- | Untyped version of 'randomPositiveBraidWord'
+_randomPositiveBraidWord 
+  :: (RandomGen g) 
+  => Int             -- ^ number of strands
+  -> Int             -- ^ length of the random word
+  -> g -> ([BrGen], g)
+_randomPositiveBraidWord n len = runRand $ replicateM len $ do
+  liftM Sigma $ randChoose (1,n-1)
 
 --------------------------------------------------------------------------------
+
diff --git a/Math/Combinat/Groups/Braid/NF.hs b/Math/Combinat/Groups/Braid/NF.hs
--- a/Math/Combinat/Groups/Braid/NF.hs
+++ b/Math/Combinat/Groups/Braid/NF.hs
@@ -18,12 +18,17 @@
       DataKinds, KindSignatures, Rank2Types #-}
 
 module Math.Combinat.Groups.Braid.NF  
-  ( BraidNF (..)
+  ( -- * Normal form
+    BraidNF (..)
   , nfReprWord
   , braidNormalForm
   , braidNormalForm'
-#ifdef QUICKCHECK
-#endif
+  , braidNormalFormNaive'
+    -- * Starting and finishing sets
+  , permWordStartingSet
+  , permWordFinishingSet    
+  , permutationStartingSet
+  , permutationFinishingSet    
   )
   where
 
@@ -53,10 +58,6 @@
 
 import Math.Combinat.Groups.Braid
 
-#ifdef QUICKCHECK
-import Test.QuickCheck
-#endif
-
 --------------------------------------------------------------------------------
 
 -- | A unique normal form for braids, called the /left-greedy normal form/.
@@ -95,11 +96,19 @@
 braidNormalForm' braid@(Braid gens) = BraidNF (dexp+pexp) perms where
   n = numberOfStrands braid
   invless = replaceInverses n gens
-  -- invless = replaceInversesNaive gens
   (dexp,posxword) = moveDeltasLeft n invless
   factors = leftGreedyFactors n $ expandPosXWord n posxword
   (pexp,perms) = normalizePermFactors n $ map (_braidPermutation n) factors
 
+-- | This one uses the naive inverse replacement method. Probably somewhat slower than 'braidNormalForm''.
+braidNormalFormNaive' :: KnownNat n => Braid n -> BraidNF n
+braidNormalFormNaive' braid@(Braid gens) = BraidNF (dexp+pexp) perms where
+  n = numberOfStrands braid
+  invless = replaceInversesNaive gens
+  (dexp,posxword) = moveDeltasLeft n invless
+  factors = leftGreedyFactors n $ expandPosXWord n posxword
+  (pexp,perms) = normalizePermFactors n $ map (_braidPermutation n) factors
+
 --------------------------------------------------------------------------------
 
 -- | Replaces groups of @sigma_i^-1@ generators by @(Delta^-1 * P)@, 
@@ -421,7 +430,7 @@
     | otherwise                = (e   , [oddTau (e+ep) p] )
 
   oddTau :: Int -> Permutation -> Permutation
-  oddTau !e p = if even e then p else permTau p
+  oddTau !e p = if even e then p else tauPerm p
 
 {-
   checkDelta :: Int -> Permutation -> [Permutation] -> (Int,[Permutation])
@@ -431,11 +440,6 @@
     | otherwise                  = let (e',rs) = worker e rest in (e', oddTau e p : rs)
 -}        
 
--- | The involution tau on permutation
-permTau :: Permutation -> Permutation
-permTau (Permutation arr) = Permutation $ listArray (1,n) [ (n+1) - arr!(n-i) | i<-[0..n-1] ] where
-  (1,n) = bounds arr
-
 -------------------------------------------------------------------------------- 
 
 -- | Given a /positive/ word, we apply left-greedy factorization of
@@ -528,17 +532,3 @@
 
 --------------------------------------------------------------------------------
 
-#ifdef QUICKCHECK
-
-prop_braidnf_reduce :: KnownNat n => Braid n -> Bool
-prop_braidnf_reduce braid = (braidNormalForm' braid == braidNormalForm braid)
-
-prop_braidnf_reprs :: KnownNat n => Braid n -> Bool
-prop_braidnf_reprs braid = (nf == nf') where
-  nf  = braidNormalForm braid 
-  nf' = braidNormalForm braid'
-  braid' = nfReprWord nf
-
-#endif
-
---------------------------------------------------------------------------------
diff --git a/Math/Combinat/Groups/Free.hs b/Math/Combinat/Groups/Free.hs
--- a/Math/Combinat/Groups/Free.hs
+++ b/Math/Combinat/Groups/Free.hs
@@ -493,7 +493,7 @@
       
 {-
 
--- some basic testing. TODO: QuickCheck tests
+-- some basic testing. TODO: real tests
 
 import Math.Combinat.Helper
 import Math.Combinat.Groups.Free
diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs
--- a/Math/Combinat/Helper.hs
+++ b/Math/Combinat/Helper.hs
@@ -1,12 +1,14 @@
 
 -- | Miscellaneous helper functions
 
-{-# LANGUAGE BangPatterns, PolyKinds #-}
+{-# LANGUAGE BangPatterns, PolyKinds, GeneralizedNewtypeDeriving #-}
 module Math.Combinat.Helper where
 
 --------------------------------------------------------------------------------
 
 import Control.Monad
+import Control.Applicative ( Applicative(..) )    -- required before AMP (before GHC 7.10)
+import Data.Functor.Identity
 
 import Data.List
 import Data.Ord
@@ -17,6 +19,9 @@
 
 import Debug.Trace
 
+import System.Random
+import Control.Monad.Trans.State
+
 --------------------------------------------------------------------------------
 -- * debugging
 
@@ -24,24 +29,6 @@
 debug x y = trace ("-- " ++ show x ++ "\n") y
 
 --------------------------------------------------------------------------------
--- * proxy
-
-proxyUndef :: Proxy a -> a
-proxyUndef _ = error "proxyUndef"
-
-proxyOf :: a -> Proxy a
-proxyOf _ = Proxy
-
-proxyOf1 :: f a -> Proxy a
-proxyOf1 _ = Proxy
-
-proxyOf2 :: g (f a) -> Proxy a
-proxyOf2 _ = Proxy
-
-asProxyTypeOf1 :: f a -> Proxy a -> f a 
-asProxyTypeOf1 y _ = y
-
---------------------------------------------------------------------------------
 -- * pairs
 
 swap :: (a,b) -> (b,a)
@@ -254,4 +241,40 @@
 -}
 
 --------------------------------------------------------------------------------
-  
+-- * random
+
+-- | A simple random monad to make life suck less
+type Rand g = RandT g Identity
+
+runRand :: Rand g a -> g -> (a,g)
+runRand action g = runIdentity (runRandT action g)
+
+flipRunRand :: Rand s a -> s -> (s,a)
+flipRunRand action g = runIdentity (flipRunRandT action g)
+
+
+-- | The Rand monad transformer
+newtype RandT g m a = RandT (StateT g m a) deriving (Functor,Applicative,Monad)
+
+runRandT :: RandT g m a -> g -> m (a,g)
+runRandT (RandT stuff) = runStateT stuff
+
+-- | This may be occasionally useful
+flipRunRandT :: Monad m => RandT s m a -> s -> m (s,a)
+flipRunRandT action ini = liftM swap $ runRandT action ini
+
+
+-- | Puts a standard-conforming random function into the monad
+rand :: (g -> (a,g)) -> Rand g a
+rand user = RandT (state user)
+
+randRoll :: (RandomGen g, Random a) => Rand g a
+randRoll = rand random
+
+randChoose :: (RandomGen g, Random a) => (a,a) -> Rand g a
+randChoose uv = rand (randomR uv)
+
+randProxy1 :: Rand g (f n) -> Proxy n -> Rand g (f n)
+randProxy1 action _ = action
+
+--------------------------------------------------------------------------------
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
@@ -20,11 +20,6 @@
 import Math.Combinat.Partitions.Integer
 import Math.Combinat.Helper
 
-#ifdef QUICKCHECK
-import System.Random
-import Test.QuickCheck
-#endif
-
 --------------------------------------------------------------------------------
 -- * Trivial series
 
@@ -55,6 +50,10 @@
 addSeries :: Num a => [a] -> [a] -> [a]
 addSeries xs ys = longZipWith 0 0 (+) xs ys
 
+sumSeries :: Num a => [[a]] -> [a]
+sumSeries [] = [0]
+sumSeries xs = foldl1' addSeries xs
+
 subSeries :: Num a => [a] -> [a] -> [a]
 subSeries xs ys = longZipWith 0 0 (-) xs ys
 
@@ -301,25 +300,6 @@
 
 -- Reciprocals of polynomials, without coefficients
 
-#ifdef QUICKCHECK
--- | Expansion of @1 / (1-x^k)@. Included for completeness only; 
--- it equals to @coinSeries [k]@, and for example
--- for @k=4@ it is simply
--- 
--- > [1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0...]
---
-pseries1 :: Int -> [Integer]
-pseries1 k1 = convolveWithPSeries1 k1 unitSeries 
-
--- | The expansion of @1 / (1-x^k_1-x^k_2)@
-pseries2 :: Int -> Int -> [Integer]
-pseries2 k1 k2 = convolveWithPSeries2 k1 k2 unitSeries 
-
--- | The expansion of @1 / (1-x^k_1-x^k_2-x^k_3)@
-pseries3 :: Int -> Int -> Int -> [Integer]
-pseries3 k1 k2 k3 = convolveWithPSeries3 k1 k2 k3 unitSeries
-#endif
-
 -- | The power series expansion of 
 --
 -- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)
@@ -327,33 +307,6 @@
 pseries :: [Int] -> [Integer]
 pseries ks = convolveWithPSeries ks unitSeries
 
-#ifdef QUICKCHECK
--- | Convolve with (the expansion of) @1 / (1-x^k1)@
-convolveWithPSeries1 :: Int -> [Integer] -> [Integer]
-convolveWithPSeries1 k1 series1 = xs where
-  series = series1 ++ repeat 0 
-  xs = zipWith (+) series ( replicate k1 0 ++ xs )
-
--- | Convolve with (the expansion of) @1 / (1-x^k1-x^k2)@
-convolveWithPSeries2 :: Int -> Int -> [Integer] -> [Integer]
-convolveWithPSeries2 k1 k2 series1 = xs where
-  series = series1 ++ repeat 0 
-  xs = zipWith3 (\x y z -> x + y + z)
-    series
-    ( replicate k1 0 ++ xs )
-    ( replicate k2 0 ++ xs )
-    
--- | Convolve with (the expansion of) @1 / (1-x^k_1-x^k_2-x^k_3)@
-convolveWithPSeries3 :: Int -> Int -> Int -> [Integer] -> [Integer]
-convolveWithPSeries3 k1 k2 k3 series1 = xs where
-  series = series1 ++ repeat 0 
-  xs = zipWith4 (\x y z w -> x + y + z + w)
-    series
-    ( replicate k1 0 ++ xs )
-    ( replicate k2 0 ++ xs )
-    ( replicate k3 0 ++ xs )
-#endif
-
 -- | Convolve with (the expansion of) 
 --
 -- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)
@@ -369,24 +322,6 @@
 --------------------------------------------------------------------------------
 --  Reciprocals of polynomials, with coefficients
 
-#ifdef QUICKCHECK
--- | @1 / (1 - a*x^k)@. 
--- For example, for @a=3@ and @k=2@ it is just
--- 
--- > [1,0,3,0,9,0,27,0,81,0,243,0,729,0,2187,0,6561,0,19683,0...]
---
-pseries1' :: Num a => (a,Int) -> [a]
-pseries1' ak1 = convolveWithPSeries1' ak1 unitSeries
-
--- | @1 / (1 - a_1*x^k_1 - a_2*x^k_2)@
-pseries2' :: Num a => (a,Int) -> (a,Int) -> [a]
-pseries2' ak1 ak2 = convolveWithPSeries2' ak1 ak2 unitSeries
-
--- | @1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3)@
-pseries3' :: Num a => (a,Int) -> (a,Int) -> (a,Int) -> [a]
-pseries3' ak1 ak2 ak3 = convolveWithPSeries3' ak1 ak2 ak3 unitSeries
-#endif
-
 -- | The expansion of 
 --
 -- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)
@@ -394,35 +329,6 @@
 pseries' :: Num a => [(a,Int)] -> [a]
 pseries' aks = convolveWithPSeries' aks unitSeries
 
-#ifdef QUICKCHECK
--- | Convolve with @1 / (1 - a*x^k)@. 
-convolveWithPSeries1' :: Num a => (a,Int) -> [a] -> [a]
-convolveWithPSeries1' (a1,k1) series1 = xs where
-  series = series1 ++ repeat 0 
-  xs = zipWith (+)
-    series
-    ( replicate k1 0 ++ map (*a1) xs )
-
--- | Convolve with @1 / (1 - a_1*x^k_1 - a_2*x^k_2)@
-convolveWithPSeries2' :: Num a => (a,Int) -> (a,Int) -> [a] -> [a]
-convolveWithPSeries2' (a1,k1) (a2,k2) series1 = xs where
-  series = series1 ++ repeat 0 
-  xs = zipWith3 (\x y z -> x + y + z)
-    series
-    ( replicate k1 0 ++ map (*a1) xs )
-    ( replicate k2 0 ++ map (*a2) xs )
-    
--- | Convolve with @1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3)@
-convolveWithPSeries3' :: Num a => (a,Int) -> (a,Int) -> (a,Int) -> [a] -> [a]
-convolveWithPSeries3' (a1,k1) (a2,k2) (a3,k3) series1 = xs where
-  series = series1 ++ repeat 0 
-  xs = zipWith4 (\x y z w -> x + y + z + w)
-    series
-    ( replicate k1 0 ++ map (*a1) xs )
-    ( replicate k2 0 ++ map (*a2) xs )
-    ( replicate k3 0 ++ map (*a3) xs )
-#endif
-
 -- | Convolve with (the expansion of) 
 --
 -- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)
@@ -467,186 +373,4 @@
      
 --------------------------------------------------------------------------------
 
-#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
-(=!=) xs1 ys1 = (take m xs == take m ys) where 
-  m = 1000
-  xs = xs1 ++ repeat 0
-  ys = ys1 ++ repeat 0
-infix 4 =!=
-
-newtype Nat = Nat { fromNat :: Int } deriving (Eq,Ord,Show,Num,Random)
-newtype Ser = Ser { fromSer :: [Integer] } deriving (Eq,Ord,Show)
-newtype Exp  = Exp  { fromExp  ::  Int  } deriving (Eq,Ord,Show,Num,Random)
-newtype Exps = Exps { fromExps :: [Int] } deriving (Eq,Ord,Show)
-newtype CoeffExp  = CoeffExp  { fromCoeffExp  ::  (Integer,Int)  } deriving (Eq,Ord,Show)
-newtype CoeffExps = CoeffExps { fromCoeffExps :: [(Integer,Int)] } deriving (Eq,Ord,Show)
-
-minSerSize = 0    :: Int
-maxSerSize = 1000 :: Int
-
-minSerValue = -10000 :: Integer
-maxSerValue =  10000 :: Integer
-
-rndList :: (RandomGen g, Random a) => Int -> (a, a) -> g -> ([a], g)
-rndList n minmax g = swap $ mapAccumL f g [1..n] where
-  f g _ = swap $ randomR minmax g 
-
-instance Arbitrary Nat where
-  arbitrary = choose (Nat 0 , Nat 750)
-
-instance Arbitrary Exp where
-  arbitrary = choose (Exp 1 , Exp 32)
-
-instance Arbitrary CoeffExp where
-  arbitrary = do
-    coeff <- choose (minSerValue, maxSerValue) :: Gen Integer
-    exp <- arbitrary :: Gen Exp
-    return $ CoeffExp (coeff,fromExp exp)
-   
-instance Random Ser where
-  random g = (Ser list, g2) where
-    (size,g1) = randomR (minSerSize,maxSerSize) g
-    (list,g2) = rndList size (minSerValue,maxSerValue) g1
-  randomR _ = random
-
-instance Random Exps where
-  random g = (Exps list, g2) where
-    (size,g1) = randomR (0,10) g
-    (list,g2) = rndList size (1,32) g1
-  randomR _ = random
-
-instance Random CoeffExps where
-  random g = (CoeffExps (zip list2 list1), g3) where
-    (size,g1) = randomR (0,10) g
-    (list1,g2) = rndList size (1,32) g1
-    (list2,g3) = rndList size (minSerValue,maxSerValue) g2
-  randomR _ = random
-  
-instance Arbitrary Ser where
-  arbitrary = choose undefined
-
-instance Arbitrary Exps where
-  arbitrary = choose undefined
-
-instance Arbitrary CoeffExps where
-  arbitrary = choose undefined
-  
--- TODO: quickcheck test properties
-
-checkAll :: IO ()
-checkAll = do
-  let f :: Testable a => a -> IO ()
-      f = quickCheck
-{- 
-  -- these are very slow, because random is slow
-  putStrLn "leftIdentity"  ; f prop_leftIdentity
-  putStrLn "rightIdentity" ; f prop_rightIdentity
-  putStrLn "commutativity" ; f prop_commutativity
-  putStrLn "associativity" ; f prop_associativity
--}
-  putStrLn "convPSeries1 vs generic" ; f prop_conv1_vs_gen
-  putStrLn "convPSeries2 vs generic" ; f prop_conv2_vs_gen
-  putStrLn "convPSeries3 vs generic" ; f prop_conv3_vs_gen
-  putStrLn "convPSeries1' vs generic" ; f prop_conv1_vs_gen'
-  putStrLn "convPSeries2' vs generic" ; f prop_conv2_vs_gen'
-  putStrLn "convPSeries3' vs generic" ; f prop_conv3_vs_gen'
-  putStrLn "convolve_pseries"  ; f prop_convolve_pseries 
-  putStrLn "convolve_pseries'" ; f prop_convolve_pseries' 
-  putStrLn "coinSeries vs pseries"  ; f prop_coin_vs_pseries
-  putStrLn "coinSeries vs pseries'" ; f prop_coin_vs_pseries'
-     
-prop_leftIdentity ser = ( xs =!= unitSeries `convolve` xs ) where 
-  xs = fromSer ser 
-
-prop_rightIdentity ser = ( unitSeries `convolve` xs =!= xs ) where 
-  xs = fromSer ser 
-
-prop_commutativity ser1 ser2 = ( xs `convolve` ys =!= ys `convolve` xs ) where 
-  xs = fromSer ser1
-  ys = fromSer ser2
-
-prop_associativity ser1 ser2 ser3 = ( one =!= two ) where
-  one = (xs `convolve` ys) `convolve` zs
-  two = xs `convolve` (ys `convolve` zs)
-  xs = fromSer ser1
-  ys = fromSer ser2
-  zs = fromSer ser3
-  
-prop_conv1_vs_gen exp1 ser = ( one =!= two ) where
-  one = convolveWithPSeries1 k1 xs 
-  two = convolveWithPSeries [k1] xs
-  k1 = fromExp exp1
-  xs = fromSer ser  
-
-prop_conv2_vs_gen exp1 exp2 ser = (one =!= two) where
-  one = convolveWithPSeries2 k1 k2 xs 
-  two = convolveWithPSeries [k2,k1] xs
-  k1 = fromExp exp1
-  k2 = fromExp exp2
-  xs = fromSer ser  
-
-prop_conv3_vs_gen exp1 exp2 exp3 ser = (one =!= two) where
-  one = convolveWithPSeries3 k1 k2 k3 xs 
-  two = convolveWithPSeries [k2,k3,k1] xs
-  k1 = fromExp exp1
-  k2 = fromExp exp2
-  k3 = fromExp exp3
-  xs = fromSer ser  
-
-prop_conv1_vs_gen' exp1 ser = ( one =!= two ) where
-  one = convolveWithPSeries1' ak1 xs 
-  two = convolveWithPSeries' [ak1] xs
-  ak1 = fromCoeffExp exp1
-  xs = fromSer ser  
-
-prop_conv2_vs_gen' exp1 exp2 ser = (one =!= two) where
-  one = convolveWithPSeries2' ak1 ak2 xs 
-  two = convolveWithPSeries' [ak2,ak1] xs
-  ak1 = fromCoeffExp exp1
-  ak2 = fromCoeffExp exp2
-  xs = fromSer ser  
-
-prop_conv3_vs_gen' exp1 exp2 exp3 ser = (one =!= two) where
-  one = convolveWithPSeries3' ak1 ak2 ak3 xs 
-  two = convolveWithPSeries' [ak2,ak3,ak1] xs
-  ak1 = fromCoeffExp exp1
-  ak2 = fromCoeffExp exp2
-  ak3 = fromCoeffExp exp3
-  xs = fromSer ser  
-
-prop_convolve_pseries exps1 ser = (one =!= two) where
-  one = convolveWithPSeries ks1 xs 
-  two = xs `convolve` pseries ks1 
-  ks1 = fromExps exps1
-  xs = fromSer ser  
-
-prop_convolve_pseries' cexps1 ser = (one =!= two) where
-  one = convolveWithPSeries' aks1 xs 
-  two = xs `convolve` pseries' aks1 
-  aks1 = fromCoeffExps cexps1
-  xs = fromSer ser  
-
-prop_coin_vs_pseries exps1 = (one =!= two) where
-  one = coinSeries ks1 
-  two = convolveMany (map pseries1 ks1)
-  ks1 = fromExps exps1
-
-prop_coin_vs_pseries' cexps1 = (one =!= two) where
-  one = coinSeries' aks1 
-  two = convolveMany (map pseries1' aks1)
-  aks1 = fromCoeffExps cexps1
-    
-#endif 
-
---------------------------------------------------------------------------------
 
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,12 +17,13 @@
 -- <<svg/ferrers.svg>>
 -- 
 
-{-# LANGUAGE CPP, BangPatterns #-}
+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables #-}
 module Math.Combinat.Partitions.Integer where
 
 --------------------------------------------------------------------------------
 
 import Data.List
+import Control.Monad ( liftM , replicateM )
 
 -- import Data.Map (Map)
 -- import qualified Data.Map as Map
@@ -32,6 +33,9 @@
 import Math.Combinat.Numbers (factorial,binomial,multinomial)
 import Math.Combinat.Helper
 
+import Data.Array
+import System.Random
+
 --------------------------------------------------------------------------------
 -- * Type and basic stuff
 
@@ -208,9 +212,69 @@
   go _  0  = [[]]
   go !h !n = [ a:as | a<-[1..min n h], as <- go a (n-a) ]
 
+-- | Number of partitions of @n@
 countPartitions :: Int -> Integer
-countPartitions d = countPartitions' (d,d) d
+countPartitions n = partitionCountList !! n
 
+-- | This uses 'countPartitions'', and thus is slow
+countPartitionsNaive :: Int -> Integer
+countPartitionsNaive d = countPartitions' (d,d) d
+
+--------------------------------------------------------------------------------
+
+-- | Infinite list of number of partitions of @0,1,2,...@
+--
+-- This uses the infinite product formula the generating function of partitions, recursively
+-- expanding it; it is quite fast.
+--
+-- > partitionCountList == map countPartitions [0..]
+--
+partitionCountList :: [Integer]
+partitionCountList = final where
+
+  final = go 1 (1:repeat 0) 
+
+  go !k (x:xs) = x : go (k+1) ys where
+    ys = zipWith (+) xs (take k final ++ ys)
+    -- explanation:
+    --   xs == drop k $ f (k-1)
+    --   ys == drop k $ f (k  )  
+
+{-
+
+Full explanation of 'partitionCountList':
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+let f k = productPSeries $ map (:[]) [1..k]
+
+f 0 = [1,0,0,0,0,0,0,0...]
+f 1 = [1,1,1,1,1,1,1,1...]
+f 2 = [1,1,2,2,3,3,4,4...]
+f 3 = [1,1,2,3,4,5,7,8...]
+
+observe: 
+
+* take (k+1) (f k) == take (k+1) partitionCountList
+* f (k+1) == zipWith (+) (f k) (replicate (k+1) 0 ++ f (k+1))
+
+now apply (drop (k+1)) to the second one : 
+
+* drop (k+1) (f (k+1)) == zipWith (+) (drop (k+1) $ f k) (f (k+1))
+* f (k+1) = take (k+1) final ++ drop (k+1) (f (k+1))
+
+-}
+
+--------------------------------------------------------------------------------
+
+-- | Naive infinite list of number of partitions of @0,1,2,...@
+--
+-- > partitionCountListNaive == map countPartitionsNaive [0..]
+--
+-- This is much slower than the power series expansion above.
+--
+partitionCountListNaive :: [Integer]
+partitionCountListNaive = map countPartitionsNaive [0..]
+
 -- | All integer partitions up to a given degree (that is, all integer partitions whose sum is less or equal to @d@)
 allPartitions :: Int -> [Partition]
 allPartitions d = concat [ partitions i | i <- [0..d] ]
@@ -266,7 +330,65 @@
 countPartitions' (h,w) d = sum
   [ countPartitions' (i,w-1) (d-i) | i <- [1..min d h] ] 
 
+
 ---------------------------------------------------------------------------------
+-- * Random partitions
+
+-- | Uniformly random partition of the given weight. 
+--
+-- NOTE: This algorithm is effective for small @n@-s (say @n@ up to a few hundred \/ one thousand it should work nicely),
+-- and the first time it is executed may be slower (as it needs to build the table 'partitionCountList' first)
+--
+-- Algorithm of Nijenhuis and Wilf (1975); see
+--
+-- * Knuth Vol 4A, pre-fascicle 3B, exercise 47;
+--
+-- * Nijenhuis and Wilf: Combinatorial Algorithms for Computers and Calculators, chapter 10
+--
+randomPartition :: RandomGen g => Int -> g -> (Partition, g)
+randomPartition n g = (p, g') where
+  ([p], g') = randomPartitions 1 n g
+
+-- | Generates several uniformly random partitions of @n@ at the same time.
+-- Should be a little bit faster then generating them individually.
+--
+randomPartitions 
+  :: forall g. RandomGen g 
+  => Int   -- ^ number of partitions to generate
+  -> Int   -- ^ the weight of the partitions
+  -> g -> ([Partition], g)
+randomPartitions howmany n = runRand $ replicateM howmany (worker n []) where
+
+  table = listArray (0,n) $ take (n+1) partitionCountList :: Array Int Integer
+  cnt k = table ! k
+ 
+  finish :: [(Int,Int)] -> Partition
+  finish = mkPartition . concatMap f where f (j,d) = replicate j d
+
+  fi :: Int -> Integer 
+  fi = fromIntegral
+
+  find_jd :: Int -> Integer -> (Int,Int)
+  find_jd m capm = go 0 [ (j,d) | j<-[1..n], d<-[1..div m j] ] where
+    go :: Integer -> [(Int,Int)] -> (Int,Int)
+    go !s []   = (1,1)       -- ??
+    go !s [jd] = jd          -- ??
+    go !s (jd@(j,d):rest) = 
+      if s' > capm 
+        then jd 
+        else go s' rest
+      where
+        s' = s + fi d * cnt (m - j*d)
+
+  worker :: Int -> [(Int,Int)] -> Rand g Partition
+  worker  0 acc = return $ finish acc
+  worker !m acc = do
+    capm <- randChoose (0, (fi m) * cnt m - 1)
+    let jd@(!j,!d) = find_jd m capm
+    worker (m - j*d) (jd:acc)
+
+
+---------------------------------------------------------------------------------
 -- * Dominance order 
 
 -- | @q \`dominates\` p@ returns @True@ if @q >= p@ in the dominance order of partitions
@@ -584,42 +706,3 @@
 
 --------------------------------------------------------------------------------
 
-#ifdef QUICKCHECK
-
--- * 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 ])
-
-#endif
-
---------------------------------------------------------------------------------
diff --git a/Math/Combinat/Partitions/Skew.hs b/Math/Combinat/Partitions/Skew.hs
--- a/Math/Combinat/Partitions/Skew.hs
+++ b/Math/Combinat/Partitions/Skew.hs
@@ -133,38 +133,3 @@
 
 --------------------------------------------------------------------------------
 
-#ifdef QUICKCHECK
-
-prop_dual_dual :: SkewPartition -> Bool
-prop_dual_dual sp = (dualSkewPartition (dualSkewPartition sp) == sp)
-
-prop_dual_from :: SkewPartition -> Bool
-prop_dual_from sp = (p==p' && q==q') where
-  (p,q)   = fromSkewPartition sp
-  sp'     = dualSkewPartition sp
-  (p',q') = fromSkewPartition sp'
-
-prop_from_to :: SkewPartition -> Bool
-prop_from_to sp = (mkSkewPartition (fromSkewPartition sp) == sp)
-
-prop_to_from :: (Partition,Partition) -> Bool
-prop_to_from (p,q) = 
-  case mb of
-    Nothing -> True
-    Just sp -> fromSkewPartition sp == (p,q)
-  where
-    mb = safeSkewPartition (p,q)
-
-prop_from_to_from :: SkewPartition -> Bool
-prop_from_to_from sp = (pq == pq') where
-  pq  = fromSkewPartition sp
-  sp' = mkSkewPartition pq
-  pq' = fromSkewPartition sp'
-
-prop_weight :: SkewPartition -> Bool
-prop_weight sp = (skewPartitionWeight sp == weight p - weight q) where
-  (p,q) = fromSkewPartition sp
-
-#endif
-
---------------------------------------------------------------------------------
diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs
--- a/Math/Combinat/Permutations.hs
+++ b/Math/Combinat/Permutations.hs
@@ -8,7 +8,7 @@
 -- are represented internally. Also now they act on the /right/ by default!
 --
 
-{-# LANGUAGE CPP, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-}
+{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, GeneralizedNewtypeDeriving, FlexibleContexts #-}
 module Math.Combinat.Permutations 
   ( -- * The Permutation type
     Permutation (..)
@@ -45,6 +45,13 @@
   , cycleLeft
   , cycleRight
   , reversePermutation
+    -- * Inversions
+  , inversions
+  , numberOfInversions
+  , numberOfInversionsNaive
+  , numberOfInversionsMerge
+  , bubbleSort2
+  , bubbleSort
     -- * Permutation groups
   , identity
   , inverse
@@ -60,7 +67,8 @@
   , asciiPermutation
   , asciiDisjointCycles
   , twoLineNotation 
-  , twoLineNotation'
+  , inverseTwoLineNotation
+  , genericTwoLineNotation
     -- * List of permutations
   , permutations
   , _permutations
@@ -78,11 +86,6 @@
   , permuteMultiset
   , countPermuteMultiset
   , fasc2B_algorithm_L
-
-#ifdef QUICKCHECK
-    -- * QuickCheck 
-  , checkAll
-#endif 
   ) 
   where
 
@@ -91,7 +94,8 @@
 import Control.Monad
 import Control.Monad.ST
 
-import Data.List hiding (permutations)
+import Data.List hiding ( permutations )
+import Data.Ord ( comparing )
 
 import Data.Array (Array)
 import Data.Array.ST
@@ -104,18 +108,15 @@
 import Math.Combinat.Classes
 import Math.Combinat.Helper
 import Math.Combinat.Sign
-import Math.Combinat.Numbers (factorial,binomial)
+import Math.Combinat.Numbers ( factorial , binomial )
 
 import System.Random
 
-#ifdef QUICKCHECK
-import Test.QuickCheck
-#endif
-
 --------------------------------------------------------------------------------
 -- * Types
 
--- | A permutation. Internally it is an (unboxed) array of the integers @[1..n]@. 
+-- | A permutation. Internally it is an (unboxed) array of the integers @[1..n]@, with 
+-- indexing range also being @(1,n)@. 
 --
 -- If this array of integers is @[p1,p2,...,pn]@, then in two-line 
 -- notations, that represents the permutation
@@ -186,7 +187,7 @@
   -- the zero index is an unidiomatic hack
   ar = (accumArray (+) 0 (0,n) $ map f xs) :: UArray Int Int
   f :: Int -> (Int,Int)
-  f j = if j<1 || j>n then (0,1) else (j,1)
+  f !j = if j<1 || j>n then (0,1) else (j,1)
 
 -- | Checks whether the input is a permutation of the numbers @[1..n]@.
 maybePermutation :: [Int] -> Maybe Permutation
@@ -233,16 +234,28 @@
 asciiDisjointCycles :: DisjointCycles -> ASCII
 asciiDisjointCycles (DisjointCycles cycles) = final where
   final = hCatWith VTop (HSepSpaces 1) boxes 
-  boxes = [ twoLineNotation' (f cyc) | cyc <- cycles ]
+  boxes = [ genericTwoLineNotation (f cyc) | cyc <- cycles ]
   f cyc = pairs (cyc ++ [head cyc])
 
 -- | The standard two-line notation, moving the element indexed by the top row into
 -- the place indexed by the corresponding element in the bottom row.
 twoLineNotation :: Permutation -> ASCII
-twoLineNotation (Permutation arr) = twoLineNotation' $ zip [1..] (elems arr)
+twoLineNotation (Permutation arr) = genericTwoLineNotation $ zip [1..] (elems arr)
 
-twoLineNotation' :: [(Int,Int)] -> ASCII
-twoLineNotation' xys = asciiFromLines [ topLine, botLine ] where
+-- | The inverse two-line notation, where the it\'s the bottom line 
+-- which is in standard order. The columns of this are a permutation
+-- of the columns 'twoLineNotation'.
+--
+-- Remark: the top row of @inverseTwoLineNotation perm@ is the same 
+-- as the bottom row of @twoLineNotation (inverse perm)@.
+--
+inverseTwoLineNotation :: Permutation -> ASCII
+inverseTwoLineNotation (Permutation arr) =
+  genericTwoLineNotation $ sortBy (comparing snd) $ zip [1..] (elems arr) 
+
+-- | Two-line notation for any set of numbers
+genericTwoLineNotation :: [(Int,Int)] -> ASCII
+genericTwoLineNotation xys = asciiFromLines [ topLine, botLine ] where
   topLine = "( " ++ intercalate " " us ++ " )"
   botLine = "( " ++ intercalate " " vs ++ " )"
   pairs   = [ (show x, show y) | (x,y) <- xys ]
@@ -389,6 +402,119 @@
     DisjointCycles cycles = permutationToDisjointCycles perm
 
 --------------------------------------------------------------------------------
+-- * Inversions
+
+-- | An /inversion/ of a permutation @sigma@ is a pair @(i,j)@ such that
+-- @i<j@ and @sigma(i) > sigma(j)@.
+--
+-- This functions returns the inversion of a permutation.
+--
+inversions :: Permutation -> [(Int,Int)]
+inversions (Permutation arr) = list where
+  (_,n) = bounds arr
+  list = [ (i,j) | i<-[1..n-1], j<-[i+1..n], arr!i > arr!j ]
+
+-- | Returns the number of inversions:
+--
+-- > numberOfInversions perm = length (inversions perm)
+--
+-- Synonym for 'numberOfInversionsMerge'
+--
+numberOfInversions :: Permutation -> Int
+numberOfInversions = numberOfInversionsMerge
+
+-- | Returns the number of inversions, using the merge-sort algorithm.
+-- This should be @O(n*log(n))@
+--
+numberOfInversionsMerge :: Permutation -> Int
+numberOfInversionsMerge (Permutation arr) = fst (sortCnt n $ elems arr) where
+  (_,n) = bounds arr
+                                        
+  -- | First argument is length of the list.
+  -- Returns also the inversion count.
+  sortCnt :: Int -> [Int] -> (Int,[Int])
+  sortCnt 0 _     = (0,[] )
+  sortCnt 1 [x]   = (0,[x])
+  sortCnt 2 [x,y] = if x>y then (1,[y,x]) else (0,[x,y])
+  sortCnt n xs    = mergeCnt (sortCnt k us) (sortCnt l vs) where
+    k = div n 2
+    l = n - k 
+    (us,vs) = splitAt k xs
+
+  mergeCnt :: (Int,[Int]) -> (Int,[Int]) -> (Int,[Int])
+  mergeCnt (!c,us) (!d,vs) = (c+d+e,ws) where
+
+    (e,ws) = go 0 us vs 
+
+    go !k xs [] = ( k*length xs , xs )
+    go _  [] ys = ( 0 , ys)
+    go !k xxs@(x:xs) yys@(y:ys) = if x < y
+      then let (a,zs) = go  k     xs yys in (a+k, x:zs)
+      else let (a,zs) = go (k+1) xxs  ys in (a  , y:zs)
+
+-- | Returns the number of inversions, using the definition, thus it's @O(n^2)@.
+--
+numberOfInversionsNaive :: Permutation -> Int
+numberOfInversionsNaive (Permutation arr) = length list where
+  (_,n) = bounds arr
+  list = [ (0::Int) | i<-[1..n-1], j<-[i+1..n], arr!i > arr!j ]
+
+-- | Bubble sorts breaks a permutation into the product of adjacent transpositions:
+--
+-- > multiplyMany' n (map (transposition n) $ bubbleSort2 perm) == perm
+--
+-- Note that while this is not unique, the number of transpositions 
+-- equals the number of inversions.
+--
+bubbleSort2 :: Permutation -> [(Int,Int)]
+bubbleSort2 = map f . bubbleSort where f i = (i,i+1)
+
+-- | Another version of bubble sort. An entry @i@ in the return sequence means
+-- the transposition @(i,i+1)@:
+--
+-- > multiplyMany' n (map (adjacentTransposition n) $ bubbleSort perm) == perm
+--
+bubbleSort :: Permutation -> [Int]
+bubbleSort perm@(Permutation tgt) = runST action where
+  (_,n)           = bounds tgt
+
+  action :: forall s. ST s [Int]
+  action = do
+    fwd <- newArray_ (1,n) :: ST s (STUArray s Int Int)
+    inv <- newArray_ (1,n) :: ST s (STUArray s Int Int)
+    forM_ [1..n] $ \i -> writeArray fwd i i
+    forM_ [1..n] $ \i -> writeArray inv i i
+
+    list <- forM [1..n] $ \x -> do
+
+      let k = tgt ! x        -- we take the number which will be at the @x@-th position at the end
+      i <- readArray inv k   -- number @k@ is at the moment at position @i@
+      let j = x              -- but the final place is at @x@      
+
+      let swaps = move i j
+      forM_ swaps $ \y -> do
+
+        a <- readArray fwd  y
+        b <- readArray fwd (y+1)
+        writeArray fwd (y+1) a
+        writeArray fwd  y    b
+
+        u <- readArray inv a
+        v <- readArray inv b
+        writeArray inv b u
+        writeArray inv a v
+
+      return swaps
+  
+    return (concat list)
+
+  move :: Int -> Int -> [Int]
+  move !i !j
+    | j == i  = []
+    | j >  i  = [i..j-1]
+    | j <  i  = [i-1,i-2..j]
+
+--------------------------------------------------------------------------------
 -- * Some concrete permutations
 
 -- | The permutation @[n,n-1,n-2,...,2,1]@. Note that it is the inverse of itself.
@@ -710,6 +836,7 @@
 --   The order is lexicographic.  
 fasc2B_algorithm_L :: (Eq a, Ord a) => [a] -> [[a]] 
 fasc2B_algorithm_L xs = unfold1 next (sort xs) where
+
   -- next :: [a] -> Maybe [a]
   next xs = case findj (reverse xs,[]) of 
     Nothing -> Nothing
@@ -725,170 +852,11 @@
   findj ( [] , _ ) = Nothing
   
   -- inc :: a -> [a] -> ([a],[a]) -> [a]
-  inc u us ( (x:xs) , yys ) = if u >= x
+  inc !u us ( (x:xs) , yys ) = if u >= x
     then inc u us ( xs , x : yys ) 
     else reverse (x:us)  ++ reverse (u:yys) ++ xs
   inc _ _ ( [] , _ ) = error "permute: should not happen"
       
 --------------------------------------------------------------------------------
 
-#ifdef QUICKCHECK
-
-minPermSize = 1
-maxPermSize = 123
-
-newtype Elem = Elem Int deriving Eq
-newtype Nat  = Nat { fromNat :: Int } deriving (Eq,Ord,Show,Num,Random)
-
-naturalSet :: Permutation -> Array Int Elem
-naturalSet perm = listArray (1,n) [ Elem i | i<-[1..n] ] where
-  n = permutationSize perm
-
-permInternalSet :: Permutation -> Array Int Elem
-permInternalSet perm@(Permutation arr) = listArray (1,n) [ Elem (arr!i) | i<-[1..n] ] where
-  n = permutationSize perm
-
-sameSize :: Permutation ->  Permutation -> Bool
-sameSize perm1 perm2 = ( permutationSize perm1 == permutationSize perm2)
-
-newtype CyclicPermutation = Cyclic { fromCyclic :: Permutation } deriving Show
-
-data SameSize = SameSize Permutation Permutation deriving Show
-
-instance Random Permutation where
-  random g = randomPermutation size g1 where
-    (size,g1) = randomR (minPermSize,maxPermSize) g
-  randomR _ = random
-
-instance Random CyclicPermutation where
-  random g = (Cyclic cycl,g2) where
-    (size,g1) = randomR (minPermSize,maxPermSize) g
-    (cycl,g2) = randomCyclicPermutation size g1
-  randomR _ = random
-
-instance Random DisjointCycles where
-  random g = (disjcyc,g2) where
-    (size,g1) = randomR (minPermSize,maxPermSize) g
-    (perm,g2) = randomPermutation size g1
-    disjcyc   = permutationToDisjointCycles perm
-  randomR _ = random
-
-instance Random SameSize where
-  random g = (SameSize prm1 prm2, g3) where
-    (size,g1) = randomR (minPermSize,maxPermSize) g
-    (prm1,g2) = randomPermutation size g1 
-    (prm2,g3) = randomPermutation size g2
-  randomR _ = random
-
-instance Arbitrary Nat where
-  arbitrary = choose (Nat 0 , Nat 50)
-     
-instance Arbitrary Permutation       where arbitrary = choose undefined
-instance Arbitrary CyclicPermutation where arbitrary = choose undefined
-instance Arbitrary DisjointCycles    where arbitrary = choose undefined
-instance Arbitrary SameSize          where arbitrary = choose undefined
-
--- | Runs all quickCheck tests
-checkAll :: IO ()
-checkAll = do
-  let f :: Testable a => a -> IO ()
-      f = quickCheck
-
-  f prop_disjcyc_1
-  f prop_disjcyc_2 
-
-  f prop_disjcyc_Mathematica
-
-  f prop_randCyclic
-  f prop_inverse
-
-  f prop_mulPerm
-  f prop_mulPermLeft
-  f prop_mulPermRight
-
-  f prop_perm
-  f prop_permLeft
-  f prop_permRight
-  f prop_permLeftRight
-
-  f prop_cycleLeft
-  f prop_cycleRight
-
-  f prop_mulSign      
-  f prop_invMul
-  f prop_cyclSign
-  f prop_permIsPerm
-  f prop_isEven
-          
-prop_disjcyc_1 perm = ( perm == disjointCyclesToPermutation n (permutationToDisjointCycles perm) )
-  where n = permutationSize perm
-
-prop_disjcyc_2 k dcyc = ( dcyc == permutationToDisjointCycles (disjointCyclesToPermutation n dcyc) )
-  where 
-    n = fromNat k + m 
-    m = case fromDisjointCycles dcyc of
-      []  -> 1
-      xxs -> maximum (concat xxs)
-
--- PermutationCycles[ { 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 } ]
--- Cycles           [ {{1, 12, 18, 19, 14, 21, 23, 13, 22}, {2, 15, 8, 9, 20, 16, 10, 3, 5}, {4, 6, 7, 17}} ]
-prop_disjcyc_Mathematica = (permutationToDisjointCycles   perm == disjcyc) 
-                        && (disjointCyclesToPermutation n disjcyc == perm)
-  where
-    n       = permutationSize perm
-    perm    = toPermutation  [ 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 ]
-    disjcyc = DisjointCycles [ [1, 12, 18, 19, 14, 21, 23, 13, 22], [2, 15, 8, 9, 20, 16, 10, 3, 5], [4, 6, 7, 17] ]
-
-xperm    = toPermutation  [ 12, 15, 5, 6, 2, 7, 17, 9, 20, 3, 11, 18, 22, 21, 8, 10, 4, 19, 14, 16, 23, 1, 13 ]
-xdisjcyc = DisjointCycles [ [1, 12, 18, 19, 14, 21, 23, 13, 22], [2, 15, 8, 9, 20, 16, 10, 3, 5], [4, 6, 7, 17] ]
-
-prop_randCyclic cycl = ( isCyclicPermutation (fromCyclic cycl) )
-
-prop_inverse perm = ( perm == inverse (inverse perm) ) 
-
-prop_mulPerm (SameSize perm1 perm2) = 
-    ( permute perm2 (permute perm1 set) == permute (perm1 `multiply` perm2) set ) 
-  where 
-    set = naturalSet perm1
-
-prop_mulPermRight (SameSize perm1 perm2) = 
-    ( permuteRight perm2 (permuteRight perm1 set) == permuteRight (perm1 `multiply` perm2) set ) 
-  where 
-    set = naturalSet perm1
-
-prop_mulPermLeft (SameSize perm1 perm2) = 
-    ( permuteLeft perm2 (permuteLeft perm1 set) == permuteLeft (perm2 `multiply` perm1) set ) 
-  where 
-    set = naturalSet perm1
-
-prop_perm          perm = permute perm (naturalSet perm) == permInternalSet perm
-prop_permLeft      perm = permuteLeft  perm (permInternalSet perm) == naturalSet perm
-prop_permRight     perm = permuteRight perm (naturalSet perm) == permInternalSet perm
-prop_permLeftRight perm = permuteLeft (inverse perm) (naturalSet perm) == permuteRight (perm) (naturalSet perm) 
-
-prop_cycleLeft  = permuteList (cycleLeft  5) "abcde" == "bcdea"
-prop_cycleRight = permuteList (cycleRight 5) "abcde" == "eabcd"
-
-prop_mulSign (SameSize perm1 perm2) = 
-    ( sgn perm1 * sgn perm2 == sgn (perm1 `multiply` perm2) ) 
-  where 
-    sgn = signValue . signOfPermutation :: Permutation -> Int
-
-prop_invMul (SameSize perm1 perm2) =   
-  ( inverse perm2 `multiply` inverse perm1 == inverse (perm1 `multiply` perm2) ) 
-
-prop_cyclSign cycl = ( isEvenPermutation perm == odd n ) where
-  perm = fromCyclic cycl
-  n = permutationSize perm
-  
-prop_permIsPerm perm = ( isPermutation (fromPermutation perm) ) 
-
-prop_isEven perm = ( isEvenPermutation perm == isEvenAlternative perm ) where
-  isEvenAlternative p = 
-    even $ sum $ map (\x->x-1) $ map length $ fromDisjointCycles $ permutationToDisjointCycles p
-
-
-#endif
-
---------------------------------------------------------------------------------
 
diff --git a/Math/Combinat/Sign.hs b/Math/Combinat/Sign.hs
--- a/Math/Combinat/Sign.hs
+++ b/Math/Combinat/Sign.hs
@@ -7,11 +7,12 @@
 --------------------------------------------------------------------------------
 
 import Data.Monoid
+import System.Random
 
 --------------------------------------------------------------------------------
 
 data Sign
-  = Plus
+  = Plus                            -- hmm, this way @Plus < Minus@, not sure about that
   | Minus
   deriving (Eq,Ord,Show,Read)
 
@@ -19,6 +20,10 @@
   mempty  = Plus
   mappend = mulSign
   mconcat = productOfSigns
+
+instance Random Sign where
+  random        g = let (b,g') = random g in (if b    then Plus else Minus, g')
+  randomR (u,v) g = let (y,g') = random g in (if u==v then u    else y    , g') 
 
 isPlus, isMinus :: Sign -> Bool
 isPlus  s = case s of { Plus  -> True ; _ -> False }
diff --git a/Math/Combinat/Tableaux/LittlewoodRichardson.hs b/Math/Combinat/Tableaux/LittlewoodRichardson.hs
--- a/Math/Combinat/Tableaux/LittlewoodRichardson.hs
+++ b/Math/Combinat/Tableaux/LittlewoodRichardson.hs
@@ -328,7 +328,7 @@
 -- | Computes the expansion of the product of Schur polynomials @s[mu]*s[nu]@ using the
 -- Littlewood-Richardson rule. Note: this is symmetric in the two arguments.
 --
--- Based on the wikipedia article <https://en.wikipedia.org/wiki/Littlewood–Richardson_rule>
+-- Based on the wikipedia article <https://en.wikipedia.org/wiki/Littlewood-Richardson_rule>
 --
 -- > lrMult mu nu == Map.fromList list where
 -- >   lamw = weight nu + weight mu
diff --git a/Math/Combinat/Tableaux/Skew.hs b/Math/Combinat/Tableaux/Skew.hs
--- a/Math/Combinat/Tableaux/Skew.hs
+++ b/Math/Combinat/Tableaux/Skew.hs
@@ -221,35 +221,3 @@
 
 --------------------------------------------------------------------------------
 
-#ifdef QUICKCHECK
-
-prop_dual_dual :: SkewTableau Int -> Bool
-prop_dual_dual st = (dualSkewTableau (dualSkewTableau st) == st)
-
-prop_rowWord :: SkewTableau Int -> Bool
-prop_rowWord st = (fillSkewPartitionWithRowWord shape content == st) where
-  shape   = skewShape st
-  content = skewTableauRowWord st
-
-prop_columnWord :: SkewTableau Int -> Bool
-prop_columnWord st = (fillSkewPartitionWithColumnWord shape content == st) where
-  shape   = skewShape st
-  content = skewTableauColumnWord st
-
-prop_fill_shape :: SkewPartition -> Bool
-prop_fill_shape shape = (shape == shape') where
-  tableau = fillSkewPartitionWithColumnWord shape [1..]
-  shape'  = skewShape tableau
-
-prop_semistandard :: SkewPartition -> Bool
-prop_semistandard shape = and 
-  [ isSemiStandardSkewTableau st 
-  | n  <- [1..nn] 
-  , st <- semiStandardSkewTableaux n shape
-  ]
-  where
-    nn = skewPartitionWeight shape
-
-#endif
-
---------------------------------------------------------------------------------
diff --git a/Math/Combinat/TypeLevel.hs b/Math/Combinat/TypeLevel.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/TypeLevel.hs
@@ -0,0 +1,117 @@
+
+-- | Type-level hackery.
+--
+-- This module is used for groups whose parameters are encoded as type-level natural numbers,
+-- for example finite cyclic groups, free groups, symmetric groups and braid groups.
+--
+
+{-# LANGUAGE PolyKinds, DataKinds, KindSignatures, ScopedTypeVariables, 
+             ExistentialQuantification, Rank2Types #-}
+
+module Math.Combinat.TypeLevel 
+  ( -- * Proxy
+    Proxy(..)
+  , proxyUndef
+  , proxyOf
+  , proxyOf1
+  , proxyOf2
+  , asProxyTypeOf   -- defined in Data.Proxy
+  , asProxyTypeOf1
+    -- * Type-level naturals as type arguments
+  , typeArg 
+  , iTypeArg
+    -- * Hiding the type parameter
+  , Some (..)
+  , withSome , withSomeM
+  , selectSome , selectSomeM
+  , withSelected , withSelectedM
+  )
+  where
+
+--------------------------------------------------------------------------------
+
+import Data.Proxy ( Proxy(..) , asProxyTypeOf )
+import GHC.TypeLits
+
+import Math.Combinat.Numbers.Primes ( isProbablyPrime )
+
+--------------------------------------------------------------------------------
+-- * Proxy
+
+proxyUndef :: Proxy a -> a
+proxyUndef _ = error "proxyUndef"
+
+proxyOf :: a -> Proxy a
+proxyOf _ = Proxy
+
+proxyOf1 :: f a -> Proxy a
+proxyOf1 _ = Proxy
+
+proxyOf2 :: g (f a) -> Proxy a
+proxyOf2 _ = Proxy
+
+asProxyTypeOf1 :: f a -> Proxy a -> f a 
+asProxyTypeOf1 y _ = y
+
+--------------------------------------------------------------------------------
+-- * Type-level naturals as type arguments
+
+typeArg :: KnownNat n => f (n :: Nat) -> Integer
+typeArg = natVal . proxyOf1
+
+iTypeArg :: KnownNat n => f (n :: Nat) -> Int
+iTypeArg = fromIntegral . typeArg
+
+--------------------------------------------------------------------------------
+-- * Hiding the type parameter
+
+-- | Hide the type parameter of a functor. Example: @Some Braid@
+data Some f = forall n. KnownNat n => Some (f n)
+
+-- | Uses the value inside a 'Some'
+withSome :: Some f -> (forall n. KnownNat n => f n -> a) -> a
+withSome some polyFun = case some of { Some stuff -> polyFun stuff }
+
+-- | Monadic version of 'withSome'
+withSomeM :: Monad m => Some f -> (forall n. KnownNat n => f n -> m a) -> m a
+withSomeM some polyAct = case some of { Some stuff -> polyAct stuff }
+
+-- | Given a polymorphic value, we select at run time the
+-- one specified by the second argument
+selectSome :: Integral int => (forall n. KnownNat n => f n) -> int -> Some f
+selectSome poly n = case someNatVal (fromIntegral n :: Integer) of
+  Nothing   -> error "selectSome: not a natural number"
+  Just snat -> case snat of
+    SomeNat pxy -> Some (asProxyTypeOf1 poly pxy)
+
+-- | Monadic version of 'selectSome'
+selectSomeM :: forall m f int. (Integral int, Monad m) => (forall n. KnownNat n => m (f n)) -> int -> m (Some f)
+selectSomeM runPoly n = case someNatVal (fromIntegral n :: Integer) of
+  Nothing   -> error "selectSomeM: not a natural number"
+  Just snat -> case snat of
+    SomeNat pxy -> do
+      poly <- runPoly 
+      return $ Some (asProxyTypeOf1 poly pxy)
+
+-- | Combination of 'selectSome' and 'withSome': we make a temporary structure
+-- of the given size, but we immediately consume it.
+withSelected 
+  :: Integral int 
+  => (forall n. KnownNat n => f n -> a) 
+  -> (forall n. KnownNat n => f n) 
+  -> int 
+  -> a
+withSelected f x n = withSome (selectSome x n) f
+
+-- | (Half-)monadic version of 'withSelected'
+withSelectedM 
+  :: forall m f int a. (Integral int, Monad m) 
+  => (forall n. KnownNat n => f n -> a) 
+  -> (forall n. KnownNat n => m (f n)) 
+  -> int 
+  -> m a
+withSelectedM f mx n = do 
+  s <- selectSomeM mx n
+  return (withSome s f)
+
+--------------------------------------------------------------------------------
diff --git a/combinat.cabal b/combinat.cabal
--- a/combinat.cabal
+++ b/combinat.cabal
@@ -1,5 +1,5 @@
 Name:                combinat
-Version:             0.2.8.0
+Version:             0.2.8.1
 Synopsis:            Generate and manipulate various combinatorial objects.
 Description:         A collection of functions to generate, count, manipulate
                      and visualize all kinds of combinatorial objects like 
@@ -22,18 +22,11 @@
 extra-source-files:  svg/*.svg
                      svg/src/gen_figures.hs                     
 
-Flag withQuickCheck
-  Description: Compile with the QuickCheck tests. 
-  default: False
-  
+ 
 Library
 
   Build-Depends:       base >= 4 && < 5, array >= 0.5, containers, random, transformers
 
-  if flag(withQuickCheck)
-    Build-Depends:       QuickCheck
-    cpp-options:         -DQUICKCHECK
-
   Exposed-Modules:     Math.Combinat
                        Math.Combinat.Classes
                        Math.Combinat.Numbers
@@ -68,6 +61,7 @@
                        Math.Combinat.LatticePaths
                        Math.Combinat.ASCII
                        Math.Combinat.Helper
+                       Math.Combinat.TypeLevel
 
   Default-Extensions:  CPP, BangPatterns
   Other-Extensions:    MultiParamTypeClasses, ScopedTypeVariables, 
@@ -80,3 +74,17 @@
 
   ghc-options:         -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports
     
+
+test-suite combinat-tests
+                      
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             TestSuite.hs
+
+  build-depends:       base >= 4 && < 5, array >= 0.5, containers, random, transformers,
+                       combinat,
+                       QuickCheck >= 2, test-framework, test-framework-quickcheck2
+
+  Default-Language:    Haskell2010
+  Default-Extensions:  CPP, BangPatterns
+
diff --git a/test/TestSuite.hs b/test/TestSuite.hs
new file mode 100644
--- /dev/null
+++ b/test/TestSuite.hs
@@ -0,0 +1,41 @@
+
+module Main where
+
+--------------------------------------------------------------------------------
+
+import Test.Framework
+import Test.Framework.Providers.QuickCheck2
+
+import Tests.Permutations       ( testgroup_Permutations      )
+import Tests.Partitions.Integer ( testgroup_IntegerPartitions )
+import Tests.Partitions.Skew    ( testgroup_SkewPartitions    )
+import Tests.Braid              ( testgroup_Braid 
+                                , testgroup_Braid_NF          )
+import Tests.Series             ( testgroup_PowerSeries       )
+import Tests.SkewTableaux       ( testgroup_SkewTableaux      )
+import Tests.Thompson           ( testgroup_ThompsonF         )
+import Tests.LatticePaths       ( testgroup_LatticePaths      )
+
+--------------------------------------------------------------------------------
+
+main :: IO ()
+main = defaultMain tests
+
+tests :: [Test]
+tests = 
+  [ testgroup_Permutations
+  , testGroup "Partitions" 
+      [ testgroup_IntegerPartitions
+      , testgroup_SkewPartitions
+      ]
+  , testgroup_SkewTableaux
+  , testgroup_ThompsonF
+  , testgroup_LatticePaths
+  , testGroup "Braids" 
+      [ testgroup_Braid 
+      , testgroup_Braid_NF 
+      ]
+  , testgroup_PowerSeries  
+  ]
+
+--------------------------------------------------------------------------------