diff --git a/Data/Permutation.hs b/Data/Permutation.hs
new file mode 100644
--- /dev/null
+++ b/Data/Permutation.hs
@@ -0,0 +1,265 @@
+{-# 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/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) Patrick Perry <patperry@stanford.edu> 2008
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,14 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> import System.Cmd
+> import System.Exit ( ExitCode(..) )
+>
+> testing _ _ _ _ = do
+>     err <- system "make -C tests"
+>     system "make -s -C tests clean"
+>     if err /= ExitSuccess
+>         then ioError $ userError $ "failed"
+>         else return ()
+>
+> main = defaultMainWithHooks defaultUserHooks
+>        {runTests=testing}
diff --git a/permutation.cabal b/permutation.cabal
new file mode 100644
--- /dev/null
+++ b/permutation.cabal
@@ -0,0 +1,26 @@
+name:            permutation
+version:         0.1
+homepage:        http://stat.stanford.edu/~patperry/code/permutation
+synopsis:        Permutations and combinations
+description:
+    A library for representing and applying permutations.
+    .
+category:        Math
+license:         BSD3
+license-file:    LICENSE
+copyright:       (c) 2008. Patrick Perry <patperry@stanford.edu>
+author:          Patrick Perry
+maintainer:      Patrick Perry <patperry@stanford.edu>
+cabal-version: >= 1.2.0
+build-type:      Custom
+tested-with:     GHC ==6.8.2
+
+extra-source-files:  tests/Unit.hs
+                     tests/Makefile
+
+library
+    exposed-modules: Data.Permutation
+
+    build-depends:   base, containers
+
+    ghc-options:     -Wall
diff --git a/tests/Makefile b/tests/Makefile
new file mode 100644
--- /dev/null
+++ b/tests/Makefile
@@ -0,0 +1,8 @@
+all:
+	ghc -O -i. -i.. Unit.hs --make
+	./Unit
+
+clean:
+	rm -f *~ Unit *.hi *.o \
+		../Data/Permutation.hi ../Data/Permutation.o \
+		../System/Random/Permutation.o
diff --git a/tests/Unit.hs b/tests/Unit.hs
new file mode 100644
--- /dev/null
+++ b/tests/Unit.hs
@@ -0,0 +1,214 @@
+
+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
+