diff --git a/Math/Combinat/Numbers.hs b/Math/Combinat/Numbers.hs
--- a/Math/Combinat/Numbers.hs
+++ b/Math/Combinat/Numbers.hs
@@ -6,6 +6,8 @@
 
 module Math.Combinat.Numbers where
 
+--------------------------------------------------------------------------------
+
 import Data.Array
 
 --------------------------------------------------------------------------------
@@ -124,24 +126,6 @@
     f k = toRational (paritySign (n+k) * factorial k * stirling2nd n k) 
         / toRational (k+1)
 
---------------------------------------------------------------------------------
--- * Power series
-
--- | Power series expansion of 
--- 
--- > @1 / ( (1-x^a_1) * (1-x^a_2) * ... * (1-x^a_n) )@
---
--- Example:
---
--- @(coinSeries [2,3,5]) !! k@ is the number of ways 
--- to pay @k@ dollars with coins of two, three and five dollars.
---
--- TODO: better name?
-coinSeries :: [Int] -> [Integer]
-coinSeries []  = 1 : repeat 0
-coinSeries (k:ks) = xs where
-  xs = zipWith (+) (coinSeries ks) (replicate k 0 ++ xs) 
-  
 --------------------------------------------------------------------------------
 
  
diff --git a/Math/Combinat/Numbers/Series.hs b/Math/Combinat/Numbers/Series.hs
new file mode 100644
--- /dev/null
+++ b/Math/Combinat/Numbers/Series.hs
@@ -0,0 +1,452 @@
+
+-- | Some basic power series expansions.
+-- This module is not re-exported by "Math.Combinat".
+--
+-- Note: the \"@convolveWithXXX@\" functions are much faster than the equivalent
+-- @(XXX \`convolve\`)@!
+-- 
+-- TODO: better names for these functions.
+--
+
+{-# LANGUAGE CPP, GeneralizedNewtypeDeriving #-}
+module Math.Combinat.Numbers.Series where
+
+--------------------------------------------------------------------------------
+
+import Data.List
+
+#ifdef QUICKCHECK
+import System.Random
+import Test.QuickCheck
+#endif
+
+--------------------------------------------------------------------------------
+
+-- | 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!
+convolve :: Num a => [a] -> [a] -> [a]
+convolve xs1 ys1 = res where
+  res = [ foldl' (+) 0 (zipWith (*) xs (reverse (take n ys)))
+        | n<-[1..] 
+        ]
+  xs = xs1 ++ repeat 0
+  ys = ys1 ++ repeat 0
+
+-- | Convolution of many series. Still slow!
+convolveMany :: Num a => [[a]] -> [a]
+convolveMany []  = 1 : repeat 0
+convolveMany xss = foldl1 convolve xss
+
+--------------------------------------------------------------------------------
+-- * \"Coin\" series
+
+-- | Power series expansion of 
+-- 
+-- > 1 / ( (1-x^k_1) * (1-x^k_2) * ... * (1-x^k_n) )
+--
+-- Example:
+--
+-- @(coinSeries [2,3,5])!!k@ is the number of ways 
+-- to pay @k@ dollars with coins of two, three and five dollars.
+--
+-- TODO: better name?
+coinSeries :: [Int] -> [Integer]
+coinSeries [] = 1 : repeat 0
+coinSeries (k:ks) = xs where
+  xs = zipWith (+) (coinSeries ks) (replicate k 0 ++ xs) 
+
+-- | Generalization of the above to include coefficients: expansion of 
+--  
+-- > 1 / ( (1-a_1*x^k_1) * (1-a_2*x^k_2) * ... * (1-a_n*x^k_n) ) 
+-- 
+coinSeries' :: Num a => [(a,Int)] -> [a]
+coinSeries' [] = 1 : repeat 0
+coinSeries' ((a,k):aks) = xs where
+  xs = zipWith (+) (coinSeries' aks) (replicate k 0 ++ map (*a) xs) 
+
+convolveWithCoinSeries :: [Int] -> [Integer] -> [Integer]
+convolveWithCoinSeries ks series1 = worker ks where
+  series = series1 ++ repeat 0
+  worker [] = series
+  worker (k:ks) = xs where
+    xs = zipWith (+) (worker ks) (replicate k 0 ++ xs)
+
+convolveWithCoinSeries' :: Num a => [(a,Int)] -> [a] -> [a]
+convolveWithCoinSeries' ks series1 = worker ks where
+  series = series1 ++ repeat 0
+  worker [] = series
+  worker ((a,k):aks) = xs where
+    xs = zipWith (+) (worker aks) (replicate k 0 ++ map (*a) xs)
+
+--------------------------------------------------------------------------------
+-- * Reciprocals of products of polynomials
+
+-- | Convolution of many 'pseries', that is, the expansion of the reciprocal
+-- of a product of polynomials
+productPSeries :: [[Int]] -> [Integer]
+productPSeries = foldl (flip convolveWithPSeries) unitSeries
+
+-- | The same, with coefficients.
+productPSeries' :: Num a => [[(a,Int)]] -> [a]
+productPSeries' = foldl (flip convolveWithPSeries') unitSeries
+
+convolveWithProductPSeries :: [[Int]] -> [Integer] -> [Integer]
+convolveWithProductPSeries kss ser = foldl (flip convolveWithPSeries) ser kss
+
+-- | This is the most general function in this module; all the others
+-- are special cases of this one.  
+convolveWithProductPSeries' :: Num a => [[(a,Int)]] -> [a] -> [a] 
+convolveWithProductPSeries' akss ser = foldl (flip convolveWithPSeries') ser akss
+  
+--------------------------------------------------------------------------------
+-- * Reciprocals of polynomials
+
+-- 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)
+--
+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)
+--
+convolveWithPSeries :: [Int] -> [Integer] -> [Integer]
+convolveWithPSeries ks series1 = ys where 
+  series = series1 ++ repeat 0 
+  ys = worker ks ys 
+  worker [] _ = series 
+  worker (k:ks) ys = xs where
+    xs = zipWith (+) (replicate k 0 ++ ys) (worker ks ys)
+
+--------------------------------------------------------------------------------
+--  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)
+--
+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)
+--
+convolveWithPSeries' :: Num a => [(a,Int)] -> [a] -> [a]
+convolveWithPSeries' aks series1 = ys where 
+  series = series1 ++ repeat 0 
+  ys = worker aks ys 
+  worker [] _ = series
+  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
+
+-- | Convolve with (the expansion of) 
+--
+-- > 1 / (1 +- x^k_1 +- x^k_2 +- x^k_3 +- ... +- x^k_n)
+--
+-- Should be faster than using `convolveWithPSeries'`.
+-- Note: 'Plus' corresponds to the coefficient @-1@ in `pseries'` (since
+-- there is a minus sign in the definition there)!
+convolveWithSignedPSeries :: [(Sign,Int)] -> [Integer] -> [Integer]
+convolveWithSignedPSeries aks series1 = ys where 
+  series = series1 ++ repeat 0 
+  ys = worker aks ys 
+  worker [] _ = series
+  worker ((a,k):aks) ys = xs where
+    xs = case a of
+      Minus -> zipWith (+) one two 
+      Plus  -> zipWith (-) one two
+    one = worker aks ys
+    two = replicate k 0 ++ ys
+     
+--------------------------------------------------------------------------------
+
+#ifdef QUICKCHECK
+
+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/Permutations.hs b/Math/Combinat/Permutations.hs
--- a/Math/Combinat/Permutations.hs
+++ b/Math/Combinat/Permutations.hs
@@ -45,6 +45,11 @@
   , permuteMultiset
   , countPermuteMultiset
   , fasc2B_algorithm_L
+
+#ifdef QUICKCHECK
+    -- * QuickCheck 
+  , checkAll
+#endif QUICKCHECK
   ) 
   where
 
@@ -440,6 +445,22 @@
 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_disjcyc1
+  f prop_disjcyc2 
+  f prop_randCyclic
+  f prop_inverse
+  f prop_mulPerm
+  f prop_mulSign      
+  f prop_invMul
+  f prop_cyclSign
+  f prop_permIsPerm
+  f prop_isEven
+          
 prop_disjcyc1 perm = ( perm == disjointCyclesToPermutation n (permutationToDisjointCycles perm) )
   where n = permutationSize perm
 prop_disjcyc2 k dcyc = ( dcyc == permutationToDisjointCycles (disjointCyclesToPermutation n dcyc) )
diff --git a/combinat.cabal b/combinat.cabal
--- a/combinat.cabal
+++ b/combinat.cabal
@@ -1,5 +1,5 @@
 Name:                combinat
-Version:             0.2.3.1
+Version:             0.2.4
 Synopsis:            Generation of various combinatorial objects.
 Description:         A collection of functions to generate combinatorial
                      objects like partitions, combinations, permutations,
@@ -43,6 +43,7 @@
 
   Exposed-Modules:     Math.Combinat, 
                        Math.Combinat.Numbers,
+                       Math.Combinat.Numbers.Series,
                        Math.Combinat.Sets,
                        Math.Combinat.Tuples, 
                        Math.Combinat.Combinations,
@@ -57,7 +58,8 @@
   
   Other-Modules:       Math.Combinat.Helper
 
-  Extensions:          MultiParamTypeClasses, ScopedTypeVariables, CPP
+  Extensions:          CPP, MultiParamTypeClasses, ScopedTypeVariables, 
+                       GeneralizedNewtypeDeriving 
 
   Hs-Source-Dirs:      .