diff --git a/Data/CLongArray.hs b/Data/CLongArray.hs
deleted file mode 100644
--- a/Data/CLongArray.hs
+++ /dev/null
@@ -1,155 +0,0 @@
-{-# LANGUAGE MagicHash, UnboxedTuples, ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- Convenience functions for dealing with arrays of 'CLong's.
-
-module Data.CLongArray
-    (
-    -- * Data type
-    CLongArray
-
-    -- * Conversions
-    , fromList
-    , toList
-    , slice
-    , unsafeSlice
-
-    -- * Accessors
-    , size
-    , at
-    , unsafeAt
-
-    -- * Map
-    , imap
-
-    -- * Low level functions
-    , unsafeNew
-    , unsafeWith
-    ) where
-
-import Data.Ord
-import Foreign
-import Foreign.C.Types
-import GHC.Base
-
-infixl 9 `at`
-infixl 9 `unsafeAt`
-
-inlinePerformIO :: IO a -> a
-inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
-{-# INLINE inlinePerformIO #-}
-
-
--- Data type
--- ---------
-
--- | An array of 'CLong's
-data CLongArray = CArr {-# UNPACK #-} !(ForeignPtr CLong) -- elements
-                       {-# UNPACK #-} !Int                -- size
-
-instance Show CLongArray where
-    show w = "fromList " ++ show (toList w)
-
-instance Eq CLongArray where
-    u == v = toList u == toList v
-
-instance Ord CLongArray where
-    compare u v =
-        case comparing size u v of
-          EQ -> comparing toList u v
-          x  -> x
-
-
--- Conversions
--- -----------
-
--- | Construct an array from a list of elements.
-fromList :: [Int] -> CLongArray
-fromList xs = CArr p (length xs)
-    where p = inlinePerformIO $ newForeignPtr finalizerFree =<< newArray (map fromIntegral xs)
-{-# INLINE fromList #-}
-
--- | The list of elements.
-toList :: CLongArray -> [Int]
-toList w = map fromIntegral . inlinePerformIO . unsafeWith w $ peekArray (size w)
-{-# INLINE toList #-}
-
--- | Slice a 'CLongArray' into contiguous segments of the given
--- sizes. Each segment size must be positive and they must sum to the
--- size of the array.
-slice :: [Int] -> CLongArray -> [CLongArray]
-slice ks w
-    | any (<=0) ks     = error "Data.CLongArray.slice: zero or negative parts"
-    | sum ks /= size w = error "Data.CLongArray.slice: parts doesn't sum to size of array"
-    | otherwise        = unsafeSlice ks w
-
--- | Like 'slice' but without range checking.
-unsafeSlice :: [Int] -> CLongArray -> [CLongArray]
-unsafeSlice parts w = inlinePerformIO . unsafeWith w $ go parts
-    where
-      go []     _ = return []
-      go (k:ks) p = do
-        vs <- go ks (advancePtr p k)
-        v  <- unsafeNew k $ \q -> copyArray q p k
-        return (v:vs)
-
-
--- Accessors
--- ---------
-
--- | The size/length of the given array.
-size :: CLongArray -> Int
-size (CArr _ n) = n
-{-# INLINE size #-}
-
--- | @w \`at\` i@ is the value of @w@ at @i@, where @i@ is in @[0..size w-1]@.
-at :: CLongArray -> Int -> Int
-at w i =
-    let n = size w
-    in if (i < 0 || i >= n)
-       then error $ "Data.CLongArray.at: " ++ show i ++ " not in [0.." ++ show (n-1) ++ "]"
-       else unsafeAt w i
-{-# INLINE at #-}
-
--- | Like 'at' but without range checking.
-unsafeAt :: CLongArray -> Int -> Int
-unsafeAt w = fromIntegral . inlinePerformIO . unsafeWith w . flip peekElemOff
-{-# INLINE unsafeAt #-}
-
-
--- Map
--- ---
-
--- | Apply a function to every element of an array and its index.
-imap :: (Int -> CLong -> CLong) -> CLongArray -> CLongArray
-imap f w = inlinePerformIO . unsafeWith w $ \p -> unsafeNew n (go 0 p)
-    where
-      n = size w
-      go i p q
-        | i >= n = return ()
-        | otherwise = do
-            x <- peek p
-            poke q (f i x)
-            go (i+1) (advancePtr p 1) (advancePtr q 1)
-
-
--- Low level functions
--- -------------------
-
--- | Create a new array of the given size that is initialized through
--- an IO action.
-unsafeNew :: Int -> (Ptr CLong -> IO ()) -> IO CLongArray
-unsafeNew n act = do
-  q <- newForeignPtr finalizerFree =<< mallocArray n
-  withForeignPtr q act
-  return $ CArr q n
-{-# INLINE unsafeNew #-}
-
--- | Pass a pointer to the array to an IO action; the array may not be
--- modified through the pointer.
-unsafeWith :: CLongArray -> (Ptr CLong -> IO a) -> IO a
-unsafeWith (CArr p _) = withForeignPtr p
-{-# INLINE unsafeWith #-}
diff --git a/Data/Perm.hs b/Data/Perm.hs
deleted file mode 100644
--- a/Data/Perm.hs
+++ /dev/null
@@ -1,85 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface, TypeSynonymInstances #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- Generating permutations: rank and unrank
-
-module Data.Perm
-    (
-      module Data.CLongArray
-    , Perm
-    , emptyperm
-    , one
-    , idperm
-    , ebb
-    , mkPerm
-    , rank
-    , unrank
-    , perms
-    ) where
-
-import Data.List
-import Data.CLongArray
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
--- | A permutation is just a 'CLongArray'. By convention a permutation
--- of size @n@ is understood to be a permutation of @[0..n-1]@.
-type Perm = CLongArray
-
--- | The unique permutation length zero.
-emptyperm :: Perm
-emptyperm = fromList []
-
--- | The unique permutation length one.
-one :: Perm
-one = fromList [0]
-
--- | The identity permutation.
-idperm :: Int -> Perm
-idperm n = fromList [0..n-1]
-
--- | The reverse of the identity permutation.
-ebb :: Int -> Perm
-ebb n = fromList [n-1,n-2..0]
-
--- | Construct a permutation from a list of elements. As opposed to
--- 'fromList' this is a safe function in the sense that the output of
--- @mkPerm xs@ is guaranteed to be a permutation of @[0..length xs-1]@.
--- E.g., @mkPerm \"baxa\" == fromList [2,0,3,1]@.
-mkPerm :: Ord a => [a] -> Perm
-mkPerm xs =
-    let sti ys = map snd . sort $ zip ys [ 0::Int .. ]
-    in fromList $ (sti . sti) xs
-
-foreign import ccall unsafe "rank.h rank" c_rank
-    :: Ptr CLong -> CLong -> IO CDouble
-
--- | The rank of the given permutation, where the rank is defined as
--- in [W. Myrvold and F. Ruskey, Ranking and Unranking Permutations in
--- Linear Time, Information Processing Letters, 79 (2001) 281-284].
-rank :: Perm -> Integer
-rank w =
-    let n = fromIntegral (size w)
-    in truncate . unsafeDupablePerformIO . unsafeWith w $ flip c_rank n
-{-# INLINE rank #-}
-
-foreign import ccall unsafe "rank.h unrank" c_unrank
-    :: Ptr CLong -> CLong -> CDouble -> IO ()
-
--- | The permutation of size @n@ whose rank is @r@, where the rank
--- is defined as in [W. Myrvold and F. Ruskey, Ranking and Unranking
--- Permutations in Linear Time, Information Processing Letters, 79
--- (2001) 281-284].
-unrank :: Int -> Integer -> Perm
-unrank n r =
-    unsafeDupablePerformIO . unsafeNew n $ \ptr ->
-        c_unrank ptr (fromIntegral n) (fromIntegral r)
-{-# INLINE unrank #-}
-
--- | All permutations of a given size.
-perms :: Int -> [Perm]
-perms n = map (unrank n) [0..nFac-1] where nFac = product [1..toInteger n]
diff --git a/Data/Perm/Internal.hs b/Data/Perm/Internal.hs
deleted file mode 100644
--- a/Data/Perm/Internal.hs
+++ /dev/null
@@ -1,89 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Data.Perm.Internal
-    (
-      Set
-    , normalize
-    , subsets
-    ) where
-
-import Data.List
-import Data.CLongArray
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
-
--- | A set is represented by an increasing array of non-negative
--- integers.
-type Set = CLongArray
-
--- Utils
--- -----
-
--- | Sort and remove duplicates.
-normalize :: Ord a => [a] -> [a]
-normalize = map head . group . sort
-
-
--- Bitmasks
--- --------
-
--- 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 -> CLongArray
-    ones m = fromList . take (popCount m) $ filter (testBit m) [0..]
-
-instance Bitmask CLong where
-    next = nextCLong
-    ones = onesCLong
-
-instance Bitmask Integer where
-    next = 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
-    | n <= 32   = map ones (bitmasks n k :: [CLong])
-    | otherwise = map ones (bitmasks n k :: [Integer])
-
-foreign import ccall unsafe "bit.h next" c_next :: CLong -> CLong
-
--- | Lexicographically, the next 'CLong' with the same Hamming weight.
-nextCLong :: CLong -> CLong
-nextCLong = c_next
-
-foreign import ccall unsafe "bit.h ones" c_ones :: Ptr CLong -> CLong -> IO ()
-
--- | @onesCLong m@ gives the indices whose bits are set in @m@.
-onesCLong :: CLong -> CLongArray
-onesCLong m = unsafeDupablePerformIO . unsafeNew (popCount m) $ flip c_ones m
-
--- | 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/Data/Permgram.hs b/Data/Permgram.hs
deleted file mode 100644
--- a/Data/Permgram.hs
+++ /dev/null
@@ -1,96 +0,0 @@
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- Permutation diagrams, or permutations as monads.
-
-module Data.Permgram
-    (
-    -- * Data types
-      Label
-    , Permgram
-
-    -- * Accessors
-    , perm
-    , label
-    , size
-
-    -- * Construct permgrams
-    , permgram
-    , inverse
-    ) where
-
-import Data.Ord
-import Data.List
-import Data.Perm (Perm)
-import qualified Data.Perm as P
-import Data.Array.Unboxed
-
--- | The purpose of this data type is to assign labels to the indices of
--- a given permutation.
-type Label a = Array Int a
-
--- | A permgram consists of a permutation together with a label for each
--- index of the permutation.
-data Permgram a = PGram {
-      -- | The underlying permutation.
-      perm  :: Perm
-      -- | The assignment of labels to indices.
-    , label :: Label a
-    }
-
-constituents :: Permgram a -> (Perm, [a])
-constituents (PGram v f) = (v, elems f)
-
-instance Show a => Show (Permgram a) where
-    show w =
-        let (v, ys) = constituents w
-        in unwords ["permgram", "(" ++ show v ++ ")", show ys]
-
-instance Eq a => Eq (Permgram a) where
-    u == v = constituents u == constituents v
-
-instance Ord a => Ord (Permgram a) where
-    compare u v =
-        case comparing size u v of
-          EQ -> comparing constituents u v
-          x  -> x
-
--- | Construct a permgram from an underlying permutation and a list of
--- labels.
-permgram :: Perm -> [a] -> Permgram a
-permgram v = PGram v . listArray (0, P.size v - 1) . cycle
-
--- | The inverse permgram. It's obtained by mirroring the permgram in
--- the /x=y/ diagonal.
-inverse :: Permgram a -> Permgram a
-inverse (PGram u f) = PGram (P.fromList v) (listArray (0,n-1) (map (f!) v))
-    where
-      v = map snd . sort $ zip (P.toList u) [0..] -- v = u^{-1}
-      n = P.size u
-
--- | The size of a permgram is the size of the underlying permutation.
-size :: Permgram a -> Int
-size = P.size . perm
-
-instance Functor Permgram where
-    fmap f w = w { label = amap f (label w) }
-
-instance Monad Permgram where
-    return x = permgram (P.fromList [0]) [x]
-    w >>= f  = joinPermgram $ fmap f w
-
-joinPermgram :: Permgram (Permgram a) -> Permgram a
-joinPermgram w@(PGram u f) = PGram (P.fromList xs) (listArray (0,m-1) ys)
-    where
-      len = amap size f
-      m = sum $ elems len
-      n = size w
-      uInverse = map snd . sort $ zip (P.toList u) [0..]
-      a :: UArray Int Int
-      a = listArray (0,n-1) . scanl (+) 0 $ map (len!) uInverse
-      (xs, ys) = unzip $ do
-        i <- [0..n-1]
-        let PGram v g = f!i
-        let d = a ! (u `P.unsafeAt` i)
-        [ (d + v `P.unsafeAt` j, g!j) | j <- [0..len!i-1] ]
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c)2012, 2013, Anders Claesson
+Copyright (c)2012-2014, Anders Claesson
 
 All rights reserved.
 
diff --git a/Math/Perm.hs b/Math/Perm.hs
deleted file mode 100644
--- a/Math/Perm.hs
+++ /dev/null
@@ -1,18 +0,0 @@
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- A meta module collecting all Perm-modules, except those that are best
--- imported \"qualified\".
-
-module Math.Perm (module P) where
-
-import Data.Perm                 as P
-import Math.Perm.Class           as P
-import Math.Perm.Component       as P
-import Math.Perm.Constructions   as P
-import Math.Perm.Bijection       as P
-import Math.Perm.Group           as P
-import Math.Perm.Pattern         as P
-import Math.Perm.Simple          as P
-import Math.Perm.Sort            as P
diff --git a/Math/Perm/Bijection.hs b/Math/Perm/Bijection.hs
deleted file mode 100644
--- a/Math/Perm/Bijection.hs
+++ /dev/null
@@ -1,39 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.Bijection
-    (
-      simionSchmidt
-    , simionSchmidt'
-    ) where
-
-import Data.Perm
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
-foreign import ccall unsafe "bij.h simion_schmidt" c_simion_schmidt
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
-foreign import ccall unsafe "bij.h simion_schmidt_inverse" c_simion_schmidt'
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
-marshal :: (Ptr CLong -> Ptr CLong -> CLong -> IO ()) -> Perm -> Perm
-marshal bij w =
-    unsafeDupablePerformIO . unsafeWith w $ \p -> do
-      let n = size w
-      unsafeNew n $ \q -> bij q p (fromIntegral n)
-{-# INLINE marshal #-}
-
--- | The Simion-Schmidt bijection from Av(123) onto Av(132).
-simionSchmidt :: Perm -> Perm
-simionSchmidt = marshal c_simion_schmidt
-
--- | The inverse of the Simion-Schmidt bijection. It is a function
--- from Av(132) to Av(123).
-simionSchmidt' :: Perm -> Perm
-simionSchmidt' = marshal c_simion_schmidt'
diff --git a/Math/Perm/Class.hs b/Math/Perm/Class.hs
deleted file mode 100644
--- a/Math/Perm/Class.hs
+++ /dev/null
@@ -1,187 +0,0 @@
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.Class
-    (
-      inc
-    , dec
-    , av1
-    , av12
-    , av21
-    , av123
-    , av132
-    , av213
-    , av231
-    , av312
-    , av321
-    , av1243
-    , av1324
-    , av2134
-    , av
-    , vee
-    , caret
-    , gt
-    , lt
-    , wedges
-    , separables
-    , kLayered
-    , layered
-    , kFibonacci
-    , fibonacci
-    ) where
-
-import Data.Perm
-import Math.Perm.Bijection
-import Math.Perm.Constructions
-import Data.Perm.Internal
-import Math.Perm.Pattern
-import qualified Math.Perm.D8 as D8
-
--- | The class of increasing permutations.
-inc :: Int -> [Perm]
-inc = av21
-
--- | The class of decreasing permutations.
-dec :: Int -> [Perm]
-dec = av12
-
--- | Av(1)
-av1 :: Int -> [Perm]
-av1 0 = [emptyperm]
-av1 _ = []
-
--- | Av(12)
-av12 :: Int -> [Perm]
-av12 n = [ebb n]
-
--- | Av(21)
-av21 :: Int -> [Perm]
-av21 n = [idperm n]
-
--- | Av(123)
-av123 :: Int -> [Perm]
-av123 = map simionSchmidt' . av132
-
--- | Av(132)
-av132 :: Int -> [Perm]
-av132 = map D8.reverse . av231
-
--- | Av(213)
-av213 :: Int -> [Perm]
-av213 = map D8.complement . av231
-
--- | Av(231); also know as the stack sortable permutations.
-av231 :: Int -> [Perm]
-av231 0 = [emptyperm]
-av231 n = do
-  k <- [0..n-1]
-  s <- streamAv231 !! k
-  t <- streamAv231 !! (n-k-1)
-  return $ s /+/ (one \-\ t)
-
-streamAv231 :: [[Perm]]
-streamAv231 = map av231 [0..]
-
--- | Av(312)
-av312 :: Int -> [Perm]
-av312 = map D8.inverse . av231
-
--- | Av(321)
-av321 :: Int -> [Perm]
-av321 = map D8.complement . av123
-
--- | Av(1243)
-av1243 :: Int -> [Perm]
-av1243 n = avoiders [fromList [0,1,3,2]] (perms n)
-
--- | Av(1324)
-av1324 :: Int -> [Perm]
-av1324 n = avoiders [fromList [0,2,1,3]] (perms n)
-
--- | Av(2134)
-av2134 :: Int -> [Perm]
-av2134 n = avoiders [fromList [1,0,2,3]] (perms n)
-
--- | Av(s) where s is a string of one or more patterns, using space as a
--- seperator.
-av :: String -> Int -> [Perm]
-av s = avoiders (map mkPerm (words s)) . perms
-
--- | The V-class is Av(132, 231). It is so named because the diagram of
--- a typical permutation in this class is shaped like a V.
-vee :: Int -> [Perm]
-vee = (streamVee !!)
-
-streamVee :: [[Perm]]
-streamVee = [emptyperm] : [one] : zipWith (++) vee_n n_vee
-    where
-      n_vee = (map.map) (one \-\) ws
-      vee_n = (map.map) (/+/ one) ws
-      ws    = tail streamVee
-
--- | The ∧-class is Av(213, 312). It is so named because the diagram of
--- a typical permutation in this class is shaped like a ∧.
-caret :: Int -> [Perm]
-caret = map D8.complement . vee
-
--- | The >-class is Av(132, 312). It is so named because the diagram of
--- a typical permutation in this class is shaped like a >.
-gt :: Int -> [Perm]
-gt = map D8.rotate . vee
-
--- | The <-class is Av(213, 231). It is so named because the diagram of
--- a typical permutation in this class is shaped like a <.
-lt :: Int -> [Perm]
-lt = map D8.reverse . gt
-
-union :: [Int -> [Perm]] -> Int -> [Perm]
-union cs n = normalize $ concat [ c n | c <- cs ]
-
--- | The union of 'vee', 'caret', 'gt' and 'lt'.
-wedges :: Int -> [Perm]
-wedges = union [vee, caret, gt, lt]
-
-compositions :: Int -> Int -> [[Int]]
-compositions 0 0 = [[]]
-compositions 0 _ = []
-compositions _ 0 = []
-compositions k n = [1..n] >>= \i -> map (i:) (compositions (k-1) (n-i))
-
-boundedCompositions :: Int -> Int -> Int -> [[Int]]
-boundedCompositions _ 0 0 = [[]]
-boundedCompositions _ 0 _ = []
-boundedCompositions _ _ 0 = []
-boundedCompositions b k n = [1..b] >>= \i -> map (i:) (boundedCompositions b (k-1) (n-i))
-
--- | The class of separable permutations; it is identical to Av(2413,3142).
-separables :: Int -> [Perm]
-separables 0 = [emptyperm]
-separables 1 = [one]
-separables n = pIndec n ++ mIndec n
-    where
-      comps  m = [2..m] >>= \k -> compositions k m
-      pIndec 0 = []
-      pIndec 1 = [one]
-      pIndec m = comps m >>= map skewSum . mapM (streamMIndec !!)
-      mIndec m = map D8.complement $ pIndec m
-      streamMIndec = map mIndec [0..]
-
--- | The class of layered permutations with /k/ layers.
-kLayered :: Int -> Int -> [Perm]
-kLayered k = map (directSum . map ebb) . compositions k
-
--- | The class of layered permutations.
-layered :: Int -> [Perm]
-layered n = [1..n] >>= flip kLayered n
-
--- | The class of Fibonacci permutations with /k/ layers. A /Fibonacci permutation/
--- is a layered permutation whose layers are all of size 1 or 2.
-kFibonacci :: Int -> Int -> [Perm]
-kFibonacci k = map (directSum . map ebb) . boundedCompositions 2 k
-
--- | The class of Fibonacci permutations. A /Fibonacci permutation/ is a
--- layered permutation whose layers are all of size 1 or 2.
-fibonacci :: Int -> [Perm]
-fibonacci n = [1..n] >>= flip kFibonacci n
diff --git a/Math/Perm/Component.hs b/Math/Perm/Component.hs
deleted file mode 100644
--- a/Math/Perm/Component.hs
+++ /dev/null
@@ -1,76 +0,0 @@
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- Components of permutations.
--- 
-
-module Math.Perm.Component
-    (
-      components
-    , skewComponents
-    , leftMaxima
-    , leftMinima
-    , rightMaxima
-    , rightMinima
-    ) where
-
-import Foreign
-import System.IO.Unsafe
-import Data.Perm
-import qualified Math.Perm.D8 as D8
-
--- Positions /i/ such that /max{ w[j] : j <= i } = i/. These positions
--- mark the boundaries of components.
-comps :: Perm -> [Int]
-comps w = unsafeDupablePerformIO . unsafeWith w $ go [] 0 0
-    where
-      n = size w
-      go ks m i p
-        | i >= n = return (reverse ks)
-        | otherwise =
-            do y <- fromIntegral `fmap` peek p
-               let p'  = advancePtr p 1
-               let i'  = i+1
-               let m'  = if y > m then y else m
-               let ks' = if m' == i then i:ks else ks
-               go ks' m' i' p'
-
--- | The list of (plus) components.
-components :: Perm -> [Perm]
-components w =
-    let ds = 0 : map (+1) (comps w)
-        ks = zipWith (-) (tail ds) ds
-        ws = unsafeSlice ks w
-    in zipWith (\d v -> imap (\_ x -> x - fromIntegral d) v) ds ws
-
--- | The list of skew components, also called minus components.
-skewComponents :: Perm -> [Perm]
-skewComponents = map D8.complement . components . D8.complement
-
-records :: (a -> a -> Bool) -> [a] -> [a]
-records _ []     = []
-records f (x:xs) = recs [x] xs
-    where
-      recs rs@(r:_) (y:ys) = recs ((if f r y then y else r):rs) ys
-      recs rs       _      = rs
-
--- | For each position, left-to-right, records the largest value seen
--- thus far.
-leftMaxima :: Perm -> [Int]
-leftMaxima w = reverse $ records (<) (toList w)
-
--- | For each position, left-to-right, records the smallest value seen
--- thus far.
-leftMinima :: Perm -> [Int]
-leftMinima w = reverse $ records (>) (toList w)
-
--- | For each position, /right-to-left/, records the largest value seen
--- thus far.
-rightMaxima :: Perm -> [Int]
-rightMaxima w = records (<) (reverse (toList w))
-
--- | For each position, /right-to-left/, records the smallest value seen
--- thus far.
-rightMinima :: Perm -> [Int]
-rightMinima w = records (>) (reverse (toList w))
diff --git a/Math/Perm/Constructions.hs b/Math/Perm/Constructions.hs
deleted file mode 100644
--- a/Math/Perm/Constructions.hs
+++ /dev/null
@@ -1,62 +0,0 @@
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- Sum, skew sum, etc
--- 
-
-module Math.Perm.Constructions
-    (
-      (/+/)
-    , (\-\)
-    , directSum
-    , skewSum
-    , inflate
-    ) where
-
-import Foreign
-import System.IO.Unsafe
-import Control.Monad
-import Data.Perm
-import qualified Data.Permgram as G
-import qualified Math.Perm.D8 as D8
-
-infixl 6 /+/
-infixl 6 \-\
-
--- | The /direct sum/ of two permutations.
-(/+/) :: Perm -> Perm -> Perm
-(/+/) u v =
-   let k  = size u
-       l  = size v
-       v' = imap (\_ x -> x + fromIntegral k) v
-   in unsafeDupablePerformIO . unsafeNew (k+l) $ \p ->
-       let q = advancePtr p k
-       in unsafeWith u  $ \uPtr ->
-          unsafeWith v' $ \vPtr -> do
-              copyArray p uPtr k
-              copyArray q vPtr l
-
--- | The direct sum of a list of permutations.
-directSum :: [Perm] -> Perm
-directSum = foldr (/+/) emptyperm
-
--- | The /skew sum/ of two permutations.
-(\-\) :: Perm -> Perm -> Perm
-(\-\) u v = D8.complement $ D8.complement u /+/ D8.complement v
-
--- | The skew sum of a list of permutations.
-skewSum :: [Perm] -> Perm
-skewSum = foldr (\-\) emptyperm
-
--- | @inflate w vs@ is the /inflation/ of @w@ by @vs@. It is the
--- permutation of length @sum (map size vs)@ obtained by replacing
--- each entry @w!i@ by an interval that is order isomorphic to @vs!i@
--- in such a way that the intervals are order isomorphic to @w@. In
--- particular,
--- 
--- > u /+/ v == inflate (mkPerm "12") [u,v]
--- > u \-\ v == inflate (mkPerm "21") [u,v]
--- 
-inflate :: Perm -> [Perm] -> Perm
-inflate w = G.perm . join . G.permgram w . map (`G.permgram` [()])
diff --git a/Math/Perm/D8.hs b/Math/Perm/D8.hs
deleted file mode 100644
--- a/Math/Perm/D8.hs
+++ /dev/null
@@ -1,156 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.D8
-    (
-    -- * The group elements
-      r0, r1, r2, r3
-    , s0, s1, s2, s3
-
-    -- * D8, the klein four-group, and orbits
-    , d8
-    , klein4
-    , orbit
-    , symmetryClasses
-    , d8Classes
-    , klein4Classes
-
-    -- * Aliases
-    , rotate
-    , complement
-    , reverse
-    , inverse
-    ) where
-
-import Data.List hiding (reverse)
-import Prelude hiding (reverse)
-import Data.Perm
-import Data.Perm.Internal
-import Foreign hiding (complement, rotate)
-import Foreign.C.Types
-import System.IO.Unsafe
-
-
--- The group elements
--- ------------------
-
--- | Ration by 0 degrees, i.e. the identity map.
-r0 :: Perm -> Perm
-r0 w = w
-
--- | Ration by 90 degrees clockwise.
-r1 :: Perm -> Perm
-r1 = s2 . s1
-
--- | Ration by 2*90 = 180 degrees clockwise.
-r2 :: Perm -> Perm
-r2 = r1 . r1
-
--- | Ration by 3*90 = 270 degrees clockwise.
-r3 :: Perm -> Perm
-r3 = r2 . r1
-
--- | Reflection through a horizontal axis (also called 'complement').
-s0 :: Perm -> Perm
-s0 = complement
-
--- | Reflection through a vertical axis (also called 'reverse').
-s1 :: Perm -> Perm
-s1 = reverse
-
--- | Reflection through the main diagonal (also called 'inverse').
-s2 :: Perm -> Perm
-s2 = inverse
-
--- | Reflection through the anti-diagonal.
-s3 :: Perm -> Perm
-s3 = s1 . r1
-
-
--- D8, the klein four-group, and orbits
--- ------------------------------------
-
--- | The dihedral group of order 8 (the symmetries of a square); that is,
--- 
--- > d8 = [r0, r1, r2, r3, s0, s1, s2, s3]
--- 
-d8 :: [Perm -> Perm]
-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 :: [Perm -> Perm]
-klein4 = [r0, r2, s0, s1]
-
--- | @orbit fs x@ is the orbit of @x@ under the /group/ of function @fs@. E.g.,
--- 
--- > orbit klein4 "2314" == ["1423","2314","3241","4132"]
--- 
-orbit :: [Perm -> Perm] -> Perm -> [Perm]
-orbit fs x = normalize [ f x | f <- fs ]
-
--- | @symmetryClasses fs xs@ is the list of equivalence classes under
--- the action of the /group/ of functions @fs@.
-symmetryClasses :: [Perm -> Perm] -> [Perm] -> [[Perm]]
-symmetryClasses _  [] = []
-symmetryClasses fs xs@(x:xt) = insert orb $ symmetryClasses fs ys
-    where
-      orb = [ w | w <- orbit fs x, w `elem` xs ]
-      ys  = [ y | y <- xt, y `notElem` orb ]
-
--- | Symmetry classes with respect to D8.
-d8Classes :: [Perm] -> [[Perm]]
-d8Classes = symmetryClasses d8
-
--- | Symmetry classes with respect to Klein4
-klein4Classes :: [Perm] -> [[Perm]]
-klein4Classes = symmetryClasses klein4
-
-
--- Aliases
--- -------
-
-marshal :: (Ptr CLong -> Ptr CLong -> CLong -> IO ()) -> Perm -> Perm
-marshal op w =
-    unsafeDupablePerformIO . unsafeWith w $ \p -> do
-      let n = size w
-      unsafeNew n $ \q -> op q p (fromIntegral n)
-{-# INLINE marshal #-}
-
-foreign import ccall unsafe "d8.h inverse" c_inverse
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
--- | The group theoretical inverse: if @v = inverse u@ then
--- @v \`at\` (u \`at\` i) = i@.
-inverse :: Perm -> Perm
-inverse = marshal c_inverse
-{-# INLINE inverse #-}
-
-foreign import ccall unsafe "d8.h reverse" c_reverse
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
--- | The reverse of the given permutation: if @v = reverse u@ then
--- @v \`at\` i = u \`at\` (n-1-i)@.
-reverse :: Perm -> Perm
-reverse = marshal c_reverse
-{-# INLINE reverse #-}
-
-foreign import ccall unsafe "d8.h complement" c_complement
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
--- | The complement of the given permutation: if @v = complement u@ then
--- @v \`at\` i = n - 1 - u \`at\` i@.
-complement :: Perm -> Perm
-complement = marshal c_complement
-{-# INLINE complement #-}
-
--- | @rotate = r1 = inverse . reverse@
-rotate :: Perm -> Perm
-rotate = r1
diff --git a/Math/Perm/Group.hs b/Math/Perm/Group.hs
deleted file mode 100644
--- a/Math/Perm/Group.hs
+++ /dev/null
@@ -1,48 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.Group
-    (
-      compose
-    , act
-    ) where
-
-import Data.Perm
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
-
-marshal :: (Ptr CLong -> Ptr CLong -> CLong -> Ptr CLong -> CLong -> IO ())
-        -> Perm -> Perm -> Perm
-marshal op u v =
-    unsafeDupablePerformIO $
-    unsafeWith u $ \u' ->
-    unsafeWith v $ \v' -> do
-      let k = size u
-      let n = size v
-      let m = max k n
-      unsafeNew m $ \p -> op p u' (fromIntegral k) v' (fromIntegral n)
-{-# INLINE marshal #-}
-
-foreign import ccall unsafe "group.h compose" c_compose
-    :: Ptr CLong -> Ptr CLong -> CLong -> Ptr CLong -> CLong -> IO ()
-
--- | The product/composition of @u@ and @v@: if @w = u `compose` v@
--- then @w `at ` i = v \`at\` (u \`at\` i)@.
-compose :: Perm -> Perm -> Perm
-compose = marshal c_compose
-{-# INLINE compose #-}
-
-foreign import ccall unsafe "group.h act" c_act
-    :: Ptr CLong -> Ptr CLong -> CLong -> Ptr CLong -> CLong -> IO ()
-
--- | The (left) group action of Perm on itself: if @w = u `act` v@
--- then @w `at ` (u \`at\` i) = v \`at\` i@.
-act :: Perm -> Perm -> Perm
-act = marshal c_act
-{-# INLINE act #-}
diff --git a/Math/Perm/Pattern.hs b/Math/Perm/Pattern.hs
deleted file mode 100644
--- a/Math/Perm/Pattern.hs
+++ /dev/null
@@ -1,107 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.Pattern
-    (
-      Pattern
-    , Set
-    , ordiso
-    , subsets
-    , copiesOf
-    , contains
-    , avoids
-    , avoidsAll
-    , avoiders
-    , minima
-    , maxima
-    , coeff
-    ) where
-
-import Data.Perm (Perm, perms)
-import Data.Perm.Internal
-import Data.CLongArray
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
--- | Pattern is just an alias for permutation.
-type Pattern = Perm
-
-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 :: Perm -> Perm -> Set -> Bool
-ordiso u v m =
-    let k = fromIntegral (size u)
-    in unsafeDupablePerformIO $
-       unsafeWith u $ \u' ->
-       unsafeWith v $ \v' ->
-       unsafeWith m $ \m' ->
-           return . toBool $ c_ordiso u' v' m' k
-{-# INLINE ordiso #-}
-
--- | @copies p w@ is the list of sets that represent copies of @p@ in @w@.
-copiesOf :: Pattern -> Perm -> [Set]
-copiesOf p w = filter (ordiso p w) $ subsets (size w) (size p)
-{-# INLINE copiesOf #-}
-
--- | @w `contains` p@ is a predicate determining if @w@ contains the pattern @p@.
-contains :: Perm -> Pattern -> Bool
-w `contains` p = not $ w `avoids` p
-
--- | @w `avoids` p@ is a predicate determining if @w@ avoids the pattern @p@.
-avoids :: Perm -> Pattern -> Bool
-w `avoids` p = null $ copiesOf p w
-
--- | @w `avoidsAll` ps@ is a predicate determining if @w@ avoids the patterns @ps@.
-avoidsAll :: Perm -> [Pattern] -> Bool
-w `avoidsAll` ps = all (w `avoids`) ps
-
--- | @avoiders ps ws@ is the list of permutations in @ws@ avoiding the
--- patterns in @ps@.
-avoiders :: [Pattern] -> [Perm] -> [Perm]
-avoiders ps ws = foldl (flip avoiders1) ws ps
-
--- @avoiders1 p ws@ is the list of permutations in @ws@ avoiding the
--- pattern @p@.
-avoiders1 :: Pattern -> [Perm] -> [Perm]
-avoiders1 _ [] = []
-avoiders1 q vs@(v:_) = filter avoids_q us ++ filter (`avoids` q) ws
-    where
-      n = size v
-      k = size q
-      (us, ws) = span (\u -> size u == n) vs
-      xs = subsets n k
-      avoids_q u = not $ any (ordiso q u) xs
-
--- | The set of minimal elements with respect to containment.
-minima :: [Pattern] -> [Pattern]
-minima [] = []
-minima ws = let (v:vs) = normalize ws
-            in v : minima (avoiders [v] vs)
-
--- | The set of maximal elements with respect to containment.
-maxima :: [Pattern] -> [Pattern]
-maxima [] = []
-maxima ws = let (v:vs) = reverse $ normalize ws
-            in v : maxima (filter (avoids v) vs)
-
--- | @coeff f v@ is the coefficient of @v@ when expanding the
--- permutation statistic @f@ as a sum of permutations/patterns. See
--- Petter Brändén and Anders Claesson: Mesh patterns and the expansion
--- of permutation statistics as sums of permutation patterns, The
--- Electronic Journal of Combinatorics 18(2) (2011),
--- <http://www.combinatorics.org/ojs/index.php/eljc/article/view/v18i2p5>.
-coeff :: (Pattern -> Int) -> Pattern -> Int
-coeff f v = f v + sum [ (-1)^(k - j) * c * f u |
-                        j <- [0 .. k-1]
-                      , u <- perms j
-                      , let c = length $ copiesOf u v
-                      , c > 0
-                      ] where k = size v
diff --git a/Math/Perm/Simple.hs b/Math/Perm/Simple.hs
deleted file mode 100644
--- a/Math/Perm/Simple.hs
+++ /dev/null
@@ -1,25 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.Simple
-    (
-     simple
-    ) where
-
-import Data.Perm
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
-foreign import ccall unsafe "simple.h simple" c_simple
-    :: Ptr CLong -> CLong -> CInt
-
--- | Is the permutation simple?
-simple :: Perm -> Bool
-simple w = toBool . unsafeDupablePerformIO $
-    let n = fromIntegral (size w)
-    in unsafeWith w $ \ptr -> return $ c_simple ptr n
diff --git a/Math/Perm/Sort.hs b/Math/Perm/Sort.hs
deleted file mode 100644
--- a/Math/Perm/Sort.hs
+++ /dev/null
@@ -1,40 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Perm.Sort
-    (
-      stackSort
-    , bubbleSort
-    ) where
-
-import Data.Perm
-import Foreign
-import Foreign.C.Types
-import System.IO.Unsafe
-
-foreign import ccall unsafe "sortop.h stacksort" c_stacksort
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
-foreign import ccall unsafe "sortop.h bubblesort" c_bubblesort
-    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
-
-marshal :: (Ptr CLong -> Ptr CLong -> CLong -> IO ()) -> Perm -> Perm
-marshal op w =
-    unsafeDupablePerformIO . unsafeWith w $ \p -> do
-      let n = size w
-      unsafeNew n $ \q -> op q p (fromIntegral n)
-{-# INLINE marshal #-}
-
--- | One pass of stack-sort.
-stackSort :: Perm -> Perm
-stackSort = marshal c_stacksort
-{-# INLINE stackSort #-}
-
--- | One pass of bubble-sort.
-bubbleSort :: Perm -> Perm
-bubbleSort = marshal c_bubblesort
-{-# INLINE bubbleSort #-}
diff --git a/Math/Perm/Stat.hs b/Math/Perm/Stat.hs
deleted file mode 100644
--- a/Math/Perm/Stat.hs
+++ /dev/null
@@ -1,275 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
--- Common permutation statistics. To avoid name clashes this module is
--- best imported @qualified@; e.g.
--- 
--- > import qualified Math.Perm.Stat as S
--- 
-
-module Math.Perm.Stat 
-    (
-      asc         -- ascents
-    , des         -- descents
-    , exc         -- excedances
-    , fp          -- fixed points
-    , sfp         -- strong fixed points
-    , cyc         -- cycles
-    , inv         -- inversions
-    , maj         -- the major index
-    , comaj       -- the co-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
-    , scomp       -- skew components
-    , ep          -- rank a la Elizalde & Pak
-    , dim         -- dimension
-    , asc0        -- small ascents
-    , des0        -- small descents
---    , shad        -- shadow
-    ) where
-
-import Prelude hiding (head, last)
-import Data.Perm
-import qualified Math.Perm.D8 as D8
-import Foreign.Ptr
-import Foreign.C.Types
-import System.IO.Unsafe
-
-marshal :: (Ptr CLong -> CLong -> CLong) -> Perm -> Int
-marshal f w =
-    fromIntegral . unsafeDupablePerformIO . unsafeWith w $ \p ->
-        return $ f p (fromIntegral (size w))
-{-# INLINE marshal #-}
-
-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 sfp" c_sfp
-    :: Ptr CLong -> CLong -> CLong
-
-foreign import ccall unsafe "stat.h cyc" c_cyc
-    :: 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 comaj" c_comaj
-    :: 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
-
-foreign import ccall unsafe "stat.h dim" c_dim
-    :: Ptr CLong -> CLong -> CLong
-
-foreign import ccall unsafe "stat.h asc0" c_asc0
-    :: Ptr CLong -> CLong -> CLong
-
-foreign import ccall unsafe "stat.h des0" c_des0
-    :: Ptr CLong -> CLong -> CLong
-
--- | The number of ascents. An /ascent/ in @w@ is an index @i@ such
--- that @w[i] \< w[i+1]@.
-asc :: Perm -> Int
-asc = marshal c_asc
-
--- | The number of descents. A /descent/ in @w@ is an index @i@ such
--- that @w[i] > w[i+1]@.
-des :: Perm -> Int
-des = marshal c_des
-
--- | The number of /excedances/: positions @i@ such that @w[i] > i@.
-exc :: Perm -> Int
-exc = marshal c_exc
-
--- | The number of /fixed points/: positions @i@ such that @w[i] == i@.
-fp :: Perm -> Int
-fp = marshal c_fp
-
--- | The number of /strong fixed points/ (also called splitters):
--- positions @i@ such that @w[j] \< i@ for @j \< i@ and @w[j] \> i@ for @j \> i@.
-sfp :: Perm -> Int
-sfp = marshal c_sfp
-
--- | The number of /cycles/:
--- orbits of the permutation when viewed as a function.
-cyc :: Perm -> Int
-cyc = marshal c_cyc
-
--- | The number of /inversions/:
--- pairs @\(i,j\)@ such that @i \< j@ and @w[i] > w[j]@.
-inv :: Perm -> Int
-inv = marshal c_inv
-
--- | /The major index/ is the sum of descents.
-maj :: Perm -> Int
-maj = marshal c_maj
-
--- | /The co-major index/ is the sum of descents.
-comaj :: Perm -> Int
-comaj = marshal c_comaj
-
--- | The number of /peaks/:
--- positions @i@ such that @w[i-1] \< w[i]@ and @w[i] \> w[i+1]@.
-peak :: Perm -> Int
-peak = marshal c_peak
-
--- | The number of /valleys/:
--- positions @i@ such that @w[i-1] \> w[i]@ and @w[i] \< w[i+1]@.
-vall :: Perm -> Int
-vall = marshal c_vall
-
--- | The number of /double ascents/:
--- positions @i@ such that @w[i-1] \<  w[i] \< w[i+1]@.
-dasc :: Perm -> Int
-dasc = marshal c_dasc
-
--- | The number of /double descents/:
--- positions @i@ such that @w[i-1] \>  w[i] \> w[i+1]@.
-ddes :: Perm -> Int
-ddes = marshal c_ddes
-
--- | The number of /left-to-right minima/:
--- positions @i@ such that @w[i] \< w[j]@ for all @j \< i@.
-lmin :: Perm -> Int
-lmin = marshal c_lmin
-
--- | The number of /left-to-right maxima/:
--- positions @i@ such that @w[i] \> w[j]@ for all @j \< i@.
-lmax :: Perm -> Int
-lmax = marshal c_lmax
-
--- | The number of /right-to-left minima/:
--- positions @i@ such that @w[i] \< w[j]@ for all @j \> i@.
-rmin :: Perm -> Int
-rmin = lmin . D8.reverse
-
--- | The number of /right-to-left maxima/:
--- positions @i@ such that @w[i] \> w[j]@ for all @j \> i@.
-rmax :: Perm -> Int
-rmax = lmax . D8.reverse
-
--- | The first (left-most) element in the standardization. E.g.,
--- @head \"231\" = head (fromList [1,2,0]) = 1@.
-head :: Perm -> Int
-head w | size w > 0 = w `unsafeAt` 0
-       | otherwise  = 0
-
--- | The last (right-most) element in the standardization. E.g.,
--- @last \"231\" = last (fromList [1,2,0]) = 0@.
-last :: Perm -> Int
-last w | size w > 0 = w `unsafeAt` (size w - 1)
-       | otherwise  = 0
-
--- | Length of the left-most increasing run: largest @i@ such that
--- @w[0] \< w[1] \< ... \< w[i-1]@.
-lir :: Perm -> Int
-lir = marshal c_lir
-
--- | Length of the left-most decreasing run: largest @i@ such that
--- @w[0] \> w[1] \> ... \> w[i-1]@.
-ldr :: Perm -> Int
-ldr = marshal c_ldr
-
--- | Length of the right-most increasing run: largest @i@ such that
--- @w[n-i] \< ... \< w[n-2] \< w[n-1]@.
-rir :: Perm -> Int
-rir = ldr . D8.reverse
-
--- | Length of the right-most decreasing run: largest @i@ such that
--- @w[n-i] \> ... \> w[n-2] \> w[n-1]@.
-rdr :: Perm -> Int
-rdr = lir . D8.reverse
-
--- | 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 -> Int
-comp = marshal c_comp
-
--- | The number of skew components. E.g., @[5,7,4,6,3,1,0,2]@ has three
--- skew components: @[5,7,4,6]@, @[3]@ and @[1,0,2]@.
-scomp :: Perm -> Int
-scomp = comp . D8.complement
-
--- | 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 -> Int
-ep = marshal c_ep
-
--- | The dimension of a permutation is defined as the largest
--- non-fixed-point, or zero if all points are fixed.
-dim :: Perm -> Int
-dim = marshal c_dim
-
--- | The number of small ascents. A /small ascent/ in @w@ is an index
--- @i@ such that @w[i] + 1 == w[i+1]@.
-asc0 :: Perm -> Int
-asc0 = marshal c_asc0
-
--- | The number of small descents. A /small descent/ in @w@ is an
--- index @i@ such that @w[i] == w[i+1] + 1@.
-des0 :: Perm -> Int
-des0 = marshal c_des0
-
--- | The size of the shadow of @w@. That is, the number of different
--- one point deletions of @w@.
--- shad :: Perm -> Int
--- shad = length . shadow . return . st
diff --git a/Math/Sym.hs b/Math/Sym.hs
deleted file mode 100644
--- a/Math/Sym.hs
+++ /dev/null
@@ -1,125 +0,0 @@
-{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
-
--- |
--- Copyright   : Anders Claesson 2013
--- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
---
-
-module Math.Sym
-    (
-      Permutation(..)
-    , perms
-    , lift
-    , lift2
-    ) where
-
-import Data.Ord
-import Data.List
-import Math.Perm (Perm)
-import qualified Math.Perm as P
-import qualified Math.Perm.D8 as D8
-
-
--- The permutation typeclass
--- -------------------------
-
--- | The class of permutations. Minimal complete definition: 'st',
--- 'act' and 'idperm'. The default implementation of 'size' can be
--- somewhat slow, so you may want to implement it as well.
-class Ord a => Permutation 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 -> Perm
-
-    -- | A (left) /group action/ of 'Perm' on @a@. As for any group
-    -- action it should hold that
-    -- 
-    -- > (u `act` v) `act` w == u `act` (v `act` w)   &&   idperm n `act` v == v
-    -- 
-    -- where @v,w::a@ and @u::Perm@ are of size @n@.
-    act :: Perm -> a -> a
-
-    -- | The size of a permutation. The default implementation derived from
-    -- 
-    -- > size == size . st
-    -- 
-    -- This is not a circular definition as 'size' on 'Perm' 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 = P.size . st
-
-    -- | The identity permutation of the given size.
-    idperm :: Int -> a
-
-    -- | The group theoretical inverse. It should hold that
-    -- 
-    -- > inverse == unst . inverse . st
-    -- 
-    -- and this is the default implementation.
-    {-# INLINE inverse #-}
-    inverse :: a -> a
-    inverse = unst . D8.inverse . st
-
-    -- | 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 (size u)
-    -- 
-    {-# INLINE ordiso #-}
-    ordiso :: Perm -> a -> Bool
-    ordiso u v = u == st v
-
-    -- | The inverse of 'st'. It should hold that
-    -- 
-    -- > unst w == w `act` idperm (P.size w)
-    -- 
-    -- and this is the default implementation.
-    unst :: Permutation a => Perm -> a
-    unst w = w `act` idperm (P.size w)
-
-instance Permutation Perm where
-    st       = id
-    act      = P.act
-    idperm   = P.idperm
-    inverse  = D8.inverse
-    ordiso   = (==)
-    unst     = id
-
--- | A String viewed as a permutation of its characters. The alphabet
--- is ordered as
--- 
--- > ['1'..'9'] ++ ['A'..'Z'] ++ ['a'..]
--- 
-instance Permutation String where
-    st       = P.mkPerm
-    act v    = map snd . sortBy (comparing fst) . zip (P.toList (D8.inverse v))
-    size     = length
-    idperm n = take n $ ['1'..'9'] ++ ['A'..'Z'] ++ ['a'..]
-
--- | The list of all permutations of the given size.
-perms :: Permutation a => Int -> [a]
-perms = map unst . P.perms
-
--- | Lifts a function on 'Perm's to one on any permutations.
-lift :: (Permutation a) => (Perm -> Perm) -> a -> a
-lift f = unst . f . st
-
--- | Like 'lift' but for functions of two variables.
-lift2 :: (Permutation a) => (Perm -> Perm -> Perm) -> a -> a -> a
-lift2 f u v = unst $ f (st u) (st v)
-
diff --git a/Sym.hs b/Sym.hs
new file mode 100644
--- /dev/null
+++ b/Sym.hs
@@ -0,0 +1,134 @@
+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym
+    (
+      Permutation(..)
+    , perms
+    , lift
+    , lift2
+    ) where
+
+import Data.Ord
+import Sym.Perm.SSYT (SSYTPair (..))
+import qualified Sym.Perm.SSYT as Y
+import Data.List
+import Sym.Perm.Meta (Perm)
+import qualified Sym.Perm.Meta as P
+import qualified Sym.Perm.D8 as D8
+
+
+-- The permutation typeclass
+-- -------------------------
+
+-- | The class of permutations. Minimal complete definition: 'st',
+-- 'act' and 'idperm'. The default implementation of 'size' can be
+-- somewhat slow, so you may want to implement it as well.
+class Permutation 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 -> Perm
+
+    -- | A (left) /group action/ of 'Perm' on @a@. As for any group
+    -- action it should hold that
+    -- 
+    -- > (u `act` v) `act` w == u `act` (v `act` w)   &&   idperm n `act` v == v
+    -- 
+    -- where @v,w::a@ and @u::Perm@ are of size @n@.
+    act :: Perm -> a -> a
+
+    -- | The size of a permutation. The default implementation derived from
+    -- 
+    -- > size == size . st
+    -- 
+    -- This is not a circular definition as 'size' on 'Perm' 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 = P.size . st
+
+    -- | The identity permutation of the given size.
+    idperm :: Int -> a
+
+    -- | The group theoretical inverse. It should hold that
+    -- 
+    -- > inverse == unst . inverse . st
+    -- 
+    -- and this is the default implementation.
+    {-# INLINE inverse #-}
+    inverse :: a -> a
+    inverse = unst . D8.inverse . st
+
+    -- | 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 (size u)
+    -- 
+    {-# INLINE ordiso #-}
+    ordiso :: Perm -> a -> Bool
+    ordiso u v = u == st v
+
+    -- | The inverse of 'st'. It should hold that
+    -- 
+    -- > unst w == w `act` idperm (P.size w)
+    -- 
+    -- and this is the default implementation.
+    unst :: Permutation a => Perm -> a
+    unst w = w `act` idperm (P.size w)
+
+instance Permutation Perm where
+    st       = id
+    act      = P.act
+    idperm   = P.idperm
+    inverse  = D8.inverse
+    ordiso   = (==)
+    unst     = id
+
+-- | A String viewed as a permutation of its characters. The alphabet
+-- is ordered as
+-- 
+-- > ['1'..'9'] ++ ['A'..'Z'] ++ ['a'..]
+-- 
+instance Permutation String where
+    st       = P.mkPerm
+    act v    = map snd . sortBy (comparing fst) . zip (P.toList (D8.inverse v))
+    size     = length
+    idperm n = take n $ ['1'..'9'] ++ ['A'..'Z'] ++ ['a'..]
+
+instance Permutation SSYTPair where
+    st = Y.toPerm
+    unst = Y.fromPerm
+    u `act` v = unst $ u `act` st v
+    size (SSYTPair p _) = sum $ map length p
+    idperm n = SSYTPair p p where p = [[0..n-1]]
+    inverse (SSYTPair p q) = SSYTPair q p
+
+-- | The list of all permutations of the given size.
+perms :: Permutation a => Int -> [a]
+perms = map unst . P.perms
+
+-- | Lifts a function on 'Perm's to one on any permutations.
+lift :: (Permutation a) => (Perm -> Perm) -> a -> a
+lift f = unst . f . st
+
+-- | Like 'lift' but for functions of two variables.
+lift2 :: (Permutation a) => (Perm -> Perm -> Perm) -> a -> a -> a
+lift2 f u v = unst $ f (st u) (st v)
diff --git a/Sym/Internal/CLongArray.hs b/Sym/Internal/CLongArray.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Internal/CLongArray.hs
@@ -0,0 +1,170 @@
+{-# LANGUAGE MagicHash, UnboxedTuples #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Convenience functions for dealing with arrays of 'CLong's.
+
+module Sym.Internal.CLongArray
+    (
+    -- * Data type
+    CLongArray
+
+    -- * Conversions
+    , fromList
+    , toList
+    , slice
+    , unsafeSlice
+
+    -- * Accessors
+    , size
+    , at
+    , unsafeAt
+
+    -- * Map
+    , imap
+    , izipWith
+
+    -- * Low level functions
+    , unsafeNew
+    , unsafeWith
+    ) where
+
+import Data.Ord
+import Sym.Internal.Size
+import Foreign
+import Foreign.C.Types
+import GHC.Base
+
+infixl 9 `at`
+infixl 9 `unsafeAt`
+
+inlinePerformIO :: IO a -> a
+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
+{-# INLINE inlinePerformIO #-}
+
+
+-- Data type
+-- ---------
+
+-- | An array of 'CLong's
+data CLongArray = CArr {-# UNPACK #-} !(ForeignPtr CLong) -- elements
+                       {-# UNPACK #-} !Int                -- size
+
+instance Show CLongArray where
+    show w = "fromList " ++ show (toList w)
+
+instance Eq CLongArray where
+    u == v = toList u == toList v
+
+instance Ord CLongArray where
+    compare u v =
+        case comparing size u v of
+          EQ -> comparing toList u v
+          x  -> x
+
+instance Size CLongArray where
+    size (CArr _ n) = n
+    {-# INLINE size #-}
+
+
+-- Conversions
+-- -----------
+
+-- | Construct an array from a list of elements.
+fromList :: [Int] -> CLongArray
+fromList xs = CArr p (length xs)
+  where
+    p = inlinePerformIO $ newForeignPtr finalizerFree =<< newArray (map fromIntegral xs)
+{-# INLINE fromList #-}
+
+-- | The list of elements.
+toList :: CLongArray -> [Int]
+toList w = map fromIntegral . inlinePerformIO . unsafeWith w $ peekArray (size w)
+{-# INLINE toList #-}
+
+-- | Slice a 'CLongArray' into contiguous segments of the given
+-- sizes. Each segment size must be positive and they must sum to the
+-- size of the array.
+slice :: [Int] -> CLongArray -> [CLongArray]
+slice ks w
+    | any (<=0) ks     = error "Sym.Internal.CLongArray.slice: zero or negative parts"
+    | sum ks /= size w = error "Sym.Internal.CLongArray.slice: parts doesn't sum to size of array"
+    | otherwise        = unsafeSlice ks w
+
+-- | Like 'slice' but without range checking.
+unsafeSlice :: [Int] -> CLongArray -> [CLongArray]
+unsafeSlice parts w = inlinePerformIO . unsafeWith w $ go parts
+  where
+    go []     _ = return []
+    go (k:ks) p = do
+      vs <- go ks (advancePtr p k)
+      v  <- unsafeNew k $ \q -> copyArray q p k
+      return (v:vs)
+
+
+-- Accessors
+-- ---------
+
+-- | @w \`at\` i@ is the value of @w@ at @i@, where @i@ is in @[0..size w-1]@.
+at :: CLongArray -> Int -> Int
+at w i =
+    let n = size w
+    in if i < 0 || i >= n
+       then error $ "Sym.Internal.CLongArray.at: " ++ show i ++ " not in [0.." ++ show (n-1) ++ "]"
+       else unsafeAt w i
+{-# INLINE at #-}
+
+-- | Like 'at' but without range checking.
+unsafeAt :: CLongArray -> Int -> Int
+unsafeAt w = fromIntegral . inlinePerformIO . unsafeWith w . flip peekElemOff
+{-# INLINE unsafeAt #-}
+
+
+-- Map and Zip
+-- -----------
+
+-- | Apply a function to every element of an array and its index.
+imap :: (Int -> CLong -> CLong) -> CLongArray -> CLongArray
+imap f w = inlinePerformIO . unsafeWith w $ \p -> unsafeNew n (go 0 p)
+  where
+    n = size w
+    go i p q
+      | i >= n = return ()
+      | otherwise = do
+          x <- peek p
+          poke q (f i x)
+          go (i+1) (advancePtr p 1) (advancePtr q 1)
+
+-- | Apply a function to corresponding pairs of elements and their (shared) index.
+izipWith :: (Int -> CLong -> CLong -> CLong) -> CLongArray -> CLongArray -> CLongArray
+izipWith f u v =
+    inlinePerformIO . unsafeWith u $ \p -> unsafeWith v $ \q -> unsafeNew n (go 0 p q)
+  where
+    n = min (size u) (size v)
+    go i p q r
+      | i >= n = return ()
+      | otherwise = do
+          x <- peek p
+          y <- peek q
+          poke r (f i x y)
+          go (i+1) (advancePtr p 1) (advancePtr q 1) (advancePtr r 1)
+
+-- Low level functions
+-- -------------------
+
+-- | Create a new array of the given size that is initialized through
+-- an IO action.
+unsafeNew :: Int -> (Ptr CLong -> IO ()) -> IO CLongArray
+unsafeNew n act = do
+  q <- newForeignPtr finalizerFree =<< mallocArray n
+  withForeignPtr q act
+  return $ CArr q n
+{-# INLINE unsafeNew #-}
+
+-- | Pass a pointer to the array to an IO action; the array may not be
+-- modified through the pointer.
+unsafeWith :: CLongArray -> (Ptr CLong -> IO a) -> IO a
+unsafeWith (CArr p _) = withForeignPtr p
+{-# INLINE unsafeWith #-}
diff --git a/Sym/Internal/Size.hs b/Sym/Internal/Size.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Internal/Size.hs
@@ -0,0 +1,16 @@
+module Sym.Internal.Size (Size (..)) where
+
+import qualified Data.Set as Set
+
+class Size a where
+    size :: a -> Int
+
+instance Size [a] where
+    size = length
+
+instance Size (Set.Set a) where
+    size = Set.size
+
+instance Size a => Size (Maybe a) where
+    size Nothing  = 0
+    size (Just x) = size x
diff --git a/Sym/Internal/SubSeq.hs b/Sym/Internal/SubSeq.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Internal/SubSeq.hs
@@ -0,0 +1,79 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Internal.SubSeq
+    (
+      module Sym.Internal.CLongArray
+    , SubSeq
+    , choose
+    ) where
+
+import Sym.Internal.CLongArray
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+-- | A SubSeq is represented by an increasing array of non-negative
+-- integers.
+type SubSeq = CLongArray
+
+-- Bitmasks
+-- --------
+
+-- 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 / subsequence of indices whose bits are
+    -- set in @m@. Default implementation:
+    -- 
+    -- > ones m = fromListN (popCount m) $ filter (testBit m) [0..]
+    -- 
+    ones :: a -> SubSeq
+    ones m = fromList . take (popCount m) $ filter (testBit m) [0..]
+
+instance Bitmask CLong where
+    next = nextCLong
+    ones = onesCLong
+
+instance Bitmask Integer where
+    next = 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']
+
+-- | @n \`choose\` k@ is the list of subsequences of @[0..n-1]@ with @k@
+-- elements.
+choose :: Int -> Int -> [SubSeq]
+choose n k
+    | n <= 32   = map ones (bitmasks n k :: [CLong])
+    | otherwise = map ones (bitmasks n k :: [Integer])
+
+foreign import ccall unsafe "bit.h next" c_next :: CLong -> CLong
+
+-- | Lexicographically, the next 'CLong' with the same Hamming weight.
+nextCLong :: CLong -> CLong
+nextCLong = c_next
+
+foreign import ccall unsafe "bit.h ones" c_ones :: Ptr CLong -> CLong -> IO ()
+
+-- | @onesCLong m@ gives the indices whose bits are set in @m@.
+onesCLong :: CLong -> CLongArray
+onesCLong m = unsafeDupablePerformIO . unsafeNew (popCount m) $ flip c_ones m
+
+-- | 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/Sym/Internal/Util.hs b/Sym/Internal/Util.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Internal/Util.hs
@@ -0,0 +1,53 @@
+-- |
+-- Copyright   : Anders Claesson 2014
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Internal.Util
+    (
+      minima
+    , maxima
+    , kSubsets
+    , powerset
+    , nubSort
+    ) where
+
+import Data.List
+import Data.Ord
+import Data.Set (Set)
+import qualified Data.Set as S
+
+-- | The set of minimal elements with respect to inclusion.
+minima :: Ord a => [Set a] -> [Set a]
+minima = minima' . sortBy (comparing S.size)
+  where
+    minima' [] = []
+    minima' (x:xs) = x : minima' [ y | y<-xs, not (x `S.isSubsetOf` y) ]
+
+-- | The set of maximal elements with respect to the given order.
+maxima :: Ord a => [Set a] -> [Set a]
+maxima = maxima' . sortBy (comparing $ \x -> -S.size x)
+  where
+    maxima' [] = []
+    maxima' (x:xs) = x : maxima' [ y | y<-xs, not (y `S.isSubsetOf` x) ]
+
+-- | A list of all k element subsets of the given set.
+kSubsets :: Ord a => Int -> Set a -> [Set a]
+kSubsets 0 _    = [ S.empty ]
+kSubsets k s 
+    | S.null s  = []
+    | otherwise = kSubsets k t ++ map (S.insert x) (kSubsets (k-1) t)
+  where
+    (x,t) = S.deleteFindMin s
+
+-- | A list of all subsets of the given set.
+powerset :: Ord a => Set a -> [Set a]
+powerset s 
+    | S.null s  = [s]
+    | otherwise = ts ++ map (S.insert x) ts
+  where
+    (x,t) = S.deleteFindMin s; ts = powerset t
+
+-- | Sort and remove duplicates.
+nubSort :: Ord a => [a] -> [a]
+nubSort = map head . group . sort
diff --git a/Sym/Perm.hs b/Sym/Perm.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm.hs
@@ -0,0 +1,85 @@
+{-# LANGUAGE ForeignFunctionInterface, TypeSynonymInstances #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Generating permutations: rank and unrank
+
+module Sym.Perm
+    (
+      module Sym.Internal.CLongArray
+    , Perm
+    , emptyperm
+    , one
+    , idperm
+    , ebb
+    , mkPerm
+    , rank
+    , unrank
+    , perms
+    ) where
+
+import Data.List
+import Sym.Internal.CLongArray
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+-- | A permutation is just a 'CLongArray'. By convention a permutation
+-- of size @n@ is understood to be a permutation of @[0..n-1]@.
+type Perm = CLongArray
+
+-- | The unique permutation length zero.
+emptyperm :: Perm
+emptyperm = fromList []
+
+-- | The unique permutation length one.
+one :: Perm
+one = fromList [0]
+
+-- | The identity permutation.
+idperm :: Int -> Perm
+idperm n = fromList [0..n-1]
+
+-- | The reverse of the identity permutation.
+ebb :: Int -> Perm
+ebb n = fromList [n-1,n-2..0]
+
+-- | Construct a permutation from a list of elements. As opposed to
+-- 'fromList' this is a safe function in the sense that the output of
+-- @mkPerm xs@ is guaranteed to be a permutation of @[0..length xs-1]@.
+-- E.g., @mkPerm \"baxa\" == fromList [2,0,3,1]@.
+mkPerm :: Ord a => [a] -> Perm
+mkPerm xs =
+    let sti ys = map snd . sort $ zip ys [ 0::Int .. ]
+    in fromList $ (sti . sti) xs
+
+foreign import ccall unsafe "rank.h rank" c_rank
+    :: Ptr CLong -> CLong -> IO CDouble
+
+-- | The rank of the given permutation, where the rank is defined as
+-- in [W. Myrvold and F. Ruskey, Ranking and Unranking Permutations in
+-- Linear Time, Information Processing Letters, 79 (2001) 281-284].
+rank :: Perm -> Integer
+rank w =
+    let n = fromIntegral (size w)
+    in truncate . unsafeDupablePerformIO . unsafeWith w $ flip c_rank n
+{-# INLINE rank #-}
+
+foreign import ccall unsafe "rank.h unrank" c_unrank
+    :: Ptr CLong -> CLong -> CDouble -> IO ()
+
+-- | The permutation of size @n@ whose rank is @r@, where the rank
+-- is defined as in [W. Myrvold and F. Ruskey, Ranking and Unranking
+-- Permutations in Linear Time, Information Processing Letters, 79
+-- (2001) 281-284].
+unrank :: Int -> Integer -> Perm
+unrank n r =
+    unsafeDupablePerformIO . unsafeNew n $ \ptr ->
+        c_unrank ptr (fromIntegral n) (fromIntegral r)
+{-# INLINE unrank #-}
+
+-- | All permutations of a given size.
+perms :: Int -> [Perm]
+perms n = map (unrank n) [0..nFac-1] where nFac = product [1..toInteger n]
diff --git a/Sym/Perm/Bijection.hs b/Sym/Perm/Bijection.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Bijection.hs
@@ -0,0 +1,39 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.Bijection
+    (
+      simionSchmidt
+    , simionSchmidt'
+    ) where
+
+import Sym.Perm
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+foreign import ccall unsafe "bij.h simion_schmidt" c_simion_schmidt
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+foreign import ccall unsafe "bij.h simion_schmidt_inverse" c_simion_schmidt'
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+marshal :: (Ptr CLong -> Ptr CLong -> CLong -> IO ()) -> Perm -> Perm
+marshal bij w =
+    unsafeDupablePerformIO . unsafeWith w $ \p -> do
+      let n = size w
+      unsafeNew n $ \q -> bij q p (fromIntegral n)
+{-# INLINE marshal #-}
+
+-- | The Simion-Schmidt bijection from Av(123) onto Av(132).
+simionSchmidt :: Perm -> Perm
+simionSchmidt = marshal c_simion_schmidt
+
+-- | The inverse of the Simion-Schmidt bijection. It is a function
+-- from Av(132) to Av(123).
+simionSchmidt' :: Perm -> Perm
+simionSchmidt' = marshal c_simion_schmidt'
diff --git a/Sym/Perm/Class.hs b/Sym/Perm/Class.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Class.hs
@@ -0,0 +1,187 @@
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.Class
+    (
+      inc
+    , dec
+    , av1
+    , av12
+    , av21
+    , av123
+    , av132
+    , av213
+    , av231
+    , av312
+    , av321
+    , av1243
+    , av1324
+    , av2134
+    , av
+    , vee
+    , caret
+    , gt
+    , lt
+    , wedges
+    , separables
+    , kLayered
+    , layered
+    , kFibonacci
+    , fibonacci
+    ) where
+
+import Sym.Internal.Util
+import Sym.Perm
+import Sym.Perm.Bijection
+import Sym.Perm.Constructions
+import Sym.Perm.Pattern
+import qualified Sym.Perm.D8 as D8
+
+-- | The class of increasing permutations.
+inc :: Int -> [Perm]
+inc = av21
+
+-- | The class of decreasing permutations.
+dec :: Int -> [Perm]
+dec = av12
+
+-- | Av(1)
+av1 :: Int -> [Perm]
+av1 0 = [emptyperm]
+av1 _ = []
+
+-- | Av(12)
+av12 :: Int -> [Perm]
+av12 n = [ebb n]
+
+-- | Av(21)
+av21 :: Int -> [Perm]
+av21 n = [idperm n]
+
+-- | Av(123)
+av123 :: Int -> [Perm]
+av123 = map simionSchmidt' . av132
+
+-- | Av(132)
+av132 :: Int -> [Perm]
+av132 = map D8.reverse . av231
+
+-- | Av(213)
+av213 :: Int -> [Perm]
+av213 = map D8.complement . av231
+
+-- | Av(231); also know as the stack sortable permutations.
+av231 :: Int -> [Perm]
+av231 0 = [emptyperm]
+av231 n = do
+  k <- [0..n-1]
+  s <- streamAv231 !! k
+  t <- streamAv231 !! (n-k-1)
+  return $ s /+/ (one \-\ t)
+
+streamAv231 :: [[Perm]]
+streamAv231 = map av231 [0..]
+
+-- | Av(312)
+av312 :: Int -> [Perm]
+av312 = map D8.inverse . av231
+
+-- | Av(321)
+av321 :: Int -> [Perm]
+av321 = map D8.complement . av123
+
+-- | Av(1243)
+av1243 :: Int -> [Perm]
+av1243 n = avoiders [fromList [0,1,3,2]] (perms n)
+
+-- | Av(1324)
+av1324 :: Int -> [Perm]
+av1324 n = avoiders [fromList [0,2,1,3]] (perms n)
+
+-- | Av(2134)
+av2134 :: Int -> [Perm]
+av2134 n = avoiders [fromList [1,0,2,3]] (perms n)
+
+-- | Av(s) where s is a string of one or more patterns, using space as a
+-- seperator.
+av :: String -> Int -> [Perm]
+av s = avoiders (map mkPerm (words s)) . perms
+
+-- | The V-class is Av(132, 231). It is so named because the diagram of
+-- a typical permutation in this class is shaped like a V.
+vee :: Int -> [Perm]
+vee = (streamVee !!)
+
+streamVee :: [[Perm]]
+streamVee = [emptyperm] : [one] : zipWith (++) vee_n n_vee
+    where
+      n_vee = (map.map) (one \-\) ws
+      vee_n = (map.map) (/+/ one) ws
+      ws    = tail streamVee
+
+-- | The ∧-class is Av(213, 312). It is so named because the diagram of
+-- a typical permutation in this class is shaped like a ∧.
+caret :: Int -> [Perm]
+caret = map D8.complement . vee
+
+-- | The >-class is Av(132, 312). It is so named because the diagram of
+-- a typical permutation in this class is shaped like a >.
+gt :: Int -> [Perm]
+gt = map D8.rotate . vee
+
+-- | The <-class is Av(213, 231). It is so named because the diagram of
+-- a typical permutation in this class is shaped like a <.
+lt :: Int -> [Perm]
+lt = map D8.reverse . gt
+
+union :: [Int -> [Perm]] -> Int -> [Perm]
+union cs n = nubSort $ concat [ c n | c <- cs ]
+
+-- | The union of 'vee', 'caret', 'gt' and 'lt'.
+wedges :: Int -> [Perm]
+wedges = union [vee, caret, gt, lt]
+
+compositions :: Int -> Int -> [[Int]]
+compositions 0 0 = [[]]
+compositions 0 _ = []
+compositions _ 0 = []
+compositions k n = [1..n] >>= \i -> map (i:) (compositions (k-1) (n-i))
+
+boundedCompositions :: Int -> Int -> Int -> [[Int]]
+boundedCompositions _ 0 0 = [[]]
+boundedCompositions _ 0 _ = []
+boundedCompositions _ _ 0 = []
+boundedCompositions b k n = [1..b] >>= \i -> map (i:) (boundedCompositions b (k-1) (n-i))
+
+-- | The class of separable permutations; it is identical to Av(2413,3142).
+separables :: Int -> [Perm]
+separables 0 = [emptyperm]
+separables 1 = [one]
+separables n = pIndec n ++ mIndec n
+    where
+      comps  m = [2..m] >>= \k -> compositions k m
+      pIndec 0 = []
+      pIndec 1 = [one]
+      pIndec m = comps m >>= map skewSum . mapM (streamMIndec !!)
+      mIndec m = map D8.complement $ pIndec m
+      streamMIndec = map mIndec [0..]
+
+-- | The class of layered permutations with /k/ layers.
+kLayered :: Int -> Int -> [Perm]
+kLayered k = map (directSum . map ebb) . compositions k
+
+-- | The class of layered permutations.
+layered :: Int -> [Perm]
+layered n = [1..n] >>= flip kLayered n
+
+-- | The class of Fibonacci permutations with /k/ layers. A /Fibonacci permutation/
+-- is a layered permutation whose layers are all of size 1 or 2.
+kFibonacci :: Int -> Int -> [Perm]
+kFibonacci k = map (directSum . map ebb) . boundedCompositions 2 k
+
+-- | The class of Fibonacci permutations. A /Fibonacci permutation/ is a
+-- layered permutation whose layers are all of size 1 or 2.
+fibonacci :: Int -> [Perm]
+fibonacci n = [1..n] >>= flip kFibonacci n
diff --git a/Sym/Perm/Component.hs b/Sym/Perm/Component.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Component.hs
@@ -0,0 +1,76 @@
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Components of permutations.
+-- 
+
+module Sym.Perm.Component
+    (
+      components
+    , skewComponents
+    , leftMaxima
+    , leftMinima
+    , rightMaxima
+    , rightMinima
+    ) where
+
+import Foreign
+import System.IO.Unsafe
+import Sym.Perm
+import qualified Sym.Perm.D8 as D8
+
+-- Positions /i/ such that /max{ w[j] : j <= i } = i/. These positions
+-- mark the boundaries of components.
+comps :: Perm -> [Int]
+comps w = unsafeDupablePerformIO . unsafeWith w $ go [] 0 0
+    where
+      n = size w
+      go ks m i p
+        | i >= n = return (reverse ks)
+        | otherwise =
+            do y <- fromIntegral `fmap` peek p
+               let p'  = advancePtr p 1
+               let i'  = i+1
+               let m'  = if y > m then y else m
+               let ks' = if m' == i then i:ks else ks
+               go ks' m' i' p'
+
+-- | The list of (plus) components.
+components :: Perm -> [Perm]
+components w =
+    let ds = 0 : map (+1) (comps w)
+        ks = zipWith (-) (tail ds) ds
+        ws = unsafeSlice ks w
+    in zipWith (\d v -> imap (\_ x -> x - fromIntegral d) v) ds ws
+
+-- | The list of skew components, also called minus components.
+skewComponents :: Perm -> [Perm]
+skewComponents = map D8.complement . components . D8.complement
+
+records :: (a -> a -> Bool) -> [a] -> [a]
+records _ []     = []
+records f (x:xs) = recs [x] xs
+    where
+      recs rs@(r:_) (y:ys) = recs ((if f r y then y else r):rs) ys
+      recs rs       _      = rs
+
+-- | For each position, left-to-right, records the largest value seen
+-- thus far.
+leftMaxima :: Perm -> [Int]
+leftMaxima w = reverse $ records (<) (toList w)
+
+-- | For each position, left-to-right, records the smallest value seen
+-- thus far.
+leftMinima :: Perm -> [Int]
+leftMinima w = reverse $ records (>) (toList w)
+
+-- | For each position, /right-to-left/, records the largest value seen
+-- thus far.
+rightMaxima :: Perm -> [Int]
+rightMaxima w = records (<) (reverse (toList w))
+
+-- | For each position, /right-to-left/, records the smallest value seen
+-- thus far.
+rightMinima :: Perm -> [Int]
+rightMinima w = records (>) (reverse (toList w))
diff --git a/Sym/Perm/Constructions.hs b/Sym/Perm/Constructions.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Constructions.hs
@@ -0,0 +1,62 @@
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Sum, skew sum, etc
+-- 
+
+module Sym.Perm.Constructions
+    (
+      (/+/)
+    , (\-\)
+    , directSum
+    , skewSum
+    , inflate
+    ) where
+
+import Foreign
+import System.IO.Unsafe
+import Control.Monad
+import Sym.Perm
+import qualified Sym.Permgram as G
+import qualified Sym.Perm.D8 as D8
+
+infixl 6 /+/
+infixl 6 \-\
+
+-- | The /direct sum/ of two permutations.
+(/+/) :: Perm -> Perm -> Perm
+(/+/) u v =
+   let k  = size u
+       l  = size v
+       v' = imap (\_ x -> x + fromIntegral k) v
+   in unsafeDupablePerformIO . unsafeNew (k+l) $ \p ->
+       let q = advancePtr p k
+       in unsafeWith u  $ \uPtr ->
+          unsafeWith v' $ \vPtr -> do
+              copyArray p uPtr k
+              copyArray q vPtr l
+
+-- | The direct sum of a list of permutations.
+directSum :: [Perm] -> Perm
+directSum = foldr (/+/) emptyperm
+
+-- | The /skew sum/ of two permutations.
+(\-\) :: Perm -> Perm -> Perm
+(\-\) u v = D8.complement $ D8.complement u /+/ D8.complement v
+
+-- | The skew sum of a list of permutations.
+skewSum :: [Perm] -> Perm
+skewSum = foldr (\-\) emptyperm
+
+-- | @inflate w vs@ is the /inflation/ of @w@ by @vs@. It is the
+-- permutation of length @sum (map size vs)@ obtained by replacing
+-- each entry @w!i@ by an interval that is order isomorphic to @vs!i@
+-- in such a way that the intervals are order isomorphic to @w@. In
+-- particular,
+-- 
+-- > u /+/ v == inflate (mkPerm "12") [u,v]
+-- > u \-\ v == inflate (mkPerm "21") [u,v]
+-- 
+inflate :: Perm -> [Perm] -> Perm
+inflate w = G.perm . join . G.permgram w . map (`G.permgram` [()])
diff --git a/Sym/Perm/D8.hs b/Sym/Perm/D8.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/D8.hs
@@ -0,0 +1,156 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.D8
+    (
+    -- * The group elements
+      r0, r1, r2, r3
+    , s0, s1, s2, s3
+
+    -- * D8, the klein four-group, and orbits
+    , d8
+    , klein4
+    , orbit
+    , symmetryClasses
+    , d8Classes
+    , klein4Classes
+
+    -- * Aliases
+    , rotate
+    , complement
+    , reverse
+    , inverse
+    ) where
+
+import Prelude           hiding (reverse)
+import Data.List         hiding (reverse)
+import Sym.Internal.Util
+import Sym.Perm
+import Foreign           hiding (complement, rotate)
+import Foreign.C.Types
+import System.IO.Unsafe
+
+
+-- The group elements
+-- ------------------
+
+-- | Ration by 0 degrees, i.e. the identity map.
+r0 :: Perm -> Perm
+r0 w = w
+
+-- | Ration by 90 degrees clockwise.
+r1 :: Perm -> Perm
+r1 = s2 . s1
+
+-- | Ration by 2*90 = 180 degrees clockwise.
+r2 :: Perm -> Perm
+r2 = r1 . r1
+
+-- | Ration by 3*90 = 270 degrees clockwise.
+r3 :: Perm -> Perm
+r3 = r2 . r1
+
+-- | Reflection through a horizontal axis (also called 'complement').
+s0 :: Perm -> Perm
+s0 = complement
+
+-- | Reflection through a vertical axis (also called 'reverse').
+s1 :: Perm -> Perm
+s1 = reverse
+
+-- | Reflection through the main diagonal (also called 'inverse').
+s2 :: Perm -> Perm
+s2 = inverse
+
+-- | Reflection through the anti-diagonal.
+s3 :: Perm -> Perm
+s3 = s1 . r1
+
+
+-- D8, the klein four-group, and orbits
+-- ------------------------------------
+
+-- | The dihedral group of order 8 (the symmetries of a square); that is,
+-- 
+-- > d8 = [r0, r1, r2, r3, s0, s1, s2, s3]
+-- 
+d8 :: [Perm -> Perm]
+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 :: [Perm -> Perm]
+klein4 = [r0, r2, s0, s1]
+
+-- | @orbit fs x@ is the orbit of @x@ under the /group/ of function @fs@. E.g.,
+-- 
+-- > orbit klein4 "2314" == ["1423","2314","3241","4132"]
+-- 
+orbit :: [Perm -> Perm] -> Perm -> [Perm]
+orbit fs x = nubSort [ f x | f <- fs ]
+
+-- | @symmetryClasses fs xs@ is the list of equivalence classes under
+-- the action of the /group/ of functions @fs@.
+symmetryClasses :: [Perm -> Perm] -> [Perm] -> [[Perm]]
+symmetryClasses _  [] = []
+symmetryClasses fs xs@(x:xt) = insert orb $ symmetryClasses fs ys
+    where
+      orb = [ w | w <- orbit fs x, w `elem` xs ]
+      ys  = [ y | y <- xt, y `notElem` orb ]
+
+-- | Symmetry classes with respect to D8.
+d8Classes :: [Perm] -> [[Perm]]
+d8Classes = symmetryClasses d8
+
+-- | Symmetry classes with respect to Klein4
+klein4Classes :: [Perm] -> [[Perm]]
+klein4Classes = symmetryClasses klein4
+
+
+-- Aliases
+-- -------
+
+marshal :: (Ptr CLong -> Ptr CLong -> CLong -> IO ()) -> Perm -> Perm
+marshal op w =
+    unsafeDupablePerformIO . unsafeWith w $ \p -> do
+      let n = size w
+      unsafeNew n $ \q -> op q p (fromIntegral n)
+{-# INLINE marshal #-}
+
+foreign import ccall unsafe "d8.h inverse" c_inverse
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+-- | The group theoretical inverse: if @v = inverse u@ then
+-- @v \`at\` (u \`at\` i) = i@.
+inverse :: Perm -> Perm
+inverse = marshal c_inverse
+{-# INLINE inverse #-}
+
+foreign import ccall unsafe "d8.h reverse" c_reverse
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+-- | The reverse of the given permutation: if @v = reverse u@ then
+-- @v \`at\` i = u \`at\` (n-1-i)@.
+reverse :: Perm -> Perm
+reverse = marshal c_reverse
+{-# INLINE reverse #-}
+
+foreign import ccall unsafe "d8.h complement" c_complement
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+-- | The complement of the given permutation: if @v = complement u@ then
+-- @v \`at\` i = n - 1 - u \`at\` i@.
+complement :: Perm -> Perm
+complement = marshal c_complement
+{-# INLINE complement #-}
+
+-- | @rotate = r1 = inverse . reverse@
+rotate :: Perm -> Perm
+rotate = r1
diff --git a/Sym/Perm/Group.hs b/Sym/Perm/Group.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Group.hs
@@ -0,0 +1,47 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.Group
+    (
+      compose
+    , act
+    ) where
+
+import Sym.Perm
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+marshal :: (Ptr CLong -> Ptr CLong -> CLong -> Ptr CLong -> CLong -> IO ())
+        -> Perm -> Perm -> Perm
+marshal op u v =
+    unsafeDupablePerformIO $
+    unsafeWith u $ \u' ->
+    unsafeWith v $ \v' -> do
+      let k = size u
+      let n = size v
+      let m = max k n
+      unsafeNew m $ \p -> op p u' (fromIntegral k) v' (fromIntegral n)
+{-# INLINE marshal #-}
+
+foreign import ccall unsafe "group.h compose" c_compose
+    :: Ptr CLong -> Ptr CLong -> CLong -> Ptr CLong -> CLong -> IO ()
+
+-- | The product/composition of @u@ and @v@: if @w = u `compose` v@
+-- then @w `at ` i = v \`at\` (u \`at\` i)@.
+compose :: Perm -> Perm -> Perm
+compose = marshal c_compose
+{-# INLINE compose #-}
+
+foreign import ccall unsafe "group.h act" c_act
+    :: Ptr CLong -> Ptr CLong -> CLong -> Ptr CLong -> CLong -> IO ()
+
+-- | The (left) group action of Perm on itself: if @w = u `act` v@
+-- then @w `at ` (u \`at\` i) = v \`at\` i@.
+act :: Perm -> Perm -> Perm
+act = marshal c_act
+{-# INLINE act #-}
diff --git a/Sym/Perm/MeshPattern.hs b/Sym/Perm/MeshPattern.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/MeshPattern.hs
@@ -0,0 +1,151 @@
+-- |
+-- Copyright   : Anders Claesson 2014
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+-- TODO: Generalize interface and share with Sym.Perm.Pattern
+
+module Sym.Perm.MeshPattern
+    ( MeshPattern (..)
+    , Mesh
+    , Box
+    , Point
+    , mkPattern
+    , pattern
+    , mesh
+    , cols
+    , rows
+    , col
+    , row
+    , box
+    , copiesOf
+    , contains
+    , avoids
+    , avoidsAll
+    , avoiders
+    , kVincular
+    , vincular
+    , bivincular
+    , meshPatterns
+    ) where
+
+import Data.List hiding (union)
+import Sym.Internal.Size
+import Sym.Perm
+import Sym.Internal.SubSeq
+import Data.Set (Set)
+import qualified Data.Set as Set
+import Sym.Internal.Util
+
+type Point = (Int, Int)
+type Box   = (Int, Int)
+type Mesh  = Set Box
+type PermTwoLine = [Point]
+
+data MeshPattern = MP
+    { getPerm :: Perm
+    , getMesh :: Mesh
+    } deriving (Show, Eq, Ord)
+
+instance Size MeshPattern where
+    size = size . getPerm
+
+mkPattern :: Ord a => [a] -> MeshPattern
+mkPattern w = MP (mkPerm w) Set.empty
+
+pattern :: Perm -> MeshPattern
+pattern w = MP w Set.empty
+
+mesh :: [Box] -> MeshPattern -> MeshPattern
+mesh r (MP w s) = MP w . Set.union s $ Set.fromList r
+
+cols :: [Int] -> MeshPattern -> MeshPattern
+cols xs p@(MP w _) = mesh [ (x,y) | y <- [0..size w], x <- xs ] p
+
+rows :: [Int] -> MeshPattern -> MeshPattern
+rows ys p@(MP w _) = mesh [ (x,y) | x <- [0..size w], y <- ys ] p
+
+col :: Int -> MeshPattern -> MeshPattern
+col y = cols [y]
+
+row :: Int -> MeshPattern -> MeshPattern
+row x = rows [x]
+
+box :: Box -> MeshPattern -> MeshPattern
+box xy = mesh [xy]
+
+kVincular :: Int -> Perm -> [MeshPattern]
+kVincular k w = (flip cols (pattern w) . toList) `fmap` ((1+size w) `choose` k)
+-- kVincular k w = (\xs -> cols (toList xs) (pattern w)) `fmap` ((1+size w) `choose` k)
+
+vincular :: Perm -> [MeshPattern]
+vincular w = [0..1+size w] >>= flip kVincular w
+
+bivincular :: Perm -> [MeshPattern]
+bivincular w =
+    [ foldr ((.) . either col row) id c $ pattern w | c <- choices ]
+  where
+    choices = powerset' $ [0..size w] >>= \z -> [Left z, Right z]
+    powerset' = fmap Set.toList . powerset . Set.fromList
+
+fullMesh :: Int -> Mesh
+fullMesh n = let zs = [0..n] in Set.fromList [ (x,y) | x <- zs, y <- zs ]
+
+meshPatterns :: Perm -> [MeshPattern]
+meshPatterns w = [ MP w r | r <- powerset (fullMesh (size w)) ]
+
+match' :: MeshPattern -> PermTwoLine -> PermTwoLine -> Bool
+match' (MP u r) v w =
+    and $ (u2==v2) : [ not $ f i j x y | (i,j) <- Set.toList r, (x,y) <- w ]
+  where
+    (v1, v2) = unzip v
+    m  = 1 + length w
+    xs = 0 : v1 ++ [m]
+    ys = 0 : sort v2 ++ [m]
+    u2 = map ((ys!!) . (+1)) (toList u)
+    f i j x y = xs!!i < x && x < xs!!(i+1) && ys!!j < y && y < ys!!(j+1)
+
+-- | @match p w m@ determines whether the subword in @w@ specified by
+-- @m@ is an occurrence of @p@.
+match :: MeshPattern -> Perm -> SubSeq -> Bool
+match p w m = match' p v w'
+  where
+    w' = twoLine w
+    v  = [ pt | pt@(x,_) <- w', x-1 `elem` toList m ]
+
+twoLine :: Perm -> PermTwoLine
+twoLine = zip [1..] . map (+1) . toList
+
+-- | @copiesOf p w@ is the list of sets that represent copies of @p@ in @w@.
+copiesOf :: MeshPattern -> Perm -> [SubSeq]
+copiesOf p w = filter (match p w) $ size w `choose` size p
+{-# INLINE copiesOf #-}
+
+-- | @w `contains` p@ is a predicate determining if @w@ contains the pattern @p@.
+contains :: Perm -> MeshPattern -> Bool
+w `contains` p = not $ w `avoids` p
+
+-- | @w `avoids` p@ is a predicate determining if @w@ avoids the pattern @p@.
+avoids :: Perm -> MeshPattern -> Bool
+w `avoids` p = null $ copiesOf p w
+
+-- | @w `avoidsAll` ps@ is a predicate determining if @w@ avoids the patterns @ps@.
+avoidsAll :: Perm -> [MeshPattern] -> Bool
+w `avoidsAll` ps = all (w `avoids`) ps
+
+-- | @avoiders ps ws@ is the list of permutations in @ws@ avoiding the
+-- patterns in @ps@.
+avoiders :: [MeshPattern] -> [Perm] -> [Perm]
+avoiders ps ws = foldl (flip avoiders1) ws ps
+
+-- @avoiders1 p ws@ is the list of permutations in @ws@ avoiding the
+-- pattern @p@.
+avoiders1 :: MeshPattern -> [Perm] -> [Perm]
+avoiders1 _ [] = []
+avoiders1 q vs@(v:_) = filter avoids_q us ++ filter (`avoids` q) ws
+    where
+      n = size v
+      k = size q
+      (us, ws) = span (\u -> size u == n) vs
+      xs = n `choose` k
+      avoids_q u = not $ any (match q u) xs
diff --git a/Sym/Perm/Meta.hs b/Sym/Perm/Meta.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Meta.hs
@@ -0,0 +1,18 @@
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- A meta module collecting all Perm-modules, except those that are best
+-- imported \"qualified\".
+
+module Sym.Perm.Meta (module P) where
+
+import Sym.Perm                 as P
+import Sym.Perm.Class           as P
+import Sym.Perm.Component       as P
+import Sym.Perm.Constructions   as P
+import Sym.Perm.Bijection       as P
+import Sym.Perm.Group           as P
+import Sym.Perm.Pattern         as P
+import Sym.Perm.Simple          as P
+import Sym.Perm.Sort            as P
diff --git a/Sym/Perm/Pattern.hs b/Sym/Perm/Pattern.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Pattern.hs
@@ -0,0 +1,109 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.Pattern
+    (
+      Pattern
+    , SubSeq
+    , ordiso
+    , choose
+    , copiesOf
+    , contains
+    , avoids
+    , avoidsAll
+    , avoiders
+    , minima
+    , maxima
+    , coeff
+    ) where
+
+import Sym.Perm (Perm, perms)
+import Sym.Internal.SubSeq
+import Sym.Internal.Util (nubSort)
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+-- | Pattern is just an alias for permutation.
+type Pattern = Perm
+
+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 :: Pattern -> Perm -> SubSeq -> Bool
+ordiso u v m =
+    let k = fromIntegral (size u)
+    in unsafeDupablePerformIO $
+       unsafeWith u $ \u' ->
+       unsafeWith v $ \v' ->
+       unsafeWith m $ \m' ->
+           return . toBool $ c_ordiso u' v' m' k
+{-# INLINE ordiso #-}
+
+-- | @copiesOf p w@ is the list of sets that represent copies of @p@ in @w@.
+copiesOf :: Pattern -> Perm -> [SubSeq]
+copiesOf p w = filter (ordiso p w) $ size w `choose` size p
+{-# INLINE copiesOf #-}
+
+-- | @w `contains` p@ is a predicate determining if @w@ contains the pattern @p@.
+contains :: Perm -> Pattern -> Bool
+w `contains` p = not $ w `avoids` p
+
+-- | @w `avoids` p@ is a predicate determining if @w@ avoids the pattern @p@.
+avoids :: Perm -> Pattern -> Bool
+w `avoids` p = null $ copiesOf p w
+
+-- | @w `avoidsAll` ps@ is a predicate determining if @w@ avoids the patterns @ps@.
+avoidsAll :: Perm -> [Pattern] -> Bool
+w `avoidsAll` ps = all (w `avoids`) ps
+
+-- | @avoiders ps ws@ is the list of permutations in @ws@ avoiding the
+-- patterns in @ps@.
+avoiders :: [Pattern] -> [Perm] -> [Perm]
+avoiders ps ws = foldl (flip avoiders1) ws ps
+
+-- @avoiders1 p ws@ is the list of permutations in @ws@ avoiding the
+-- pattern @p@.
+avoiders1 :: Pattern -> [Perm] -> [Perm]
+avoiders1 _ [] = []
+avoiders1 q vs@(v:_) = filter avoids_q us ++ filter (`avoids` q) ws
+    where
+      n = size v
+      k = size q
+      (us, ws) = span (\u -> size u == n) vs
+      xs = n `choose` k
+      avoids_q u = not $ any (ordiso q u) xs
+
+-- | The set of minimal elements with respect to containment.  FIX: Poor
+-- implementation
+minima :: [Pattern] -> [Pattern]
+minima [] = []
+minima ws = let (v:vs) = nubSort ws
+            in v : minima (avoiders [v] vs)
+
+-- | The set of maximal elements with respect to containment. FIX: Poor
+-- implementation
+maxima :: [Pattern] -> [Pattern]
+maxima [] = []
+maxima ws = let (v:vs) = reverse $ nubSort ws
+            in v : maxima (filter (avoids v) vs)
+
+-- | @coeff f v@ is the coefficient of @v@ when expanding the
+-- permutation statistic @f@ as a sum of permutations/patterns. See
+-- Petter Brändén and Anders Claesson: Mesh patterns and the expansion
+-- of permutation statistics as sums of permutation patterns, The
+-- Electronic Journal of Combinatorics 18(2) (2011),
+-- <http://www.combinatorics.org/ojs/index.php/eljc/article/view/v18i2p5>.
+coeff :: (Pattern -> Int) -> Pattern -> Int
+coeff f v = f v + sum [ (-1)^(k - j) * c * f u |
+                        j <- [0 .. k-1]
+                      , u <- perms j
+                      , let c = length $ copiesOf u v
+                      , c > 0
+                      ] where k = size v
diff --git a/Sym/Perm/SSYT.hs b/Sym/Perm/SSYT.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/SSYT.hs
@@ -0,0 +1,156 @@
+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Data types for Semistandard Young Tableaux (SSYT) and functions for
+-- converting between (generalized) permutataions and SSYT. In other
+-- words, this module implements the Robinson-Schensted-Knuth (RSK)
+-- correspondence.
+
+module Sym.Perm.SSYT
+    (
+      GeneralizedPerm
+    , Entry
+    , SSYT
+    , SSYTPair (..)
+    , Shape (..)
+    , empty
+    , null
+    , display
+    , fromPerm
+    , fromGeneralizedPerm
+    , toPerm
+    , toGeneralizedPerm
+    ) where
+
+import Prelude    hiding (null)
+import Data.List  hiding (null)
+import Sym.Perm
+
+type Row = Int
+
+-- | An entry is a non-negative integer
+type Entry = Int
+
+-- | A /Generalized Permutation/ is a lexicographically sorted list of
+-- pairs of non-negative integers.
+type GeneralizedPerm = [(Int, Int)]
+
+-- | A /Semistandard Young Tableau (SSYT)/: the entries weakly increase
+-- along each row and strictly increase down each column.
+type SSYT = [[Entry]]
+
+-- | A pair of Semistandard Young Tableaux.
+data SSYTPair = SSYTPair { insertionTableau :: SSYT
+                         , recordingTableau :: SSYT
+                         } deriving Eq
+
+class Shape a where
+    shape :: a -> [Int]
+
+instance Shape SSYT where
+    shape = map length
+
+instance Shape SSYTPair where
+    shape = shape . recordingTableau
+
+-- | A pair of empty Young tableaux.
+empty :: SSYTPair
+empty = SSYTPair [] []
+
+-- | Check if a given pair of Young tableaux are empty.
+null :: SSYTPair -> Bool
+null pq = pq == empty
+
+instance Show SSYTPair where
+    show (SSYTPair p q) = unwords ["SSYTPair", show p, show q]
+
+-- | Produce a string for pretty printing SSYT pairs.
+display :: SSYTPair -> String
+display pq@(SSYTPair p q)
+    | null pq   = "[] []"
+    | otherwise = intercalate "\n" $ zipWith (++) (pad p') q'
+       where
+         p'@(r:_) = map show p
+         q'       = map show q
+         pad      = map $ \s -> take (1+length r) (s ++ repeat ' ')
+
+-- Inserts Entry into a given Tableau returning the resulting Tableau
+-- and the row where the Tableau was extended with a new box.
+insertP :: SSYT -> Entry -> (SSYT, Row)
+insertP []     k = ([[k]], 1)
+insertP (r:rs) k =
+    let (smaller, larger) = span (<=k) r
+    in case larger of
+         []   -> ((r++[k]):rs, 1)
+         c:cs -> let (rs', i) = insertP rs c
+                 in ((smaller ++ k:cs) : rs', i+1)
+
+-- Given (i,j), inserts j at the end of row i in the given Tableau.
+insertQ :: SSYT -> Row -> Entry -> SSYT
+insertQ []     _ j = [[j]]
+insertQ (r:rs) 1 j = (r ++ [j]) : rs
+insertQ (r:rs) i j = r : insertQ rs (i-1) j
+
+-- Given (i,j) and pair of tableaux (p,q) of the same shape, inserts i
+-- into p and j into q so that the resulting pair of tableaux (p',q')
+-- still have the same shape.
+insertPQ :: SSYTPair -> (Entry, Entry) -> SSYTPair
+insertPQ (SSYTPair p q) (i,j) =
+    let (p',k) = insertP p j in SSYTPair p' (insertQ q k i)
+
+trim :: SSYT -> SSYT
+trim = takeWhile (/=[])
+
+-- The inverse of insertP
+removeP :: SSYT -> Row -> (SSYT, Entry)
+removeP p k = (trim $ reverse vs ++ [init t] ++ p2, e)
+    where
+      (p1, p2) = splitAt (k+1) p
+      (t : ts) = reverse p1 -- t is the k-th row (counting from 0)
+      (vs,  e) = unbump (last t) ts
+      unbump x [] = ([], x)
+      unbump x (r:rs) =
+          let (r1, r2) = span (<x) r
+              (us,  y) = unbump (last r1) rs
+          in ((init r1 ++ x:r2) : us, y)
+
+-- The inverse of insertQ
+removeQ :: SSYT -> (SSYT, Row, Entry)
+removeQ q = (trim q', k, e)
+    where
+      -- The last element and the length of a given row:
+      f = foldl (\(_,n) x -> (x,n+1)) (0,0) :: [Int] -> (Int, Int)
+      -- Equal elements of Q are inserted left-to-right, allowing us to
+      -- know which element, e, was the last to be inserted:
+      ((e, _), k) = maximum $ zip (map f q) [0..]
+      -- Remove e from Q:
+      q' = [ if i == k then init r else r | (r,i) <- zip q [0..] ]
+
+-- The inverse of insertPQ
+removePQ :: SSYTPair -> (SSYTPair, (Entry, Entry))
+removePQ (SSYTPair p q) = (SSYTPair p' q', (e1, e2))
+    where
+      (q', k, e1) = removeQ q
+      (p', e2)    = removeP p k
+
+-- | The Robinson-Schensted-Knuth (RSK) algorithm.
+fromGeneralizedPerm :: GeneralizedPerm -> SSYTPair
+fromGeneralizedPerm = foldl insertPQ empty
+
+-- | The Robinson-Schensted algorithm.
+fromPerm :: Perm -> SSYTPair
+fromPerm = fromGeneralizedPerm . zip [0..] . toList
+
+-- | The inverse of the Robinson-Schensted-Knuth algorithm.
+toGeneralizedPerm :: SSYTPair -> GeneralizedPerm
+toGeneralizedPerm = go []
+    where
+      go ijs pq | null pq   = ijs
+                | otherwise = let (rs,ij) = removePQ pq in go (ij:ijs) rs
+
+-- | The inverse of the Robinson-Schensted algorithm.
+toPerm :: SSYTPair -> Perm
+toPerm = fromList . map snd . toGeneralizedPerm
diff --git a/Sym/Perm/Simple.hs b/Sym/Perm/Simple.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Simple.hs
@@ -0,0 +1,25 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.Simple
+    (
+     simple
+    ) where
+
+import Sym.Perm
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+foreign import ccall unsafe "simple.h simple" c_simple
+    :: Ptr CLong -> CLong -> CInt
+
+-- | Is the permutation simple?
+simple :: Perm -> Bool
+simple w = toBool . unsafeDupablePerformIO $
+    let n = fromIntegral (size w)
+    in unsafeWith w $ \ptr -> return $ c_simple ptr n
diff --git a/Sym/Perm/Sort.hs b/Sym/Perm/Sort.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Sort.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+
+module Sym.Perm.Sort
+    (
+      stackSort
+    , bubbleSort
+    ) where
+
+import Sym.Perm
+import Foreign
+import Foreign.C.Types
+import System.IO.Unsafe
+
+foreign import ccall unsafe "sortop.h stacksort" c_stacksort
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+foreign import ccall unsafe "sortop.h bubblesort" c_bubblesort
+    :: Ptr CLong -> Ptr CLong -> CLong -> IO ()
+
+marshal :: (Ptr CLong -> Ptr CLong -> CLong -> IO ()) -> Perm -> Perm
+marshal op w =
+    unsafeDupablePerformIO . unsafeWith w $ \p -> do
+      let n = size w
+      unsafeNew n $ \q -> op q p (fromIntegral n)
+{-# INLINE marshal #-}
+
+-- | One pass of stack-sort.
+stackSort :: Perm -> Perm
+stackSort = marshal c_stacksort
+{-# INLINE stackSort #-}
+
+-- | One pass of bubble-sort.
+bubbleSort :: Perm -> Perm
+bubbleSort = marshal c_bubblesort
+{-# INLINE bubbleSort #-}
diff --git a/Sym/Perm/Stat.hs b/Sym/Perm/Stat.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Perm/Stat.hs
@@ -0,0 +1,289 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Common permutation statistics. To avoid name clashes this module is
+-- best imported @qualified@; e.g.
+-- 
+-- > import qualified Sym.Perm.Stat as S
+-- 
+
+module Sym.Perm.Stat 
+    (
+      asc         -- ascents
+    , des         -- descents
+    , exc         -- excedances
+    , fp          -- fixed points
+    , sfp         -- strong fixed points
+    , cyc         -- cycles
+    , inv         -- inversions
+    , maj         -- the major index
+    , comaj       -- the co-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
+    , scomp       -- skew components
+    , ep          -- rank a la Elizalde & Pak
+    , dim         -- dimension
+    , asc0        -- small ascents
+    , des0        -- small descents
+    , lis         -- longest increasing subsequence
+    , lds         -- longest decreasing subsequence
+--    , shad        -- shadow
+    ) where
+
+import Prelude hiding (head, last)
+import qualified Prelude
+import Sym.Perm
+import qualified Sym.Perm.SSYT as Y
+import qualified Sym.Perm.D8 as D8
+import Foreign.Ptr
+import Foreign.C.Types
+import System.IO.Unsafe
+
+marshal :: (Ptr CLong -> CLong -> CLong) -> Perm -> Int
+marshal f w =
+    fromIntegral . unsafeDupablePerformIO . unsafeWith w $ \p ->
+        return $ f p (fromIntegral (size w))
+{-# INLINE marshal #-}
+
+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 sfp" c_sfp
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h cyc" c_cyc
+    :: 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 comaj" c_comaj
+    :: 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
+
+foreign import ccall unsafe "stat.h dim" c_dim
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h asc0" c_asc0
+    :: Ptr CLong -> CLong -> CLong
+
+foreign import ccall unsafe "stat.h des0" c_des0
+    :: Ptr CLong -> CLong -> CLong
+
+-- | The number of ascents. An /ascent/ in @w@ is an index @i@ such
+-- that @w[i] \< w[i+1]@.
+asc :: Perm -> Int
+asc = marshal c_asc
+
+-- | The number of descents. A /descent/ in @w@ is an index @i@ such
+-- that @w[i] > w[i+1]@.
+des :: Perm -> Int
+des = marshal c_des
+
+-- | The number of /excedances/: positions @i@ such that @w[i] > i@.
+exc :: Perm -> Int
+exc = marshal c_exc
+
+-- | The number of /fixed points/: positions @i@ such that @w[i] == i@.
+fp :: Perm -> Int
+fp = marshal c_fp
+
+-- | The number of /strong fixed points/ (also called splitters):
+-- positions @i@ such that @w[j] \< i@ for @j \< i@ and @w[j] \> i@ for @j \> i@.
+sfp :: Perm -> Int
+sfp = marshal c_sfp
+
+-- | The number of /cycles/:
+-- orbits of the permutation when viewed as a function.
+cyc :: Perm -> Int
+cyc = marshal c_cyc
+
+-- | The number of /inversions/:
+-- pairs @\(i,j\)@ such that @i \< j@ and @w[i] > w[j]@.
+inv :: Perm -> Int
+inv = marshal c_inv
+
+-- | /The major index/ is the sum of descents.
+maj :: Perm -> Int
+maj = marshal c_maj
+
+-- | /The co-major index/ is the sum of descents.
+comaj :: Perm -> Int
+comaj = marshal c_comaj
+
+-- | The number of /peaks/:
+-- positions @i@ such that @w[i-1] \< w[i]@ and @w[i] \> w[i+1]@.
+peak :: Perm -> Int
+peak = marshal c_peak
+
+-- | The number of /valleys/:
+-- positions @i@ such that @w[i-1] \> w[i]@ and @w[i] \< w[i+1]@.
+vall :: Perm -> Int
+vall = marshal c_vall
+
+-- | The number of /double ascents/:
+-- positions @i@ such that @w[i-1] \<  w[i] \< w[i+1]@.
+dasc :: Perm -> Int
+dasc = marshal c_dasc
+
+-- | The number of /double descents/:
+-- positions @i@ such that @w[i-1] \>  w[i] \> w[i+1]@.
+ddes :: Perm -> Int
+ddes = marshal c_ddes
+
+-- | The number of /left-to-right minima/:
+-- positions @i@ such that @w[i] \< w[j]@ for all @j \< i@.
+lmin :: Perm -> Int
+lmin = marshal c_lmin
+
+-- | The number of /left-to-right maxima/:
+-- positions @i@ such that @w[i] \> w[j]@ for all @j \< i@.
+lmax :: Perm -> Int
+lmax = marshal c_lmax
+
+-- | The number of /right-to-left minima/:
+-- positions @i@ such that @w[i] \< w[j]@ for all @j \> i@.
+rmin :: Perm -> Int
+rmin = lmin . D8.reverse
+
+-- | The number of /right-to-left maxima/:
+-- positions @i@ such that @w[i] \> w[j]@ for all @j \> i@.
+rmax :: Perm -> Int
+rmax = lmax . D8.reverse
+
+-- | The first (left-most) element in the standardization. E.g.,
+-- @head \"231\" = head (fromList [1,2,0]) = 1@.
+head :: Perm -> Int
+head w | size w > 0 = w `unsafeAt` 0
+       | otherwise  = 0
+
+-- | The last (right-most) element in the standardization. E.g.,
+-- @last \"231\" = last (fromList [1,2,0]) = 0@.
+last :: Perm -> Int
+last w | size w > 0 = w `unsafeAt` (size w - 1)
+       | otherwise  = 0
+
+-- | Length of the left-most increasing run: largest @i@ such that
+-- @w[0] \< w[1] \< ... \< w[i-1]@.
+lir :: Perm -> Int
+lir = marshal c_lir
+
+-- | Length of the left-most decreasing run: largest @i@ such that
+-- @w[0] \> w[1] \> ... \> w[i-1]@.
+ldr :: Perm -> Int
+ldr = marshal c_ldr
+
+-- | Length of the right-most increasing run: largest @i@ such that
+-- @w[n-i] \< ... \< w[n-2] \< w[n-1]@.
+rir :: Perm -> Int
+rir = ldr . D8.reverse
+
+-- | Length of the right-most decreasing run: largest @i@ such that
+-- @w[n-i] \> ... \> w[n-2] \> w[n-1]@.
+rdr :: Perm -> Int
+rdr = lir . D8.reverse
+
+-- | 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 -> Int
+comp = marshal c_comp
+
+-- | The number of skew components. E.g., @[5,7,4,6,3,1,0,2]@ has three
+-- skew components: @[5,7,4,6]@, @[3]@ and @[1,0,2]@.
+scomp :: Perm -> Int
+scomp = comp . D8.complement
+
+-- | 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 -> Int
+ep = marshal c_ep
+
+-- | The dimension of a permutation is defined as the largest
+-- non-fixed-point, or zero if all points are fixed.
+dim :: Perm -> Int
+dim = marshal c_dim
+
+-- | The number of small ascents. A /small ascent/ in @w@ is an index
+-- @i@ such that @w[i] + 1 == w[i+1]@.
+asc0 :: Perm -> Int
+asc0 = marshal c_asc0
+
+-- | The number of small descents. A /small descent/ in @w@ is an
+-- index @i@ such that @w[i] == w[i+1] + 1@.
+des0 :: Perm -> Int
+des0 = marshal c_des0
+
+-- | The longest increasing subsequence.
+lis :: Perm -> Int
+lis w = case Y.shape (Y.fromPerm w) of
+          []    -> 0
+          (x:_) -> x
+
+-- | The longest decreasing subsequence.
+lds :: Perm -> Int
+lds = length . Y.recordingTableau . Y.fromPerm
+
+-- | The size of the shadow of @w@. That is, the number of different
+-- one point deletions of @w@.
+-- shad :: Perm -> Int
+-- shad = length . shadow . return . st
diff --git a/Sym/Permgram.hs b/Sym/Permgram.hs
new file mode 100644
--- /dev/null
+++ b/Sym/Permgram.hs
@@ -0,0 +1,96 @@
+-- |
+-- Copyright   : Anders Claesson 2013
+-- Maintainer  : Anders Claesson <anders.claesson@gmail.com>
+--
+-- Permutation diagrams, or permutations as monads.
+
+module Sym.Permgram
+    (
+    -- * Data types
+      Label
+    , Permgram
+
+    -- * Accessors
+    , perm
+    , label
+    , size
+
+    -- * Construct permgrams
+    , permgram
+    , inverse
+    ) where
+
+import Data.Ord
+import Data.List
+import Sym.Perm (Perm)
+import qualified Sym.Perm as P
+import Data.Array.Unboxed
+
+-- | The purpose of this data type is to assign labels to the indices of
+-- a given permutation.
+type Label a = Array Int a
+
+-- | A permgram consists of a permutation together with a label for each
+-- index of the permutation.
+data Permgram a = PGram {
+      -- | The underlying permutation.
+      perm  :: Perm
+      -- | The assignment of labels to indices.
+    , label :: Label a
+    }
+
+constituents :: Permgram a -> (Perm, [a])
+constituents (PGram v f) = (v, elems f)
+
+instance Show a => Show (Permgram a) where
+    show w =
+        let (v, ys) = constituents w
+        in unwords ["permgram", "(" ++ show v ++ ")", show ys]
+
+instance Eq a => Eq (Permgram a) where
+    u == v = constituents u == constituents v
+
+instance Ord a => Ord (Permgram a) where
+    compare u v =
+        case comparing size u v of
+          EQ -> comparing constituents u v
+          x  -> x
+
+-- | Construct a permgram from an underlying permutation and a list of
+-- labels.
+permgram :: Perm -> [a] -> Permgram a
+permgram v = PGram v . listArray (0, P.size v - 1) . cycle
+
+-- | The inverse permgram. It's obtained by mirroring the permgram in
+-- the /x=y/ diagonal.
+inverse :: Permgram a -> Permgram a
+inverse (PGram u f) = PGram (P.fromList v) (listArray (0,n-1) (map (f!) v))
+    where
+      v = map snd . sort $ zip (P.toList u) [0..] -- v = u^{-1}
+      n = P.size u
+
+-- | The size of a permgram is the size of the underlying permutation.
+size :: Permgram a -> Int
+size = P.size . perm
+
+instance Functor Permgram where
+    fmap f w = w { label = amap f (label w) }
+
+instance Monad Permgram where
+    return x = permgram (P.fromList [0]) [x]
+    w >>= f  = joinPermgram $ fmap f w
+
+joinPermgram :: Permgram (Permgram a) -> Permgram a
+joinPermgram w@(PGram u f) = PGram (P.fromList xs) (listArray (0,m-1) ys)
+    where
+      len = amap size f
+      m = sum $ elems len
+      n = size w
+      uInverse = map snd . sort $ zip (P.toList u) [0..]
+      a :: UArray Int Int
+      a = listArray (0,n-1) . scanl (+) 0 $ map (len!) uInverse
+      (xs, ys) = unzip $ do
+        i <- [0..n-1]
+        let PGram v g = f!i
+        let d = a ! (u `P.unsafeAt` i)
+        [ (d + v `P.unsafeAt` j, g!j) | j <- [0..len!i-1] ]
diff --git a/include/bij.h b/include/bij.h
deleted file mode 100644
--- a/include/bij.h
+++ /dev/null
@@ -1,4 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-void simion_schmidt(long *, const long *, long);
-void simion_schmidt_inverse(long *, const long *, long);
diff --git a/include/bit.h b/include/bit.h
deleted file mode 100644
--- a/include/bit.h
+++ /dev/null
@@ -1,4 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-unsigned int next(unsigned int);
-void ones(unsigned int *, const unsigned int);
diff --git a/include/d8.h b/include/d8.h
deleted file mode 100644
--- a/include/d8.h
+++ /dev/null
@@ -1,5 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-void inverse(long, const long *, long);
-void reverse(long, const long *, long);
-void complement(long, const long *, long);
diff --git a/include/group.h b/include/group.h
deleted file mode 100644
--- a/include/group.h
+++ /dev/null
@@ -1,4 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-void compose (long *, const long *, const long *, long);
-void act     (long *, const long *, const long *, long);
diff --git a/include/ordiso.h b/include/ordiso.h
deleted file mode 100644
--- a/include/ordiso.h
+++ /dev/null
@@ -1,3 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-int ordiso(const long *, const long *, const long *, long);
diff --git a/include/rank.h b/include/rank.h
deleted file mode 100644
--- a/include/rank.h
+++ /dev/null
@@ -1,4 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-void unrank(long *, long, double);
-double rank(long *, long);
diff --git a/include/simple.h b/include/simple.h
deleted file mode 100644
--- a/include/simple.h
+++ /dev/null
@@ -1,3 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-int simple(const long *, long);
diff --git a/include/sortop.h b/include/sortop.h
deleted file mode 100644
--- a/include/sortop.h
+++ /dev/null
@@ -1,4 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-void stacksort (long *, long *, long);
-void bubblesort(long *, long *, long);
diff --git a/include/stat.h b/include/stat.h
deleted file mode 100644
--- a/include/stat.h
+++ /dev/null
@@ -1,24 +0,0 @@
-/* (c) Anders Claesson 2013 */
-
-long asc  (const long *, long);
-long des  (const long *, long);
-long exc  (const long *, long);
-long fp   (const long *, long);
-long sfp  (const long *, long);
-long cyc  (const long *, long);
-long inv  (const long *, long);
-long maj  (const long *, long);
-long comaj(const long *, long);
-long peak (const long *, long);
-long vall (const long *, long);
-long dasc (const long *, long);
-long ddes (const long *, long);
-long lmin (const long *, long);
-long lmax (const long *, long);
-long lir  (const long *, long);
-long ldr  (const long *, long);
-long comp (const long *, long);
-long ep   (const long *, long);
-long dim  (const long *, long);
-long asc0 (const long *, long);
-long des0 (const long *, long);
diff --git a/sym.cabal b/sym.cabal
--- a/sym.cabal
+++ b/sym.cabal
@@ -1,5 +1,5 @@
 name:                sym
-version:             0.9
+version:             0.11
 synopsis:            Permutations, patterns, and statistics
 description:
   Definitions for permutations with an emphasis on permutation
@@ -10,7 +10,7 @@
 license-file:        LICENSE
 author:              Anders Claesson
 maintainer:          anders.claesson@gmail.com
-category:            Data
+category:            Math
 build-type:          Simple
 cabal-version:       >=1.8
 
@@ -19,25 +19,32 @@
   location:            git://github.com/akc/sym.git
 
 library
-  exposed-modules:     Data.CLongArray
-                       Data.Perm
-                       Data.Permgram
-                       Math.Perm
-                       Math.Perm.Component
-                       Math.Perm.Constructions
-                       Math.Perm.D8
-                       Math.Perm.Group
-                       Math.Perm.Bijection
-                       Math.Perm.Stat
-                       Math.Perm.Sort
-                       Math.Perm.Simple
-                       Math.Perm.Pattern
-                       Math.Perm.Class
-                       Math.Sym
-
-  other-modules:       Data.Perm.Internal
+  exposed-modules:     Sym
+                       Sym.Perm
+                       Sym.Perm.Meta
+                       Sym.Perm.SSYT
+                       Sym.Perm.Component
+                       Sym.Perm.Constructions
+                       Sym.Perm.D8
+                       Sym.Perm.Group
+                       Sym.Perm.Bijection
+                       Sym.Perm.Stat
+                       Sym.Perm.Sort
+                       Sym.Perm.Simple
+                       Sym.Perm.Pattern
+                       Sym.Perm.MeshPattern
+                       Sym.Perm.Class
+                       Sym.Permgram
+                       Sym.Internal.SubSeq
+                       Sym.Internal.CLongArray
+                       Sym.Internal.Size
+                       Sym.Internal.Util
 
-  build-depends:       base >= 3 && <= 4.7, array >=0.4, hashable >=1.1, QuickCheck >=2.5
+  build-depends:       base >= 3 && <= 4.7,
+                       array >=0.4,
+                       hashable >=1.1,
+                       containers,
+                       QuickCheck >=2.5
 
   ghc-prof-options:    -auto-all
   ghc-options:         -Wall
@@ -52,7 +59,3 @@
                        cbits/bit.c
                        cbits/simple.c
                        cbits/sortop.c
-
-  include-dirs:        include
-  includes:            rank.h, stat.h, d8.h, group.h, bij.h, ordiso.h, bit.h, simple.h, sortop.h
-  install-includes:    rank.h, stat.h, d8.h, group.h, bij.h, ordiso.h, bit.h, simple.h, sortop.h