diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2012, Anders Claesson
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Anders Claesson nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Math/Sym.hs b/Math/Sym.hs
new file mode 100644
--- /dev/null
+++ b/Math/Sym.hs
@@ -0,0 +1,342 @@
+-- |
+-- Module      : Math.Sym
+-- Copyright   : (c) Anders Claesson 2012
+-- License     : BSD-style
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+-- 
+-- Provides an efficient definition of standard permutations,
+-- 'StPerm', together with an abstract class, 'Perm', whose
+-- functionality is largely inherited from 'StPerm' using a group
+-- action and the standardization map.
+
+module Math.Sym
+    (
+    -- * Standard permutations
+      StPerm
+    , toVector        -- :: StPerm -> Vector Int
+    , fromVector      -- :: Vector Int -> StPerm
+    , toList          -- :: StPerm -> [Int]
+    , fromList        -- :: [Int] -> StPerm
+    , (/-/)           -- :: StPerm -> StPerm -> StPerm
+    , unrankStPerm    -- :: Int -> Integer -> StPerm
+    , sym             -- :: Int -> [StPerm]
+
+    -- * The permutation typeclass
+    , Perm (..)
+
+    -- * Generalize
+    , generalize      -- :: Perm a => (StPerm -> StPerm) -> a -> a
+
+    -- * Generating permutations
+    , unrankPerm      -- :: Perm a => a -> Integer -> a
+    , randomPerm      -- :: Perm a => a -> IO a
+    , perms           -- :: Perm a => a -> [a]
+
+    -- * Sorting operators
+    , stackSort       -- :: Perm a => a -> a
+    , bubbleSort      -- :: Perm a => a -> a
+
+    -- * Permutation patterns
+    , copies          -- :: Perm a => StPerm -> a -> [Set]
+    , avoids          -- :: Perm a => [StPerm] -> a -> Bool
+    , avoiders        -- :: Perm a => [StPerm] -> [a] -> [a]
+    , av              -- :: [StPerm] -> Int -> [StPerm]
+
+    -- * Subsets
+    , Set
+    , subsets         -- :: Int -> Int -> [Set]
+    ) where
+
+import Control.Monad (liftM)
+import Data.Ord (comparing)
+import Data.Monoid (Monoid(..))
+import Data.Bits (Bits, bitSize, testBit, popCount, shiftL)
+import Data.List (sort, sortBy)
+import Data.Vector.Storable (Vector)
+import qualified Data.Vector.Storable as SV (Vector, toList, fromList, fromListN, empty, map, (++))
+import qualified Math.Sym.Internal as I
+import Foreign.C.Types (CUInt(..))
+
+
+-- Standard permutations
+-- ---------------------
+
+-- | By a /standard permutation/ we shall mean a permutations of
+-- @[0..k-1]@.
+newtype StPerm = StPerm { perm0 :: I.Perm0 } deriving (Eq, Ord)
+
+instance Show StPerm where
+    show = show . toVector
+
+instance Monoid StPerm where
+    mempty = fromVector SV.empty
+    mappend u v = fromVector $ (SV.++) u' v'
+        where
+          u' = toVector u
+          v' = SV.map ( + size u) $ toVector v
+
+-- | Convert a standard permutation to a vector.
+toVector :: StPerm -> Vector Int
+toVector = perm0
+
+-- | Convert a vector to a standard permutation. The vector should a
+-- permutation of the elements @[0..k-1]@ for some positive @k@. No
+-- checks for this are done.
+fromVector :: Vector Int -> StPerm
+fromVector = StPerm
+
+-- | Convert a standard permutation to a list.
+toList :: StPerm -> [Int]
+toList = SV.toList . toVector
+
+-- | Convert a list to a standard permutation. The list should a
+-- permutation of the elements @[0..k-1]@ for some positive @k@. No
+-- checks for this are done.
+fromList :: [Int] -> StPerm
+fromList = fromVector . SV.fromList
+
+infixl 6 /-/
+
+-- | The /skew sum/ of two permutations. (A definition of the
+-- /direct sum/ is provided by the Monoid instance.)
+(/-/) :: StPerm -> StPerm -> StPerm
+u /-/ v = fromVector $ (SV.++) u' v'
+    where
+      u' = SV.map ( + size v) $ toVector u
+      v' = toVector v
+
+-- | @unrankStPerm n rank@ is the @rank@-th (Myrvold & Ruskey)
+-- permutation of @[0..n-1]@. E.g.,
+-- 
+-- > unrankStPerm 16 19028390 == fromList [6,15,4,11,7,8,9,2,5,0,10,3,12,13,14,1]
+-- 
+unrankStPerm :: Int -> Integer -> StPerm
+unrankStPerm n = fromVector . I.unrankPerm n
+
+-- | The list of standard permutations of the given size (the symmetric group). E.g.,
+-- 
+-- > sym 2 == [fromList [0,1], fromList [1,0]]
+-- 
+sym :: Int -> [StPerm]
+sym n = map (unrankStPerm n) [0 .. product [1 .. toInteger n] - 1]
+
+
+-- The permutation typeclass
+-- -------------------------
+
+-- | The class of permutations. Minimal complete definition: 'st' and
+-- 'act'. The default implementations of 'size' and 'idperm' can be
+-- somewhat slow, so you may want to consider implementing them as
+-- well.
+class Perm a where
+
+    -- | The standardization map. If there is an underlying linear
+    -- order on @a@ then @st@ is determined by the unique order
+    -- preserving map from @[0..]@ to that order. In any case, the
+    -- standardization map should be equivariant with respect to the
+    -- group action defined below; i.e., it should hold that
+    -- 
+    -- > st (u `act` v) == u `act` st v
+    -- 
+    st :: a -> StPerm
+
+    -- | A (left) /group action/ of 'StPerm' on @a@. As for any group
+    -- action it should hold that
+    -- 
+    -- > (u `act` v) `act` w == u `act` (v `act` w)   &&   idperm u `act` v == v
+    -- 
+    act :: StPerm -> a -> a
+
+    -- | The size of a permutation. The default implementation derived from
+    -- 
+    -- > size == size . st
+    -- 
+    -- This is not a circular definition as 'size' on 'StPerm' is
+    -- implemented independently. If the implementation of 'st' is
+    -- slow, then it can be worth while to override the standard
+    -- definiton; any implementation should, however, satisfy the
+    -- identity above.
+    {-# INLINE size #-}
+    size :: a -> Int
+    size = size . st
+
+    -- | The identity permutation on the same underlying set as the
+    -- given permutation. It should hold that
+    -- 
+    -- > st (idperm u) == idperm (st u)
+    -- 
+    -- Group theoretically, it should also hold that @u . inverse u ==
+    -- idperm u@. In terms of the group action this means
+    -- 
+    -- > idperm u == inverse (st u) `act` u
+    -- 
+    -- and this is the default implementation.
+    {-# INLINE idperm #-}
+    idperm :: a -> a
+    idperm u = inverse (st u) `act` u
+
+    -- | The group theoretical inverse. It should hold that
+    -- 
+    -- > inverse u == inverse (st u) `act` idperm u
+    -- 
+    -- and this is the default implementation.
+    {-# INLINE inverse #-}
+    inverse :: a -> a
+    inverse u = inverse (st u) `act` idperm u
+
+    -- | Predicate determining if two permutations are
+    -- order-isomorphic. The default implementation uses
+    -- 
+    -- > u `ordiso` v  ==  u == st v
+    -- 
+    -- Equivalently, one could use
+    -- 
+    -- > u `ordiso` v  ==  inverse u `act` v == idperm v
+    -- 
+    {-# INLINE ordiso #-}
+    ordiso :: StPerm -> a -> Bool
+    ordiso u v = u == st v
+
+instance Perm StPerm where
+    st         = id
+    act u v    = fromVector $ I.act (toVector u) (toVector v)
+    size       = I.size . toVector
+    idperm     = fromVector . I.idperm . size
+    inverse    = fromVector . I.inverse . toVector
+    ordiso     = (==)
+
+-- Auxiliary function: @w = act' u v@ iff @w[u[i]] = v[i]@.
+-- Caveat: @act'@ is not a proper group action.
+act' :: Ord a => [a] -> [b] -> [b]
+act' u = map snd . sortBy (comparing fst) . zip u
+
+instance (Enum a, Ord a) => Perm [a] where
+    st         = fromVector . I.st . I.fromList . map fromEnum
+    act u      = act' $ toList (inverse u)
+    inverse v  = act' v (idperm v)
+    size       = length
+    idperm     = sort
+
+
+-- Generalize
+-- ----------
+
+-- | Generalize a function on 'StPerm' to a function on any permutations:
+-- 
+-- > generalize f v = f (st v) `act` idperm v
+-- 
+-- Note that this will only work as intended if @f@ is size preserving.
+generalize :: Perm a => (StPerm -> StPerm) -> a -> a
+generalize f v = f (st v) `act` idperm v
+
+
+-- Generating permutations
+-- -----------------------
+
+-- | @unrankPerm u rank@ is the @rank@-th (Myrvold & Ruskey)
+-- permutation of @u@. E.g.,
+-- 
+-- > unrankPerm ['1'..'9'] 88888 == "561297843"
+-- 
+unrankPerm :: Perm a => a -> Integer -> a
+unrankPerm u = (`act` u) . fromVector . I.unrankPerm (size u)
+
+-- | @randomPerm u@ is a random permutation of @u@.
+randomPerm :: Perm a => a -> IO a
+randomPerm u = ((`act` u) . fromVector . I.fromLehmercode) `liftM` I.randomLehmercode (size u)
+
+-- | All permutations of a given permutation. E.g.,
+-- 
+-- > perms "123" == ["123","213","321","132","231","312"]
+-- 
+perms :: Perm a => a -> [a]
+perms u = map (`act` u) $ sym (size u)
+
+
+-- Sorting operators
+-- -----------------
+
+-- | One pass of stack-sort.
+stackSort :: Perm a => a -> a
+stackSort = generalize (fromVector . I.stackSort . toVector)
+
+-- | One pass of bubble-sort.
+bubbleSort :: Perm a => a -> a
+bubbleSort = generalize (fromVector . I.bubbleSort . toVector)
+
+
+-- Permutation patterns
+-- --------------------
+
+-- | @copies p w@ is the list of (indices of) copies of the pattern
+-- @p@ in the permutation @w@. E.g.,
+-- 
+-- > copies (st "21") "2431" == [fromList [1,2],fromList [0,3],fromList [1,3],fromList [2,3]]
+-- 
+copies :: Perm a => StPerm -> a -> [Set]
+copies p w = I.copies subsets (toVector p) (toVector $ st w)
+
+-- | @avoids ps w@ is a predicate determining if @w@ avoids the patterns @ps@.
+avoids :: Perm a => [StPerm] -> a -> Bool
+avoids ps w = all null [ copies p w | p <- ps ]
+
+-- | @avoiders ps v@ is the list of permutations of @v@ avoiding the
+-- patterns @ps@. This is equivalent to the definition
+-- 
+-- > avoiders ps = filter (avoids ps)
+-- 
+-- but is usually much faster.
+avoiders :: Perm a => [StPerm] -> [a] -> [a]
+avoiders ps = I.avoiders subsets (toVector . st) (map toVector ps)
+
+-- | @av ps n@ is the list of permutations of @[0..n-1]@ avoiding the
+-- patterns @ps@. E.g.,
+-- 
+-- > map (length . av [st "132", st "321"]) [1..8] == [1,2,4,7,11,16,22,29]
+-- 
+av :: [StPerm] -> Int -> [StPerm]
+av ps = avoiders ps . sym
+
+
+-- Subsets
+-- -------
+
+-- | A set is represented by an increasing vector of non-negative
+-- integers.
+type Set = SV.Vector Int
+
+-- A sub-class of 'Bits' used internally. Minimal complete definiton: 'next'.
+class (Bits a, Integral a) => Bitmask a where
+    -- | Lexicographically, the next bitmask with the same Hamming weight.
+    next :: a -> a
+
+    -- | @ones k m@ is the set of indices whose bits are set in
+    -- @m@. Default implementation:
+    -- 
+    -- > ones m = fromListN (popCount m) $ filter (testBit m) [0..]
+    -- 
+    ones :: a -> Set
+    ones m = SV.fromListN (popCount m) $ filter (testBit m) [0..]
+
+instance Bitmask CUInt where
+    next = I.nextCUInt
+    ones = I.onesCUInt
+
+instance Bitmask Integer where
+    next = I.nextIntegral
+
+-- @bitmasks n k@ is the list of bitmasks with Hamming weight @k@ and
+-- size less than @2^n@.
+bitmasks :: Bitmask a => Int -> Int -> [a]
+bitmasks n k = take binomial (iterate next ((1 `shiftL` k) - 1))
+    where
+      n' = toInteger n
+      k' = toInteger k
+      binomial = fromIntegral $ product [n', n'-1 .. n'-k'+1] `div` product [1..k']
+
+-- | @subsets n k@ is the list of subsets of @[0..n-1]@ with @k@
+-- elements.
+subsets :: Int -> Int -> [Set]
+subsets n k = if n <= bitSize (0 :: CUInt)
+              then map ones (bitmasks n k :: [CUInt])
+              else map ones (bitmasks n k :: [Integer])
diff --git a/Math/Sym/D8.hs b/Math/Sym/D8.hs
new file mode 100644
--- /dev/null
+++ b/Math/Sym/D8.hs
@@ -0,0 +1,118 @@
+-- |
+-- Module      : Math.Sym.D8
+-- Copyright   : (c) Anders Claesson 2012
+-- License     : BSD-style
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+-- 
+-- The dihedral group of order 8 acting on permutations.
+-- 
+-- To avoid name clashes this module is best imported @qualified@;
+-- e.g.
+-- 
+-- > import qualified Math.Sym.D8 as D8
+-- 
+
+module Math.Sym.D8
+    (
+    -- * The group elements
+      r0            -- :: Perm a => a -> a
+    , r1            -- :: Perm a => a -> a
+    , r2            -- :: Perm a => a -> a
+    , r3            -- :: Perm a => a -> a
+    , s0            -- :: Perm a => a -> a
+    , s1            -- :: Perm a => a -> a
+    , s2            -- :: Perm a => a -> a
+    , s3            -- :: Perm a => a -> a
+
+    -- * D8, the klein four-group, and orbits
+    , d8            -- :: Perm a => [a -> a]
+    , klein4        -- :: Perm a => [a -> a]
+    , orbit         -- :: Ord a => Perm a => [a -> a] -> a -> [a]
+
+    -- * Aliases
+    , id            -- :: Perm a => a -> a
+    , rotate        -- :: Perm a => a -> a
+    , complement    -- :: Perm a => a -> a
+    , reverse       -- :: Perm a => a -> a
+    , inverse       -- :: Perm a => a -> a
+    ) where
+
+import Prelude hiding (reverse, id)
+import Data.List (group, sort)
+import Math.Sym (Perm (size), fromVector, act)
+import qualified Math.Sym (inverse)
+import Math.Sym.Internal (revIdperm)
+
+-- The group elements
+-- ------------------
+
+r0, r1, r2, r3, s0, s1, s2, s3 :: Perm a => a -> a
+
+-- | Ration by 0 degrees, i.e. the identity map.
+r0 w = w
+
+-- | Ration by 90 degrees clockwise.
+r1 = s2 . s1
+
+-- | Ration by 2*90 = 180 degrees clockwise.
+r2 = r1 . r1
+
+-- | Ration by 3*90 = 270 degrees clockwise.
+r3 = r2 . r1
+
+-- | Reflection through a horizontal axis (also called 'complement').
+s0 = r1 . s2
+
+-- | Reflection through a vertical axis (also called 'reverse').
+s1 w = (fromVector . revIdperm . size) w `act` w
+
+-- | Reflection through the main diagonal (also called 'inverse').
+s2 = Math.Sym.inverse
+
+-- | Reflection through the anti-diagonal.
+s3 = s1 . r1
+
+-- D8, the klein four-group, and orbits
+-- ------------------------------------
+
+d8, klein4 :: Perm a => [a -> a]
+
+-- | The dihedral group of order 8 (the symmetries of a square); that is,
+-- 
+-- > d8 = [r0, r1, r2, r3, s0, s1, s2, s3]
+-- 
+d8 = [r0, r1, r2, r3, s0, s1, s2, s3]
+
+-- | The Klein four-group (the symmetries of a non-equilateral
+-- rectangle); that is,
+-- 
+-- > klein4 = [r0, r2, s0, s1]
+-- 
+klein4 = [r0, r2, s0, s1]
+
+-- | @orbit fs x@ is the orbit of @x@ under the functions in @fs@. E.g.,
+-- 
+-- > orbit klein4 "2314" == ["1423","2314","3241","4132"]
+-- 
+orbit :: Ord a => Perm a => [a -> a] -> a -> [a]
+orbit fs x = map head . group $ sort [ f x | f <- fs ]
+
+-- Aliases
+-- -------
+
+id, rotate, complement, reverse, inverse :: Perm a => a -> a
+
+-- | @id = r0@
+id = r0
+
+-- | @rotate = r1@
+rotate = r1
+
+-- | @complement = s0@
+complement = s0
+
+-- | @reverse = s1@
+reverse = s1
+
+-- | @inverse = s2@
+inverse = s2
diff --git a/Math/Sym/Internal.hs b/Math/Sym/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Math/Sym/Internal.hs
@@ -0,0 +1,481 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Module      : Math.Sym.Internal
+-- Copyright   : (c) Anders Claesson 2012
+-- License     : BSD-style
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+-- 
+-- An internal module used by the sym package.
+-- 
+-- A Lehmercode is a vector of integers @w@ such @w!i <= length w - 1 - i@
+-- for each @i@ in @[0..length w - 1]@; such a vector encodes a permutation.
+-- This module implements /O(n)/ algorithms for unranking Lehmercodes and
+-- permutations; the algorithms are due to W. Myrvold and F. Ruskey
+-- [Ranking and Unranking Permutations in Linear Time, Information Processing
+-- Letters, 79 (2001) 281-284].
+-- 
+-- In addition, this module implements sorting operators, the
+-- symmetries in D8 acting on permutations, as well as most of the
+-- common permutation statistics.
+
+module Math.Sym.Internal
+    (
+      Lehmercode
+    , Perm0
+
+    -- * Lehmercodes
+    , unrankLehmercode
+    , fromLehmercode
+    , randomLehmercode
+    , lehmercodes
+
+    -- * Permutations
+    , size
+    , toList
+    , fromList
+    , act
+    , unrankPerm
+    , randomPerm
+    , sym
+    , idperm
+    , revIdperm
+    , sti
+    , st
+    , ordiso
+    , copies
+    , avoiders
+
+    -- * Permutation symmetries
+    , reverse
+    , complement
+    , inverse
+    , rotate
+
+    -- * Permutation statistics
+    , asc     -- ascents
+    , des     -- descents
+    , exc     -- excedances
+    , fp      -- fixed points
+    , inv     -- inversions
+    , maj     -- the major index
+    , peak    -- peaks
+    , vall    -- valleys
+    , dasc    -- double ascents
+    , ddes    -- double descents
+    , lmin    -- left-to-right minima
+    , lmax    -- left-to-right maxima
+    , rmin    -- right-to-left minima
+    , rmax    -- right-to-left maxima
+    , head    -- the first element
+    , last    -- the last element
+    , lir     -- left-most increasing run
+    , ldr     -- left-most decreasing run
+    , rir     -- right-most increasing run
+    , rdr     -- right-most decreasing run
+    , comp    -- components
+    , ep      -- rank a la Elizalde & Pak
+
+    -- * Sorting operators
+    , stackSort
+    , bubbleSort
+
+    -- * Bitmasks
+    , onesCUInt
+    , nextCUInt
+    , nextIntegral
+    ) where
+
+import Prelude hiding (reverse, head, last)
+import qualified Prelude (head)
+import System.Random (getStdRandom, randomR)
+import Control.Monad (forM_, liftM)
+import Control.Monad.ST (runST)
+import Data.List (group)
+import Data.Bits (Bits, shiftR, (.|.), (.&.), popCount)
+import qualified Data.Vector.Storable as SV
+    ( Vector, toList, fromList, length, (!), thaw
+    , unsafeFreeze, unsafeWith, enumFromN, enumFromStepN
+    , head, last, filter, maximum, minimum, null, reverse, map
+    )
+import qualified Data.Vector.Storable.Mutable as MV
+    ( unsafeNew, unsafeWrite, unsafeWith, swap, replicate
+    )
+import Foreign (Ptr, castPtr)
+import System.IO.Unsafe (unsafePerformIO)
+import Foreign.C.Types (CLong(..), CInt(..), CUInt(..))
+import Foreign.Marshal.Utils (toBool)
+
+-- | A Lehmercode is a vector of integers @w@ such @w!i <= length w - 1 - i@
+-- for each @i@ in @[0..length w - 1]@.
+type Lehmercode = SV.Vector Int
+
+-- | By convention, a member of @Perm0@ is a permutation of some
+-- finite subset of @[0..]@.
+type Perm0 = SV.Vector Int
+
+
+-- Lehmercodes
+-- -----------
+
+-- | @unrankLehmercode n rank@ is the @rank@-th Lehmercode of length @n@.
+unrankLehmercode :: Int -> Integer -> Lehmercode
+unrankLehmercode n rank = runST $ do
+  v <- MV.unsafeNew n
+  iter v n rank (toInteger n)
+  SV.unsafeFreeze v
+    where
+      {-# INLINE iter #-}
+      iter _ 0 _ _ = return ()
+      iter v i r m = do
+        let (r',j) = quotRem r m
+        MV.unsafeWrite v (n-i) (fromIntegral j)
+        iter v (i-1) r' (m-1)
+
+-- | Build a permutation from its Lehmercode.
+fromLehmercode :: Lehmercode -> Perm0
+fromLehmercode code = runST $ do
+  let n = SV.length code
+  v <- MV.unsafeNew n
+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v i i
+  forM_ [0..n-1] $ \i -> MV.swap v i (i + (SV.!) code i)
+  SV.unsafeFreeze v
+
+-- | A random Lehmercode of the given length.
+randomLehmercode :: Int -> IO Lehmercode
+randomLehmercode n = unrankLehmercode n `liftM` getStdRandom (randomR (0, factorial n - 1))
+
+-- | The list of Lehmercodes of a given length.
+lehmercodes :: Int -> [Lehmercode]
+lehmercodes n = map (unrankLehmercode n) [0 .. factorial n - 1]
+
+
+-- Permutations
+-- ------------
+
+-- | The size of a permutation; the number of elements.
+size :: Perm0 -> Int
+size = SV.length
+
+-- | The list of images of a permutation.
+toList :: Perm0 -> [Int]
+toList = SV.toList
+
+-- | Make a permutation from a list of images.
+fromList :: [Int] -> Perm0
+fromList = SV.fromList
+
+-- | @act u v@ is the permutation /w/ defined by /w(u(i)) = v(i)/.
+act :: Perm0 -> Perm0 -> Perm0
+act u v = runST $ do
+  let n = SV.length u
+  w <- MV.unsafeNew n
+  forM_ [0..n-1] $ \i -> MV.unsafeWrite w i ((SV.!) v ((SV.!) u i))
+  SV.unsafeFreeze w
+
+factorial :: Integral a => a -> Integer
+factorial = product . enumFromTo 1 . toInteger 
+
+-- | @unrankPerm n rank@ is the @rank@-th (Myrvold & Ruskey) permutation of length @n@.
+unrankPerm :: Int -> Integer -> Perm0
+unrankPerm n = fromLehmercode . unrankLehmercode n
+
+-- | A random permutation of the given length.
+randomPerm :: Int -> IO Perm0
+randomPerm n = fromLehmercode `liftM` randomLehmercode n
+
+-- | @sym n@ is the list of permutations of @[0..n-1]@ (the symmetric group).
+sym :: Int -> [Perm0]
+sym n = map (unrankPerm n) [0 .. factorial n - 1]
+
+-- | The identity permutation of the given length.
+idperm :: Int -> Perm0
+idperm = SV.enumFromN 0
+
+-- | The reverse of the identity permutation.
+revIdperm :: Int -> Perm0
+revIdperm n = SV.enumFromStepN (n-1) (-1) n
+
+-- | @sti w@ is the inverse of the standardization of @w@ (a
+-- permutation on @[0..length w-1]@). E.g., @sti \<4,9,2\> ==
+-- \<2,0,1\>@.
+sti :: Perm0 -> Perm0
+sti w = runST $ do
+  let a = if SV.null w then 0 else SV.minimum w
+  let b = if SV.null w then 0 else SV.maximum w
+  let n = SV.length w
+  v <- MV.replicate (1 + b - a) (-1)
+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v ((SV.!) w i - a) i
+  SV.filter (>=0) `liftM` SV.unsafeFreeze v
+
+-- | The standardization map.
+st :: Perm0 -> Perm0
+st = inverse . sti
+
+foreign import ccall unsafe "ordiso.h ordiso" c_ordiso
+    :: Ptr CLong -> Ptr CLong -> Ptr CLong -> CLong -> CInt
+
+-- | @ordiso u v m@ determines whether the subword in @v@ specified by
+-- @m@ is order isomorphic to @u@.
+ordiso :: Perm0 -> Perm0 -> SV.Vector Int -> Bool
+ordiso u v m =
+    let k = fromIntegral (SV.length u)
+    in  unsafePerformIO $
+        SV.unsafeWith u $ \u' ->
+        SV.unsafeWith v $ \v' ->
+        SV.unsafeWith m $ \m' ->
+        return . toBool $ c_ordiso (castPtr u') (castPtr v') (castPtr m') k
+
+-- | @copies subsets p w@ is the list of bitmasks that represent copies of @p@ in @w@.
+copies :: (Int -> Int -> [SV.Vector Int]) -> Perm0 -> Perm0 -> [SV.Vector Int]
+copies subsets p w = filter (ordiso p w) $ subsets n k
+    where
+      n = SV.length w
+      k = SV.length p
+
+avoiders1 :: (Int -> Int -> [SV.Vector Int]) -> (a -> Perm0) -> Perm0 -> [a] -> [a]
+avoiders1 subsets f p ws =
+    let ws0 = map f ws
+        ws2 = zip ws0 ws
+    in case group (map SV.length ws0) of
+         []  -> []
+         [_] -> let k = SV.length p
+                    n = SV.length (Prelude.head ws0)
+                in  [ v | (v0,v) <- ws2,  not $ any (ordiso p v0) (subsets n k) ]
+         _   ->     [ v | (v0,v) <- ws2, null $ copies subsets p v0 ] 
+
+-- | @avoiders subsets st ps ws@ is the list of permutations in @ws@
+-- avoiding the patterns in @ps@.
+avoiders :: (Int -> Int -> [SV.Vector Int]) -> (a -> Perm0) -> [Perm0] -> [a] -> [a]
+avoiders _       _   []   ws = ws
+avoiders subsets f (p:ps) ws = avoiders subsets f ps $ avoiders1 subsets f p ws
+
+
+-- Permutation symmetries
+-- ----------------------
+
+-- | @reverse \<a_1,...,a_n\> == \<a_n,,...,a_1\>@. E.g., @reverse
+-- \<9,3,7,2\> == \<2,7,3,9\>@.
+reverse :: Perm0 -> Perm0
+reverse = SV.reverse
+
+-- | @complement \<a_1,...,a_n\> == \<b_1,,...,b_n\>@, where @b_i = n - a_i - 1@.
+-- E.g., @complement \<3,4,0,1,2\> == \<1,0,4,3,2\>@.
+complement :: Perm0 -> Perm0
+complement w = SV.map (\x -> SV.length w - x - 1) w
+
+-- | @inverse w@ is the group theoretical inverse of @w@. E.g.,
+-- @inverse \<1,2,0\> == \<2,0,1\>@.
+inverse :: Perm0 -> Perm0
+inverse w = runST $ do
+  let n = SV.length w
+  v <- MV.unsafeNew n
+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v ((SV.!) w i) i
+  SV.unsafeFreeze v
+
+-- | The clockwise rotatation through 90 degrees. E.g.,
+-- @rotate \<1,0,2\> == \<1,2,0\>@.
+rotate :: Perm0 -> Perm0
+rotate w = runST $ do
+  let n = SV.length w
+  v <- MV.unsafeNew n
+  forM_ [0..n-1] $ \i -> MV.unsafeWrite v ((SV.!) w (n-1-i)) i
+  SV.unsafeFreeze v
+
+
+-- Permutation statistics
+-- ----------------------
+
+foreign import ccall unsafe "stat.h asc" c_asc
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h des" c_des
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h exc" c_exc
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h fp" c_fp
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h inv" c_inv
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h maj" c_maj
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h peak" c_peak
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h vall" c_vall
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h dasc" c_dasc
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h ddes" c_ddes
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h lmin" c_lmin
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h lmax" c_lmax
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h lir" c_lir
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h ldr" c_ldr
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h comp" c_comp
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h ep" c_ep
+    :: Ptr CLong -> CLong -> CLong
+
+-- Marshal a permutation statistic defined in C to on in Haskell.
+stat :: (Ptr CLong -> CLong -> CLong) -> Perm0 -> Int
+stat f w = unsafePerformIO $
+           SV.unsafeWith w $ \ptr ->
+               return . fromIntegral $ f (castPtr ptr) (fromIntegral (SV.length w))
+
+-- | First (left-most) value of a permutation.
+head :: Perm0 -> Int
+head = SV.head
+
+-- | Last (right-most) value of a permutation.
+last :: Perm0 -> Int
+last = SV.last
+
+-- | The number of right-to-left minima.
+rmin :: Perm0 -> Int
+rmin = lmin . SV.reverse
+
+-- | The number of right-to-left maxima.
+rmax :: Perm0 -> Int
+rmax = lmax . SV.reverse
+
+-- | The right-most increasing run.
+rir :: Perm0 -> Int
+rir = ldr . SV.reverse
+
+-- | The right-most decreasing run.
+rdr :: Perm0 -> Int
+rdr = lir . SV.reverse
+
+-- | The number of ascents.
+asc :: Perm0 -> Int
+asc = stat c_asc
+
+-- | The number of descents.
+des :: Perm0 -> Int
+des = stat c_des
+
+-- | The number of inversions.
+inv :: Perm0 -> Int
+inv = stat c_inv
+
+-- | The major index.
+maj :: Perm0 -> Int
+maj = stat c_maj
+
+-- | The number of peaks.
+peak :: Perm0 -> Int
+peak = stat c_peak
+
+-- | The number of valleys.
+vall :: Perm0 -> Int
+vall = stat c_vall
+
+-- | The number of double ascents.
+dasc :: Perm0 -> Int
+dasc = stat c_dasc
+
+-- | The number of double descents.
+ddes :: Perm0 -> Int
+ddes = stat c_ddes
+
+-- | The number of left-to-right minima.
+lmin :: Perm0 -> Int
+lmin = stat c_lmin
+
+-- | The number of left-to-right maxima.
+lmax :: Perm0 -> Int
+lmax = stat c_lmax
+
+-- | The left-most increasing run.
+lir :: Perm0 -> Int
+lir = stat c_lir
+
+-- | The left-most decreasing run.
+ldr :: Perm0 -> Int
+ldr = stat c_ldr
+
+-- | The number of excedances.
+exc :: Perm0 -> Int
+exc = stat c_exc
+
+-- | The number of fixed points.
+fp :: Perm0 -> Int
+fp = stat c_fp
+
+-- | The number of components.
+comp :: Perm0 -> Int
+comp = stat c_comp
+
+-- | Rank as defined by Elizalde & Pak.
+ep :: Perm0 -> Int
+ep = stat c_ep
+
+
+-- Sorting operators
+-- -----------------
+
+foreign import ccall unsafe "sortop.h stacksort" c_stacksort
+    :: Ptr CLong -> CLong -> IO ()
+
+foreign import ccall unsafe "sortop.h bubblesort" c_bubblesort
+    :: Ptr CLong -> CLong -> IO ()
+
+-- Marshal a sorting operator defined in C to on in Haskell.
+sortop :: (Ptr CLong -> CLong -> IO ()) -> Perm0 -> Perm0
+sortop f w = unsafePerformIO $ do
+               v <- SV.thaw w
+               MV.unsafeWith v $ \ptr -> do
+                 f (castPtr ptr) (fromIntegral (SV.length w))
+                 SV.unsafeFreeze v
+
+-- | One pass of stack-sort.
+stackSort :: Perm0 -> Perm0
+stackSort = sortop c_stacksort
+
+-- | One pass of bubble-sort.
+bubbleSort :: Perm0 -> Perm0
+bubbleSort = sortop c_bubblesort
+
+
+-- Bitmasks
+-- --------
+
+foreign import ccall unsafe "bit.h next" c_next :: CUInt -> CUInt
+
+-- | Lexicographically, the next 'CUInt' with the same Hamming weight.
+nextCUInt :: CUInt -> CUInt
+nextCUInt = c_next
+
+foreign import ccall unsafe "bit.h ones" c_ones :: Ptr CUInt -> CUInt -> IO ()
+
+-- | @onesCUInt k m@ gives the @k@ smallest indices whose bits are set in @m@.
+onesCUInt :: CUInt -> SV.Vector Int
+onesCUInt m = SV.map fromIntegral . unsafePerformIO $ do
+                v <- MV.unsafeNew (popCount m)
+                MV.unsafeWith v $ \ptr -> do
+                  c_ones ptr m
+                  SV.unsafeFreeze v
+
+-- | Lexicographically, the next integral number with the same Hamming weight.
+nextIntegral :: (Integral a, Bits a) => a -> a
+nextIntegral a =
+    let b = (a .|. (a - 1)) + 1
+    in  b .|. ((((b .&. (-b)) `div` (a .&. (-a))) `shiftR` 1) - 1)
diff --git a/Math/Sym/Stat.hs b/Math/Sym/Stat.hs
new file mode 100644
--- /dev/null
+++ b/Math/Sym/Stat.hs
@@ -0,0 +1,154 @@
+-- |
+-- Module      : Math.Sym.Stat
+-- Copyright   : (c) Anders Claesson 2012
+-- License     : BSD-style
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+-- 
+-- Common permutation statistics. Please contact the maintainer if you
+-- feel that there is a statistic that is missing; even better, send a
+-- patch or make a pull request.
+-- 
+-- To avoid name clashes this module is best imported @qualified@;
+-- e.g.
+-- 
+-- > import qualified Math.Sym.Stat as S
+-- 
+-- For any permutation statistic @f@, below, it holds that @f w == f
+-- (st w)@, and therefore the explanations below will be done on
+-- standard permutations for convenience.
+
+module Math.Sym.Stat 
+    (
+      asc     -- ascents
+    , des     -- descents
+    , exc     -- excedances
+    , fp      -- fixed points
+    , inv     -- inversions
+    , maj     -- the major index
+    , peak    -- peaks
+    , vall    -- valleys
+    , dasc    -- double ascents
+    , ddes    -- double descents
+    , lmin    -- left-to-right minima
+    , lmax    -- left-to-right maxima
+    , rmin    -- right-to-left minima
+    , rmax    -- right-to-left maxima
+    , head    -- the first element
+    , last    -- the last element
+    , lir     -- left-most increasing run
+    , ldr     -- left-most decreasing run
+    , rir     -- right-most increasing run
+    , rdr     -- right-most decreasing run
+    , comp    -- components
+    , ep      -- rank a la Elizalde & Pak
+    ) where
+
+import Prelude hiding (head, last)
+import Math.Sym (Perm, toVector, st)
+import Math.Sym.Internal (Perm0)
+import qualified Math.Sym.Internal as I 
+    ( asc, des, exc, fp, inv, maj, peak, vall, dasc, ddes, lmin, lmax, rmin, rmax
+    , head, last, lir, ldr, rir, rdr, comp, ep
+    )
+
+generalize :: Perm a => (Perm0 -> Int) -> a -> Int
+generalize f = f . toVector . st
+
+-- | The number of ascents. An /ascent/ in @w@ is an index @i@ such
+-- that @w[i] \< w[i+1]@.
+asc :: Perm a => a -> Int
+asc = generalize I.asc
+
+-- | The number of descents. A /descent/ in @w@ is an index @i@ such
+-- that @w[i] > w[i+1]@.
+des :: Perm a => a -> Int
+des = generalize I.des
+
+-- | The number of /excedances/: positions @i@ such that @w[i] > i@;
+exc :: Perm a => a -> Int
+exc = generalize I.exc
+
+-- | The number of /fixed points/: positions @i@ such that @w[i] == i@;
+fp :: Perm a => a -> Int
+fp = generalize I.fp
+
+-- | The number of /inversions/: pairs @\(i,j\)@ such that @i \< j@ and @w[i] > w[j]@
+inv :: Perm a => a -> Int
+inv = generalize I.inv
+
+-- | /The major index/ is the sum of descents.
+maj :: Perm a => a -> Int
+maj = generalize I.maj
+
+-- | The number of /peaks/: positions @i@ such that @w[i-1] \< w[i]@ and @w[i] \> w[i+1]@.
+peak :: Perm a => a -> Int
+peak = generalize I.peak
+
+-- | The number of /valleys/: positions @i@ such that @w[i-1] \> w[i]@ and @w[i] \< w[i+1]@.
+vall :: Perm a => a -> Int
+vall = generalize I.vall
+
+-- | The number of /double ascents/: positions @i@ such that @w[i-1] \<  w[i] \< w[i+1]@.
+dasc :: Perm a => a -> Int
+dasc = generalize I.dasc
+
+-- | The number of /double descents/: positions @i@ such that @w[i-1] \>  w[i] \> w[i+1]@.
+ddes :: Perm a => a -> Int
+ddes = generalize I.ddes
+
+-- | The number of /left-to-right minima/: positions @i@ such that @w[i] \< w[j]@ for all @j \< i@.
+lmin :: Perm a => a -> Int
+lmin = generalize I.lmin
+
+-- | The number of /left-to-right maxima/: positions @i@ such that @w[i] \> w[j]@ for all @j \< i@.
+lmax :: Perm a => a -> Int
+lmax = generalize I.lmax
+
+-- | The number of /right-to-left minima/: positions @i@ such that @w[i] \< w[j]@ for all @j \> i@.
+rmin :: Perm a => a -> Int
+rmin = generalize I.rmin
+
+-- | The number of /right-to-left maxima/: positions @i@ such that @w[i] \> w[j]@ for all @j \> i@.
+rmax :: Perm a => a -> Int
+rmax = generalize I.rmax
+
+-- | The first (left-most) element in the standardization. E.g., @head \"231\" = head (fromList [1,2,0]) = 1@.
+head :: Perm a => a -> Int
+head = generalize I.head
+
+-- | The last (right-most) element in the standardization. E.g., @last \"231\" = last (fromList [1,2,0]) = 0@.
+last :: Perm a => a -> Int
+last = generalize I.last
+
+-- | Length of the left-most increasing run: largest @i@ such that
+-- @w[0] \< w[1] \< ... \< w[i-1]@.
+lir :: Perm a => a -> Int
+lir = generalize I.lir
+
+-- | Length of the left-most decreasing run: largest @i@ such that
+-- @w[0] \> w[1] \> ... \> w[i-1]@.
+ldr :: Perm a => a -> Int
+ldr = generalize I.ldr
+
+-- | Length of the right-most increasing run: largest @i@ such that
+-- @w[n-i] \< ... \< w[n-2] \< w[n-1]@.
+rir :: Perm a => a -> Int
+rir = generalize I.rir
+
+-- | Length of the right-most decreasing run: largest @i@ such that
+-- @w[n-i] \> ... \> w[n-2] \> w[n-1]@.
+rdr :: Perm a => a -> Int
+rdr = generalize I.rdr
+
+-- | The number of components. E.g., @[2,0,3,1,4,6,7,5]@ has three
+-- components: @[2,0,3,1]@, @[4]@ and @[6,7,5]@.
+comp :: Perm a => a -> Int
+comp = generalize I.comp
+
+-- | The rank as defined by Elizalde and Pak [Bijections for
+-- refined restricted permutations, /J. Comb. Theory, Ser. A/, 2004]:
+-- 
+-- > maximum [ k | k <- [0..n-1], w[i] >= k for all i < k ]
+-- 
+ep :: Perm a => a -> Int
+ep = generalize I.ep
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/cbits/bit.c b/cbits/bit.c
new file mode 100644
--- /dev/null
+++ b/cbits/bit.c
@@ -0,0 +1,19 @@
+ #include <strings.h>
+
+/* Lexicographically, the next bitmask with the same Hamming weight */
+unsigned int
+next(const unsigned int v)
+{
+	unsigned int t = v | (v - 1);
+	return ((t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1)));
+}
+
+/* Positions of bits set */
+void
+ones(unsigned int *u, const unsigned int a)
+{
+	unsigned int b;
+
+	for (b = a; b; b &= b-1)
+		*u++ = ffs(b) - 1;
+}
diff --git a/cbits/ordiso.c b/cbits/ordiso.c
new file mode 100644
--- /dev/null
+++ b/cbits/ordiso.c
@@ -0,0 +1,27 @@
+#include <stdlib.h>
+
+/*
+ * Determines whether the subword in v specified by m is order
+ * isomorphic to u; len is the length of u.
+ */
+int
+ordiso(const long *u, const long *v, const long *m, long len)
+{
+	register int i;
+	long *w = malloc(len*sizeof(*w));
+	long *w0 = w;
+
+        /* Let w = v.m.u^{-1} */
+	for (i = 0; i < len; i++, u++, m++)
+		w[(int)*u] = v[(int)*m];
+
+        /* Return 1 if w is increasing, 0 otherwise */
+        for (; len > 1; len--, w++) {
+                if (*w > *(w+1)) {
+			free(w0);
+			return 0;
+		}
+        }
+	free(w0);
+	return 1;
+}
diff --git a/cbits/sortop.c b/cbits/sortop.c
new file mode 100644
--- /dev/null
+++ b/cbits/sortop.c
@@ -0,0 +1,35 @@
+
+/* One pass of stack-sort implemented a la Petter Br\"and\'en [Actions
+ * on permutations and unimodality of descent polynomials, European
+ * J. Combin. 29 (2008)]
+ */
+void
+stacksort(long *w, long len) {
+        int i = 0;
+        int j = 0;
+        int y;
+        while (i < len) {
+                j = i;
+                y = w[j];
+                while (y > w[j+1] && j+1 < len) {
+                        w[j] = w[j+1];
+                        j++;
+                }
+                w[j] = y;
+                if (j == i)
+			i++;
+        }
+}
+
+/* On pass of bubble-sort */
+void
+bubblesort(long *w, long len) {
+        int tmp;
+	for (; len > 1; len--, w++) {
+		if (*w > *(w+1)) {
+			tmp    = *w;
+			*w     = *(w+1);
+			*(w+1) = tmp;
+		}
+	}
+}
diff --git a/cbits/stat.c b/cbits/stat.c
new file mode 100644
--- /dev/null
+++ b/cbits/stat.c
@@ -0,0 +1,244 @@
+#include <string.h>
+
+/* The number of ascents */
+long
+asc(const long *w, long len)
+{
+	long acc = 0;
+
+	for (; len > 1; len--, w++) {
+		if (*w < *(w+1))
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of descents */
+long
+des(const long *w, long len)
+{
+	long acc = 0;
+
+	for (; len > 1; len--, w++) {
+		if (*w > *(w+1))
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of excedances */
+long
+exc(const long *w, long len)
+{
+	long i, acc = 0;
+
+	for (i = 0; i < len; i++, w++) {
+		if (*w > i)
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of fixed points */
+long
+fp(const long *w, long len)
+{
+	long i, acc = 0;
+
+	for (i = 0; i < len; i++, w++) {
+		if (*w == i)
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of inversions */
+long
+inv(const long *w, long len)
+{
+	long i, j;
+	long acc = 0;
+	long *v;
+
+	for (i = 0; i < len; i++, w++) {
+		for (j = i+1, v = (long*)w+1; j < len; j++, v++) {
+			if (*w > *v)
+				acc++;
+		}
+	}
+	return acc;
+}
+
+
+/* The major index */
+long
+maj(const long *w, long len)
+{
+	long i, sum = 0;
+
+	for (i = 1; i < len; i++, w++) {
+		if (*w > *(w+1))
+			sum += i;
+	}
+	return sum;
+}
+
+
+/* The number of peaks */
+long
+peak(const long *w, long len)
+{
+	long acc = 0;
+
+	for (; len > 2; len--, w++) {
+		if (*w < *(w+1) && *(w+1) > *(w+2))
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of valleys */
+long
+vall(const long *w, long len)
+{
+	long acc = 0;
+
+	for (; len > 2; len--, w++) {
+		if (*w > *(w+1) && *(w+1) < *(w+2))
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of double ascents */
+long
+dasc(const long *w, long len)
+{
+	long acc = 0;
+
+	for (; len > 2; len--, w++) {
+		if (*w < *(w+1) && *(w+1) < *(w+2))
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of double descents */
+long
+ddes(const long *w, long len)
+{
+	long acc = 0;
+
+	for (; len > 2; len--, w++) {
+		if (*w > *(w+1) && *(w+1) > *(w+2))
+			acc++;
+	}
+	return acc;
+}
+
+
+/* The number of left-to-right minima */
+long
+lmin(const long *w, long len)
+{
+	long m = *w + 1;
+	long acc = 0;
+
+	for (; len > 0; len--, w++) {
+		if (*w < m) {
+			m = *w;
+			acc++;
+		}
+	}
+	return acc;
+}
+
+
+/* The number of left-to-right maxima */
+long
+lmax(const long *w, long len)
+{
+	long m = *w - 1;
+	long acc = 0;
+
+	for (; len > 0; len--, w++) {
+		if (*w > m) {
+			m = *w;
+			acc++;
+		}
+	}
+	return acc;
+}
+
+/* The left-most increasing run */
+long
+lir(const long *w, long len)
+{
+	long acc;
+
+	if (len == 0)
+		return 0;
+
+	for (acc = 1; len > 1 && *w < *(w+1); len--, w++)
+		acc++;
+
+	return acc;
+}
+
+/* The left-most decreasing run */
+long
+ldr(const long *w, long len)
+{
+	long acc;
+
+	if (len == 0)
+		return 0;
+
+	for (acc = 1; len > 1 && *w > *(w+1); len--, w++)
+		acc++;
+
+	return acc;
+}
+
+/* The number of components; O(len) */
+long
+comp(const long *w, long len)
+{
+	long i;
+	long m = *w - 1;
+	long acc = 0;
+
+	for (i = 0; i < len; i++, w++) {
+		if (*w > m)
+			m = *w;
+		if (m == i)
+			acc++;
+	}
+	return acc;
+}
+
+/* rank as defined by Elizalde & Pak */
+long
+ep(const long *w, long len)
+{
+	long i;
+	long m = *w;
+
+	if (len == 0)
+		return 0;
+
+	for (i = 0; i < len; i++, w++) {
+		if (*w <= m)
+			m = *w;
+		if (m <= i)
+			return i;
+	}
+	return len;
+}
diff --git a/include/bit.h b/include/bit.h
new file mode 100644
--- /dev/null
+++ b/include/bit.h
@@ -0,0 +1,2 @@
+unsigned int next(unsigned int);
+void ones(unsigned int *, const unsigned int);
diff --git a/include/ordiso.h b/include/ordiso.h
new file mode 100644
--- /dev/null
+++ b/include/ordiso.h
@@ -0,0 +1,1 @@
+int ordiso(const long *, const long *, const long *, long);
diff --git a/include/sortop.h b/include/sortop.h
new file mode 100644
--- /dev/null
+++ b/include/sortop.h
@@ -0,0 +1,1 @@
+void stacksort(long *, long);
diff --git a/include/stat.h b/include/stat.h
new file mode 100644
--- /dev/null
+++ b/include/stat.h
@@ -0,0 +1,16 @@
+long asc  (const long *, long); /* ascents */
+long des  (const long *, long); /* descents */
+long exc  (const long *, long); /* excedances */
+long fp   (const long *, long); /* fixed points */
+long inv  (const long *, long); /* inversions */
+long maj  (const long *, long); /* major index */
+long peak (const long *, long); /* peaks */
+long vall (const long *, long); /* valleys */
+long dasc (const long *, long); /* double ascents */
+long ddes (const long *, long); /* double descents */
+long lmin (const long *, long); /* left-to-right minima */
+long lmax (const long *, long); /* left-to-right maxima */
+long lir  (const long *, long); /* left-most increasing run */
+long ldr  (const long *, long); /* left-most decreasing run */
+long comp (const long *, long); /* components */
+long ep   (const long *, long); /* rank a la Elizalde & Pak */
diff --git a/sym.cabal b/sym.cabal
new file mode 100644
--- /dev/null
+++ b/sym.cabal
@@ -0,0 +1,50 @@
+Name:                sym
+Version:             0.1
+Synopsis:            Permutations, patterns, and statistics
+Description:         
+  Definitions for permutations with an emphasis on permutation
+  patterns and statistics.
+  .
+  ["Math.Sym"] Provides an efficient definition of standard
+  permutations (@StPerm@) together with a typeclass (@Perm@) whose
+  functionality is largely inherited from @StPerm@ using a group
+  action and the standardization map.
+  .
+  ["Math.Sym.D8"] The dihedral group of order 8 acting on permutations.
+  .
+  ["Math.Sym.Stat"] Common permutation statistics, such as @des@,
+  @inv@, @exc@, @maj@, @fp@, @comp@, @lmin@, @lmax@, ...
+
+Homepage:            http://github.com/akc/sym
+
+License:             BSD3
+License-file:        LICENSE
+Author:              Anders Claesson
+Maintainer:          anders.claesson@gmail.com
+Category:            Math
+Build-type:          Simple
+
+Extra-source-files:  tests/Properties.hs
+
+Cabal-version:       >=1.6
+
+source-repository head
+  type:                git
+  location:            git://github.com/akc/sym.git
+
+Library
+  Exposed-modules:     Math.Sym
+                       Math.Sym.D8
+                       Math.Sym.Stat
+                       Math.Sym.Internal
+
+  Build-depends:       base >= 3 && < 5, random, vector
+  
+  ghc-prof-options:    -auto-all -caf-all
+  ghc-options:         -Wall -O2
+  cc-options:          -Wall
+
+  c-sources:           cbits/stat.c, cbits/sortop.c, cbits/ordiso.c, cbits/bit.c
+  include-dirs:        include
+  includes:            stat.h, sortop.h, ordiso.h, bit.h
+  install-includes:    stat.h, sortop.h, ordiso.h, bit.h
diff --git a/tests/Properties.hs b/tests/Properties.hs
new file mode 100644
--- /dev/null
+++ b/tests/Properties.hs
@@ -0,0 +1,508 @@
+-- |
+-- Copyright   : (c) Anders Claesson 2012
+-- License     : BSD-style
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+
+import Data.List
+import Data.Monoid
+import Control.Monad
+import qualified Math.Sym as Sym
+import qualified Math.Sym.D8 as D8
+import qualified Math.Sym.Stat as S
+import qualified Math.Sym.Internal as I
+import qualified Data.Vector.Storable as SV
+import Test.QuickCheck
+
+check :: Testable prop => prop -> IO ()
+check = quickCheck
+
+---------------------------------------------------------------------------------
+-- Generators
+---------------------------------------------------------------------------------
+
+rank :: Int -> Gen Integer
+rank n = choose (0, product [1..fromIntegral n] - 1)
+
+lenRank :: Gen (Int, Integer)
+lenRank = sized $ \m -> do
+            n <- choose (0, m)
+            r <- rank n
+            return (n, r)
+
+lenRank2 :: Gen (Int, Integer, Integer)
+lenRank2 = do (n, r1) <- lenRank
+              r2 <- rank n
+              return (n, r1, r2)
+
+moreThan :: Int -> Gen Int
+moreThan x = (\d -> x + abs d) `liftM` choose (1, 100)
+
+vecFrom :: Int -> Int -> Gen [Int]
+vecFrom 0 _ = return []
+vecFrom n x = moreThan x >>= liftM (x:) . vecFrom (n-1)
+
+incVec :: Int -> Gen [Int]
+incVec n = arbitrary >>= vecFrom n
+
+-- The sub-permutation determined by a set of indices.
+subperm :: Sym.Set -> Sym.StPerm -> Sym.StPerm
+subperm m w = Sym.fromVector . I.st $ SV.map ((SV.!) (Sym.toVector w)) m
+
+subperms :: Int -> Sym.StPerm -> [Sym.StPerm]
+subperms k w = [ subperm m w | m <- Sym.subsets (Sym.size w) k ]
+
+instance Arbitrary Sym.StPerm where
+    arbitrary = uncurry Sym.unrankStPerm `liftM` lenRank
+    shrink w = nub $ [0 .. Sym.size w - 1] >>= \k -> subperms k w
+
+perm2 :: Gen (Sym.StPerm, [Int])
+perm2 = do u <- arbitrary
+           v <- incVec (Sym.size u)
+           return (u, v)
+
+perm3 :: Gen (Sym.StPerm, Sym.StPerm, [Int])
+perm3 = do (n,r1,r2) <- lenRank2
+           let u = Sym.unrankStPerm n r1
+           let v = Sym.unrankStPerm n r2
+           w <- incVec n
+           return (u, v, w)
+
+perm :: Gen [Int]
+perm = liftM (uncurry Sym.act) perm2
+
+newtype Symmetry = Symmetry (Sym.StPerm -> Sym.StPerm, String)
+
+d8Symmetries :: [Symmetry]
+d8Symmetries = [ Symmetry (D8.r0, "r0")
+               , Symmetry (D8.r1, "r1")
+               , Symmetry (D8.r2, "r2")
+               , Symmetry (D8.r3, "r3")
+               , Symmetry (D8.s0, "s0")
+               , Symmetry (D8.s1, "s1")
+               , Symmetry (D8.s2, "s2")
+               , Symmetry (D8.s3, "s3")
+               ]
+
+instance Show Symmetry where
+    show (Symmetry (_,s)) = s
+
+instance Arbitrary Symmetry where
+    arbitrary = liftM (d8Symmetries !!) $ choose (0, length d8Symmetries - 1)
+
+
+---------------------------------------------------------------------------------
+-- Properties for Math.Sym
+---------------------------------------------------------------------------------
+
+prop_monoid_mempty1 w = mempty <> w == (w :: Sym.StPerm)
+prop_monoid_mempty2 w = w <> mempty == (w :: Sym.StPerm)
+prop_monoid_associative u v w = u <> (v <> w) == (u <> v) <> (w :: Sym.StPerm)
+
+newtype S = S {unS :: Sym.StPerm} deriving (Eq, Show)
+
+instance Arbitrary S where
+    arbitrary = liftM S arbitrary
+
+prop_monoid_mempty1_S w = mempty <> w == (w :: S)
+prop_monoid_mempty2_S w = w <> mempty == (w :: S)
+prop_monoid_associative_S u v w = u <> (v <> w) == (u <> v) <> (w :: S)
+
+instance Monoid S where
+    mempty = S $ Sym.fromVector SV.empty
+    mappend u v = S $ (Sym./-/) (unS u) (unS v)
+
+prop_unrankStPerm_distinct =
+    forAll lenRank $ \(n, r) ->
+        let w = Sym.toList (Sym.unrankStPerm n r) in nub w == w
+
+prop_unrankStPerm_injective =
+    forAll lenRank2 $ \(n, r1, r2) ->
+        (Sym.unrankStPerm n r1 :: Sym.StPerm) /= Sym.unrankStPerm n r2 || r1 == r2
+
+prop_sym = and [ sort (Sym.sym n) == sort (sym' n) | n<-[0..6] ]
+    where
+      sym' n = map Sym.fromList $ Data.List.permutations [0..fromIntegral n - 1]
+
+prop_perm =
+    and [ sort (Sym.perms [1..n]) == sort (permutations [1..n]) | n<-[0..6] ]
+
+prop_st =
+    forAll perm2 $ \(u,v) -> Sym.st (u `Sym.act` v) == u `Sym.act` Sym.st v
+
+prop_act_def =
+    forAll perm2 $ \(u,v) -> u `Sym.act` v == map (v!!) (Sym.toList u)
+
+prop_act_id =
+    forAll perm2 $ \(u,v) -> Sym.idperm u `Sym.act` v == v
+
+prop_act_associative =
+    forAll perm3 $ \(u,v,w) -> (u `Sym.act` v) `Sym.act` w == u `Sym.act` (v `Sym.act` w)
+
+prop_size =
+    forAll perm $ \v -> Sym.size v == Sym.size (Sym.st v)
+
+prop_idperm =
+    forAll perm2 $ \(u,v) -> Sym.idperm u == Sym.inverse (Sym.st u) `Sym.act` u
+
+prop_inverse =
+    forAll perm $ \v -> Sym.inverse v == Sym.inverse (Sym.st v) `Sym.act` Sym.idperm v
+
+prop_ordiso1 =
+    forAll perm2 $ \(u,v) -> u `Sym.ordiso` v  ==  (u == Sym.st v)
+
+prop_ordiso2 =
+    forAll perm2 $ \(u,v) -> u `Sym.ordiso` v  ==  (Sym.inverse u `Sym.act` v == Sym.idperm v)
+
+prop_unrankPerm =
+    forAll perm $ \w ->
+    forAll (choose (0, product [1..fromIntegral (length w) - 1])) $ \r ->
+        Sym.st (Sym.unrankPerm (sort w) r) == Sym.unrankStPerm (length w) r
+
+prop_stackSort = forAll perm $ \v -> Sym.stackSort v == stack v
+
+prop_stackSort_231 =
+    forAll perm $ \v -> (Sym.stackSort v == Sym.idperm v) == (Sym.avoids [Sym.st "231"] v)
+
+prop_bubbleSort = forAll perm $ \v -> Sym.bubbleSort v == bubble v
+
+prop_bubbleSort_231_321 =
+    forAll perm $ \v -> (Sym.bubbleSort v == Sym.idperm v) == (Sym.avoids [Sym.st "231", Sym.st "321"] v)
+
+prop_subperm_copies p =
+    forAll (resize 21 perm) $ \w -> and [ subperm m (Sym.st w) == p | m <- Sym.copies p w ]
+
+prop_copies =
+    forAll (resize  6 arbitrary) $ \p ->
+    forAll (resize 12 perm)      $ \w ->
+        sort (Sym.copies p w) == sort (map I.fromList $ copies (Sym.toList p) w)
+
+prop_copies_self =
+    forAll perm $ \v -> Sym.copies (Sym.st v) v == [SV.fromList [0 .. length v - 1]]
+
+prop_copies_d8 (Symmetry (f,_)) =
+    forAll (resize  6 arbitrary) $ \p ->
+    forAll (resize 20 perm)      $ \w ->
+        let p' = f p
+            w' = Sym.generalize f w
+        in length (Sym.copies p w) == length (Sym.copies p' w')
+
+prop_avoiders_avoid =
+    forAll (resize 20 arbitrary) $ \ws ->
+    forAll (resize  6 arbitrary) $ \ps ->
+        all (Sym.avoids ps) $ Sym.avoiders ps (ws :: [Sym.StPerm])
+
+prop_avoiders_idempotent =
+    forAll (resize 18 arbitrary) $ \vs ->
+    forAll (resize  5 arbitrary) $ \ps ->
+        let ws = Sym.avoiders ps (vs :: [Sym.StPerm]) in ws == Sym.avoiders ps ws
+
+prop_avoiders_d8 (Symmetry (f,_)) =
+    forAll (choose (0, 5))      $ \n ->
+    forAll (resize 5 arbitrary) $ \p ->
+        let ws = Sym.sym n in sort (map f $ Sym.avoiders [p] ws) == sort (Sym.avoiders [f p] ws)
+
+prop_avoiders_d8' (Symmetry (f,_)) =
+    forAll (choose (0, 5))      $ \n ->
+    forAll (resize 5 arbitrary) $ \ps ->
+        let ws = Sym.sym n in sort (map f $ Sym.avoiders ps ws) == sort (Sym.avoiders (map f ps) (map f ws))
+
+prop_avoiders_d8'' (Symmetry (f,_)) =
+    forAll (resize 18 arbitrary) $ \ws ->
+    forAll (resize  5 arbitrary) $ \ps ->
+        sort (map f $ Sym.avoiders ps ws) == sort (Sym.avoiders (map f ps) (map f ws :: [Sym.StPerm]))
+
+prop_av_cardinality =
+    forAll (resize 3 arbitrary) $ \p ->
+        let spec = [ length $ Sym.av [p] n | n<-[0..6] ]
+        in case Sym.size p of
+             0 -> spec == [0,0,0,0,0,0,0]
+             1 -> spec == [1,0,0,0,0,0,0]
+             2 -> spec == [1,1,1,1,1,1,1]
+             3 -> spec == [1,1,2,5,14,42,132]
+             _ -> True
+
+binomial n k = fromIntegral $ product [n', n'-1 .. n'-k'+1] `div` product [1..k']
+    where
+      n' = toInteger n
+      k' = toInteger k
+
+kSubsequences :: Int -> [a] -> [[a]]
+kSubsequences 0 _      = [[]]
+kSubsequences _ []     = []
+kSubsequences k (x:xs) = map (x:) (kSubsequences (k-1) xs) ++ kSubsequences k xs
+
+copies :: [Int] -> [Int] -> [[Int]]
+copies p w = [ is | js <- u, let (is, q) = unzip (f js (zip [0..] w)), st q == p ]
+    where
+      k = length p
+      n = length w
+      u = kSubsequences k [0..n-1]
+      f s@(j:t) ((i,x):v) = if i == j then (i,x) : f t v else f s v
+      f _       _         = []
+
+prop_subsets1 =
+    forAll (choose (0,14)) $ \n ->
+    forAll (choose (0,14)) $ \k ->
+        sort (kSubsequences k [0..n-1]) == sort (map SV.toList $ Sym.subsets n k)
+
+prop_subsets2 =
+    forAll (choose (0,35)) $ \n ->
+    forAll (choose (0,3))  $ \k ->
+        sort (kSubsequences k [0..n-1]) == sort (map SV.toList $ Sym.subsets n k)
+
+prop_subsets_singleton =
+    forAll (choose (0,500)) $ \n ->
+        let [v] = Sym.subsets n n in SV.toList v == [0..n-1]
+
+prop_subsets_cardinality1 =
+    forAll (choose (0,20)) $ \n ->
+    forAll (choose (0,20)) $ \k ->
+        length (Sym.subsets n k) == binomial n k
+
+prop_subsets_cardinality2 =
+    forAll (choose (0,20)) $ \n ->
+    forAll (choose (0,20)) $ \k ->
+        let cs = map (SV.length) (Sym.subsets n k) in ((k > n) && null cs) || ([k] == nub cs)
+
+testsPerm =
+    [ ("monoid/mempty/1",                check prop_monoid_mempty1)
+    , ("monoid/mempty/2",                check prop_monoid_mempty2)
+    , ("monoid/mempty/associative",      check prop_monoid_associative)
+    , ("monoid/mempty/1/skew",           check prop_monoid_mempty1_S)
+    , ("monoid/mempty/2/skew",           check prop_monoid_mempty2_S)
+    , ("monoid/mempty/associative/skew", check prop_monoid_associative_S)
+    , ("unrankStPerm/distinct",          check prop_unrankStPerm_distinct)
+    , ("unrankStPerm/injective",         check prop_unrankStPerm_injective)
+    , ("sym",                            check prop_sym)
+    , ("perm",                           check prop_perm)
+    , ("st",                             check prop_st)
+    , ("act/def",                        check prop_act_def)
+    , ("act/id",                         check prop_act_id)
+    , ("act/associative",                check prop_act_associative)
+    , ("size",                           check prop_size)
+    , ("idperm",                         check prop_idperm)
+    , ("inverse",                        check prop_inverse)
+    , ("ordiso/1",                       check prop_ordiso1)
+    , ("ordiso/2",                       check prop_ordiso2)
+    , ("unrankPerm",                     check prop_unrankPerm)
+    , ("stackSort",                      check prop_stackSort)
+    , ("stackSort/231",                  check prop_stackSort_231)
+    , ("bubbleSort",                     check prop_bubbleSort)
+    , ("bubbleSort/231&321",             check prop_bubbleSort_231_321)
+    , ("subperm/copies",                 check prop_subperm_copies)
+    , ("copies",                         check prop_copies)
+    , ("copies/self",                    check prop_copies_self)
+    , ("copies/D8",                      check prop_copies_d8)
+    , ("avoiders/avoid",                 check prop_avoiders_avoid)
+    , ("avoiders/idempotent",            check prop_avoiders_idempotent)
+    , ("avoiders/D8/0",                  check prop_avoiders_d8)
+    , ("avoiders/D8/1",                  check prop_avoiders_d8')
+    , ("avoiders/D8/2",                  check prop_avoiders_d8'')
+    , ("av/cardinality",                 check prop_av_cardinality)
+    , ("subsets/1",                      check prop_subsets1)
+    , ("subsets/2",                      check prop_subsets2)
+    , ("subsets/singleton",              check prop_subsets_singleton)
+    , ("subsets/cardinality/1",          check prop_subsets_cardinality1)
+    , ("subsets/cardinality/2",          check prop_subsets_cardinality2)
+    ]
+
+---------------------------------------------------------------------------------
+-- Properties for Math.Sym.D8
+---------------------------------------------------------------------------------
+
+prop_D8_orbit fs w = all (`elem` orbD8) $ D8.orbit (map fn fs) w
+    where
+      orbD8 = D8.orbit D8.d8 w
+      fn (Symmetry (f,_)) = f
+
+prop_D8_reverse w    = I.reverse    (Sym.toVector w) == Sym.toVector (D8.reverse w)
+prop_D8_complement w = I.complement (Sym.toVector w) == Sym.toVector (D8.complement w)
+prop_D8_inverse w    = I.inverse    (Sym.toVector w) == Sym.toVector (D8.inverse w)
+prop_D8_rotate w     = I.rotate     (Sym.toVector w) == Sym.toVector (D8.rotate w)
+
+testsD8 =
+    [ ("D8/orbit",       check prop_D8_orbit)
+    , ("D8/reverse",     check prop_D8_reverse)
+    , ("D8/complement",  check prop_D8_complement)
+    , ("D8/inverse",     check prop_D8_inverse)
+    , ("D8/rotate",      check prop_D8_rotate)
+    ]
+
+---------------------------------------------------------------------------------
+-- Properties for Math.Sym.Stat
+---------------------------------------------------------------------------------
+
+-- the group theoretical inverse of w
+inverse :: (Ord a) => [a] -> [Int]
+inverse w = map snd . sort $ zip w [0..]
+
+-- the standardization of w
+st :: (Ord a) => [a] -> [Int]
+st = inverse . inverse
+
+ascents, descents :: (Ord a) => [a] -> [(a, a)]
+ascents  w = filter (uncurry (<)) $ zip w (tail w)
+descents w = filter (uncurry (>)) $ zip w (tail w)
+
+peaks          w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x < y, y > z ]
+valleys        w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x > y, y < z ]
+doubleAscents  w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x < y, y < z ]
+doubleDescents w = [ v | v@(x,y,z) <- zip3 w (tail w) (tail (tail w)), x > y, y > z ]
+
+inversions :: (Ord a) => [a] -> [(a, a)]
+inversions w = init (tails w) >>= \(x:xs) -> [ (x,y) | y<-xs, x > y ]
+
+records :: (a -> a -> Bool) -> [a] -> [a]
+records f []     = []
+records f (x:xs) = records' f [x] xs where
+    records' f recs       []     = recs
+    records' f recs@(r:_) (x:xs) = records' f (if f r x then x:recs else recs) xs
+
+lMinima, lMaxima, rMinima, rMaxima :: (Ord a) => [a] -> [a]
+
+lMinima = reverse . records (>)
+lMaxima = reverse . records (<)
+rMinima = records (>) . reverse
+rMaxima = records (<) . reverse
+
+excedances  xs = map fst . filter (\(i,a)->i <  fromIntegral a) $ zip [0..] xs
+fixedpoints xs = map fst . filter (\(i,a)->i == fromIntegral a) $ zip [0..] xs
+
+exc, fp :: [Int] -> Int
+exc = length . excedances . st
+fp  = length . fixedpoints . st
+
+runs :: Ord a => (a -> a -> Bool) -> [a] -> [a] -> [[a]]
+runs _ [] [] = []
+runs _ rs [] = [rs]
+runs f [] (x:xs) = runs f [x] xs
+runs f u@(r:_) v@(x:xs) | f r x = runs f (x:u) xs
+                        | otherwise = u : runs f [x] xs
+
+decruns :: Ord a => [a] -> [[a]]
+decruns = runs (>) []
+
+incruns :: Ord a => [a] -> [[a]]
+incruns = runs (<) []
+
+ldr, rdr, lir, rir :: (Ord a) => [a] -> Int
+
+ldr [] = 0
+ldr xs = length . head $ decruns xs
+
+rdr [] = 0
+rdr xs = length . last $ decruns xs
+
+lir [] = 0
+lir xs = length . head $ incruns xs
+
+rir [] = 0
+rir xs = length . last $ incruns xs
+
+-- The stack-sort operator
+stack [] = []
+stack xs = stack left ++ stack right ++ [n]
+    where
+      (left, n:right) = span ( < maximum xs) xs
+
+-- The bubble-sort operator; i.e. one pass of the classical bubble
+-- sort algorithm
+bubble :: Ord a => [a] -> [a]
+bubble = bub []
+    where
+      bub xs []       = reverse xs
+      bub [] (y:ys)   = bub [y] ys
+      bub (x:xs) (y:ys)
+          | x < y     = bub (y:x:xs) ys
+          | otherwise = bub (x:y:xs) ys
+
+-- Like Data.List.intersect, but by assuming that the lists are sorted
+-- uses a faster algorithm
+cap :: Ord a => [a] -> [a] -> [a]
+cap [] ys = []
+cap xs [] = []
+cap xs@(x:xt) ys@(y:yt) = case compare x y of
+                            EQ -> x : cap xt yt
+                            LT -> cap xt ys
+                            GT -> cap xs yt
+
+-- The number of components in a permutation
+comp w = length $ lMaxima w `cap` rMinima (bubble w)
+
+-- rank a la Elizalde
+ep = fst . last . filter (\(k,ys) -> all (k<=) ys) . zip [0..] . inits . st
+
+des, asc, inv, lmin, lmax, rmin, rmax, peak, vall :: [Int] -> Int
+dasc, ddes, maj, comp, ep :: [Int] -> Int
+
+maj w = sum [ i | (i,x,y) <- zip3 [1..] w (tail w), x > y ]
+des  = length . descents
+asc  = length . ascents
+inv  = length . inversions
+lmin = length . lMinima
+lmax = length . lMaxima
+rmin = length . rMinima
+rmax = length . rMaxima
+peak = length . peaks
+vall = length . valleys
+dasc = length . doubleAscents
+ddes = length . doubleDescents
+
+prop_asc  = forAll perm $ \w -> asc  w == S.asc  w
+prop_des  = forAll perm $ \w -> des  w == S.des  w
+prop_exc  = forAll perm $ \w -> exc  w == S.exc  w
+prop_fp   = forAll perm $ \w -> fp   w == S.fp   w
+prop_inv  = forAll perm $ \w -> inv  w == S.inv  w
+prop_maj  = forAll perm $ \w -> maj  w == S.maj  w
+prop_lmin = forAll perm $ \w -> lmin w == S.lmin w
+prop_lmax = forAll perm $ \w -> lmax w == S.lmax w
+prop_rmin = forAll perm $ \w -> rmin w == S.rmin w
+prop_rmax = forAll perm $ \w -> rmax w == S.rmax w
+prop_head = forAll perm $ \w -> not (null w) ==> head (st w) == S.head w
+prop_last = forAll perm $ \w -> not (null w) ==> last (st w) == S.last w
+prop_peak = forAll perm $ \w -> peak w == S.peak w
+prop_vall = forAll perm $ \w -> vall w == S.vall w
+prop_dasc = forAll perm $ \w -> dasc w == S.dasc w
+prop_ddes = forAll perm $ \w -> ddes w == S.ddes w
+prop_ep   = forAll perm $ \w -> ep  w == S.ep  w
+prop_lir  = forAll perm $ \w -> lir  w == S.lir  w
+prop_ldr  = forAll perm $ \w -> ldr  w == S.ldr  w
+prop_rir  = forAll perm $ \w -> rir  w == S.rir  w
+prop_rdr  = forAll perm $ \w -> rdr  w == S.rdr  w
+prop_comp = forAll perm $ \w -> comp w == S.comp w
+
+prop_inv_21 = forAll perm $ \w -> S.inv w == length (Sym.copies (Sym.st "21") w)
+
+testsStat =
+    [ ("asc",    check prop_asc)
+    , ("des",    check prop_des)
+    , ("exc",    check prop_exc)
+    , ("fp",     check prop_fp)
+    , ("inv",    check prop_inv)
+    , ("maj",    check prop_maj)
+    , ("lmin",   check prop_lmin)
+    , ("lmax",   check prop_lmax)
+    , ("rmin",   check prop_rmin)
+    , ("rmax",   check prop_rmax)
+    , ("head",   check prop_head)
+    , ("last",   check prop_last)
+    , ("peak",   check prop_peak)
+    , ("vall",   check prop_vall)
+    , ("dasc",   check prop_dasc)
+    , ("ddes",   check prop_ddes)
+    , ("ep",     check prop_ep)
+    , ("lir",    check prop_lir)
+    , ("ldr",    check prop_ldr)
+    , ("rir",    check prop_rir)
+    , ("rdr",    check prop_rdr)
+    , ("comp",   check prop_comp)
+    , ("inv/21", check prop_inv_21)
+    ]
+
+---------------------------------------------------------------------------------
+-- Main
+---------------------------------------------------------------------------------
+
+tests = testsPerm ++ testsD8 ++ testsStat
+
+main = mapM_ (\(name, t) -> putStr (name ++ ":\t") >> t) tests