diff --git a/Data/Permutation.hs b/Data/Permutation.hs
deleted file mode 100644
--- a/Data/Permutation.hs
+++ /dev/null
@@ -1,265 +0,0 @@
-{-# LANGUAGE CPP #-}
-{-# OPTIONS_GHC -fglasgow-exts #-}
------------------------------------------------------------------------------
--- |
--- Module     : Data.Permutation
--- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
--- License    : BSD3
--- Maintainer : Patrick Perry <patperry@stanford.edu>
--- Stability  : experimental
---
-
-module Data.Permutation (
-    -- * The Permutation type
-    Permutation,
-
-    -- * Creating permutations
-    permutation,
-    identity,
-    inverse,
-    
-    -- * Permutation properties
-    size,
-    apply,
-
-    -- * Applying permutations
-    applyWith,
-    invertWith,
-
-    -- * Converstion to/from other types
-    -- ** @ForeignPtr@s
-    fromForeignPtr,
-    toForeignPtr,
-    
-    -- ** Lists
-    toList,
-    fromList,
-
-    -- * Unsafe operations
-    withPermutationPtr,
-    unsafePermutation,
-    unsafeApply,
-    
-    ) where
-        
-import Control.Monad         ( foldM, liftM )
-import Data.IntSet           ( IntSet )
-import qualified Data.IntSet as IntSet
-import Foreign               ( Ptr, ForeignPtr, mallocForeignPtrArray, 
-                               withForeignPtr, pokeArray, peekArray, 
-                               advancePtr, peek, peekElemOff, pokeElemOff ) 
-import System.IO.Unsafe      ( unsafePerformIO )
-
-#if defined(__GLASGOW_HASKELL__)
-import GHC.Base                 ( realWorld# )
-import GHC.IOBase               ( IO(IO) )
-#endif
-
-inlinePerformIO :: IO a -> a
-#if defined(__GLASGOW_HASKELL__)
-inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
-#else
-inlinePerformIO = unsafePerformIO
-#endif
-{-# INLINE inlinePerformIO #-}
-
-
--- | Represents a permutation of the integers @[0..n)@.        
-data Permutation =
-        Perm {-# UNPACK #-} !Int             
-             {-# UNPACK #-} !(ForeignPtr Int)
-             {-# UNPACK #-} !Int
-
--- | Get the size, raw data array and offset
-toForeignPtr :: Permutation -> (Int, ForeignPtr Int, Int)
-toForeignPtr (Perm n f o) = (n, f, o)
-
--- | Convert size and raw array to a permutation.  No validation is
--- performed on the arguments.
-fromForeignPtr :: Int -> ForeignPtr Int -> Int -> Permutation
-fromForeignPtr = Perm
-
--- | Get the size of the permutation.
-size :: Permutation -> Int
-size (Perm n _ _) = n
-{-# INLINE size #-}
-
--- | Perform an operation, given a pointer to the start of the
--- permutation data
-withPermutationPtr :: Permutation -> (Ptr Int -> IO a) -> IO a
-withPermutationPtr (Perm _ fptr off) f =
-    withForeignPtr fptr $ \ptr ->
-        f (ptr `advancePtr` off)
-
--- | Apply a permutation to an integer.  The integer must be in the range
---   @[0..n)@.
-apply :: Permutation -> Int -> Int
-apply p@(Perm n _ _) i 
-    | i < 0 || i >= n =
-        error $ 
-            "applyPerm: Tried to apply permutation of size `" ++ show n ++
-            "' to the value `" ++ show i ++ "'."
-    | otherwise =
-        unsafeApply p i
-{-# INLINE apply #-}
-
--- | Same as 'apply' but does not range-check the argument.
-unsafeApply :: Permutation -> Int -> Int
-unsafeApply p i =
-    inlinePerformIO $ do
-        withPermutationPtr p $ flip peekElemOff i
-{-# INLINE unsafeApply #-}
-
--- | Create a permutation from a list of values.  The list must be of length
---   @n@ and contain each integer in @{0, 1, ..., (n-1) }@ exactly once.
---   The permutation that is returned will send the integer @i@ to its index
---   in the list.
-permutation :: Int -> [Int] -> Permutation
-permutation n is = 
-    let p = unsafePermutation n is
-    in case isValid p of
-           False -> error $ "Not a valid permutation."
-           True  -> p
-
--- | Same as 'permutation', but does not check that the inputs are valid.
-unsafePermutation :: Int -> [Int] -> Permutation
-unsafePermutation n is = 
-    unsafePerformIO $ do
-        fptr <- mallocForeignPtrArray n
-        withForeignPtr fptr $ \ptr -> pokeArray ptr is
-        return $ fromForeignPtr n fptr 0
-{-# NOINLINE unsafePermutation #-}
-
-
-fromList :: [Int] -> Permutation
-fromList is = permutation (length is) is
-
-toList :: Permutation -> [Int]
-toList p = unsafePerformIO $ 
-               withPermutationPtr p $ peekArray (size p)
-
--- | Create an identity permutation of the given size.
-identity :: Int -> Permutation
-identity n =
-    unsafePerformIO $ do
-        fptr <- mallocForeignPtrArray n
-        withForeignPtr fptr $ \ptr -> pokeArray ptr [0..(n-1)]
-        return $ fromForeignPtr n fptr 0
-
-
--- | Get the inverse of a permutation.
-inverse :: Permutation -> Permutation
-inverse p =
-    let n = size p
-    in
-        unsafePerformIO $ do
-            fptr <- mallocForeignPtrArray n
-            withForeignPtr fptr $ \ptr -> do
-                pokeArray ptr [0..(n-1)]
-                invertWith (swap ptr) p
-            return $ fromForeignPtr n fptr 0
-    where
-        swap :: Ptr Int -> Int -> Int -> IO ()
-        swap ptr i j = do
-            x <- peekElemOff ptr i
-            y <- peekElemOff ptr j
-            pokeElemOff ptr i y
-            pokeElemOff ptr j x
-{-# NOINLINE inverse #-}
-
-
-
-
-isValid :: Permutation -> Bool
-isValid p@(Perm n _ _) =
-    unsafePerformIO $
-        withPermutationPtr p $ \ptr -> do
-            liftM and $
-                mapM (\i -> peekElemOff ptr i 
-                            >>= \x -> isValidI ptr x i)
-                     [0..(n-1)]
-             
-    where
-        isValidI :: Ptr Int -> Int -> Int -> IO Bool
-        isValidI ptr x i =
-            liftM and $ 
-                sequence [ inRange x, isUnique x ptr i ]
-    
-        inRange :: Int -> IO Bool
-        inRange x =
-            return $ x >= 0 && x < n
-        
-        isUnique :: Int -> Ptr Int -> Int -> IO Bool
-        isUnique x ptr' n'
-            | n' == 0 =
-                return True
-            | otherwise = do
-                x' <- peek ptr'
-                if x' == x 
-                    then return False
-                    else isUnique x (ptr' `advancePtr` 1) (n'-1)
-    
--- | @applyWith swap perm@ applies the permutation as a sequence of swaps.  After 
--- this function is applied, @OUT[i] = IN[P[i]]@
-applyWith :: (Monad m) => (Int -> Int -> m ()) -> Permutation -> m ()
-applyWith swap p =
-    let n = size p
-    in foldM (flip $ doCycle swap) IntSet.empty [0..(n-1)] >> return ()
-    
-    where
-        doCycle :: (Monad m) => 
-            (Int -> Int -> m ()) -> Int -> IntSet -> m (IntSet)
-        doCycle swp i visited = 
-            if i `IntSet.member` visited 
-                then return visited
-                else let visited' = IntSet.insert i visited
-                         next     = unsafeApply p i
-                     in doCycle' swp i i next visited'
-            
-        doCycle' :: (Monad m) => 
-            (Int -> Int -> m ()) -> Int -> Int -> Int -> IntSet -> m (IntSet)
-        doCycle' swp start cur next visited
-            | next == start =
-                return visited 
-            | otherwise = 
-                let visited' = IntSet.insert next visited
-                    next'    = unsafeApply p next
-                in do
-                    swp cur next
-                    doCycle' swp start next next' visited'
-
--- | @invertWith swap p@ applies the inverse of the permutation as a 
--- sequence of swaps.  After this function is applied, @OUT[P[i]] = IN[i]@
-invertWith :: (Monad m) => (Int -> Int -> m ()) -> Permutation -> m ()
-invertWith swap p =
-    let n = size p
-    in foldM (flip $ doCycle swap) IntSet.empty [0..(n-1)] >> return ()
-    
-    where
-        doCycle :: Monad m => 
-            (Int -> Int -> m ()) -> Int -> IntSet -> m (IntSet)
-        doCycle swp i visited = 
-            if i `IntSet.member` visited 
-                then return visited
-                else let visited' = IntSet.insert i visited
-                         cur      = unsafeApply p i
-                     in doCycle' swp i cur visited'
-            
-        doCycle' :: Monad m => (Int -> Int -> m ()) -> Int -> Int -> IntSet -> m (IntSet)
-        doCycle' swp start cur visited
-            | cur == start =
-                return visited 
-            | otherwise = 
-                let visited' = IntSet.insert cur visited
-                    cur'     = unsafeApply p cur
-                in do
-                    swp start cur
-                    doCycle' swp start cur' visited'
-
-
-instance Show Permutation where
-    show p = "permutation "     ++ show (size p) ++ " " ++ show (toList p)
-    
-instance Eq Permutation where
-    (==) p q = (size p == size q) && (toList p == toList q)
-    
diff --git a/Setup.lhs b/Setup.lhs
--- a/Setup.lhs
+++ b/Setup.lhs
@@ -10,5 +10,5 @@
 >         then ioError $ userError $ "failed"
 >         else return ()
 >
-> main = defaultMainWithHooks defaultUserHooks
+> main = defaultMainWithHooks simpleUserHooks
 >        {runTests=testing}
diff --git a/examples/Enumerate.hs b/examples/Enumerate.hs
new file mode 100644
--- /dev/null
+++ b/examples/Enumerate.hs
@@ -0,0 +1,81 @@
+module Main where
+    
+import Control.Monad
+import Control.Monad.ST
+import Data.List( permutations )
+import Data.STRef
+import System.Environment
+
+import Data.Permute
+import Data.Permute.MPermute
+    
+-- | Execute an action on every permutation of a given order.  This function
+-- is unsafe, because it only allocates space for a single permutation.  The
+-- action @f@ should not retain any references to the passed-in @Permute@
+-- object, otherwise bad things will happen.  For instance, running
+-- >
+-- >    forAllPermutes 2 id
+-- >
+-- in ghci yields @[listPermute 2 [1,0],listPermute 2 [1,0]]@.
+--
+forAllPermutes :: Int -> (Permute -> a) -> [a]
+forAllPermutes n f = runST $ do
+    -- Allocate a mutable permutation initialized to the identity
+    p  <- newPermute n
+    
+    -- Run the action on all successors of p
+    runOnSuccessors p
+    
+  where
+    runOnSuccessors p = do
+        -- Cast the mutable permutation to an immutable one
+        -- and the action on the immutable
+        a <- liftM f (unsafeFreeze p)
+
+        -- Set the permutation to be equal to its successor
+        hasNext <- setNext p
+        
+        -- If a successor exists, recurse, otherwise stop
+        as <- unsafeInterleaveST $
+            if hasNext then runOnSuccessors p 
+                       else return []
+        return (a:as)
+
+
+forAllPermutesM_ :: (MPermute p m) => Int -> (Permute -> m a) -> m () 
+forAllPermutesM_ n f = sequence_ $ forAllPermutes n f
+{-# INLINE forAllPermutesM_ #-}
+
+-- | Count the number of permutations of a given order
+countAllPermutes :: Int -> Int
+countAllPermutes n = length $ forAllPermutes n id
+
+-- | Another version of the same function.  This one is slightly slower.
+countAllPermutes2 :: Int -> Int
+countAllPermutes2 n = runST $ do
+    count <- newSTRef 0
+    forAllPermutesM_ n $ (const $ modifySTRef' count (+1))
+    readSTRef count
+  where
+    modifySTRef' var f = do
+        old <- readSTRef var
+        writeSTRef var $! f old
+
+-- | Yet another version, this time using 'permutations' from Data.List.
+-- This version is faster but uses more memory.
+countAllPermutes3 :: Int -> Int
+countAllPermutes3 n = length $ permutations [0 .. n-1]
+
+
+-- | Print all permutations of a given order.
+printAllPermutes :: Int -> IO ()
+printAllPermutes n =
+    forAllPermutesM_ n (putStrLn . show . elems)
+
+
+main = do
+    n  <- fmap (read . head) getArgs
+    let count = countAllPermutes n
+    putStrLn $ 
+        "There are " ++ show count ++ " permutations of order " ++ show n ++ "."
+    
diff --git a/lib/Data/IntArray.hs b/lib/Data/IntArray.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/IntArray.hs
@@ -0,0 +1,151 @@
+{-# LANGUAGE MagicHash, UnboxedTuples #-}
+{-# OPTIONS_GHC -funbox-strict-fields #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.IntArray
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.IntArray (
+    IntArray,
+    STIntArray,
+
+    numElements,
+    unsafeAt,
+    elems,
+
+    newArray_,
+    
+    numElementsSTIntArray,
+    getNumElements,
+    unsafeRead,
+    unsafeWrite,
+    unsafeSwap,
+    
+    readElems,
+    writeElems,
+    
+    unsafeFreeze,
+    unsafeThaw,
+    
+    ) where
+
+import GHC.Base( Int(..) )
+import GHC.Prim
+import GHC.ST
+import Foreign( sizeOf )
+
+-----------------------------  Immutable arrays -----------------------------
+
+data IntArray = IntArray !Int (ByteArray#)
+
+{-# INLINE numElements #-}
+numElements :: IntArray -> Int
+numElements (IntArray n _) = n
+
+{-# INLINE unsafeAt #-}
+unsafeAt :: IntArray -> Int -> Int
+unsafeAt (IntArray _ arr#) (I# i#) =
+    case indexIntArray# arr# i# of { e# ->
+    I# e# }
+
+{-# INLINE elems #-}
+elems :: IntArray -> [Int]
+elems arr@(IntArray n _) =
+    [ unsafeAt arr i | i <- [0 .. n-1]]
+
+
+------------------------------  Mutable arrays ------------------------------
+
+data STIntArray s = STIntArray !Int (MutableByteArray# s)
+
+{-# INLINE newArray_ #-}
+newArray_ :: Int -> ST s (STIntArray s)
+newArray_ n@(I# n#) =
+    ST $ \s1# ->
+        case newByteArray# (n# *# sizeOfInt) s1# of { (# s2#, marr# #) ->
+        (# s2#, STIntArray n marr# #) }
+  where
+    sizeOfInt = case sizeOf (0 :: Int) of (I# s#) -> s#
+    
+{-# INLINE numElementsSTIntArray #-}
+numElementsSTIntArray :: STIntArray s -> Int
+numElementsSTIntArray (STIntArray n _) = n
+
+{-# INLINE getNumElements #-}
+getNumElements :: STIntArray s -> ST s Int
+getNumElements arr = return $! numElementsSTIntArray arr
+
+{-# INLINE unsafeRead #-}
+unsafeRead :: STIntArray s -> Int -> ST s Int
+unsafeRead (STIntArray _ marr#) (I# i#) =
+    ST $ \s1# -> 
+        case readIntArray# marr# i# s1# of { (# s2#, e# #) ->
+        let e = I# e# in
+        (# s2#, e #) }
+
+{-# INLINE unsafeWrite #-}
+unsafeWrite :: STIntArray s -> Int -> Int -> ST s ()
+unsafeWrite (STIntArray _ marr#) (I# i#) (I# e#) =
+    ST $ \s1# -> 
+        case writeIntArray# marr# i# e# s1# of { s2# ->
+        (# s2#, () #) }
+
+{-# INLINE unsafeSwap #-}
+unsafeSwap :: STIntArray s -> Int -> Int -> ST s ()
+unsafeSwap (STIntArray _ marr#) (I# i#) (I# j#) =
+    ST $ \s1# ->
+        let doSwap =
+                case readIntArray#  marr# i# s1# of { (# s2#, e# #) ->
+                case readIntArray#  marr# j# s2# of { (# s3#, f# #) ->
+                case writeIntArray# marr# i# f# s3# of { s4# ->
+                     writeIntArray# marr# j# e# s4# }}} in
+        if i# ==# j# then (# s1#, () #)
+                     else case doSwap of { s2# ->
+                          (# s2#, () #) }
+
+{-# INLINE readElems #-}
+readElems :: STIntArray s -> ST s [Int]
+readElems (STIntArray (I# n#) marr#) =
+    ST $ \s1# ->
+        let inlineReadList i# | i# ==# n# = []
+                              | otherwise = 
+                case readIntArray# marr# i# s1# of { (# _, e# #) ->
+                let e  = I# e#
+                    es = inlineReadList (i# +# 1#) in
+                (e:es) } in
+        case inlineReadList 0# of { es ->
+        (# s1#, es #)}
+        
+{-# INLINE writeElems #-}
+writeElems :: STIntArray s -> [Int] -> ST s ()
+writeElems (STIntArray (I# n#) marr#) es =
+    ST $ \s1# ->
+        let fillFromList i# xs s2# | i# ==# n# = s2#
+                                   | otherwise = case xs of
+                []         -> s2#
+                (I# y#):ys -> case writeIntArray# marr# i# y# s2# of { s3# ->
+                              fillFromList (i# +# 1#) ys s3# } in
+        case fillFromList 0# es s1# of { s2# ->
+        (# s2#, () #) }
+
+{-# INLINE unsafeFreeze #-}
+unsafeFreeze :: STIntArray s -> ST s IntArray
+unsafeFreeze (STIntArray n marr#) =
+    ST $ \s1# ->
+        case unsafeFreezeByteArray# marr# s1# of { (# s2#, arr# #) ->
+        let arr = IntArray n arr# in
+        (# s2#, arr #)}
+
+{-# INLINE unsafeThaw #-}
+unsafeThaw :: IntArray -> ST s (STIntArray s)
+unsafeThaw (IntArray n arr#) =
+    ST $ \s1# ->
+        let coerceArray :: State# s -> MutableByteArray# s
+            coerceArray _ = unsafeCoerce# arr#
+            marr# = coerceArray s1#
+            marr  = STIntArray n marr# in
+        (# s1#, marr #)
diff --git a/lib/Data/Permute.hs b/lib/Data/Permute.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Permute.hs
@@ -0,0 +1,168 @@
+{-# LANGUAGE Rank2Types #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Permute
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+-- Immutable permutations.
+
+module Data.Permute (
+    -- * Permutations
+    Permute,
+    
+    -- * Creating permutations
+    permute,
+    listPermute,
+    swapsPermute,
+
+    -- * Accessing permutation elements
+    apply,
+    unsafeApply,
+
+    -- * Permutation properties
+    size,
+    elems,
+    
+    -- * Permutation functions
+    inverse,
+    next,
+    prev,
+    
+    -- * Applying permutations
+    swaps,
+    invSwaps,
+    
+    -- * Sorting
+    sort,
+    sortBy,
+    order,
+    orderBy,
+    rank,
+    rankBy,
+    
+    ) where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Permute.Base
+import Data.Permute.ST
+
+-- | Construct an identity permutation of the given size.
+permute :: Int -> Permute
+permute n = runST $
+    unsafeFreeze =<< newPermute n
+
+-- | Construct a permutation from a list of elements.  
+-- @listPermute n is@ creates a permuation of size @n@ with
+-- the @i@th element equal to @is !! i@.  For the permutation to be valid,
+-- the list @is@ must have length @n@ and contain the indices @0..(n-1)@ 
+-- exactly once each.
+listPermute :: Int -> [Int] -> Permute
+listPermute n is = runST $
+    unsafeFreeze =<< newListPermute n is
+
+-- | Construct a permutation from a list of swaps.
+-- @swapsPermute n ss@ creats a permuation of size @n@ given by a
+-- sequence of swaps.
+-- If @ss@ is @[(i0,j0), (i1,j1), ..., (ik,jk)]@, the
+-- sequence of swaps is
+-- @i0 \<-> j0@, then 
+-- @i1 \<-> j1@, and so on until
+-- @ik \<-> jk@.
+swapsPermute :: Int -> [(Int,Int)] -> Permute
+swapsPermute n ss = runST $
+    unsafeFreeze =<< newSwapsPermute n ss
+
+-- | @apply p i@ gets the value of the @i@th element of the permutation
+-- @p@.  The index @i@ must be in the range @0..(n-1)@, where @n@ is the
+-- size of the permutation.
+apply :: Permute -> Int -> Int
+apply p i
+    | i >= 0 && i < size p = 
+        unsafeApply p i
+    | otherwise =
+        error "Invalid index"
+{-# INLINE apply #-}
+
+-- | Get the inverse of a permutation
+inverse :: Permute -> Permute
+inverse p = runST $ 
+    unsafeFreeze =<< getInverse =<< unsafeThaw p
+
+-- | Return the next permutation in lexicographic order, or @Nothing@ if
+-- there are no further permutations.  Starting with the identity permutation
+-- and repeatedly calling this function will iterate through all permutations
+-- of a given order.
+next :: Permute -> Maybe Permute
+next = nextPrevHelp setNext
+
+-- | Return the previous permutation in lexicographic order, or @Nothing@
+-- if there is no such permutation.
+prev :: Permute -> Maybe Permute
+prev = nextPrevHelp setPrev
+
+nextPrevHelp :: (forall s. STPermute s -> ST s Bool) 
+             -> Permute -> Maybe Permute
+nextPrevHelp set p = runST $ do
+    p' <- thaw p
+    set p' >>= \valid ->
+        if valid
+            then liftM Just $ unsafeFreeze p'
+            else return Nothing
+
+-- | Get a list of swaps equivalent to the permutation.  A result of
+-- @[ (i0,j0), (i1,j1), ..., (ik,jk) ]@ means swap @i0 \<-> j0@, 
+-- then @i1 \<-> j1@, and so on until @ik \<-> jk@.
+swaps :: Permute -> [(Int,Int)]
+swaps p = runST $
+    getSwaps =<< unsafeThaw p
+
+-- | Get a list of swaps equivalent to the inverse of permutation.
+invSwaps :: Permute -> [(Int,Int)]
+invSwaps p = runST $
+    getInvSwaps =<< unsafeThaw p
+
+
+-- | @sort n xs@ sorts the first @n@ elements of @xs@ and returns a 
+-- permutation which transforms @xs@ into sorted order.  The results are
+-- undefined if @n@ is greater than the length of @xs@.  This is a special 
+-- case of 'sortBy'.
+sort :: (Ord a) => Int -> [a] -> ([a], Permute)
+sort n xs = runST $ do
+    (xs',mp) <- getSort n xs
+    p <- unsafeFreeze mp
+    return (xs',p)
+    
+sortBy :: (a -> a -> Ordering) -> Int -> [a] -> ([a], Permute)
+sortBy cmp n xs = runST $ do
+    (xs',mp) <- getSortBy cmp n xs
+    p <- unsafeFreeze mp
+    return (xs',p)
+
+
+-- | @order n xs@ returns a permutation which rearranges the first @n@
+-- elements of @xs@ into ascending order. The results are undefined if @n@ is 
+-- greater than the length of @xs@.  This is a special case of 'orderBy'.
+order :: (Ord a) => Int -> [a] -> Permute
+order n xs = runST $ 
+    unsafeFreeze =<< getOrder n xs
+
+orderBy :: (a -> a -> Ordering) -> Int -> [a] -> Permute
+orderBy cmp n xs = runST $
+    unsafeFreeze =<< getOrderBy cmp n xs
+
+-- | @rank n xs@ eturns a permutation, the inverse of which rearranges the 
+-- first @n@ elements of @xs@ into ascending order. The returned permutation, 
+-- @p@, has the property that @p[i]@ is the rank of the @i@th element of @xs@. 
+-- The results are undefined if @n@ is greater than the length of @xs@.
+-- This is a special case of 'rankBy'.  
+rank :: (Ord a) => Int -> [a] -> Permute
+rank n xs = runST $
+    unsafeFreeze =<< getRank n xs
+
+rankBy :: (a -> a -> Ordering) -> Int -> [a] -> Permute
+rankBy cmp n xs = runST $
+    unsafeFreeze =<< getRankBy cmp n xs
diff --git a/lib/Data/Permute/Base.hs b/lib/Data/Permute/Base.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Permute/Base.hs
@@ -0,0 +1,110 @@
+{-# LANGUAGE Rank2Types #-}
+{-# OPTIONS_GHC -XMagicHash -XUnboxedTuples #-}
+{-# OPTIONS_HADDOCK hide #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Permute.Base
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Permute.Base 
+    where
+
+import Control.Monad
+import Control.Monad.ST
+import Foreign
+
+import Data.IntArray ( IntArray, STIntArray )
+import qualified Data.IntArray as Arr
+import qualified Data.IntArray as ArrST
+
+
+--------------------------------- Permute ---------------------------------
+
+-- | The immutable permutation data type.
+-- Internally, a permutation of size @n@ is stored as an
+-- @0@-based array of @n@ 'Int's.  The permutation represents a reordering of
+-- the integers @0, ..., (n-1)@.  The @i@th element of the array stores
+-- the value @p[i]@. 
+newtype Permute = Permute IntArray
+
+unsafeApply :: Permute -> Int -> Int
+unsafeApply (Permute p) i = Arr.unsafeAt p i
+{-# INLINE unsafeApply #-}
+
+-- | Get the size of the permutation.
+size :: Permute -> Int
+size (Permute p) = Arr.numElements p
+{-# INLINE size #-}
+
+-- | Get a list of the permutation elements.
+elems :: Permute -> [Int]
+elems (Permute p) = Arr.elems p
+{-# INLINE elems #-}
+
+instance Show Permute where
+    show p = "listPermute " ++ show (size p) ++ " " ++ show (elems p)
+    
+instance Eq Permute where
+    (==) p q = (size p == size q) && (elems p == elems q)
+
+
+--------------------------------- STPermute --------------------------------
+
+-- | A mutable permutation that can be manipulated in the 'ST' monad. The
+-- type argument @s@ is the state variable argument for the 'ST' type.
+newtype STPermute s = STPermute (STIntArray s)
+
+getSizeSTPermute :: STPermute s -> ST s Int
+getSizeSTPermute (STPermute marr) = ArrST.getNumElements marr
+{-# INLINE getSizeSTPermute #-}
+
+sizeSTPermute :: STPermute s -> Int
+sizeSTPermute (STPermute marr) = ArrST.numElementsSTIntArray marr
+{-# INLINE sizeSTPermute #-}
+
+newSTPermute :: Int -> ST s (STPermute s)
+newSTPermute n = do
+    p@(STPermute marr) <- newSTPermute_ n
+    ArrST.writeElems marr [0 .. n-1]
+    return $! p
+{-# INLINE newSTPermute #-}
+
+newSTPermute_ :: Int -> ST s (STPermute s)
+newSTPermute_ n = do
+    when (n < 0) $ fail "invalid size"
+    liftM STPermute $ ArrST.newArray_ n
+{-# INLINE newSTPermute_ #-}
+
+unsafeGetElemSTPermute :: STPermute s -> Int -> ST s Int
+unsafeGetElemSTPermute (STPermute marr) i = ArrST.unsafeRead marr i
+{-# INLINE unsafeGetElemSTPermute #-}
+
+unsafeSetElemSTPermute :: STPermute s -> Int -> Int -> ST s ()
+unsafeSetElemSTPermute (STPermute marr) i x = ArrST.unsafeWrite marr i x
+{-# INLINE unsafeSetElemSTPermute #-}
+
+unsafeSwapElemsSTPermute :: STPermute s -> Int -> Int -> ST s ()
+unsafeSwapElemsSTPermute (STPermute marr) i j = ArrST.unsafeSwap marr i j
+{-# INLINE unsafeSwapElemsSTPermute #-}
+
+getElemsSTPermute :: STPermute s -> ST s [Int]
+getElemsSTPermute (STPermute marr) = ArrST.readElems marr
+{-# INLINE getElemsSTPermute #-}
+
+setElemsSTPermute :: STPermute s -> [Int] -> ST s ()
+setElemsSTPermute (STPermute marr) is = ArrST.writeElems marr is
+{-# INLINE setElemsSTPermute #-}
+
+unsafeFreezeSTPermute :: STPermute s -> ST s Permute
+unsafeFreezeSTPermute (STPermute marr) = 
+    (liftM Permute . ArrST.unsafeFreeze) marr
+{-# INLINE unsafeFreezeSTPermute #-}
+
+unsafeThawSTPermute :: Permute -> ST s (STPermute s)
+unsafeThawSTPermute (Permute arr) =
+    (liftM STPermute . ArrST.unsafeThaw) arr
+{-# INLINE unsafeThawSTPermute #-}
diff --git a/lib/Data/Permute/IO.hs b/lib/Data/Permute/IO.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Permute/IO.hs
@@ -0,0 +1,20 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Permute.IO
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+-- Mutable permutations in the 'IO' monad.
+
+module Data.Permute.IO (
+    -- * Permutations
+    IOPermute,
+    
+    -- * Overloaded mutable permutation interface
+    module Data.Permute.MPermute
+    ) where
+
+import Data.Permute.IOBase( IOPermute )
+import Data.Permute.MPermute
diff --git a/lib/Data/Permute/IOBase.hs b/lib/Data/Permute/IOBase.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Permute/IOBase.hs
@@ -0,0 +1,64 @@
+{-# OPTIONS_HADDOCK hide #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Permute.IOBase
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Permute.IOBase
+    where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Permute.Base
+
+-- | A mutable permutation that can be manipulated in the 'IO' monad. The
+-- type argument @s@ is the state variable argument for the 'IO' type.
+newtype IOPermute = IOPermute (STPermute RealWorld)
+
+newIOPermute :: Int -> IO (IOPermute)
+newIOPermute n = liftM IOPermute $ stToIO (newSTPermute n)
+{-# INLINE newIOPermute #-}
+
+newIOPermute_ :: Int -> IO (IOPermute)
+newIOPermute_ n = liftM IOPermute $ stToIO (newSTPermute_ n)
+{-# INLINE newIOPermute_ #-}
+
+getSizeIOPermute :: IOPermute -> IO Int
+getSizeIOPermute (IOPermute p) = stToIO $ getSizeSTPermute p
+{-# INLINE getSizeIOPermute #-}
+
+sizeIOPermute :: IOPermute -> Int
+sizeIOPermute (IOPermute p) = sizeSTPermute p
+{-# INLINE sizeIOPermute #-}
+        
+unsafeGetElemIOPermute :: IOPermute -> Int -> IO Int
+unsafeGetElemIOPermute (IOPermute p) i = stToIO $ unsafeGetElemSTPermute p i
+{-# INLINE unsafeGetElemIOPermute #-}
+
+unsafeSetElemIOPermute :: IOPermute -> Int -> Int -> IO ()
+unsafeSetElemIOPermute (IOPermute p) i x = stToIO $ unsafeSetElemSTPermute p i x
+{-# INLINE unsafeSetElemIOPermute #-}
+
+unsafeSwapElemsIOPermute :: IOPermute -> Int -> Int -> IO ()
+unsafeSwapElemsIOPermute (IOPermute p) i j = stToIO $ unsafeSwapElemsSTPermute p i j
+{-# INLINE unsafeSwapElemsIOPermute #-}
+
+getElemsIOPermute :: IOPermute -> IO [Int]
+getElemsIOPermute (IOPermute p) = stToIO $ getElemsSTPermute p
+{-# INLINE getElemsIOPermute #-}
+
+setElemsIOPermute :: IOPermute -> [Int] -> IO ()
+setElemsIOPermute (IOPermute p) is = stToIO $ setElemsSTPermute p is
+{-# INLINE setElemsIOPermute #-}
+
+unsafeFreezeIOPermute :: IOPermute -> IO Permute
+unsafeFreezeIOPermute (IOPermute p) = stToIO $ unsafeFreezeSTPermute p
+{-# INLINE unsafeFreezeIOPermute #-}
+
+unsafeThawIOPermute :: Permute -> IO (IOPermute)
+unsafeThawIOPermute p = liftM IOPermute $ stToIO (unsafeThawSTPermute p)
+{-# INLINE unsafeThawIOPermute #-}
diff --git a/lib/Data/Permute/MPermute.hs b/lib/Data/Permute/MPermute.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Permute/MPermute.hs
@@ -0,0 +1,474 @@
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, 
+        FlexibleContexts #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Permute.MPermute
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+-- An overloaded interface to mutable permutations. For permutation types which 
+-- can be used with this interface, see "Data.Permute.IO" and "Data.Permute.ST".
+--
+
+module Data.Permute.MPermute (
+    -- * Class of mutable permutation types
+    MPermute, 
+
+    -- * Constructing mutable permutations
+    newPermute,
+    newPermute_,
+    newListPermute,
+    newSwapsPermute,
+    newCopyPermute,
+    copyPermute,
+    setIdentity,
+    
+    -- * Accessing permutation elements
+    getElem,
+    setElem,
+    swapElems,
+    
+    -- * Permutation properties
+    getSize,
+    getElems,
+    setElems,
+    isValid,
+
+    -- * Permutation functions
+    getInverse,
+    copyInverse,
+    setNext,
+    setPrev,
+    
+    -- * Applying permutations
+    getSwaps,
+    getInvSwaps,
+    
+    -- * Sorting
+    getSort,
+    getSortBy,
+    getOrder,
+    getOrderBy,
+    getRank,
+    getRankBy,
+    
+    -- * Converstions between mutable and immutable permutations
+    freeze,
+    unsafeFreeze,
+    thaw,
+    unsafeThaw,
+    
+    -- * Unsafe operations
+    unsafeNewListPermute,
+    unsafeNewSwapsPermute,
+    unsafeGetElem,
+    unsafeSetElem,
+    unsafeSwapElems,
+
+    ) where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Function( on )
+import qualified Data.List as List
+import System.IO.Unsafe( unsafeInterleaveIO )
+
+import Data.Permute.Base
+import Data.Permute.IOBase
+
+
+--------------------------------- MPermute --------------------------------
+
+-- | Class for representing a mutable permutation.  The type is parameterized
+-- over the type of the monad, @m@, in which the mutable permutation will be
+-- manipulated.
+class (Monad m) => MPermute p m | p -> m, m -> p where
+    -- | Get the size of a permutation.
+    getSize :: p -> m Int
+
+    -- | Create a new permutation initialized to be the identity.
+    newPermute :: Int -> m p
+
+    -- | Allocate a new permutation but do not initialize it.
+    newPermute_ :: Int -> m p
+
+    unsafeGetElem :: p -> Int -> m Int
+    unsafeSetElem :: p -> Int -> Int -> m ()
+    unsafeSwapElems :: p -> Int -> Int -> m ()    
+
+    -- | Get a lazy list of the permutation elements.  The laziness makes this
+    -- function slightly dangerous if you are modifying the permutation.
+    getElems :: p -> m [Int]
+
+    -- | Set all the values of a permutation from a list of elements.
+    setElems :: p -> [Int] -> m ()
+
+    unsafeFreeze :: p -> m Permute
+    unsafeThaw :: Permute -> m p
+    
+    unsafeInterleaveM :: m a -> m a
+    
+        
+-- | Construct a permutation from a list of elements.  
+-- @newListPermute n is@ creates a permuation of size @n@ with
+-- the @i@th element equal to @is !! i@.  For the permutation to be valid,
+-- the list @is@ must have length @n@ and contain the indices @0..(n-1)@ 
+-- exactly once each.
+newListPermute :: (MPermute p m) => Int -> [Int] -> m p
+newListPermute n is = do
+    p <- unsafeNewListPermute n is
+    valid <- isValid p
+    when (not valid) $ fail "invalid permutation"
+    return $! p
+{-# INLINE newListPermute #-}
+
+unsafeNewListPermute :: (MPermute p m) => Int -> [Int] -> m p
+unsafeNewListPermute n is = do
+    p <- newPermute_ n
+    setElems p is
+    return $! p
+{-# INLINE unsafeNewListPermute #-}
+
+-- | Construct a permutation from a list of swaps.
+-- @newSwapsPermute n ss@ creates a permuation of size @n@ given a
+-- sequence of swaps.
+-- If @ss@ is @[(i0,j0), (i1,j1), ..., (ik,jk)]@, the
+-- sequence of swaps is
+-- @i0 \<-> j0@, then 
+-- @i1 \<-> j1@, and so on until
+-- @ik \<-> jk@.
+newSwapsPermute :: (MPermute p m) => Int -> [(Int,Int)] -> m p
+newSwapsPermute = newSwapsPermuteHelp swapElems
+{-# INLINE newSwapsPermute #-}
+
+unsafeNewSwapsPermute :: (MPermute p m) => Int -> [(Int,Int)] -> m p
+unsafeNewSwapsPermute = newSwapsPermuteHelp unsafeSwapElems
+{-# INLINE unsafeNewSwapsPermute #-}
+
+newSwapsPermuteHelp :: (MPermute p m) => (p -> Int -> Int -> m ())
+                       -> Int -> [(Int,Int)] -> m p
+newSwapsPermuteHelp swap n ss = do
+    p <- newPermute n
+    mapM_ (uncurry $ swap p) ss
+    return $! p
+{-# INLINE newSwapsPermuteHelp #-}
+
+-- | Construct a new permutation by copying another.
+newCopyPermute :: (MPermute p m) => p -> m p
+newCopyPermute p = do
+    n  <- getSize p
+    p' <- newPermute_ n
+    copyPermute p' p
+    return $! p'
+{-# INLINE newCopyPermute #-}
+
+-- | @copyPermute dst src@ copies the elements of the permutation @src@
+-- into the permtuation @dst@.  The two permutations must have the same
+-- size.
+copyPermute :: (MPermute p m) => p -> p -> m ()
+copyPermute dst src =
+    getElems src >>= setElems dst
+{-# INLINE copyPermute #-}
+
+-- | Set a permutation to the identity.
+setIdentity :: (MPermute p m) => p -> m ()
+setIdentity p = do
+    n <- getSize p
+    setElems p [0 .. n-1]
+{-# INLINE setIdentity #-}
+
+-- | @getElem p i@ gets the value of the @i@th element of the permutation
+-- @p@.  The index @i@ must be in the range @0..(n-1)@, where @n@ is the
+-- size of the permutation.
+getElem :: (MPermute p m) => p -> Int -> m Int
+getElem p i = do
+    n <- getSize p
+    when (i < 0 || i >= n) $ fail "getElem: invalid index"
+    unsafeGetElem p i
+{-# INLINE getElem #-}
+
+-- | @setElem p i x@ sets the value of the @i@th element of the permutation
+-- @p@.  The index @i@ must be in the range @0..(n-1)@, where @n@ is the
+-- size of the permutation.
+setElem :: (MPermute p m) => p -> Int -> Int -> m ()
+setElem p i x = do
+    n <- getSize p
+    when (i < 0 || i >= n) $ fail "getElem: invalid index"
+    unsafeSetElem p i x
+{-# INLINE setElem #-}
+
+-- | @swapElems p i j@ exchanges the @i@th and @j@th elements of the 
+-- permutation @p@.
+swapElems :: (MPermute p m) => p -> Int -> Int -> m ()
+swapElems p i j = do
+    n <- getSize p
+    when (i < 0 || i >= n || j < 0 || j >= n) $ fail "swapElems: invalid index"
+    unsafeSwapElems p i j
+{-# INLINE swapElems #-}
+
+-- | Returns whether or not the permutation is valid.  For it to be valid,
+-- the numbers @0,...,(n-1)@ must all appear exactly once in the stored
+-- values @p[0],...,p[n-1]@.
+isValid :: (MPermute p m) => p -> m Bool
+isValid p = do
+    n <- getSize p
+    liftM and $ validIndices n
+  where
+    j `existsIn` i = do
+        seen <- liftM (take i) $ getElems p
+        return $ (any (==j)) seen
+        
+    isValidIndex n i = do
+        i' <- unsafeGetElem p i
+        valid  <- return $ i' >= 0 && i' < n
+        unique <- liftM not (i' `existsIn` i)
+        return $ valid && unique
+
+    validIndices n = validIndicesHelp n 0
+
+    validIndicesHelp n i
+        | i == n = return []
+        | otherwise = do
+            a  <- isValidIndex n i
+            as <- unsafeInterleaveM $ validIndicesHelp n (i+1)
+            return (a:as)
+{-# INLINE isValid #-}
+
+-- | Compute the inverse of a permutation.  
+getInverse :: (MPermute p m) => p -> m p
+getInverse p = do
+    n <- getSize p
+    q <- newPermute_ n
+    copyInverse q p
+    return $! q
+{-# INLINE getInverse #-}
+
+-- | Set one permutation to be the inverse of another.  
+-- @copyInverse inv p@ computes the inverse of @p@ and stores it in @inv@.
+-- The two permutations must have the same size.
+copyInverse :: (MPermute p m) => p -> p -> m ()
+copyInverse dst src = do
+    n  <- getSize src
+    n' <- getSize dst
+    when (n /= n') $ fail "permutation size mismatch"
+    forM_ [0 .. n-1] $ \i -> do
+        i' <- unsafeGetElem src i
+        unsafeSetElem dst i' i
+{-# INLINE copyInverse #-}
+    
+-- | Advance a permutation to the next permutation in lexicogrphic order and
+-- return @True@.  If no further permutaitons are available, return @False@ and
+-- leave the permutation unmodified.  Starting with the idendity permutation 
+-- and repeatedly calling @setNext@ will iterate through all permutations of a 
+-- given size.
+setNext :: (MPermute p m) => p -> m Bool
+setNext = setNextBy compare
+{-# INLINE setNext #-}
+
+-- | Step backwards to the previous permutation in lexicographic order and
+-- return @True@.  If there is no previous permutation, return @False@ and
+-- leave the permutation unmodified.
+setPrev :: (MPermute p m) => p -> m Bool
+setPrev = setNextBy (flip compare)
+{-# INLINE setPrev #-}
+
+setNextBy :: (MPermute p m) => (Int -> Int -> Ordering) -> p -> m Bool
+setNextBy cmp p = do
+    n <- getSize p
+    if n > 1
+        then do
+            findLastAscent (n-2) >>=
+                maybe (return False) (\i -> do
+                    i'     <- unsafeGetElem p i
+                    i1'    <- unsafeGetElem p (i+1)
+                    (k,k') <- findSmallestLargerThan n i' (i+2) (i+1) i1'
+                    
+                    -- swap i and k
+                    unsafeSetElem p i k'
+                    unsafeSetElem p k i'
+                    
+                    reverseElems (i+1) (n-1)
+                    
+                    return True
+                )
+        else 
+            return False
+        
+  where
+    i `lt` j = cmp i j == LT
+    i `gt` j = cmp i j == GT
+    
+    findLastAscent i = do
+        ascent <- isAscent i
+        if ascent then return (Just i) else recurse
+      where
+        recurse = if i /= 0 then findLastAscent (i-1) else return Nothing 
+    
+    findSmallestLargerThan n i' j k k'
+        | j < n = do
+            j' <- unsafeGetElem p j
+            if j' `gt` i' && j' `lt` k'
+                then findSmallestLargerThan n i' (j+1) j j'
+                else findSmallestLargerThan n i' (j+1) k k'
+        | otherwise =
+            return (k,k')
+            
+    isAscent i = liftM2 lt (unsafeGetElem p i) (unsafeGetElem p (i+1))
+    
+    reverseElems i j
+        | i >= j = return ()
+        | otherwise = do
+            unsafeSwapElems p i j
+            reverseElems (i+1) (j-1)
+{-# INLINE setNextBy #-}  
+
+
+-- | Get a lazy list of swaps equivalent to the permutation.  A result of
+-- @[ (i0,j0), (i1,j1), ..., (ik,jk) ]@ means swap @i0 \<-> j0@, 
+-- then @i1 \<-> j1@, and so on until @ik \<-> jk@.  The laziness makes this
+-- function slightly dangerous if you are modifying the permutation.
+getSwaps :: (MPermute p m) => p -> m [(Int,Int)]
+getSwaps = getSwapsHelp False
+{-# INLINE getSwaps #-}
+
+-- | Get a lazy list of swaps equivalent to the inverse of a permutation.
+getInvSwaps :: (MPermute p m) => p -> m [(Int,Int)]
+getInvSwaps = getSwapsHelp True
+{-# INLINE getInvSwaps #-}
+
+getSwapsHelp :: (MPermute p m) => Bool -> p -> m [(Int,Int)]
+getSwapsHelp inv p = do
+    n <- getSize p
+    liftM concat $
+        forM [0..(n-1)] $ \i -> do
+            k <- unsafeGetElem p i
+            least <- isLeast i k
+            if least 
+                then do
+                    i' <- unsafeGetElem p i
+                    unsafeInterleaveM $ doCycle i i i'
+                else
+                    return []
+  where
+    isLeast i k 
+        | k > i = do
+            k' <- unsafeGetElem p k
+            isLeast i k'
+        | k < i     = return False
+        | otherwise = return True
+        
+    doCycle start i i'
+        | i' == start = return []
+        | otherwise = do
+            i'' <- unsafeGetElem p i'
+            let s = if inv then (start,i') else (i,i')
+            ss <- unsafeInterleaveM $ doCycle start i' i''
+            return (s:ss)
+{-# INLINE getSwapsHelp #-}
+
+-- | Convert a mutable permutation to an immutable one.
+freeze :: (MPermute p m) => p -> m Permute
+freeze p = unsafeFreeze =<< newCopyPermute p
+{-# INLINE freeze #-}
+
+-- | Convert an immutable permutation to a mutable one.
+thaw :: (MPermute p m) => Permute -> m p
+thaw p = newCopyPermute =<< unsafeThaw p
+{-# INLINE thaw #-}
+
+-- | @getSort n xs@ sorts the first @n@ elements of @xs@ and returns a 
+-- permutation which transforms @xs@ into sorted order.  The results are
+-- undefined if @n@ is greater than the length of @xs@.  This is a special 
+-- case of 'getSortBy'.
+getSort :: (Ord a, MPermute p m) => Int -> [a] -> m ([a], p)
+getSort = getSortBy compare
+{-# INLINE getSort #-}
+    
+getSortBy :: (MPermute p m) => (a -> a -> Ordering) -> Int -> [a] -> m ([a], p)
+getSortBy cmp n xs =
+    let ys       = take n xs
+        (is,ys') = (unzip . List.sortBy (cmp `on` snd) . zip [0..]) ys
+    in liftM ((,) ys') $ unsafeNewListPermute n is
+{-# INLINE getSortBy #-}
+
+-- | @getOrder n xs@ returns a permutation which rearranges the first @n@
+-- elements of @xs@ into ascending order. The results are undefined if @n@ is 
+-- greater than the length of @xs@.  This is a special case of 'getOrderBy'.
+getOrder :: (Ord a, MPermute p m) => Int -> [a] -> m p
+getOrder = getOrderBy compare
+{-# INLINE getOrder #-}
+
+getOrderBy :: (MPermute p m) => (a -> a -> Ordering) -> Int -> [a] -> m p
+getOrderBy cmp n xs =
+    liftM snd $ getSortBy cmp n xs
+{-# INLINE getOrderBy #-}
+
+-- | @getRank n xs@ eturns a permutation, the inverse of which rearranges the 
+-- first @n@ elements of @xs@ into ascending order. The returned permutation, 
+-- @p@, has the property that @p[i]@ is the rank of the @i@th element of @xs@. 
+-- The results are undefined if @n@ is greater than the length of @xs@.
+-- This is a special case of 'getRankBy'.  
+getRank :: (Ord a, MPermute p m) => Int -> [a] -> m p
+getRank = getRankBy compare
+{-# INLINE getRank #-}
+
+getRankBy :: (MPermute p m) => (a -> a -> Ordering) -> Int -> [a] -> m p
+getRankBy cmp n xs = do
+    p <- getOrderBy cmp n xs
+    getInverse p
+{-# INLINE getRankBy #-}
+
+
+--------------------------------- Instances ---------------------------------
+
+instance MPermute (STPermute s) (ST s) where
+    getSize = getSizeSTPermute
+    {-# INLINE getSize #-}
+    newPermute = newSTPermute
+    {-# INLINE newPermute #-}
+    newPermute_ = newSTPermute_
+    {-# INLINE newPermute_ #-}
+    unsafeGetElem = unsafeGetElemSTPermute
+    {-# INLINE unsafeGetElem #-}
+    unsafeSetElem = unsafeSetElemSTPermute
+    {-# INLINE unsafeSetElem #-}
+    unsafeSwapElems = unsafeSwapElemsSTPermute
+    {-# INLINE unsafeSwapElems #-}
+    getElems = getElemsSTPermute
+    {-# INLINE getElems #-}
+    setElems = setElemsSTPermute
+    {-# INLINE setElems #-}
+    unsafeFreeze = unsafeFreezeSTPermute
+    {-# INLINE unsafeFreeze #-}
+    unsafeThaw = unsafeThawSTPermute
+    {-# INLINE unsafeThaw #-}
+    unsafeInterleaveM = unsafeInterleaveST
+    {-# INLINE unsafeInterleaveM #-}
+
+
+instance MPermute IOPermute IO where
+    getSize = getSizeIOPermute
+    {-# INLINE getSize #-}
+    newPermute = newIOPermute
+    {-# INLINE newPermute #-}
+    newPermute_ = newIOPermute_
+    {-# INLINE newPermute_ #-}
+    unsafeGetElem = unsafeGetElemIOPermute
+    {-# INLINE unsafeGetElem #-}
+    unsafeSetElem = unsafeSetElemIOPermute
+    {-# INLINE unsafeSetElem #-}
+    unsafeSwapElems = unsafeSwapElemsIOPermute
+    {-# INLINE unsafeSwapElems #-}
+    getElems = getElemsIOPermute
+    {-# INLINE getElems #-}
+    setElems = setElemsIOPermute
+    {-# INLINE setElems #-}
+    unsafeFreeze = unsafeFreezeIOPermute
+    {-# INLINE unsafeFreeze #-}
+    unsafeThaw = unsafeThawIOPermute
+    {-# INLINE unsafeThaw #-}
+    unsafeInterleaveM = unsafeInterleaveIO
+    {-# INLINE unsafeInterleaveM #-}
+
diff --git a/lib/Data/Permute/ST.hs b/lib/Data/Permute/ST.hs
new file mode 100644
--- /dev/null
+++ b/lib/Data/Permute/ST.hs
@@ -0,0 +1,32 @@
+{-# LANGUAGE Rank2Types #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Permute.ST
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+-- Mutable permutations in the 'ST' monad.
+
+module Data.Permute.ST (
+    -- * Permutations
+    STPermute,
+    runSTPermute,
+    
+    -- * Overloaded mutable permutation interface
+    module Data.Permute.MPermute
+    ) where
+
+import Control.Monad.ST
+
+import Data.Permute.Base( Permute, STPermute, unsafeFreezeSTPermute )
+import Data.Permute.MPermute
+
+-- | A safe way to create and work with a mutable permutation before returning 
+-- an immutable permutation for later perusal. This function avoids copying the 
+-- permutation before returning it - it uses unsafeFreeze internally, but this 
+-- wrapper is a safe interface to that function. 
+runSTPermute :: (forall s. ST s (STPermute s)) -> Permute
+runSTPermute p = runST (p >>= unsafeFreezeSTPermute)
+{-# INLINE runSTPermute #-}
diff --git a/permutation.cabal b/permutation.cabal
--- a/permutation.cabal
+++ b/permutation.cabal
@@ -1,11 +1,19 @@
 name:            permutation
-version:         0.1
+version:         0.2
 homepage:        http://stat.stanford.edu/~patperry/code/permutation
-synopsis:        Permutations and combinations
+synopsis:        A library for representing and applying permutations.
 description:
-    A library for representing and applying permutations.
+    This library includes data types for storing permutations.  It
+    implements pure and impure types, the latter which can be modified
+    in-place.  The main utility of the library is converting between
+    the linear representation of a permutation to a sequence of swaps.
+    This allows, for instance, applying a permutation or its inverse
+    to an array with O(1) memory use.
     .
-category:        Math
+    Much of the interface for the library is based on the permutation 
+    functions in the GNU Scientific Library (GSL).
+    .
+category:        Data Structures, Math
 license:         BSD3
 license-file:    LICENSE
 copyright:       (c) 2008. Patrick Perry <patperry@stanford.edu>
@@ -13,14 +21,32 @@
 maintainer:      Patrick Perry <patperry@stanford.edu>
 cabal-version: >= 1.2.0
 build-type:      Custom
-tested-with:     GHC ==6.8.2
+tested-with:     GHC ==6.8.2, GHC ==6.10.1
 
-extra-source-files:  tests/Unit.hs
+extra-source-files:  examples/Enumerate.hs
+                     tests/Test/Permute.hs
+                     tests/Driver.hs
+                     tests/Main.hs
+                     tests/Pure.hs
+                     tests/ST.hs
                      tests/Makefile
 
 library
-    exposed-modules: Data.Permutation
+    hs-source-dirs:  lib
+    exposed-modules: Data.Permute
+                     Data.Permute.MPermute
+                     Data.Permute.IO
+                     Data.Permute.ST
+                     
+    other-modules:   Data.IntArray
+                     Data.Permute.Base
+                     Data.Permute.IOBase
 
-    build-depends:   base, containers
+    build-depends:   base
+    extensions:      MultiParamTypeClasses, FunctionalDependencies, 
+                     FlexibleContexts, Rank2Types, MagicHash, UnboxedTuples
 
     ghc-options:     -Wall
+
+    if impl(ghc >= 6.9)
+      build-depends: ghc-prim
diff --git a/tests/Driver.hs b/tests/Driver.hs
new file mode 100644
--- /dev/null
+++ b/tests/Driver.hs
@@ -0,0 +1,176 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Driver
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Driver (
+    Natural(..),
+    Index(..),
+    ListPermute(..),
+    SwapsPermute(..),
+    Sort(..),
+    SortBy(..),
+    Swap(..),
+    
+    mytest,
+    mycheck,
+    mytests,
+    done,
+
+    ) where
+
+import Control.Monad
+import Data.List
+import Data.Ord
+import System.IO
+import System.Random
+import Test.QuickCheck
+import Text.Printf
+import Text.Show.Functions
+
+
+newtype Natural = Nat Int deriving (Eq,Show)
+instance Arbitrary Natural where
+    arbitrary = do
+        n <- arbitrary
+        return $ Nat (abs n)
+    
+    coarbitrary = undefined
+
+data Index = Index Int Int deriving (Eq,Show)
+instance Arbitrary Index where
+    arbitrary = do
+        (Nat n) <- arbitrary
+        i <- choose (0, n)
+        return $ Index (n + 1) i
+
+    coarbitrary = undefined
+
+data ListPermute = ListPermute Int [Int] deriving (Eq,Show)
+instance Arbitrary ListPermute where
+    arbitrary = do
+        (Nat n) <- arbitrary
+        xs <- vector n :: Gen [Int]
+        return . ListPermute n $ 
+            (snd . unzip) $ sortBy (comparing fst) $ zip xs [0..]
+
+    coarbitrary = undefined
+
+data SwapsPermute = SwapsPermute Int [(Int,Int)] deriving (Eq,Show)
+instance Arbitrary SwapsPermute where
+    arbitrary = do
+        (Nat n) <- arbitrary
+        let n' = n + 1
+        (Nat k) <- arbitrary
+        ss <- replicateM k (swap n')
+        return $ SwapsPermute n' ss
+
+    coarbitrary = undefined
+
+swap n = do
+    i <- choose (0,n-1)
+    j <- choose (0,n-1)
+    return (i,j)
+
+
+data Swap = Swap Int Int Int deriving (Eq,Show)
+instance Arbitrary Swap where
+    arbitrary = do
+        (Index n i) <- arbitrary
+        j <- choose (0,n-1)
+        return $ Swap n i j
+
+    coarbitrary = undefined
+
+instance Arbitrary Ordering where
+    arbitrary   = elements [ LT, GT, EQ ]
+    coarbitrary = coarbitrary . fromEnum
+
+data Sort = Sort Int [Int] deriving (Eq,Show)
+instance Arbitrary Sort where
+    arbitrary = do
+        (Index n i) <- arbitrary
+        xs <- vector n
+        return $ Sort i xs
+        
+    coarbitrary = undefined
+
+data SortBy = SortBy (Int -> Int -> Ordering) Int [Int] deriving (Show)
+instance Arbitrary SortBy where
+    arbitrary = do
+        cmp <- arbitrary
+        (Sort n xs) <- arbitrary
+        return $ SortBy cmp n xs
+    
+    coarbitrary = undefined
+
+
+------------------------------------------------------------------------
+--
+-- QC driver ( taken from xmonad-0.6 )
+--
+
+debug = False
+
+mytest :: Testable a => a -> Int -> IO (Bool, Int)
+mytest a n = mycheck defaultConfig
+    { configMaxTest=n
+    , configEvery   = \n args -> let s = show n in s ++ [ '\b' | _ <- s ] } a
+ -- , configEvery= \n args -> if debug then show n ++ ":\n" ++ unlines args else [] } a
+
+mycheck :: Testable a => Config -> a -> IO (Bool, Int)
+mycheck config a = do
+    rnd <- newStdGen
+    mytests config (evaluate a) rnd 0 0 []
+
+mytests :: Config -> Gen Result -> StdGen -> Int -> Int -> [[String]] -> IO (Bool, Int)
+mytests config gen rnd0 ntest nfail stamps
+    | ntest == configMaxTest config = done "OK," ntest stamps >> return (True, ntest)
+    | nfail == configMaxFail config = done "Arguments exhausted after" ntest stamps >> return (True, ntest)
+    | otherwise               =
+      do putStr (configEvery config ntest (arguments result)) >> hFlush stdout
+         case ok result of
+           Nothing    ->
+             mytests config gen rnd1 ntest (nfail+1) stamps
+           Just True  ->
+             mytests config gen rnd1 (ntest+1) nfail (stamp result:stamps)
+           Just False ->
+             putStr ( "Falsifiable after "
+                   ++ show ntest
+                   ++ " tests:\n"
+                   ++ unlines (arguments result)
+                    ) >> hFlush stdout >> return (False, ntest)
+     where
+      result      = generate (configSize config ntest) rnd2 gen
+      (rnd1,rnd2) = split rnd0
+
+done :: String -> Int -> [[String]] -> IO ()
+done mesg ntest stamps = putStr ( mesg ++ " " ++ show ntest ++ " tests" ++ table )
+  where
+    table = display
+            . map entry
+            . reverse
+            . sort
+            . map pairLength
+            . group
+            . sort
+            . filter (not . null)
+            $ stamps
+
+    display []  = ".\n"
+    display [x] = " (" ++ x ++ ").\n"
+    display xs  = ".\n" ++ unlines (map (++ ".") xs)
+
+    pairLength xss@(xs:_) = (length xss, xs)
+    entry (n, xs)         = percentage n ntest
+                       ++ " "
+                       ++ concat (intersperse ", " xs)
+
+    percentage n m        = show ((100 * n) `div` m) ++ "%"
+
+------------------------------------------------------------------------
+
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,31 @@
+
+import Control.Monad
+import System.Environment
+import Text.Printf
+
+import Driver
+import Pure
+import ST
+
+main :: IO ()
+main = do
+    args <- getArgs
+    let n = if null args then 100 else read (head args)
+    
+    (results, passed) <- liftM unzip $
+        foldM ( \prev (name,subtests) -> do
+                     printf "\n%s\n" name
+                     printf "%s\n" $ replicate (length name) '-'
+                     cur <- mapM (\(s,a) -> printf "%-30s: " s >> a n) subtests
+                     return (prev ++ cur)
+              )
+              []
+              tests
+               
+    printf "\nPassed %d tests!\n\n" (sum passed)
+    when (not . and $ results) $ fail "\nNot all tests passed!"
+ where
+
+    tests = [ ("Permute"   , tests_Permute)
+            , ("STPermute" , tests_STPermute)
+            ]
diff --git a/tests/Makefile b/tests/Makefile
--- a/tests/Makefile
+++ b/tests/Makefile
@@ -1,8 +1,16 @@
 all:
-	ghc -O -i. -i.. Unit.hs --make
-	./Unit
+	ghc -O -i. -i../lib Main.hs --make -o test-permute
+	./test-permute
 
+hpc:
+	ghc -fforce-recomp -i. -i../lib -fhpc --make Main.hs -o test-permute
+	rm -f test-permute.tix
+	./test-permute
+	hpc markup test-permute
+
 clean:
-	rm -f *~ Unit *.hi *.o \
-		../Data/Permutation.hi ../Data/Permutation.o \
-		../System/Random/Permutation.o
+	find ../lib . -name '*.hi' | xargs rm -f
+	find ../lib . -name '*.o'  | xargs rm -f
+	find . -name '*.html' | xargs rm -f
+	rm -f test-permute test-permute.tix
+	
diff --git a/tests/Pure.hs b/tests/Pure.hs
new file mode 100644
--- /dev/null
+++ b/tests/Pure.hs
@@ -0,0 +1,173 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+module Pure (
+    tests_Permute
+    ) where
+    
+import Control.Monad.ST
+import Data.Array.ST
+import Data.List( foldl' )
+import qualified Data.List as List
+import Data.Maybe( fromJust )
+
+import Data.Permute
+
+import Driver
+import Test.QuickCheck
+    
+import Test.Permute()
+import qualified Test.Permute as Test
+
+
+prop_size_permute (Nat n) =
+    size (permute n) == n
+prop_elems_permute (Nat n) =
+    elems (permute n) == [0..(n-1)]
+
+prop_size_listPermute (ListPermute n is) =
+    size (listPermute n is) == n
+prop_elems_listPermute (ListPermute n is) =
+    elems (listPermute n is) == is
+
+prop_size_swapsPermute (SwapsPermute n ss) =
+    size (swapsPermute n ss) == n
+prop_elems_swapsPermute (SwapsPermute n ss) =
+    elems (swapsPermute n ss) == map at [0..(n-1)]
+  where
+    at i = foldl' doSwap i $ reverse ss
+    doSwap k (i,j) | k == i    = j
+                   | k == j    = i
+                   | otherwise = k
+
+prop_apply       = prop_apply_help apply
+prop_unsafeApply = prop_apply_help unsafeApply
+prop_apply_help a =
+    forAll arbitrary $ \(Index n i) ->
+    forAll (Test.permute n) $ \p ->
+        a p i == (elems p) !! i
+
+prop_size_inverse (p :: Permute) =
+    size (inverse p) == size p
+prop_elems_inverse (p :: Permute) =
+    all (\i -> is' !! (apply p i) == i) [0..(n-1)]
+  where
+    n   = size p
+    is' = elems (inverse p)
+
+prop_swaps (Nat n) =
+    forAll (Test.permute n) $ \p ->
+    forAll (vector n) $ \xs ->
+        let xs' = applySwaps (swaps p) xs
+        in all (\i -> xs' !! i == xs !! (apply p i)) [0..(n-1)]
+
+prop_invSwaps (Nat n) =
+    forAll (Test.permute n) $ \p ->
+    forAll (vector n) $ \xs ->
+        let xs' = applySwaps (invSwaps p) xs
+        in all (\i -> xs' !! (apply p i) == xs !! i) [0..(n-1)]
+
+prop_swaps_inverse (Nat n) =
+    forAll (Test.permute n) $ \p ->
+    forAll (vector n) $ \xs ->
+        applySwaps (swaps $ inverse p) xs == (applySwaps (invSwaps p) xs)
+    
+prop_invSwaps_inverse (Nat n) =
+    forAll (Test.permute n) $ \p ->
+    forAll (vector n) $ \xs ->
+        applySwaps (invSwaps $ inverse p) xs == (applySwaps (swaps p) xs)
+
+prop_prev_permute (Nat n) =
+    prev (permute n) == Nothing
+prop_next_last (Nat n) =
+    next (listPermute n $ reverse [0..(n-1)]) == Nothing
+
+prop_next_prev (p :: Permute) =
+    case prev p of
+        Just p' -> p == (fromJust $ next p')
+        Nothing -> p == permute n
+  where
+    n = size p
+    
+prop_prev_next (p :: Permute) =
+    case next p of
+        Just p' -> p == (fromJust $ prev p')
+        Nothing -> p == (listPermute n $ reverse [0..(n-1)])
+  where
+    n = size p
+
+prop_fst_sort (Sort n xs) = let
+    ys = take n xs
+    in (fst . sort n) xs == (List.sort ys)
+prop_snd_sort (Sort n xs) = let
+    ys = take n xs
+    in applySwaps (swaps $ snd $ sort n xs) ys == (List.sort ys)
+
+prop_fst_sortBy (SortBy cmp n xs) = let 
+    ys = take n xs
+    in (fst . sortBy cmp n) xs == (List.sortBy cmp ys)
+prop_snd_sortBy (SortBy cmp n xs) = let 
+    ys = take n xs
+    in applySwaps (swaps $ snd $ sortBy cmp n xs) ys == (List.sortBy cmp ys)
+
+prop_order (Sort n xs) = let 
+    ys = take n xs
+    in applySwaps (swaps $ order n xs) ys == (List.sort ys)
+
+prop_orderBy (SortBy cmp n xs) = let 
+    ys = take n xs
+    in applySwaps (swaps $ orderBy cmp n xs) ys == (List.sortBy cmp ys)
+
+prop_rank (Sort n xs) = let
+    ys = take n xs
+    in applySwaps (invSwaps $ rank n xs) ys == (List.sort ys)
+
+prop_rankBy (SortBy cmp n xs) = let
+    ys = take n xs
+    in applySwaps (invSwaps $ rankBy cmp n xs) ys == (List.sortBy cmp ys)
+
+prop_swapsPermute_swaps (p :: Permute) =
+    swapsPermute (size p) (swaps p) == p
+
+
+tests_Permute = 
+    [ ("size . permute"          , mytest prop_size_permute)
+    , ("elems . permute"         , mytest prop_elems_permute)
+    , ("size . listPermute"      , mytest prop_size_listPermute)
+    , ("elems . listPermute"     , mytest prop_elems_listPermute)
+    , ("size . swapsPermute"     , mytest prop_size_swapsPermute)
+    , ("elems . swapsPermute"    , mytest prop_elems_swapsPermute)
+    , ("apply"                   , mytest prop_apply)
+    , ("unsafeApply"             , mytest prop_unsafeApply)
+    , ("size . inverse"          , mytest prop_size_inverse)
+    , ("elems . inverse"         , mytest prop_elems_inverse)
+    , ("swaps"                   , mytest prop_swaps)
+    , ("invSwaps"                , mytest prop_invSwaps)
+    , ("swaps . inverse"         , mytest prop_swaps_inverse)
+    , ("invSwaps . inverse"      , mytest prop_invSwaps_inverse)
+    , ("prev . permute"          , mytest prop_prev_permute)
+    , ("next (last permutation)" , mytest prop_next_last)
+    , ("next . prev"             , mytest prop_next_prev)
+    , ("prev . next"             , mytest prop_prev_next)
+    , ("fst . sort"              , mytest prop_fst_sort)
+    , ("snd . sort"              , mytest prop_snd_sort)
+    , ("fst . sortBy"            , mytest prop_fst_sortBy)
+    , ("snd . sortBy"            , mytest prop_snd_sortBy)
+    , ("order"                   , mytest prop_order)
+    , ("orderBy"                 , mytest prop_orderBy)
+    , ("rank"                    , mytest prop_rank)
+    , ("rankBy"                  , mytest prop_rankBy)
+    , ("swapsPermute . swaps"    , mytest prop_swapsPermute_swaps)
+    ]
+
+
+applySwaps :: [(Int,Int)] -> [Int] -> [Int]
+applySwaps ss xs = runST $ do
+    arr <- newListArray (0,length xs - 1) xs :: ST s (STUArray s Int Int) 
+    mapM_ (swap arr) ss
+    getElems arr
+  where
+    swap arr (i,j) = do
+        i' <- readArray arr i
+        j' <- readArray arr j
+        writeArray arr j i'
+        writeArray arr i j'
+
diff --git a/tests/ST.hs b/tests/ST.hs
new file mode 100644
--- /dev/null
+++ b/tests/ST.hs
@@ -0,0 +1,207 @@
+{-# LANGUAGE Rank2Types #-}
+module ST (
+    tests_STPermute
+    ) where
+    
+import Control.Monad
+import Control.Monad.ST
+
+import Data.Permute
+import Data.Permute.ST
+
+
+import Driver
+import Debug.Trace
+import Test.QuickCheck
+import Text.Printf
+    
+import Test.Permute()
+import qualified Test.Permute as Test
+
+newPermute_S n = permute n
+prop_NewPermute (Nat n) = 
+    newPermute n `equivalent` newPermute_S n
+
+newListPermute_S n is = listPermute n is
+prop_NewListPermute (ListPermute n is) =
+    newListPermute n is `equivalent` newListPermute_S n is
+
+newSwapsPermute_S n ss = swapsPermute n ss
+prop_NewSwapsPermute (SwapsPermute n ss) =
+    newSwapsPermute n ss `equivalent` newSwapsPermute_S n ss
+prop_UnsafeNewSwapsPermute (SwapsPermute n ss) =
+    unsafeNewSwapsPermute n ss `equivalent` newSwapsPermute_S n ss
+
+
+
+newCopyPermute_S p = (p, p)
+prop_NewCopyPermute =
+    implements
+        (\p -> newCopyPermute p >>= unsafeFreeze)
+        (\p -> newCopyPermute_S p)
+
+copyPermute_S p q = ((), q, q)
+prop_CopyPermute =
+    copyPermute `implements2` copyPermute_S
+
+setIdentity_S p = ((), permute (size p))
+prop_SetIdentity =
+    setIdentity `implements` setIdentity_S 
+
+getElem_S p i = ((elems p) !! i, p)
+prop_GetElem (Index n i) =
+    implementsFor n
+        (\p -> getElem   p i)
+        (\p -> getElem_S p i)
+
+prop_UnsafeGetElem (Index n i) =
+    implementsFor n
+        (\p -> unsafeGetElem p i)
+        (\p -> getElem_S p i)
+
+swapElems_S p i j = ((), p')
+  where
+    (n,is) = (size p, elems p)
+      
+    at k | k == i    = is !! j
+         | k == j    = is !! i
+         | otherwise = is !! k
+    p'   = listPermute n $ map at [0..(n-1)]
+
+prop_SwapElems (Swap n i j) =
+    implementsFor n
+        (\p -> swapElems p i j) 
+        (\p -> swapElems_S p i j)
+
+prop_UnsafeSwapElems (Swap n i j) =
+    implementsFor n
+        (\p -> unsafeSwapElems p i j) 
+        (\p -> swapElems_S p i j)
+
+
+getSize_S p = (length (elems p), p)
+prop_GetSize = getSize `implements` getSize_S
+      
+getElems_S p = (elems p, p)
+prop_GetElems = getElems `implements` getElems_S
+
+
+
+
+tests_STPermute = 
+    [ ("newPermute"               , mytest prop_NewPermute)
+    , ("newListPermute"           , mytest prop_NewListPermute)
+    , ("newSwapsPermute"          , mytest prop_NewSwapsPermute)
+    , ("unsafeNewSwapsPermute"    , mytest prop_UnsafeNewSwapsPermute)
+    , ("newCopyPermute"           , mytest prop_NewCopyPermute)
+    , ("copyPermute"              , mytest prop_CopyPermute)
+    , ("setIdentity"              , mytest prop_SetIdentity)
+    , ("getElem"                  , mytest prop_GetElem)
+    , ("unsafeGetElem"            , mytest prop_UnsafeGetElem)
+    , ("swapElems"                , mytest prop_SwapElems)
+    , ("unsafeSwapElems"          , mytest prop_UnsafeSwapElems)
+    , ("getSize"                  , mytest prop_GetSize)
+    , ("getElems"                 , mytest prop_GetElems)
+    ]
+
+
+------------------------------------------------------------------------
+-- 
+-- The specification language
+--
+    
+abstract :: STPermute s -> ST s Permute
+abstract = freeze
+
+commutes :: (Eq a, Show a) =>
+    STPermute s -> (STPermute s -> ST s a) ->
+        (Permute -> (a,Permute)) -> ST s Bool
+commutes p a f = do
+    old <- abstract p
+    r   <- a p
+    new <- abstract p
+    let s      = f old
+        s'     = (r,new)
+        passed = s == s'
+        
+    when (not passed) $
+        trace (printf ("expected `%s' but got `%s'") (show s) (show s'))
+              return ()
+              
+    return passed
+
+equivalent :: (forall s . ST s (STPermute s)) -> Permute -> Bool
+equivalent p s = runST $ do
+    p' <- (p >>= abstract)
+    when (not $ p' == s) $
+        trace (printf ("expected `%s' but got `%s'") (show s) (show p'))
+            return ()
+    return (p' == s)
+    
+implements :: (Eq a, Show a) =>
+    (forall s . STPermute s -> ST s a) ->
+    (Permute -> (a,Permute)) -> 
+        Property
+a `implements` f =
+    forAll arbitrary $ \(Nat n) ->
+        implementsFor n a f
+
+implementsFor :: (Eq a, Show a) =>
+    Int ->
+    (forall s . STPermute s -> ST s a) ->
+    (Permute -> (a,Permute)) -> 
+        Property
+implementsFor n a f =
+    forAll (Test.permute n) $ \p ->
+        runST $ do
+            p' <- unsafeThaw p
+            commutes p' a f
+
+implementsIf :: (Eq a, Show a) =>
+    (forall s . STPermute s -> ST s Bool) ->
+    (forall s . STPermute s -> ST s a) ->
+    (Permute -> (a, Permute)) -> 
+        Property
+implementsIf pre a f =
+    forAll arbitrary $ \p ->
+        runST ( do
+            p' <- thaw p
+            pre p') ==>
+        runST ( do
+            p' <- unsafeThaw p
+            commutes p' a f )
+
+
+commutes2 :: (Eq a, Show a) =>
+    STPermute s -> STPermute s -> (STPermute s -> STPermute s -> ST s a) ->
+        (Permute -> Permute -> (a,Permute,Permute)) -> ST s Bool
+commutes2 p q a f = do
+    oldp <- abstract p
+    oldq <- abstract q
+    r   <- a p q
+    newp <- abstract p
+    newq <- abstract q
+    let s      = f oldp oldq
+        s'     = (r,newp,newq)
+        passed = s == s'
+        
+    when (not passed) $
+        trace (printf ("expected `%s' but got `%s'") (show s) (show s'))
+              return ()
+              
+    return passed
+
+
+implements2 :: (Eq a, Show a) =>
+    (forall s . STPermute s -> STPermute s -> ST s a) ->
+    (Permute -> Permute -> (a,Permute,Permute)) -> 
+        Property
+implements2 a f =
+    forAll arbitrary $ \(Nat n) ->
+    forAll (Test.permute n) $ \p ->
+    forAll (Test.permute n) $ \q ->
+        runST $ do
+            p' <- unsafeThaw p
+            q' <- unsafeThaw q
+            commutes2 p' q' a f
+
diff --git a/tests/Test/Permute.hs b/tests/Test/Permute.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/Permute.hs
@@ -0,0 +1,40 @@
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.Permute
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+-- An 'Arbitrary' instance and functions for generating random permutations.
+--
+
+module Test.Permute (
+    permute
+    ) where
+
+import Data.List ( sortBy )
+import Data.Ord ( comparing )
+
+import Test.QuickCheck
+
+import Data.Permute ( Permute )
+import qualified Data.Permute as P
+
+-- | Generate a random permutation of the given size.
+permute :: Int -> Gen Permute
+permute n = do
+    xs <- vector n :: Gen [Int]
+    let is = (snd . unzip) $ sortBy (comparing fst) $ zip xs [0..]
+    return $ P.listPermute n is
+
+
+instance Arbitrary Permute where
+    arbitrary = do
+        n <- arbitrary >>= return . abs
+        p <- permute n
+        return $ p
+        
+    coarbitrary p =
+        coarbitrary $ P.elems p
diff --git a/tests/Unit.hs b/tests/Unit.hs
deleted file mode 100644
--- a/tests/Unit.hs
+++ /dev/null
@@ -1,214 +0,0 @@
-
-import Control.Monad ( forM_ )
-import Data.List ( sortBy )
-import Data.Ord  ( comparing )
-import Foreign
-import System.Random
-import Test.HUnit
-
-import Data.Permutation
-
-
-
-swapElems :: Storable a => Ptr a -> Int -> Int -> IO ()
-swapElems ptr i j =  do
-    x <- peekElemOff ptr i
-    y <- peekElemOff ptr j
-    pokeElemOff ptr i y
-    pokeElemOff ptr j x
-
-newRandomArray :: (Random a, Storable a) => Int -> IO (Ptr a)
-newRandomArray n = do
-    xs <- sequence $ replicate n randomIO
-    newArray xs
-
-newArrayCopy :: Storable a => Int -> Ptr a -> IO (Ptr a)
-newArrayCopy n x = do
-    y <- mallocArray n
-    copyArray y x n
-    return y
-
-newRandomPerm :: Int -> IO (Permutation)
-newRandomPerm n = do
-    xs <- sequence $ replicate n randomIO :: IO [Int]
-    let ixs  = zip [0..] xs
-        ixs' = sortBy (comparing snd) ixs
-        is   = (fst . unzip) ixs'
-        p    = permutation n is
-    return p
-    
-
-testIdentity :: Test
-testIdentity = TestCase $ do
-    let n = 30
-        p = identity n
-    
-    x <- newRandomArray n :: IO (Ptr Double)
-    y <- newArrayCopy n x
-    applyWith (swapElems y) p
-    
-    xs <- peekArray n x
-    ys <- peekArray n y
-    
-    assertEqual ""
-        xs
-        ys
-        
-        
-testSmallCycle :: Test
-testSmallCycle = TestCase $ do
-    let n = 3
-        p = permutation 3  [ 2, 0, 1 ]
-
-    x <- newArray [ 6, 7, 8 ] :: IO (Ptr Int)
-    applyWith (swapElems x) p
-    
-    xs <- peekArray n x
-    assertEqual ""
-        [ 8, 6, 7 ]        
-        xs
-
-testSmallCycleInv :: Test
-testSmallCycleInv = TestCase $ do
-    let n = 3
-        p = permutation 3  [ 2, 0, 1 ]
-
-    x <- newArray [ 6, 7, 8 ] :: IO (Ptr Int)
-    invertWith (swapElems x) p
-    
-    xs <- peekArray n x
-    assertEqual ""
-        [ 7, 8, 6 ]
-        xs
-
-
-testSingleCycleInv :: Test
-testSingleCycleInv = TestCase $ do
-    let n = 5
-        p = permutation n [ 3, 0, 1, 4, 2 ]
-        
-    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)
-    invertWith (swapElems x) p
-    
-    xs <- peekArray n x
-    assertEqual ""
-        [ 0.382, 0.818, 0.037, 0.237, -0.413 ]
-        xs
-
-testSingleCycle :: Test
-testSingleCycle = TestCase $ do
-    let n = 5
-        p = permutation n [ 3, 0, 1, 4, 2 ]
-
-    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)
-    applyWith (swapElems x) p
-    
-    xs <- peekArray n x
-    assertEqual ""
-        [ -0.413, 0.237, 0.382, 0.037, 0.818 ]
-        xs
-        
-testMultiCycleInv :: Test
-testMultiCycleInv = TestCase $ do
-    let n = 5
-        p = permutation n [ 2, 3, 0, 4, 1 ]
-        
-    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)
-    invertWith (swapElems x) p
-    
-    xs <- peekArray n x
-    assertEqual ""
-        [ 0.818, 0.037, 0.237, 0.382, -0.413 ]
-        xs
-
-testMultiCycle :: Test
-testMultiCycle = TestCase $ do
-    let n = 5
-        p = permutation n [ 2, 3, 0, 4, 1 ]
-    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)
-    applyWith (swapElems x) p
-    
-    xs <- peekArray n x
-    assertEqual ""
-        [ 0.818, -0.413, 0.237, 0.037, 0.382 ]
-        xs
-
-testApply :: Int -> Test
-testApply n = TestCase $ do
-    x <- newRandomArray n :: IO (Ptr Double)
-    y <- newArrayCopy n x
-    p <- newRandomPerm n
-    applyWith  (swapElems y) p
-    
-    xs <- peekArray n x
-    ys <- peekArray n y
-
-    forM_ [0..(n-1)] $ \i ->
-        assertEqual ("position " ++ show i) (xs !! (apply p i)) (ys !! i)
-
-testInvert :: Int -> Test
-testInvert n = TestCase $ do
-    x <- newRandomArray n :: IO (Ptr Double)
-    y <- newArrayCopy n x
-    p <- newRandomPerm n
-    invertWith  (swapElems y) p
-    
-    xs <- peekArray n x
-    ys <- peekArray n y
-
-    forM_ [0..(n-1)] $ \i ->
-        assertEqual ("position " ++ show i) (xs !! i) (ys !! (apply p i))
-
-
-testApplyThenInv :: Int -> Test
-testApplyThenInv n = TestCase $ do
-    x <- newRandomArray n :: IO (Ptr Double)
-    y <- newArrayCopy n x
-    p <- newRandomPerm n
-    applyWith  (swapElems y) p
-    invertWith (swapElems y) p
-    
-    xs <- peekArray n x
-    ys <- peekArray n y
-    
-    assertEqual "" xs ys
-
-testInvThenApply :: Int -> Test
-testInvThenApply n = TestCase $ do
-    x <- newRandomArray n :: IO (Ptr Double)
-    y <- newArrayCopy n x
-    p <- newRandomPerm n
-    invertWith (swapElems y) p
-    applyWith  (swapElems y) p
-    
-    xs <- peekArray n x
-    ys <- peekArray n y
-    
-    assertEqual "" xs ys
-
-
-
- 
-smallTests = [ TestLabel "identity"             testIdentity
-             , TestLabel "small cycle"          testSmallCycle
-             , TestLabel "small cycle inverse"  testSmallCycleInv
-             , TestLabel "single cycle"         testSingleCycle
-             , TestLabel "single cycle inverse" testSingleCycleInv
-             , TestLabel "multi cycle"          testMultiCycle
-             , TestLabel "multi cycle inverse"  testMultiCycleInv
-             ]
-
- 
-main =
-    let ns      = [ 0, 1, 2, 3, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 ]
-        aTests  = map (\n -> TestLabel ("apply (" ++ show n ++ ")") 
-                                       (testApply n)) ns
-        iTests  = map (\n -> TestLabel ("apply (" ++ show n ++ ")") 
-                                       (testInvert n)) ns
-        aiTests = map (\n -> TestLabel ("apply then inverse (" ++ show n ++ ")") 
-                                       (testApplyThenInv n)) ns
-        iaTests = map (\n -> TestLabel ("inverse then apply (" ++ show n ++ ")") 
-                                       (testInvThenApply n)) ns
-        tests   = TestList $ smallTests ++ aTests ++ aiTests ++ iaTests
-    in runTestTT tests
-