diff --git a/Data/Vector/Algorithms/AmericanFlag.hs b/Data/Vector/Algorithms/AmericanFlag.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/AmericanFlag.hs
+++ /dev/null
@@ -1,337 +0,0 @@
-{-# LANGUAGE FlexibleContexts, ScopedTypeVariables #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.AmericanFlag
--- Copyright   : (c) 2011 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (FlexibleContexts, ScopedTypeVariables)
---
--- This module implements American flag sort: an in-place, unstable, bucket
--- sort. Also in contrast to radix sort, the values are inspected in a big
--- endian order, and buckets are sorted via recursive splitting. This,
--- however, makes it sensible for sorting strings in lexicographic order
--- (provided indexing is fast).
---
--- The algorithm works as follows: at each stage, the array is looped over,
--- counting the number of elements for each bucket. Then, starting at the
--- beginning of the array, elements are permuted in place to reside in the
--- proper bucket, following chains until they reach back to the current
--- base index. Finally, each bucket is sorted recursively. This lends itself
--- well to the aforementioned variable-length strings, and so the algorithm
--- takes a stopping predicate, which is given a representative of the stripe,
--- rather than running for a set number of iterations.
-
-module Data.Vector.Algorithms.AmericanFlag ( sort
-                                           , sortBy
-                                           , Lexicographic(..)
-                                           ) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad
-import Control.Monad.Primitive
-
-import Data.Word
-import Data.Int
-import Data.Bits
-
-import qualified Data.ByteString as B
-
-import Data.Vector.Generic.Mutable
-import qualified Data.Vector.Primitive.Mutable as PV
-
-import qualified Data.Vector.Unboxed.Mutable as U
-
-import Data.Vector.Algorithms.Common
-
-import qualified Data.Vector.Algorithms.Insertion as I
-
--- | The methods of this class specify the information necessary to sort
--- arrays using the default ordering. The name 'Lexicographic' is meant
--- to convey that index should return results in a similar way to indexing
--- into a string.
-class Lexicographic e where
-  -- | Given a representative of a stripe and an index number, this
-  -- function should determine whether to stop sorting.
-  terminate :: e -> Int -> Bool
-  -- | The size of the bucket array necessary for sorting es
-  size      :: e -> Int
-  -- | Determines which bucket a given element should inhabit for a
-  -- particular iteration.
-  index     :: Int -> e -> Int
-
-instance Lexicographic Word8 where
-  terminate _ n = n > 0
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index _ n = fromIntegral n
-  {-# INLINE index #-}
-
-instance Lexicographic Word16 where
-  terminate _ n = n > 1
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ (n `shiftR`  8) .&. 255
-  index 1 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Word32 where
-  terminate _ n = n > 3
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ (n `shiftR` 24) .&. 255
-  index 1 n = fromIntegral $ (n `shiftR` 16) .&. 255
-  index 2 n = fromIntegral $ (n `shiftR`  8) .&. 255
-  index 3 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Word64 where
-  terminate _ n = n > 7
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ (n `shiftR` 56) .&. 255
-  index 1 n = fromIntegral $ (n `shiftR` 48) .&. 255
-  index 2 n = fromIntegral $ (n `shiftR` 40) .&. 255
-  index 3 n = fromIntegral $ (n `shiftR` 32) .&. 255
-  index 4 n = fromIntegral $ (n `shiftR` 24) .&. 255
-  index 5 n = fromIntegral $ (n `shiftR` 16) .&. 255
-  index 6 n = fromIntegral $ (n `shiftR`  8) .&. 255
-  index 7 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Word where
-  terminate _ n = n > 7
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ (n `shiftR` 56) .&. 255
-  index 1 n = fromIntegral $ (n `shiftR` 48) .&. 255
-  index 2 n = fromIntegral $ (n `shiftR` 40) .&. 255
-  index 3 n = fromIntegral $ (n `shiftR` 32) .&. 255
-  index 4 n = fromIntegral $ (n `shiftR` 24) .&. 255
-  index 5 n = fromIntegral $ (n `shiftR` 16) .&. 255
-  index 6 n = fromIntegral $ (n `shiftR`  8) .&. 255
-  index 7 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Int8 where
-  terminate _ n = n > 0
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index _ n = 255 .&. fromIntegral n `xor` 128
-  {-# INLINE index #-}
-
-instance Lexicographic Int16 where
-  terminate _ n = n > 1
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ ((n `xor` minBound) `shiftR` 8) .&. 255
-  index 1 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Int32 where
-  terminate _ n = n > 3
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ ((n `xor` minBound) `shiftR` 24) .&. 255
-  index 1 n = fromIntegral $ (n `shiftR` 16) .&. 255
-  index 2 n = fromIntegral $ (n `shiftR`  8) .&. 255
-  index 3 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Int64 where
-  terminate _ n = n > 7
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = fromIntegral $ ((n `xor` minBound) `shiftR` 56) .&. 255
-  index 1 n = fromIntegral $ (n `shiftR` 48) .&. 255
-  index 2 n = fromIntegral $ (n `shiftR` 40) .&. 255
-  index 3 n = fromIntegral $ (n `shiftR` 32) .&. 255
-  index 4 n = fromIntegral $ (n `shiftR` 24) .&. 255
-  index 5 n = fromIntegral $ (n `shiftR` 16) .&. 255
-  index 6 n = fromIntegral $ (n `shiftR`  8) .&. 255
-  index 7 n = fromIntegral $ n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic Int where
-  terminate _ n = n > 7
-  {-# INLINE terminate #-}
-  size _ = 256
-  {-# INLINE size #-}
-  index 0 n = ((n `xor` minBound) `shiftR` 56) .&. 255
-  index 1 n = (n `shiftR` 48) .&. 255
-  index 2 n = (n `shiftR` 40) .&. 255
-  index 3 n = (n `shiftR` 32) .&. 255
-  index 4 n = (n `shiftR` 24) .&. 255
-  index 5 n = (n `shiftR` 16) .&. 255
-  index 6 n = (n `shiftR`  8) .&. 255
-  index 7 n = n .&. 255
-  index _ _ = 0
-  {-# INLINE index #-}
-
-instance Lexicographic B.ByteString where
-  terminate b i = i >= B.length b
-  {-# INLINE terminate #-}
-  size _ = 257
-  {-# INLINE size #-}
-  index i b
-    | i >= B.length b = 0
-    | otherwise       = fromIntegral (B.index b i) + 1
-  {-# INLINE index #-}
-
--- | Sorts an array using the default ordering. Both Lexicographic and
--- Ord are necessary because the algorithm falls back to insertion sort
--- for sufficiently small arrays.
-sort :: forall e m v. (PrimMonad m, MVector v e, Lexicographic e, Ord e)
-     => v (PrimState m) e -> m ()
-sort v = sortBy compare terminate (size e) index v
- where e :: e
-       e = undefined
-{-# INLINABLE sort #-}
-
--- | A fully parameterized version of the sorting algorithm. Again, this
--- function takes both radix information and a comparison, because the
--- algorithms falls back to insertion sort for small arrays.
-sortBy :: (PrimMonad m, MVector v e)
-       => Comparison e       -- ^ a comparison for the insertion sort flalback
-       -> (e -> Int -> Bool) -- ^ determines whether a stripe is complete
-       -> Int                -- ^ the number of buckets necessary
-       -> (Int -> e -> Int)  -- ^ the big-endian radix function
-       -> v (PrimState m) e  -- ^ the array to be sorted
-       -> m ()
-sortBy cmp stop buckets radix v
-  | length v == 0 = return ()
-  | otherwise     = do count <- new buckets
-                       pile <- new buckets
-                       countLoop (radix 0) v count
-                       flagLoop cmp stop radix count pile v
-{-# INLINE sortBy #-}
-
-flagLoop :: (PrimMonad m, MVector v e)
-         => Comparison e
-         -> (e -> Int -> Bool)           -- number of passes
-         -> (Int -> e -> Int)            -- radix function
-         -> PV.MVector (PrimState m) Int -- auxiliary count array
-         -> PV.MVector (PrimState m) Int -- auxiliary pile array
-         -> v (PrimState m) e            -- source array
-         -> m ()
-flagLoop cmp stop radix count pile v = go 0 v
- where
-
- go pass v = do e <- unsafeRead v 0
-                unless (stop e $ pass - 1) $ go' pass v
-
- go' pass v
-   | len < threshold = I.sortByBounds cmp v 0 len
-   | otherwise       = do accumulate count pile
-                          permute (radix pass) count pile v
-                          recurse 0
-  where
-  len = length v
-  ppass = pass + 1
-
-  recurse i
-    | i < len   = do j <- countStripe (radix ppass) (radix pass) count v i
-                     go ppass (unsafeSlice i (j - i) v)
-                     recurse j
-    | otherwise = return ()
-{-# INLINE flagLoop #-}
-
-accumulate :: (PrimMonad m)
-           => PV.MVector (PrimState m) Int
-           -> PV.MVector (PrimState m) Int
-           -> m ()
-accumulate count pile = loop 0 0
- where
- len = length count
-
- loop i acc
-   | i < len = do ci <- unsafeRead count i
-                  let acc' = acc + ci
-                  unsafeWrite pile i acc
-                  unsafeWrite count i acc'
-                  loop (i+1) acc'
-   | otherwise    = return ()
-{-# INLINE accumulate #-}
-
-permute :: (PrimMonad m, MVector v e)
-        => (e -> Int)                       -- radix function
-        -> PV.MVector (PrimState m) Int     -- count array
-        -> PV.MVector (PrimState m) Int     -- pile array
-        -> v (PrimState m) e                -- source array
-        -> m ()
-permute rdx count pile v = go 0
- where
- len = length v
-
- go i
-   | i < len   = do e <- unsafeRead v i
-                    let r = rdx e
-                    p <- unsafeRead pile r
-                    m <- if r > 0
-                            then unsafeRead count (r-1)
-                            else return 0
-                    case () of
-                      -- if the current element is alunsafeReady in the right pile,
-                      -- go to the end of the pile
-                      _ | m <= i && i < p  -> go p
-                      -- if the current element happens to be in the right
-                      -- pile, bump the pile counter and go to the next element
-                        | i == p           -> unsafeWrite pile r (p+1) >> go (i+1)
-                      -- otherwise follow the chain
-                        | otherwise        -> follow i e p >> go (i+1)
-   | otherwise = return ()
- 
- follow i e j = do en <- unsafeRead v j
-                   let r = rdx en
-                   p <- inc pile r
-                   if p == j
-                      -- if the target happens to be in the right pile, don't move it.
-                      then follow i e (j+1)
-                      else unsafeWrite v j e >> if i == p
-                                             then unsafeWrite v i en
-                                             else follow i en p
-{-# INLINE permute #-}
-
-countStripe :: (PrimMonad m, MVector v e)
-            => (e -> Int)                   -- radix function
-            -> (e -> Int)                   -- stripe function
-            -> PV.MVector (PrimState m) Int -- count array
-            -> v (PrimState m) e            -- source array
-            -> Int                          -- starting position
-            -> m Int                        -- end of stripe: [lo,hi)
-countStripe rdx str count v lo = do set count 0
-                                    e <- unsafeRead v lo
-                                    go (str e) e (lo+1)
- where
- len = length v
-
- go !s e i = inc count (rdx e) >>
-            if i < len
-               then do en <- unsafeRead v i
-                       if str en == s
-                          then go s en (i+1)
-                          else return i
-                else return len
-{-# INLINE countStripe #-}
-
-threshold :: Int
-threshold = 25
-
diff --git a/Data/Vector/Algorithms/Combinators.hs b/Data/Vector/Algorithms/Combinators.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Combinators.hs
+++ /dev/null
@@ -1,71 +0,0 @@
-{-# LANGUAGE Rank2Types, TypeOperators #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Combinators
--- Copyright   : (c) 2008-2010 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (rank-2 types)
---
--- The purpose of this module is to supply various combinators for commonly
--- used idioms for the algorithms in this package. Examples at the time of
--- this writing include running an algorithm keyed on some function of the
--- elements (but only computing said function once per element), and safely
--- applying the algorithms on mutable arrays to immutable arrays.
-
-module Data.Vector.Algorithms.Combinators
-       (
---       , usingKeys
---       , usingIxKeys
-       ) where
-
-import Prelude hiding (length)
-
-import Control.Monad.ST
-
-import Data.Ord
-
-import Data.Vector.Generic
-
-import qualified Data.Vector.Generic.Mutable as M
-import qualified Data.Vector.Generic.New     as N
-
-{-
--- | Uses a function to compute a key for each element which the
--- algorithm should use in lieu of the actual element. For instance:
---
--- > usingKeys sortBy f arr
---
--- should produce the same results as:
---
--- > sortBy (comparing f) arr
---
--- the difference being that usingKeys computes each key only once
--- which can be more efficient for expensive key functions.
-usingKeys :: (UA e, UA k, Ord k)
-          => (forall e'. (UA e') => Comparison e' -> MUArr e' s -> ST s ())
-          -> (e -> k)
-          -> MUArr e s
-          -> ST s ()
-usingKeys algo f arr = usingIxKeys algo (const f) arr
-{-# INLINE usingKeys #-}
-
--- | As usingKeys, only the key function has access to the array index
--- at which each element is stored.
-usingIxKeys :: (UA e, UA k, Ord k)
-            => (forall e'. (UA e') => Comparison e' -> MUArr e' s -> ST s ())
-            -> (Int -> e -> k)
-            -> MUArr e s
-            -> ST s ()
-usingIxKeys algo f arr = do
-  keys <- newMU (lengthMU arr)
-  fill len keys
-  algo (comparing fstS) (unsafeZipMU keys arr)
- where
- len = lengthMU arr
- fill k keys
-   | k < 0     = return ()
-   | otherwise = readMU arr k >>= writeMU keys k . f k >> fill (k-1) keys
-{-# INLINE usingIxKeys #-}
--}
diff --git a/Data/Vector/Algorithms/Common.hs b/Data/Vector/Algorithms/Common.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Common.hs
+++ /dev/null
@@ -1,47 +0,0 @@
-{-# LANGUAGE FlexibleContexts #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Common
--- Copyright   : (c) 2008-2011 Dan Doel
--- Maintainer  : Dan Doel
--- Stability   : Experimental
--- Portability : Portable
---
--- Common operations and utility functions for all sorts
-
-module Data.Vector.Algorithms.Common where
-
-import Prelude hiding (read, length)
-
-import Control.Monad.Primitive
-
-import Data.Vector.Generic.Mutable
-
-import qualified Data.Vector.Primitive.Mutable as PV
-
--- | A type of comparisons between two values of a given type.
-type Comparison e = e -> e -> Ordering
-
-copyOffset :: (PrimMonad m, MVector v e)
-           => v (PrimState m) e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-copyOffset from to iFrom iTo len =
-  unsafeCopy (unsafeSlice iTo len to) (unsafeSlice iFrom len from)
-{-# INLINE copyOffset #-}
-
-inc :: (PrimMonad m, MVector v Int) => v (PrimState m) Int -> Int -> m Int
-inc arr i = unsafeRead arr i >>= \e -> unsafeWrite arr i (e+1) >> return e
-{-# INLINE inc #-}
-
--- shared bucket sorting stuff
-countLoop :: (PrimMonad m, MVector v e)
-          => (e -> Int)
-          -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m ()
-countLoop rdx src count = set count 0 >> go 0
- where
- len = length src
- go i
-   | i < len    = unsafeRead src i >>= inc count . rdx >> go (i+1)
-   | otherwise  = return ()
-{-# INLINE countLoop #-}
-
diff --git a/Data/Vector/Algorithms/Heap.hs b/Data/Vector/Algorithms/Heap.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Heap.hs
+++ /dev/null
@@ -1,240 +0,0 @@
-{-# LANGUAGE TypeOperators #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Heap
--- Copyright   : (c) 2008-2011 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (type operators)
---
--- This module implements operations for working with a quaternary heap stored
--- in an unboxed array. Most heapsorts are defined in terms of a binary heap,
--- in which each internal node has at most two children. By contrast, a
--- quaternary heap has internal nodes with up to four children. This reduces
--- the number of comparisons in a heapsort slightly, and improves locality
--- (again, slightly) by flattening out the heap.
-
-module Data.Vector.Algorithms.Heap
-       ( -- * Sorting
-         sort
-       , sortBy
-       , sortByBounds
-         -- * Selection
-       , select
-       , selectBy
-       , selectByBounds
-         -- * Partial sorts
-       , partialSort
-       , partialSortBy
-       , partialSortByBounds
-         -- * Heap operations
-       , heapify
-       , pop
-       , popTo
-       , sortHeap
-       , Comparison
-       ) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad
-import Control.Monad.Primitive
-
-import Data.Bits
-
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common (Comparison)
-
-import qualified Data.Vector.Algorithms.Optimal as O
-
--- | Sorts an entire array using the default ordering.
-sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
-sort = sortBy compare
-{-# INLINABLE sort #-}
-
--- | Sorts an entire array using a custom ordering.
-sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
-sortBy cmp a = sortByBounds cmp a 0 (length a)
-{-# INLINE sortBy #-}
-
--- | Sorts a portion of an array [l,u) using a custom ordering
-sortByBounds :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
-sortByBounds cmp a l u
-  | len < 2   = return ()
-  | len == 2  = O.sort2ByOffset cmp a l
-  | len == 3  = O.sort3ByOffset cmp a l
-  | len == 4  = O.sort4ByOffset cmp a l
-  | otherwise = heapify cmp a l u >> sortHeap cmp a l (l+4) u >> O.sort4ByOffset cmp a l
- where len = u - l
-{-# INLINE sortByBounds #-}
-
--- | Moves the lowest k elements to the front of the array.
--- The elements will be in no particular order.
-select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
-select = selectBy compare
-{-# INLINE select #-}
-
--- | Moves the 'lowest' (as defined by the comparison) k elements
--- to the front of the array. The elements will be in no particular
--- order.
-selectBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> Int -> m ()
-selectBy cmp a k = selectByBounds cmp a k 0 (length a)
-{-# INLINE selectBy #-}
-
--- | Moves the 'lowest' k elements in the portion [l,u) of the
--- array into the positions [l,k+l). The elements will be in
--- no particular order.
-selectByBounds :: (PrimMonad m, MVector v e)
-               => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-selectByBounds cmp a k l u
-  | l + k <= u = heapify cmp a l (l + k) >> go l (l + k) (u - 1)
-  | otherwise  = return ()
- where
- go l m u
-   | u < m      = return ()
-   | otherwise  = do el <- unsafeRead a l
-                     eu <- unsafeRead a u
-                     case cmp eu el of
-                       LT -> popTo cmp a l m u
-                       _  -> return ()
-                     go l m (u - 1)
-{-# INLINE selectByBounds #-}
-
--- | Moves the lowest k elements to the front of the array, sorted.
-partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
-partialSort = partialSortBy compare
-{-# INLINE partialSort #-}
-
--- | Moves the lowest k elements (as defined by the comparison) to
--- the front of the array, sorted.
-partialSortBy :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> m ()
-partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a)
-{-# INLINE partialSortBy #-}
-
--- | Moves the lowest k elements in the portion [l,u) of the array
--- into positions [l,k+l), sorted.
-partialSortByBounds :: (PrimMonad m, MVector v e)
-                    => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-partialSortByBounds cmp a k l u
-  -- this potentially does more work than absolutely required,
-  -- but using a heap to find the least 2 of 4 elements
-  -- seems unlikely to be better than just sorting all of them
-  -- with an optimal sort, and the latter is obviously index
-  -- correct.
-  | len <  2   = return ()
-  | len == 2   = O.sort2ByOffset cmp a l
-  | len == 3   = O.sort3ByOffset cmp a l
-  | len == 4   = O.sort4ByOffset cmp a l
-  | u <= l + k = sortByBounds cmp a l u
-  | otherwise  = do selectByBounds cmp a k l u
-                    sortHeap cmp a l (l + 4) (l + k)
-                    O.sort4ByOffset cmp a l
- where
- len = u - l
-{-# INLINE partialSortByBounds #-}
-
--- | Constructs a heap in a portion of an array [l, u)
-heapify :: (PrimMonad m, MVector v e)
-        => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
-heapify cmp a l u = loop $ (len - 1) `shiftR` 2
-  where
- len = u - l
- loop k
-   | k < 0     = return ()
-   | otherwise = unsafeRead a (l+k) >>= \e ->
-                   siftByOffset cmp a e l k len >> loop (k - 1)
-{-# INLINE heapify #-}
-
--- | Given a heap stored in a portion of an array [l,u), swaps the
--- top of the heap with the element at u and rebuilds the heap.
-pop :: (PrimMonad m, MVector v e)
-    => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
-pop cmp a l u = popTo cmp a l u u
-{-# INLINE pop #-}
-
--- | Given a heap stored in a portion of an array [l,u) swaps the top
--- of the heap with the element at position t, and rebuilds the heap.
-popTo :: (PrimMonad m, MVector v e)
-      => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-popTo cmp a l u t = do al <- unsafeRead a l
-                       at <- unsafeRead a t
-                       unsafeWrite a t al
-                       siftByOffset cmp a at l 0 (u - l)
-{-# INLINE popTo #-}
-
--- | Given a heap stored in a portion of an array [l,u), sorts the
--- highest values into [m,u). The elements in [l,m) are not in any
--- particular order.
-sortHeap :: (PrimMonad m, MVector v e)
-         => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-sortHeap cmp a l m u = loop (u-1) >> unsafeSwap a l m
- where
- loop k
-   | m < k     = pop cmp a l k >> loop (k-1)
-   | otherwise = return ()
-{-# INLINE sortHeap #-}
-
--- Rebuilds a heap with a hole in it from start downwards. Afterward,
--- the heap property should apply for [start + off, len + off). val
--- is the new value to be put in the hole.
-siftByOffset :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> Int -> m ()
-siftByOffset cmp a val off start len = sift val start len
- where
- sift val root len
-   | child < len = do (child', ac) <- maximumChild cmp a off child len
-                      case cmp val ac of
-                        LT -> unsafeWrite a (root + off) ac >> sift val child' len
-                        _  -> unsafeWrite a (root + off) val
-   | otherwise = unsafeWrite a (root + off) val
-  where child = root `shiftL` 2 + 1
-{-# INLINE siftByOffset #-}
-
--- Finds the maximum child of a heap node, given the indx of the first child.
-maximumChild :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m (Int,  e)
-maximumChild cmp a off child1 len
-  | child4 < len = do ac1 <- unsafeRead a (child1 + off)
-                      ac2 <- unsafeRead a (child2 + off)
-                      ac3 <- unsafeRead a (child3 + off)
-                      ac4 <- unsafeRead a (child4 + off)
-                      return $ case cmp ac1 ac2 of
-                                 LT -> case cmp ac2 ac3 of
-                                         LT -> case cmp ac3 ac4 of
-                                                 LT -> (child4, ac4)
-                                                 _  -> (child3, ac3)
-                                         _  -> case cmp ac2 ac4 of
-                                                 LT -> (child4, ac4)
-                                                 _  -> (child2, ac2)
-                                 _  -> case cmp ac1 ac3 of
-                                         LT -> case cmp ac3 ac4 of
-                                                 LT -> (child4, ac4)
-                                                 _  -> (child3, ac3)
-                                         _  -> case cmp ac1 ac4 of
-                                                 LT -> (child4, ac4)
-                                                 _  -> (child1, ac1)
-  | child3 < len = do ac1 <- unsafeRead a (child1 + off)
-                      ac2 <- unsafeRead a (child2 + off)
-                      ac3 <- unsafeRead a (child3 + off)
-                      return $ case cmp ac1 ac2 of
-                                 LT -> case cmp ac2 ac3 of
-                                         LT -> (child3, ac3)
-                                         _  -> (child2, ac2)
-                                 _  -> case cmp ac1 ac3 of
-                                         LT -> (child3, ac3)
-                                         _  -> (child1, ac1)
-  | child2 < len = do ac1 <- unsafeRead a (child1 + off)
-                      ac2 <- unsafeRead a (child2 + off)
-                      return $ case cmp ac1 ac2 of
-                                 LT -> (child2, ac2)
-                                 _  -> (child1, ac1)
-  | otherwise    = do ac1 <- unsafeRead a (child1 + off) ; return (child1, ac1)
- where
- child2 = child1 + 1
- child3 = child1 + 2
- child4 = child1 + 3
-{-# INLINE maximumChild #-}
diff --git a/Data/Vector/Algorithms/Insertion.hs b/Data/Vector/Algorithms/Insertion.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Insertion.hs
+++ /dev/null
@@ -1,81 +0,0 @@
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Insertion
--- Copyright   : (c) 2008-2010 Dan Doel
--- Maintainer  : Dan Doel
--- Stability   : Experimental
--- Portability : Portable
---
--- A simple insertion sort. Though it's O(n^2), its iterative nature can be
--- beneficial for small arrays. It is used to sort small segments of an array
--- by some of the more heavy-duty, recursive algorithms.
-
-module Data.Vector.Algorithms.Insertion
-       ( sort
-       , sortBy
-       , sortByBounds
-       , sortByBounds'
-       , Comparison
-       ) where
-
-
-import Prelude hiding (read, length)
-
-import Control.Monad.Primitive
-
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common (Comparison)
-
-import qualified Data.Vector.Algorithms.Optimal as O
-
--- | Sorts an entire array using the default comparison for the type
-sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
-sort = sortBy compare
-{-# INLINABLE sort #-}
-
--- | Sorts an entire array using a given comparison
-sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
-sortBy cmp a = sortByBounds cmp a 0 (length a)
-{-# INLINE sortBy #-}
-
--- | Sorts the portion of an array delimited by [l,u)
-sortByBounds :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
-sortByBounds cmp a l u
-  | len < 2   = return ()
-  | len == 2  = O.sort2ByOffset cmp a l
-  | len == 3  = O.sort3ByOffset cmp a l
-  | len == 4  = O.sort4ByOffset cmp a l
-  | otherwise = O.sort4ByOffset cmp a l >> sortByBounds' cmp a l (l + 4) u
- where
- len = u - l
-{-# INLINE sortByBounds #-}
-
--- | Sorts the portion of the array delimited by [l,u) under the assumption
--- that [l,m) is already sorted.
-sortByBounds' :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-sortByBounds' cmp a l m u = sort m
- where
- sort i
-   | i < u     = do v <- unsafeRead a i
-                    insert cmp a l v i
-                    sort (i+1)
-   | otherwise = return ()
-{-# INLINE sortByBounds' #-}
-
--- Given a sorted array in [l,u), inserts val into its proper position,
--- yielding a sorted [l,u]
-insert :: (PrimMonad m, MVector v e)
-       => Comparison e -> v (PrimState m) e -> Int -> e -> Int -> m ()
-insert cmp a l = loop
- where
- loop val j
-   | j <= l    = unsafeWrite a l val
-   | otherwise = do e <- unsafeRead a (j - 1)
-                    case cmp val e of
-                      LT -> unsafeWrite a j e >> loop val (j - 1)
-                      _  -> unsafeWrite a j val
-{-# INLINE insert #-}
diff --git a/Data/Vector/Algorithms/Intro.hs b/Data/Vector/Algorithms/Intro.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Intro.hs
+++ /dev/null
@@ -1,211 +0,0 @@
-{-# LANGUAGE TypeOperators, BangPatterns, ScopedTypeVariables #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Intro
--- Copyright   : (c) 2008-2011 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (type operators, bang patterns)
---
--- This module implements various algorithms based on the introsort algorithm,
--- originally described by David R. Musser in the paper /Introspective Sorting
--- and Selection Algorithms/. It is also in widespread practical use, as the
--- standard unstable sort used in the C++ Standard Template Library.
---
--- Introsort is at its core a quicksort. The version implemented here has the
--- following optimizations that make it perform better in practice:
---
---   * Small segments of the array are left unsorted until a final insertion
---     sort pass. This is faster than recursing all the way down to
---     one-element arrays.
---
---   * The pivot for segment [l,u) is chosen as the median of the elements at
---     l, u-1 and (u+l)/2. This yields good behavior on mostly sorted (or
---     reverse-sorted) arrays.
---
---   * The algorithm tracks its recursion depth, and if it decides it is
---     taking too long (depth greater than 2 * lg n), it switches to a heap
---     sort to maintain O(n lg n) worst case behavior. (This is what makes the
---     algorithm introsort).
-
-module Data.Vector.Algorithms.Intro
-       ( -- * Sorting
-         sort
-       , sortBy
-       , sortByBounds
-         -- * Selecting
-       , select
-       , selectBy
-       , selectByBounds
-         -- * Partial sorting
-       , partialSort
-       , partialSortBy
-       , partialSortByBounds
-       , Comparison
-       ) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad
-import Control.Monad.Primitive
-
-import Data.Bits
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common (Comparison)
-
-import qualified Data.Vector.Algorithms.Insertion as I
-import qualified Data.Vector.Algorithms.Optimal   as O
-import qualified Data.Vector.Algorithms.Heap      as H
-
--- | Sorts an entire array using the default ordering.
-sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
-sort = sortBy compare
-{-# INLINABLE sort #-}
-
--- | Sorts an entire array using a custom ordering.
-sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
-sortBy cmp a = sortByBounds cmp a 0 (length a)
-{-# INLINE sortBy #-}
-
--- | Sorts a portion of an array [l,u) using a custom ordering
-sortByBounds :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
-sortByBounds cmp a l u
-  | len < 2   = return ()
-  | len == 2  = O.sort2ByOffset cmp a l
-  | len == 3  = O.sort3ByOffset cmp a l
-  | len == 4  = O.sort4ByOffset cmp a l
-  | otherwise = introsort cmp a (ilg len) l u
- where len = u - l
-{-# INLINE sortByBounds #-}
-
--- Internal version of the introsort loop which allows partial
--- sort functions to call with a specified bound on iterations.
-introsort :: (PrimMonad m, MVector v e)
-          => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-introsort cmp a i l u = sort i l u >> I.sortByBounds cmp a l u
- where
- sort 0 l u = H.sortByBounds cmp a l u
- sort d l u
-   | len < threshold = return ()
-   | otherwise = do O.sort3ByIndex cmp a c l (u-1) -- sort the median into the lowest position
-                    p <- unsafeRead a l
-                    mid <- partitionBy cmp a p (l+1) u
-                    unsafeSwap a l (mid - 1)
-                    sort (d-1) mid u
-                    sort (d-1) l   (mid - 1)
-  where
-  len = u - l
-  c   = (u + l) `div` 2
-{-# INLINE introsort #-}
-
--- | Moves the least k elements to the front of the array in
--- no particular order.
-select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
-select = selectBy compare
-{-# INLINE select #-}
-
--- | Moves the least k elements (as defined by the comparison) to
--- the front of the array in no particular order.
-selectBy :: (PrimMonad m, MVector v e)
-         => Comparison e -> v (PrimState m) e -> Int -> m ()
-selectBy cmp a k = selectByBounds cmp a k 0 (length a)
-{-# INLINE selectBy #-}
-
--- | Moves the least k elements in the interval [l,u) to the positions
--- [l,k+l) in no particular order.
-selectByBounds :: (PrimMonad m, MVector v e)
-               => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-selectByBounds cmp a k l u
-  | l >= u    = return ()
-  | otherwise = go (ilg len) l (l + k) u
- where
- len = u - l
- go 0 l m u = H.selectByBounds cmp a (m - l) l u
- go n l m u = do O.sort3ByIndex cmp a c l (u-1)
-                 p <- unsafeRead a l
-                 mid <- partitionBy cmp a p (l+1) u
-                 unsafeSwap a l (mid - 1)
-                 if m > mid
-                   then go (n-1) mid m u
-                   else if m < mid - 1
-                        then go (n-1) l m (mid - 1)
-                        else return ()
-  where c = (u + l) `div` 2
-{-# INLINE selectByBounds #-}
-
--- | Moves the least k elements to the front of the array, sorted.
-partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
-partialSort = partialSortBy compare
-{-# INLINE partialSort #-}
-
--- | Moves the least k elements (as defined by the comparison) to
--- the front of the array, sorted.
-partialSortBy :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> m ()
-partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a)
-{-# INLINE partialSortBy #-}
-
--- | Moves the least k elements in the interval [l,u) to the positions
--- [l,k+l), sorted.
-partialSortByBounds :: (PrimMonad m, MVector v e)
-                    => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-partialSortByBounds cmp a k l u
-  | l >= u    = return ()
-  | otherwise = go (ilg len) l (l + k) u
- where
- isort = introsort cmp a
- {-# INLINE [1] isort #-}
- len = u - l
- go 0 l m n = H.partialSortByBounds cmp a (m - l) l u
- go n l m u
-   | l == m    = return ()
-   | otherwise = do O.sort3ByIndex cmp a c l (u-1)
-                    p <- unsafeRead a l
-                    mid <- partitionBy cmp a p (l+1) u
-                    unsafeSwap a l (mid - 1)
-                    case compare m mid of
-                      GT -> do isort (n-1) l (mid - 1)
-                               go (n-1) mid m u
-                      EQ -> isort (n-1) l m
-                      LT -> go n l m (mid - 1)
-  where c = (u + l) `div` 2
-{-# INLINE partialSortByBounds #-}
-
-partitionBy :: forall m v e. (PrimMonad m, MVector v e)
-            => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
-partitionBy cmp a = partUp
- where
- -- 6.10 panics without the signatures for partUp and partDown, 6.12 and later
- -- versions don't need them
- partUp :: e -> Int -> Int -> m Int
- partUp p l u
-   | l < u = do e <- unsafeRead a l
-                case cmp e p of
-                  LT -> partUp p (l+1) u
-                  _  -> partDown p l (u-1)
-   | otherwise = return l
-
- partDown :: e -> Int -> Int -> m Int
- partDown p l u
-   | l < u = do e <- unsafeRead a u
-                case cmp p e of
-                  LT -> partDown p l (u-1)
-                  _  -> unsafeSwap a l u >> partUp p (l+1) u
-   | otherwise = return l
-{-# INLINE partitionBy #-}
-
--- computes the number of recursive calls after which heapsort should
--- be invoked given the lower and upper indices of the array to be sorted
-ilg :: Int -> Int
-ilg m = 2 * loop m 0
- where
- loop 0 !k = k - 1
- loop n !k = loop (n `shiftR` 1) (k+1)
-
--- the size of array at which the introsort algorithm switches to insertion sort
-threshold :: Int
-threshold = 18
-{-# INLINE threshold #-}
diff --git a/Data/Vector/Algorithms/Merge.hs b/Data/Vector/Algorithms/Merge.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Merge.hs
+++ /dev/null
@@ -1,95 +0,0 @@
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Merge
--- Copyright   : (c) 2008-2011 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Portable
---
--- This module implements a simple top-down merge sort. The temporary buffer
--- is preallocated to 1/2 the size of the input array, and shared through
--- the entire sorting process to ease the amount of allocation performed in
--- total. This is a stable sort.
-
-module Data.Vector.Algorithms.Merge
-       ( sort
-       , sortBy
-       , Comparison
-       ) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad.Primitive
-
-import Data.Bits
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common (Comparison, copyOffset)
-
-import qualified Data.Vector.Algorithms.Optimal   as O
-import qualified Data.Vector.Algorithms.Insertion as I
-
--- | Sorts an array using the default comparison.
-sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
-sort = sortBy compare
-{-# INLINABLE sort #-}
-
--- | Sorts an array using a custom comparison.
-sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
-sortBy cmp vec
-  | len <= 1  = return ()
-  | len == 2  = O.sort2ByOffset cmp vec 0
-  | len == 3  = O.sort3ByOffset cmp vec 0
-  | len == 4  = O.sort4ByOffset cmp vec 0
-  | otherwise = do buf <- new len
-                   mergeSortWithBuf cmp vec buf
- where
- len = length vec
-{-# INLINE sortBy #-}
-
-mergeSortWithBuf :: (PrimMonad m, MVector v e)
-                 => Comparison e -> v (PrimState m) e -> v (PrimState m) e -> m ()
-mergeSortWithBuf cmp src buf = loop 0 (length src)
- where
- loop l u
-   | len < threshold = I.sortByBounds cmp src l u
-   | otherwise       = do loop l mid
-                          loop mid u
-                          merge cmp (unsafeSlice l len src) buf (mid - l)
-  where len = u - l
-        mid = (u + l) `shiftR` 1
-{-# INLINE mergeSortWithBuf #-}
-
-merge :: (PrimMonad m, MVector v e)
-      => Comparison e -> v (PrimState m) e -> v (PrimState m) e
-      -> Int -> m ()
-merge cmp src buf mid = do unsafeCopy tmp lower
-                           eTmp <- unsafeRead tmp 0
-                           eUpp <- unsafeRead upper 0
-                           loop tmp 0 eTmp upper 0 eUpp 0
- where
- lower = unsafeSlice 0   mid                src
- upper = unsafeSlice mid (length src - mid) src
- tmp   = unsafeSlice 0   mid                buf
-
- wroteHigh low iLow eLow high iHigh iIns
-   | iHigh >= length high = unsafeCopy (unsafeSlice iIns (length low - iLow) src)
-                                       (unsafeSlice iLow (length low - iLow) low)
-   | otherwise            = do eHigh <- unsafeRead high iHigh
-                               loop low iLow eLow high iHigh eHigh iIns
-
- wroteLow low iLow high iHigh eHigh iIns
-   | iLow  >= length low  = return ()
-   | otherwise            = do eLow <- unsafeRead low iLow
-                               loop low iLow eLow high iHigh eHigh iIns
-
- loop !low !iLow !eLow !high !iHigh !eHigh !iIns = case cmp eHigh eLow of
-     LT -> do unsafeWrite src iIns eHigh
-              wroteHigh low iLow eLow high (iHigh + 1) (iIns + 1)
-     _  -> do unsafeWrite src iIns eLow
-              wroteLow low (iLow + 1) high iHigh eHigh (iIns + 1)
-{-# INLINE merge #-}
-
-threshold :: Int
-threshold = 25
-{-# INLINE threshold #-}
diff --git a/Data/Vector/Algorithms/Optimal.hs b/Data/Vector/Algorithms/Optimal.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Optimal.hs
+++ /dev/null
@@ -1,244 +0,0 @@
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Optimal
--- Copyright   : (c) 2008-2010 Dan Doel
--- Maintainer  : Dan Doel
--- Stability   : Experimental
--- Portability : Portable
---
--- Optimal sorts for very small array sizes, or for small numbers of
--- particular indices in a larger array (to be used, for instance, for
--- sorting a median of 3 values into the lowest position in an array
--- for a median-of-3 quicksort).
-
--- The code herein was adapted from a C algorithm for optimal sorts
--- of small arrays. The original code was produced for the article
--- /Sorting Revisited/ by Paul Hsieh, available here:
---
---   http://www.azillionmonkeys.com/qed/sort.html
---
--- The LICENSE file contains the relevant copyright information for
--- the reference C code.
-
-module Data.Vector.Algorithms.Optimal
-       ( sort2ByIndex
-       , sort2ByOffset
-       , sort3ByIndex
-       , sort3ByOffset
-       , sort4ByIndex
-       , sort4ByOffset
-       , Comparison
-       ) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad.Primitive
-
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common (Comparison)
-
-#include "vector.h"
-
--- | Sorts the elements at the positions 'off' and 'off + 1' in the given
--- array using the comparison.
-sort2ByOffset :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> m ()
-sort2ByOffset cmp a off = sort2ByIndex cmp a off (off + 1)
-{-# INLINABLE sort2ByOffset #-}
-
--- | Sorts the elements at the two given indices using the comparison. This
--- is essentially a compare-and-swap, although the first index is assumed to
--- be the 'lower' of the two.
-sort2ByIndex :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
-sort2ByIndex cmp a i j = UNSAFE_CHECK(checkIndex) "sort2ByIndex" i (length a)
-                       $ UNSAFE_CHECK(checkIndex) "sort2ByIndex" j (length a) $  do
-  a0 <- unsafeRead a i
-  a1 <- unsafeRead a j
-  case cmp a0 a1 of
-    GT -> unsafeWrite a i a1 >> unsafeWrite a j a0
-    _  -> return ()
-{-# INLINABLE sort2ByIndex #-}
-
--- | Sorts the three elements starting at the given offset in the array.
-sort3ByOffset :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> m ()
-sort3ByOffset cmp a off = sort3ByIndex cmp a off (off + 1) (off + 2)
-{-# INLINABLE sort3ByOffset #-}
-
--- | Sorts the elements at the three given indices. The indices are assumed
--- to be given from lowest to highest, so if 'l < m < u' then
--- 'sort3ByIndex cmp a m l u' essentially sorts the median of three into the
--- lowest position in the array.
-sort3ByIndex :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
-sort3ByIndex cmp a i j k = UNSAFE_CHECK(checkIndex) "sort3ByIndex" i (length a)
-                         $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" j (length a)
-                         $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" k (length a) $ do
-  a0 <- unsafeRead a i
-  a1 <- unsafeRead a j
-  a2 <- unsafeRead a k
-  case cmp a0 a1 of
-    GT -> case cmp a0 a2 of
-            GT -> case cmp a2 a1 of
-                    LT -> do unsafeWrite a i a2
-                             unsafeWrite a k a0
-                    _  -> do unsafeWrite a i a1
-                             unsafeWrite a j a2
-                             unsafeWrite a k a0
-            _  -> do unsafeWrite a i a1
-                     unsafeWrite a j a0
-    _  -> case cmp a1 a2 of
-            GT -> case cmp a0 a2 of
-                    GT -> do unsafeWrite a i a2
-                             unsafeWrite a j a0
-                             unsafeWrite a k a1
-                    _  -> do unsafeWrite a j a2
-                             unsafeWrite a k a1
-            _  -> return ()
-{-# INLINABLE sort3ByIndex #-}
-
--- | Sorts the four elements beginning at the offset.
-sort4ByOffset :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> m ()
-sort4ByOffset cmp a off = sort4ByIndex cmp a off (off + 1) (off + 2) (off + 3)
-{-# INLINABLE sort4ByOffset #-}
-
--- The horror...
-
--- | Sorts the elements at the four given indices. Like the 2 and 3 element
--- versions, this assumes that the indices are given in increasing order, so
--- it can be used to sort medians into particular positions and so on.
-sort4ByIndex :: (PrimMonad m, MVector v e)
-             => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> Int -> m ()
-sort4ByIndex cmp a i j k l = UNSAFE_CHECK(checkIndex) "sort4ByIndex" i (length a)
-                           $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" j (length a)
-                           $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" k (length a)
-                           $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" l (length a) $ do
-  a0 <- unsafeRead a i
-  a1 <- unsafeRead a j
-  a2 <- unsafeRead a k
-  a3 <- unsafeRead a l
-  case cmp a0 a1 of
-    GT -> case cmp a0 a2 of
-            GT -> case cmp a1 a2 of
-                    GT -> case cmp a1 a3 of
-                            GT -> case cmp a2 a3 of
-                                    GT -> do unsafeWrite a i a3
-                                             unsafeWrite a j a2
-                                             unsafeWrite a k a1
-                                             unsafeWrite a l a0
-                                    _  -> do unsafeWrite a i a2
-                                             unsafeWrite a j a3
-                                             unsafeWrite a k a1
-                                             unsafeWrite a l a0
-                            _  -> case cmp a0 a3 of
-                                    GT -> do unsafeWrite a i a2
-                                             unsafeWrite a j a1
-                                             unsafeWrite a k a3
-                                             unsafeWrite a l a0
-                                    _  -> do unsafeWrite a i a2
-                                             unsafeWrite a j a1
-                                             unsafeWrite a k a0
-                                             unsafeWrite a l a3
-                    _ -> case cmp a2 a3 of
-                           GT -> case cmp a1 a3 of
-                                   GT -> do unsafeWrite a i a3
-                                            unsafeWrite a j a1
-                                            unsafeWrite a k a2
-                                            unsafeWrite a l a0
-                                   _  -> do unsafeWrite a i a1
-                                            unsafeWrite a j a3
-                                            unsafeWrite a k a2
-                                            unsafeWrite a l a0
-                           _  -> case cmp a0 a3 of
-                                   GT -> do unsafeWrite a i a1
-                                            unsafeWrite a j a2
-                                            unsafeWrite a k a3
-                                            unsafeWrite a l a0
-                                   _  -> do unsafeWrite a i a1
-                                            unsafeWrite a j a2
-                                            unsafeWrite a k a0
-                                            -- unsafeWrite a l a3
-            _  -> case cmp a0 a3 of
-                    GT -> case cmp a1 a3 of
-                            GT -> do unsafeWrite a i a3
-                                     -- unsafeWrite a j a1
-                                     unsafeWrite a k a0
-                                     unsafeWrite a l a2
-                            _  -> do unsafeWrite a i a1
-                                     unsafeWrite a j a3
-                                     unsafeWrite a k a0
-                                     unsafeWrite a l a2
-                    _  -> case cmp a2 a3 of
-                            GT -> do unsafeWrite a i a1
-                                     unsafeWrite a j a0
-                                     unsafeWrite a k a3
-                                     unsafeWrite a l a2
-                            _  -> do unsafeWrite a i a1
-                                     unsafeWrite a j a0
-                                     -- unsafeWrite a k a2
-                                     -- unsafeWrite a l a3
-    _  -> case cmp a1 a2 of
-            GT -> case cmp a0 a2 of
-                    GT -> case cmp a0 a3 of
-                            GT -> case cmp a2 a3 of
-                                    GT -> do unsafeWrite a i a3
-                                             unsafeWrite a j a2
-                                             unsafeWrite a k a0
-                                             unsafeWrite a l a1
-                                    _  -> do unsafeWrite a i a2
-                                             unsafeWrite a j a3
-                                             unsafeWrite a k a0
-                                             unsafeWrite a l a1
-                            _  -> case cmp a1 a3 of
-                                    GT -> do unsafeWrite a i a2
-                                             unsafeWrite a j a0
-                                             unsafeWrite a k a3
-                                             unsafeWrite a l a1
-                                    _  -> do unsafeWrite a i a2
-                                             unsafeWrite a j a0
-                                             unsafeWrite a k a1
-                                             -- unsafeWrite a l a3
-                    _  -> case cmp a2 a3 of
-                            GT -> case cmp a0 a3 of
-                                    GT -> do unsafeWrite a i a3
-                                             unsafeWrite a j a0
-                                             -- unsafeWrite a k a2
-                                             unsafeWrite a l a1
-                                    _  -> do -- unsafeWrite a i a0
-                                             unsafeWrite a j a3
-                                             -- unsafeWrite a k a2
-                                             unsafeWrite a l a1
-                            _  -> case cmp a1 a3 of
-                                    GT -> do -- unsafeWrite a i a0
-                                             unsafeWrite a j a2
-                                             unsafeWrite a k a3
-                                             unsafeWrite a l a1
-                                    _  -> do -- unsafeWrite a i a0
-                                             unsafeWrite a j a2
-                                             unsafeWrite a k a1
-                                             -- unsafeWrite a l a3
-            _  -> case cmp a1 a3 of
-                    GT -> case cmp a0 a3 of
-                            GT -> do unsafeWrite a i a3
-                                     unsafeWrite a j a0
-                                     unsafeWrite a k a1
-                                     unsafeWrite a l a2
-                            _  -> do -- unsafeWrite a i a0
-                                     unsafeWrite a j a3
-                                     unsafeWrite a k a1
-                                     unsafeWrite a l a2
-                    _  -> case cmp a2 a3 of
-                            GT -> do -- unsafeWrite a i a0
-                                     -- unsafeWrite a j a1
-                                     unsafeWrite a k a3
-                                     unsafeWrite a l a2
-                            _  -> do -- unsafeWrite a i a0
-                                     -- unsafeWrite a j a1
-                                     -- unsafeWrite a k a2
-                                     -- unsafeWrite a l a3
-                                     return ()
-{-# INLINABLE sort4ByIndex #-}
diff --git a/Data/Vector/Algorithms/Radix.hs b/Data/Vector/Algorithms/Radix.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Radix.hs
+++ /dev/null
@@ -1,261 +0,0 @@
-{-# LANGUAGE ScopedTypeVariables, BangPatterns, TypeOperators #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Radix
--- Copyright   : (c) 2008-2011 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (scoped type variables, bang patterns)
---
--- This module provides a radix sort for a subclass of unboxed arrays. The
--- radix class gives information on
---   * the number of passes needed for the data type
---
---   * the size of the auxiliary arrays
---
---   * how to compute the pass-k radix of a value
---
--- Radix sort is not a comparison sort, so it is able to achieve O(n) run
--- time, though it also uses O(n) auxiliary space. In addition, there is a
--- constant space overhead of 2*size*sizeOf(Int) for the sort, so it is not
--- advisable to use this sort for large numbers of very small arrays.
---
--- A standard example (upon which one could base their own Radix instance)
--- is Word32:
---
---   * We choose to sort on r = 8 bits at a time
---
---   * A Word32 has b = 32 bits total
---
---   Thus, b/r = 4 passes are required, 2^r = 256 elements are needed in an
---   auxiliary array, and the radix function is:
---
---    > radix k e = (e `shiftR` (k*8)) .&. 256
-
-module Data.Vector.Algorithms.Radix (sort, sortBy, Radix(..)) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad
-import Control.Monad.Primitive
-
-import qualified Data.Vector.Primitive.Mutable as PV
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common
-
-import Data.Bits
-import Data.Int
-import Data.Word
-
-
-import Foreign.Storable
-
-class Radix e where
-  -- | The number of passes necessary to sort an array of es
-  passes :: e -> Int
-  -- | The size of an auxiliary array
-  size   :: e -> Int
-  -- | The radix function parameterized by the current pass
-  radix  :: Int -> e -> Int
-
-instance Radix Int where
-  passes _ = sizeOf (undefined :: Int)
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix 0 e = e .&. 255
-  radix i e
-    | i == passes e - 1 = radix' (e `xor` minBound)
-    | otherwise         = radix' e
-   where radix' e = (e `shiftR` (i `shiftL` 3)) .&. 255
-  {-# INLINE radix #-}
-
-instance Radix Int8 where
-  passes _ = 1
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix _ e = 255 .&. fromIntegral e `xor` 128
-  {-# INLINE radix #-}
-
-instance Radix Int16 where
-  passes _ = 2
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix 1 e = fromIntegral (((e `xor` minBound) `shiftR` 8) .&. 255)
-  {-# INLINE radix #-}
-
-instance Radix Int32 where
-  passes _ = 4
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
-  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
-  radix 3 e = fromIntegral (((e `xor` minBound) `shiftR` 24) .&. 255)
-  {-# INLINE radix #-}
-
-instance Radix Int64 where
-  passes _ = 8
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
-  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
-  radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)
-  radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255)
-  radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255)
-  radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255)
-  radix 7 e = fromIntegral (((e `xor` minBound) `shiftR` 56) .&. 255)
-  {-# INLINE radix #-}
-
-instance Radix Word where
-  passes _ = sizeOf (undefined :: Word)
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix i e = fromIntegral ((e `shiftR` (i `shiftL` 3)) .&. 255)
-  {-# INLINE radix #-}
-
-instance Radix Word8 where
-  passes _ = 1
-  {-# INLINE passes #-}
-  size _ = 256
-  {-# INLINE size #-}
-  radix _ = fromIntegral
-  {-# INLINE radix #-}
-
-instance Radix Word16 where
-  passes _ = 2
-  {-# INLINE passes #-}
-  size   _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
-  {-# INLINE radix #-}
-
-instance Radix Word32 where
-  passes _ = 4
-  {-# INLINE passes #-}
-  size   _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
-  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
-  radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)
-  {-# INLINE radix #-}
-
-instance Radix Word64 where
-  passes _ = 8
-  {-# INLINE passes #-}
-  size   _ = 256
-  {-# INLINE size #-}
-  radix 0 e = fromIntegral (e .&. 255)
-  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
-  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
-  radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)
-  radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255)
-  radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255)
-  radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255)
-  radix 7 e = fromIntegral ((e `shiftR` 56) .&. 255)
-  {-# INLINE radix #-}
-
-instance (Radix i, Radix j) => Radix (i, j) where
-  passes ~(i, j) = passes i + passes j
-  {-# INLINE passes #-}
-  size   ~(i, j) = size i `max` size j
-  {-# INLINE size #-}
-  radix k ~(i, j) | k < passes j = radix k j
-                     | otherwise    = radix (k - passes j) i
-  {-# INLINE radix #-}
-
--- | Sorts an array based on the Radix instance.
-sort :: forall e m v. (PrimMonad m, MVector v e, Radix e)
-     => v (PrimState m) e -> m ()
-sort arr = sortBy (passes e) (size e) radix arr
- where
- e :: e
- e = undefined
-{-# INLINABLE sort #-}
-
--- | Radix sorts an array using custom radix information
--- requires the number of passes to fully sort the array,
--- the size of of auxiliary arrays necessary (should be
--- one greater than the maximum value returned by the radix
--- function), and a radix function, which takes the pass
--- and an element, and returns the relevant radix.
-sortBy :: (PrimMonad m, MVector v e)
-       => Int               -- ^ the number of passes
-       -> Int               -- ^ the size of auxiliary arrays
-       -> (Int -> e -> Int) -- ^ the radix function
-       -> v (PrimState m) e -- ^ the array to be sorted
-       -> m ()
-sortBy passes size rdx arr = do
-  tmp    <- new (length arr)
-  count  <- new size
-  radixLoop passes rdx arr tmp count
-{-# INLINE sortBy #-}
-
-radixLoop :: (PrimMonad m, MVector v e)
-          => Int                          -- passes
-          -> (Int -> e -> Int)            -- radix function
-          -> v (PrimState m) e            -- array to sort
-          -> v (PrimState m) e            -- temporary array
-          -> PV.MVector (PrimState m) Int -- radix count array
-          -> m ()
-radixLoop passes rdx src dst count = go False 0
- where
- len = length src
- go swap k
-   | k < passes = if swap
-                    then body rdx dst src count k >> go (not swap) (k+1)
-                    else body rdx src dst count k >> go (not swap) (k+1)
-   | otherwise  = when swap (unsafeCopy src dst)
-{-# INLINE radixLoop #-}
-
-body :: (PrimMonad m, MVector v e)
-     => (Int -> e -> Int)            -- radix function
-     -> v (PrimState m) e            -- source array
-     -> v (PrimState m) e            -- destination array
-     -> PV.MVector (PrimState m) Int -- radix count
-     -> Int                          -- current pass
-     -> m ()
-body rdx src dst count k = do
-  countLoop (rdx k) src count
-  accumulate count
-  moveLoop k rdx src dst count
-{-# INLINE body #-}
-
-accumulate :: (PrimMonad m)
-           => PV.MVector (PrimState m) Int -> m ()
-accumulate count = go 0 0
- where
- len = length count
- go i acc
-   | i < len   = do ci <- unsafeRead count i
-                    unsafeWrite count i acc
-                    go (i+1) (acc + ci)
-   | otherwise = return ()
-{-# INLINE accumulate #-}
-
-moveLoop :: (PrimMonad m, MVector v e)
-         => Int -> (Int -> e -> Int) -> v (PrimState m) e
-         -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m ()
-moveLoop k rdx src dst prefix = go 0
- where
- len = length src
- go i
-   | i < len    = do srci <- unsafeRead src i
-                     pf   <- inc prefix (rdx k srci)
-                     unsafeWrite dst pf srci
-                     go (i+1)
-   | otherwise  = return ()
-{-# INLINE moveLoop #-}
-
diff --git a/Data/Vector/Algorithms/Search.hs b/Data/Vector/Algorithms/Search.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/Search.hs
+++ /dev/null
@@ -1,127 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.Search
--- Copyright   : (c) 2009-2010 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (bang patterns)
---
--- This module implements several methods of searching for indicies to insert
--- elements into a sorted vector.
-
-module Data.Vector.Algorithms.Search
-       ( binarySearch
-       , binarySearchBy
-       , binarySearchByBounds
-       , binarySearchL
-       , binarySearchLBy
-       , binarySearchLByBounds
-       , binarySearchR
-       , binarySearchRBy
-       , binarySearchRByBounds
-       , Comparison
-       ) where
-
-import Prelude hiding (read, length)
-
-import Control.Monad.Primitive
-
-import Data.Bits
-
-import Data.Vector.Generic.Mutable
-
-import Data.Vector.Algorithms.Common (Comparison)
-
--- | Finds an index in a given sorted vector at which the given element could
--- be inserted while maintaining the sortedness of the vector.
-binarySearch :: (PrimMonad m, MVector v e, Ord e)
-             => v (PrimState m) e -> e -> m Int
-binarySearch = binarySearchBy compare
-{-# INLINE binarySearch #-}
-
--- | Finds an index in a given vector, which must be sorted with respect to the
--- given comparison function, at which the given element could be inserted while
--- preserving the vector's sortedness.
-binarySearchBy :: (PrimMonad m, MVector v e)
-               => Comparison e -> v (PrimState m) e -> e -> m Int
-binarySearchBy cmp vec e = binarySearchByBounds cmp vec e 0 (length vec)
-{-# INLINE binarySearchBy #-}
-
--- | Given a vector sorted with respect to a given comparison function in indices
--- in [l,u), finds an index in [l,u] at which the given element could be inserted
--- while preserving sortedness.
-binarySearchByBounds :: (PrimMonad m, MVector v e)
-                     => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
-binarySearchByBounds cmp vec e = loop
- where
- loop !l !u
-   | u <= l    = return l
-   | otherwise = do e' <- unsafeRead vec k
-                    case cmp e' e of
-                      LT -> loop (k+1) u
-                      EQ -> return k
-                      GT -> loop l     k
-  where k = (u + l) `shiftR` 1
-{-# INLINE binarySearchByBounds #-}
-
--- | Finds the lowest index in a given sorted vector at which the given element
--- could be inserted while maintaining the sortedness.
-binarySearchL :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> e -> m Int
-binarySearchL = binarySearchLBy compare
-{-# INLINE binarySearchL #-}
-
--- | Finds the lowest index in a given vector, which must be sorted with respect to
--- the given comparison function, at which the given element could be inserted
--- while preserving the sortedness.
-binarySearchLBy :: (PrimMonad m, MVector v e)
-                => Comparison e -> v (PrimState m) e -> e -> m Int
-binarySearchLBy cmp vec e = binarySearchLByBounds cmp vec e 0 (length vec)
-{-# INLINE binarySearchLBy #-}
-
--- | Given a vector sorted with respect to a given comparison function on indices
--- in [l,u), finds the lowest index in [l,u] at which the given element could be
--- inserted while preserving sortedness.
-binarySearchLByBounds :: (PrimMonad m, MVector v e)
-                      => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
-binarySearchLByBounds cmp vec e = loop
- where
- loop !l !u
-   | u <= l    = return l
-   | otherwise = do e' <- unsafeRead vec k
-                    case cmp e' e of
-                      LT -> loop (k+1) u
-                      _  -> loop l     k
-  where k = (u + l) `shiftR` 1
-{-# INLINE binarySearchLByBounds #-}
-
--- | Finds the greatest index in a given sorted vector at which the given element
--- could be inserted while maintaining sortedness.
-binarySearchR :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> e -> m Int
-binarySearchR = binarySearchRBy compare
-{-# INLINE binarySearchR #-}
-
--- | Finds the greatest index in a given vector, which must be sorted with respect to
--- the given comparison function, at which the given element could be inserted
--- while preserving the sortedness.
-binarySearchRBy :: (PrimMonad m, MVector v e)
-                => Comparison e -> v (PrimState m) e -> e -> m Int
-binarySearchRBy cmp vec e = binarySearchRByBounds cmp vec e 0 (length vec)
-{-# INLINE binarySearchRBy #-}
-
--- | Given a vector sorted with respect to the given comparison function on indices
--- in [l,u), finds the greatest index in [l,u] at which the given element could be
--- inserted while preserving sortedness.
-binarySearchRByBounds :: (PrimMonad m, MVector v e)
-                      => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
-binarySearchRByBounds cmp vec e = loop
- where
- loop !l !u
-   | u <= l    = return l
-   | otherwise = do e' <- unsafeRead vec k
-                    case cmp e' e of
-                      GT -> loop l     k
-                      _  -> loop (k+1) u
-  where k = (u + l) `shiftR` 1
-{-# INLINE binarySearchRByBounds #-}
diff --git a/bench/Blocks.hs b/bench/Blocks.hs
--- a/bench/Blocks.hs
+++ b/bench/Blocks.hs
@@ -9,7 +9,7 @@
 
 import System.CPUTime
 
-import System.Random.Mersenne
+import System.Random.MWC (GenIO, Variate(..))
 
 -- Some conveniences for doing evil stuff in the ST monad.
 -- All the tests get run in IO, but uvector stuff happens
@@ -18,8 +18,8 @@
 clock = getCPUTime
 
 -- Strategies for filling the initial arrays
-rand :: (MTRandom e) => MTGen -> Int -> IO e
-rand g _ = random g
+rand :: Variate e => GenIO -> Int -> IO e
+rand g _ = uniform g
 
 ascend :: Num e => Int -> IO e
 ascend = return . fromIntegral
diff --git a/bench/LICENSE b/bench/LICENSE
deleted file mode 100644
--- a/bench/LICENSE
+++ /dev/null
@@ -1,30 +0,0 @@
-Copyright (c) 2009 Dan Doel
-
-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/bench/Main.hs b/bench/Main.hs
--- a/bench/Main.hs
+++ b/bench/Main.hs
@@ -23,7 +23,7 @@
 
 import System.Environment
 import System.Console.GetOpt
-import System.Random.Mersenne
+import System.Random.MWC
 
 import Blocks
 
@@ -44,7 +44,7 @@
 run :: String -> IO Integer -> IO ()
 run s t = t >>= displayTime s
 
-sortSuite :: String -> MTGen -> Int -> (MVector RealWorld Int -> IO ()) -> IO ()
+sortSuite :: String -> GenIO -> Int -> (MVector RealWorld Int -> IO ()) -> IO ()
 sortSuite str g n sort = do
   putStrLn $ "Testing: " ++ str
   run "Random            " $ speedTest n (rand g >=> modulo n) sort
@@ -54,7 +54,7 @@
   let m = 4 * (n `div` 4)
   run "Median killer     " $ speedTest m (medianKiller m) sort
 
-partialSortSuite :: String -> MTGen -> Int -> Int
+partialSortSuite :: String -> GenIO -> Int -> Int
                  -> (MVector RealWorld Int -> Int -> IO ()) -> IO ()
 partialSortSuite str g n k sort = sortSuite str g n (\a -> sort a k)
 
@@ -125,7 +125,7 @@
   [(x,t)] | all isSpace t -> Right x
   _                       -> Left s
 
-runTest :: MTGen -> Int -> Int -> Algorithm -> IO ()
+runTest :: GenIO -> Int -> Int -> Algorithm -> IO ()
 runTest g n k alg = case alg of
   DoNothing          -> sortSuite        "no algorithm"          g n   noalgo
   Allocate           -> sortSuite        "allocate"              g n   alloc
@@ -182,14 +182,13 @@
 {-# NOINLINE flagSort #-}
 
 main :: IO ()
-main = do args <- getArgs
-          gen  <- getStdGen
-          case getOpt Permute options args of
-            (fs, _, []) -> case foldl (>>=) (Right defaultOptions) fs of
-              Left err   -> putStrLn $ usageInfo err options
-              Right opts | not (usage opts) ->
-                mapM_ (runTest gen (elems opts) (portion opts)) (algos opts)
-                         | otherwise -> putStrLn $ usageInfo "uvector-algorithms-bench" options
-            (_, _, errs) -> putStrLn $ usageInfo (concat errs) options
+main = getArgs >>= \args -> withSystemRandom $ \gen ->
+  case getOpt Permute options args of
+    (fs, _, []) -> case foldl (>>=) (Right defaultOptions) fs of
+      Left err   -> putStrLn $ usageInfo err options
+      Right opts | not (usage opts) ->
+        mapM_ (runTest gen (elems opts) (portion opts)) (algos opts)
+                 | otherwise -> putStrLn $ usageInfo "uvector-algorithms-bench" options
+    (_, _, errs) -> putStrLn $ usageInfo (concat errs) options
 
 
diff --git a/bench/RadSieve.hs b/bench/RadSieve.hs
deleted file mode 100644
--- a/bench/RadSieve.hs
+++ /dev/null
@@ -1,97 +0,0 @@
--- ------------------------------------------------------------------
---
--- Module        : RadSieve
--- Copyright     : (c) 2009 Dan Doel
---
--- ------------------------------------------------------------------
--- An implementation of a radical sieve, inspired by solving Project
--- Euler problem #124.
---
--- Reproduction fo the problem text:
---
--- The radical of n, rad(n), is the product of distinct prime factors
--- of n. For example, 504 = 23 × 32 × 7, so rad(504) = 2 × 3 × 7 = 42.
---
--- If we calculate rad(n) for 1 ≤ n ≤ 10, then sort them on rad(n),
--- and sorting on n if the radical values are equal, we get:
---
---   Unsorted                 Sorted
---   n  rad(n)             n  rad(n)  k
---   1    1                1    1     1
---   2    2                2    2     2
---   3    3                4    2     3
---   4    2                8    2     4
---   5    5                3    3     5
---   6    6                9    3     6
---   7    7                5    5     7
---   8    2                6    6     8
---   9    3                7    7     9
---  10   10               10   10    10
---
--- Let E(k) be the kth element in the sorted n column; for example,
--- E(4) = 8 and E(6) = 9.
---
--- If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000).
-
-module RadSieve where
-
-import Control.Monad
-import Control.Monad.ST
-
-import Data.Array.Vector
-
--- Radicals can be sieved as follows:
---   set a[1,n] = 1
---   for i from 2 to n
---     if a[i] == 1     -- i must be prime
---      then a[j*i] *= i for positive integers j, j*i <= n
---      else do nothing -- i is composite, so its prime factors
---                      -- have been accounted for
---
--- This sieves for radicals up to the given integer.
-radSieve :: Int -> ST s (MUArr Int s)
-radSieve n = do arr <- newMU (n + 1)
-                fill arr n
-                sieve arr 1
-                return arr
- where
- fill arr i   | i < 0     = return ()
-              | otherwise = writeMU arr i 1 >> fill arr (i-1)
- sieve arr i  | n < i     = return ()
-              | otherwise = do e <- readMU arr i
-                               when (e == 1) $ mark arr i i
-                               sieve arr (i+1)
- mark arr p j | n < j     = return ()
-              | otherwise =  readMU arr j >>= writeMU arr j . (*p)
-                          >> mark arr p (j+p)
-
--- Computes the answer to the above Project Euler problem. The correct
--- answer is only generated for a stable sorting function.
-stableSortedRad :: Int -> Int
-                -> (forall s e. UA e => Comparison e -> MUArr e s -> ST s ()) 
-                -> Int
-stableSortedRad n k sortBy = runST (do rads <- radSieve n
-                                       index <- newMU (n + 1)
-                                       fillUp index n
-                                       sortBy (comparing fstS)
-                                              (unsafeZipMU rads index)
-                                       readMU k index)
- where
- fillUp arr k | k < 0     = return ()
-              | otherwise = writeMU arr k k >> fillUp arr (k-1)
-
--- Computes the answer to the above Project Euler problem. This version
--- will generate the correct answer even for unstable sorts, but may be
--- marginally slower.
-unstableSortedRad :: Int -> Int
-                  -> (forall s e. UA e => Comparison e -> MUArr e s -> ST s ()) 
-                  -> Int
-unstableSortedRad n k sortBy = runST (do rads <- radSieve n
-                                       index <- newMU (n + 1)
-                                       fillUp index n
-                                       sortBy compare (unsafeZipMU rads index)
-                                       readMU k index)
- where
- fillUp arr k | k < 0     = return ()
-              | otherwise = writeMU arr k k >> fillUp arr (k-1)
-
diff --git a/bench/vector-algorithms-bench.cabal b/bench/vector-algorithms-bench.cabal
deleted file mode 100644
--- a/bench/vector-algorithms-bench.cabal
+++ /dev/null
@@ -1,22 +0,0 @@
-name:                   vector-algorithms-bench
-version:                0.3
-license:                BSD3
-license-file:           LICENSE
-author:                 Dan Doel
-maintainer:             Dan Doel <dan.doel@gmail.com>
-homepage:               http://code.haskell.org/~doio/
-category:               Benchmark
-synopsis:               Benchmarks for vector-algorithms
-description:            A suite of various benchmarks for verifying the
-                        performance of the algorithms in vector-algorithms.
-build-type:             Simple
-cabal-version:          >= 1.2
-
-executable vec-bench
-  build-depends:        base, mersenne-random, vector, vector-algorithms, mtl
-
-  ghc-options:          -Wall -Odph
-  main-is:              Main.hs
-
-  extensions:
-      Rank2Types
diff --git a/src/Data/Vector/Algorithms/AmericanFlag.hs b/src/Data/Vector/Algorithms/AmericanFlag.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/AmericanFlag.hs
@@ -0,0 +1,339 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# lANGUAGE ScopedTypeVariables #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.AmericanFlag
+-- Copyright   : (c) 2011 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (FlexibleContexts, ScopedTypeVariables)
+--
+-- This module implements American flag sort: an in-place, unstable, bucket
+-- sort. Also in contrast to radix sort, the values are inspected in a big
+-- endian order, and buckets are sorted via recursive splitting. This,
+-- however, makes it sensible for sorting strings in lexicographic order
+-- (provided indexing is fast).
+--
+-- The algorithm works as follows: at each stage, the array is looped over,
+-- counting the number of elements for each bucket. Then, starting at the
+-- beginning of the array, elements are permuted in place to reside in the
+-- proper bucket, following chains until they reach back to the current
+-- base index. Finally, each bucket is sorted recursively. This lends itself
+-- well to the aforementioned variable-length strings, and so the algorithm
+-- takes a stopping predicate, which is given a representative of the stripe,
+-- rather than running for a set number of iterations.
+
+module Data.Vector.Algorithms.AmericanFlag ( sort
+                                           , sortBy
+                                           , Lexicographic(..)
+                                           ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import Data.Word
+import Data.Int
+import Data.Bits
+
+import qualified Data.ByteString as B
+
+import Data.Vector.Generic.Mutable
+import qualified Data.Vector.Primitive.Mutable as PV
+
+import qualified Data.Vector.Unboxed.Mutable as U
+
+import Data.Vector.Algorithms.Common
+
+import qualified Data.Vector.Algorithms.Insertion as I
+
+-- | The methods of this class specify the information necessary to sort
+-- arrays using the default ordering. The name 'Lexicographic' is meant
+-- to convey that index should return results in a similar way to indexing
+-- into a string.
+class Lexicographic e where
+  -- | Given a representative of a stripe and an index number, this
+  -- function should determine whether to stop sorting.
+  terminate :: e -> Int -> Bool
+  -- | The size of the bucket array necessary for sorting es
+  size      :: e -> Int
+  -- | Determines which bucket a given element should inhabit for a
+  -- particular iteration.
+  index     :: Int -> e -> Int
+
+instance Lexicographic Word8 where
+  terminate _ n = n > 0
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index _ n = fromIntegral n
+  {-# INLINE index #-}
+
+instance Lexicographic Word16 where
+  terminate _ n = n > 1
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ (n `shiftR`  8) .&. 255
+  index 1 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Word32 where
+  terminate _ n = n > 3
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ (n `shiftR` 24) .&. 255
+  index 1 n = fromIntegral $ (n `shiftR` 16) .&. 255
+  index 2 n = fromIntegral $ (n `shiftR`  8) .&. 255
+  index 3 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Word64 where
+  terminate _ n = n > 7
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ (n `shiftR` 56) .&. 255
+  index 1 n = fromIntegral $ (n `shiftR` 48) .&. 255
+  index 2 n = fromIntegral $ (n `shiftR` 40) .&. 255
+  index 3 n = fromIntegral $ (n `shiftR` 32) .&. 255
+  index 4 n = fromIntegral $ (n `shiftR` 24) .&. 255
+  index 5 n = fromIntegral $ (n `shiftR` 16) .&. 255
+  index 6 n = fromIntegral $ (n `shiftR`  8) .&. 255
+  index 7 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Word where
+  terminate _ n = n > 7
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ (n `shiftR` 56) .&. 255
+  index 1 n = fromIntegral $ (n `shiftR` 48) .&. 255
+  index 2 n = fromIntegral $ (n `shiftR` 40) .&. 255
+  index 3 n = fromIntegral $ (n `shiftR` 32) .&. 255
+  index 4 n = fromIntegral $ (n `shiftR` 24) .&. 255
+  index 5 n = fromIntegral $ (n `shiftR` 16) .&. 255
+  index 6 n = fromIntegral $ (n `shiftR`  8) .&. 255
+  index 7 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Int8 where
+  terminate _ n = n > 0
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index _ n = 255 .&. fromIntegral n `xor` 128
+  {-# INLINE index #-}
+
+instance Lexicographic Int16 where
+  terminate _ n = n > 1
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ ((n `xor` minBound) `shiftR` 8) .&. 255
+  index 1 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Int32 where
+  terminate _ n = n > 3
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ ((n `xor` minBound) `shiftR` 24) .&. 255
+  index 1 n = fromIntegral $ (n `shiftR` 16) .&. 255
+  index 2 n = fromIntegral $ (n `shiftR`  8) .&. 255
+  index 3 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Int64 where
+  terminate _ n = n > 7
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = fromIntegral $ ((n `xor` minBound) `shiftR` 56) .&. 255
+  index 1 n = fromIntegral $ (n `shiftR` 48) .&. 255
+  index 2 n = fromIntegral $ (n `shiftR` 40) .&. 255
+  index 3 n = fromIntegral $ (n `shiftR` 32) .&. 255
+  index 4 n = fromIntegral $ (n `shiftR` 24) .&. 255
+  index 5 n = fromIntegral $ (n `shiftR` 16) .&. 255
+  index 6 n = fromIntegral $ (n `shiftR`  8) .&. 255
+  index 7 n = fromIntegral $ n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic Int where
+  terminate _ n = n > 7
+  {-# INLINE terminate #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index 0 n = ((n `xor` minBound) `shiftR` 56) .&. 255
+  index 1 n = (n `shiftR` 48) .&. 255
+  index 2 n = (n `shiftR` 40) .&. 255
+  index 3 n = (n `shiftR` 32) .&. 255
+  index 4 n = (n `shiftR` 24) .&. 255
+  index 5 n = (n `shiftR` 16) .&. 255
+  index 6 n = (n `shiftR`  8) .&. 255
+  index 7 n = n .&. 255
+  index _ _ = 0
+  {-# INLINE index #-}
+
+instance Lexicographic B.ByteString where
+  terminate b i = i >= B.length b
+  {-# INLINE terminate #-}
+  size _ = 257
+  {-# INLINE size #-}
+  index i b
+    | i >= B.length b = 0
+    | otherwise       = fromIntegral (B.index b i) + 1
+  {-# INLINE index #-}
+
+-- | Sorts an array using the default ordering. Both Lexicographic and
+-- Ord are necessary because the algorithm falls back to insertion sort
+-- for sufficiently small arrays.
+sort :: forall e m v. (PrimMonad m, MVector v e, Lexicographic e, Ord e)
+     => v (PrimState m) e -> m ()
+sort v = sortBy compare terminate (size e) index v
+ where e :: e
+       e = undefined
+{-# INLINABLE sort #-}
+
+-- | A fully parameterized version of the sorting algorithm. Again, this
+-- function takes both radix information and a comparison, because the
+-- algorithms falls back to insertion sort for small arrays.
+sortBy :: (PrimMonad m, MVector v e)
+       => Comparison e       -- ^ a comparison for the insertion sort flalback
+       -> (e -> Int -> Bool) -- ^ determines whether a stripe is complete
+       -> Int                -- ^ the number of buckets necessary
+       -> (Int -> e -> Int)  -- ^ the big-endian radix function
+       -> v (PrimState m) e  -- ^ the array to be sorted
+       -> m ()
+sortBy cmp stop buckets radix v
+  | length v == 0 = return ()
+  | otherwise     = do count <- new buckets
+                       pile <- new buckets
+                       countLoop (radix 0) v count
+                       flagLoop cmp stop radix count pile v
+{-# INLINE sortBy #-}
+
+flagLoop :: (PrimMonad m, MVector v e)
+         => Comparison e
+         -> (e -> Int -> Bool)           -- number of passes
+         -> (Int -> e -> Int)            -- radix function
+         -> PV.MVector (PrimState m) Int -- auxiliary count array
+         -> PV.MVector (PrimState m) Int -- auxiliary pile array
+         -> v (PrimState m) e            -- source array
+         -> m ()
+flagLoop cmp stop radix count pile v = go 0 v
+ where
+
+ go pass v = do e <- unsafeRead v 0
+                unless (stop e $ pass - 1) $ go' pass v
+
+ go' pass v
+   | len < threshold = I.sortByBounds cmp v 0 len
+   | otherwise       = do accumulate count pile
+                          permute (radix pass) count pile v
+                          recurse 0
+  where
+  len = length v
+  ppass = pass + 1
+
+  recurse i
+    | i < len   = do j <- countStripe (radix ppass) (radix pass) count v i
+                     go ppass (unsafeSlice i (j - i) v)
+                     recurse j
+    | otherwise = return ()
+{-# INLINE flagLoop #-}
+
+accumulate :: (PrimMonad m)
+           => PV.MVector (PrimState m) Int
+           -> PV.MVector (PrimState m) Int
+           -> m ()
+accumulate count pile = loop 0 0
+ where
+ len = length count
+
+ loop i acc
+   | i < len = do ci <- unsafeRead count i
+                  let acc' = acc + ci
+                  unsafeWrite pile i acc
+                  unsafeWrite count i acc'
+                  loop (i+1) acc'
+   | otherwise    = return ()
+{-# INLINE accumulate #-}
+
+permute :: (PrimMonad m, MVector v e)
+        => (e -> Int)                       -- radix function
+        -> PV.MVector (PrimState m) Int     -- count array
+        -> PV.MVector (PrimState m) Int     -- pile array
+        -> v (PrimState m) e                -- source array
+        -> m ()
+permute rdx count pile v = go 0
+ where
+ len = length v
+
+ go i
+   | i < len   = do e <- unsafeRead v i
+                    let r = rdx e
+                    p <- unsafeRead pile r
+                    m <- if r > 0
+                            then unsafeRead count (r-1)
+                            else return 0
+                    case () of
+                      -- if the current element is alunsafeReady in the right pile,
+                      -- go to the end of the pile
+                      _ | m <= i && i < p  -> go p
+                      -- if the current element happens to be in the right
+                      -- pile, bump the pile counter and go to the next element
+                        | i == p           -> unsafeWrite pile r (p+1) >> go (i+1)
+                      -- otherwise follow the chain
+                        | otherwise        -> follow i e p >> go (i+1)
+   | otherwise = return ()
+ 
+ follow i e j = do en <- unsafeRead v j
+                   let r = rdx en
+                   p <- inc pile r
+                   if p == j
+                      -- if the target happens to be in the right pile, don't move it.
+                      then follow i e (j+1)
+                      else unsafeWrite v j e >> if i == p
+                                             then unsafeWrite v i en
+                                             else follow i en p
+{-# INLINE permute #-}
+
+countStripe :: (PrimMonad m, MVector v e)
+            => (e -> Int)                   -- radix function
+            -> (e -> Int)                   -- stripe function
+            -> PV.MVector (PrimState m) Int -- count array
+            -> v (PrimState m) e            -- source array
+            -> Int                          -- starting position
+            -> m Int                        -- end of stripe: [lo,hi)
+countStripe rdx str count v lo = do set count 0
+                                    e <- unsafeRead v lo
+                                    go (str e) e (lo+1)
+ where
+ len = length v
+
+ go !s e i = inc count (rdx e) >>
+            if i < len
+               then do en <- unsafeRead v i
+                       if str en == s
+                          then go s en (i+1)
+                          else return i
+                else return len
+{-# INLINE countStripe #-}
+
+threshold :: Int
+threshold = 25
+
diff --git a/src/Data/Vector/Algorithms/Common.hs b/src/Data/Vector/Algorithms/Common.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Common.hs
@@ -0,0 +1,47 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Common
+-- Copyright   : (c) 2008-2011 Dan Doel
+-- Maintainer  : Dan Doel
+-- Stability   : Experimental
+-- Portability : Portable
+--
+-- Common operations and utility functions for all sorts
+
+module Data.Vector.Algorithms.Common where
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Vector.Generic.Mutable
+
+import qualified Data.Vector.Primitive.Mutable as PV
+
+-- | A type of comparisons between two values of a given type.
+type Comparison e = e -> e -> Ordering
+
+copyOffset :: (PrimMonad m, MVector v e)
+           => v (PrimState m) e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+copyOffset from to iFrom iTo len =
+  unsafeCopy (unsafeSlice iTo len to) (unsafeSlice iFrom len from)
+{-# INLINE copyOffset #-}
+
+inc :: (PrimMonad m, MVector v Int) => v (PrimState m) Int -> Int -> m Int
+inc arr i = unsafeRead arr i >>= \e -> unsafeWrite arr i (e+1) >> return e
+{-# INLINE inc #-}
+
+-- shared bucket sorting stuff
+countLoop :: (PrimMonad m, MVector v e)
+          => (e -> Int)
+          -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m ()
+countLoop rdx src count = set count 0 >> go 0
+ where
+ len = length src
+ go i
+   | i < len    = unsafeRead src i >>= inc count . rdx >> go (i+1)
+   | otherwise  = return ()
+{-# INLINE countLoop #-}
+
diff --git a/src/Data/Vector/Algorithms/Heap.hs b/src/Data/Vector/Algorithms/Heap.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Heap.hs
@@ -0,0 +1,240 @@
+{-# LANGUAGE TypeOperators #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Heap
+-- Copyright   : (c) 2008-2011 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (type operators)
+--
+-- This module implements operations for working with a quaternary heap stored
+-- in an unboxed array. Most heapsorts are defined in terms of a binary heap,
+-- in which each internal node has at most two children. By contrast, a
+-- quaternary heap has internal nodes with up to four children. This reduces
+-- the number of comparisons in a heapsort slightly, and improves locality
+-- (again, slightly) by flattening out the heap.
+
+module Data.Vector.Algorithms.Heap
+       ( -- * Sorting
+         sort
+       , sortBy
+       , sortByBounds
+         -- * Selection
+       , select
+       , selectBy
+       , selectByBounds
+         -- * Partial sorts
+       , partialSort
+       , partialSortBy
+       , partialSortByBounds
+         -- * Heap operations
+       , heapify
+       , pop
+       , popTo
+       , sortHeap
+       , Comparison
+       ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import Data.Bits
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison)
+
+import qualified Data.Vector.Algorithms.Optimal as O
+
+-- | Sorts an entire array using the default ordering.
+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
+sort = sortBy compare
+{-# INLINABLE sort #-}
+
+-- | Sorts an entire array using a custom ordering.
+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
+sortBy cmp a = sortByBounds cmp a 0 (length a)
+{-# INLINE sortBy #-}
+
+-- | Sorts a portion of an array [l,u) using a custom ordering
+sortByBounds :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
+sortByBounds cmp a l u
+  | len < 2   = return ()
+  | len == 2  = O.sort2ByOffset cmp a l
+  | len == 3  = O.sort3ByOffset cmp a l
+  | len == 4  = O.sort4ByOffset cmp a l
+  | otherwise = heapify cmp a l u >> sortHeap cmp a l (l+4) u >> O.sort4ByOffset cmp a l
+ where len = u - l
+{-# INLINE sortByBounds #-}
+
+-- | Moves the lowest k elements to the front of the array.
+-- The elements will be in no particular order.
+select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
+select = selectBy compare
+{-# INLINE select #-}
+
+-- | Moves the 'lowest' (as defined by the comparison) k elements
+-- to the front of the array. The elements will be in no particular
+-- order.
+selectBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> Int -> m ()
+selectBy cmp a k = selectByBounds cmp a k 0 (length a)
+{-# INLINE selectBy #-}
+
+-- | Moves the 'lowest' k elements in the portion [l,u) of the
+-- array into the positions [l,k+l). The elements will be in
+-- no particular order.
+selectByBounds :: (PrimMonad m, MVector v e)
+               => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+selectByBounds cmp a k l u
+  | l + k <= u = heapify cmp a l (l + k) >> go l (l + k) (u - 1)
+  | otherwise  = return ()
+ where
+ go l m u
+   | u < m      = return ()
+   | otherwise  = do el <- unsafeRead a l
+                     eu <- unsafeRead a u
+                     case cmp eu el of
+                       LT -> popTo cmp a l m u
+                       _  -> return ()
+                     go l m (u - 1)
+{-# INLINE selectByBounds #-}
+
+-- | Moves the lowest k elements to the front of the array, sorted.
+partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
+partialSort = partialSortBy compare
+{-# INLINE partialSort #-}
+
+-- | Moves the lowest k elements (as defined by the comparison) to
+-- the front of the array, sorted.
+partialSortBy :: (PrimMonad m, MVector v e)
+              => Comparison e -> v (PrimState m) e -> Int -> m ()
+partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a)
+{-# INLINE partialSortBy #-}
+
+-- | Moves the lowest k elements in the portion [l,u) of the array
+-- into positions [l,k+l), sorted.
+partialSortByBounds :: (PrimMonad m, MVector v e)
+                    => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+partialSortByBounds cmp a k l u
+  -- this potentially does more work than absolutely required,
+  -- but using a heap to find the least 2 of 4 elements
+  -- seems unlikely to be better than just sorting all of them
+  -- with an optimal sort, and the latter is obviously index
+  -- correct.
+  | len <  2   = return ()
+  | len == 2   = O.sort2ByOffset cmp a l
+  | len == 3   = O.sort3ByOffset cmp a l
+  | len == 4   = O.sort4ByOffset cmp a l
+  | u <= l + k = sortByBounds cmp a l u
+  | otherwise  = do selectByBounds cmp a k l u
+                    sortHeap cmp a l (l + 4) (l + k)
+                    O.sort4ByOffset cmp a l
+ where
+ len = u - l
+{-# INLINE partialSortByBounds #-}
+
+-- | Constructs a heap in a portion of an array [l, u)
+heapify :: (PrimMonad m, MVector v e)
+        => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
+heapify cmp a l u = loop $ (len - 1) `shiftR` 2
+  where
+ len = u - l
+ loop k
+   | k < 0     = return ()
+   | otherwise = unsafeRead a (l+k) >>= \e ->
+                   siftByOffset cmp a e l k len >> loop (k - 1)
+{-# INLINE heapify #-}
+
+-- | Given a heap stored in a portion of an array [l,u), swaps the
+-- top of the heap with the element at u and rebuilds the heap.
+pop :: (PrimMonad m, MVector v e)
+    => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
+pop cmp a l u = popTo cmp a l u u
+{-# INLINE pop #-}
+
+-- | Given a heap stored in a portion of an array [l,u) swaps the top
+-- of the heap with the element at position t, and rebuilds the heap.
+popTo :: (PrimMonad m, MVector v e)
+      => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+popTo cmp a l u t = do al <- unsafeRead a l
+                       at <- unsafeRead a t
+                       unsafeWrite a t al
+                       siftByOffset cmp a at l 0 (u - l)
+{-# INLINE popTo #-}
+
+-- | Given a heap stored in a portion of an array [l,u), sorts the
+-- highest values into [m,u). The elements in [l,m) are not in any
+-- particular order.
+sortHeap :: (PrimMonad m, MVector v e)
+         => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+sortHeap cmp a l m u = loop (u-1) >> unsafeSwap a l m
+ where
+ loop k
+   | m < k     = pop cmp a l k >> loop (k-1)
+   | otherwise = return ()
+{-# INLINE sortHeap #-}
+
+-- Rebuilds a heap with a hole in it from start downwards. Afterward,
+-- the heap property should apply for [start + off, len + off). val
+-- is the new value to be put in the hole.
+siftByOffset :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> Int -> m ()
+siftByOffset cmp a val off start len = sift val start len
+ where
+ sift val root len
+   | child < len = do (child', ac) <- maximumChild cmp a off child len
+                      case cmp val ac of
+                        LT -> unsafeWrite a (root + off) ac >> sift val child' len
+                        _  -> unsafeWrite a (root + off) val
+   | otherwise = unsafeWrite a (root + off) val
+  where child = root `shiftL` 2 + 1
+{-# INLINE siftByOffset #-}
+
+-- Finds the maximum child of a heap node, given the indx of the first child.
+maximumChild :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m (Int,  e)
+maximumChild cmp a off child1 len
+  | child4 < len = do ac1 <- unsafeRead a (child1 + off)
+                      ac2 <- unsafeRead a (child2 + off)
+                      ac3 <- unsafeRead a (child3 + off)
+                      ac4 <- unsafeRead a (child4 + off)
+                      return $ case cmp ac1 ac2 of
+                                 LT -> case cmp ac2 ac3 of
+                                         LT -> case cmp ac3 ac4 of
+                                                 LT -> (child4, ac4)
+                                                 _  -> (child3, ac3)
+                                         _  -> case cmp ac2 ac4 of
+                                                 LT -> (child4, ac4)
+                                                 _  -> (child2, ac2)
+                                 _  -> case cmp ac1 ac3 of
+                                         LT -> case cmp ac3 ac4 of
+                                                 LT -> (child4, ac4)
+                                                 _  -> (child3, ac3)
+                                         _  -> case cmp ac1 ac4 of
+                                                 LT -> (child4, ac4)
+                                                 _  -> (child1, ac1)
+  | child3 < len = do ac1 <- unsafeRead a (child1 + off)
+                      ac2 <- unsafeRead a (child2 + off)
+                      ac3 <- unsafeRead a (child3 + off)
+                      return $ case cmp ac1 ac2 of
+                                 LT -> case cmp ac2 ac3 of
+                                         LT -> (child3, ac3)
+                                         _  -> (child2, ac2)
+                                 _  -> case cmp ac1 ac3 of
+                                         LT -> (child3, ac3)
+                                         _  -> (child1, ac1)
+  | child2 < len = do ac1 <- unsafeRead a (child1 + off)
+                      ac2 <- unsafeRead a (child2 + off)
+                      return $ case cmp ac1 ac2 of
+                                 LT -> (child2, ac2)
+                                 _  -> (child1, ac1)
+  | otherwise    = do ac1 <- unsafeRead a (child1 + off) ; return (child1, ac1)
+ where
+ child2 = child1 + 1
+ child3 = child1 + 2
+ child4 = child1 + 3
+{-# INLINE maximumChild #-}
diff --git a/src/Data/Vector/Algorithms/Insertion.hs b/src/Data/Vector/Algorithms/Insertion.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Insertion.hs
@@ -0,0 +1,81 @@
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Insertion
+-- Copyright   : (c) 2008-2010 Dan Doel
+-- Maintainer  : Dan Doel
+-- Stability   : Experimental
+-- Portability : Portable
+--
+-- A simple insertion sort. Though it's O(n^2), its iterative nature can be
+-- beneficial for small arrays. It is used to sort small segments of an array
+-- by some of the more heavy-duty, recursive algorithms.
+
+module Data.Vector.Algorithms.Insertion
+       ( sort
+       , sortBy
+       , sortByBounds
+       , sortByBounds'
+       , Comparison
+       ) where
+
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison)
+
+import qualified Data.Vector.Algorithms.Optimal as O
+
+-- | Sorts an entire array using the default comparison for the type
+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
+sort = sortBy compare
+{-# INLINABLE sort #-}
+
+-- | Sorts an entire array using a given comparison
+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
+sortBy cmp a = sortByBounds cmp a 0 (length a)
+{-# INLINE sortBy #-}
+
+-- | Sorts the portion of an array delimited by [l,u)
+sortByBounds :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
+sortByBounds cmp a l u
+  | len < 2   = return ()
+  | len == 2  = O.sort2ByOffset cmp a l
+  | len == 3  = O.sort3ByOffset cmp a l
+  | len == 4  = O.sort4ByOffset cmp a l
+  | otherwise = O.sort4ByOffset cmp a l >> sortByBounds' cmp a l (l + 4) u
+ where
+ len = u - l
+{-# INLINE sortByBounds #-}
+
+-- | Sorts the portion of the array delimited by [l,u) under the assumption
+-- that [l,m) is already sorted.
+sortByBounds' :: (PrimMonad m, MVector v e)
+              => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+sortByBounds' cmp a l m u = sort m
+ where
+ sort i
+   | i < u     = do v <- unsafeRead a i
+                    insert cmp a l v i
+                    sort (i+1)
+   | otherwise = return ()
+{-# INLINE sortByBounds' #-}
+
+-- Given a sorted array in [l,u), inserts val into its proper position,
+-- yielding a sorted [l,u]
+insert :: (PrimMonad m, MVector v e)
+       => Comparison e -> v (PrimState m) e -> Int -> e -> Int -> m ()
+insert cmp a l = loop
+ where
+ loop val j
+   | j <= l    = unsafeWrite a l val
+   | otherwise = do e <- unsafeRead a (j - 1)
+                    case cmp val e of
+                      LT -> unsafeWrite a j e >> loop val (j - 1)
+                      _  -> unsafeWrite a j val
+{-# INLINE insert #-}
diff --git a/src/Data/Vector/Algorithms/Intro.hs b/src/Data/Vector/Algorithms/Intro.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Intro.hs
@@ -0,0 +1,211 @@
+{-# LANGUAGE TypeOperators, BangPatterns, ScopedTypeVariables #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Intro
+-- Copyright   : (c) 2008-2011 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (type operators, bang patterns)
+--
+-- This module implements various algorithms based on the introsort algorithm,
+-- originally described by David R. Musser in the paper /Introspective Sorting
+-- and Selection Algorithms/. It is also in widespread practical use, as the
+-- standard unstable sort used in the C++ Standard Template Library.
+--
+-- Introsort is at its core a quicksort. The version implemented here has the
+-- following optimizations that make it perform better in practice:
+--
+--   * Small segments of the array are left unsorted until a final insertion
+--     sort pass. This is faster than recursing all the way down to
+--     one-element arrays.
+--
+--   * The pivot for segment [l,u) is chosen as the median of the elements at
+--     l, u-1 and (u+l)/2. This yields good behavior on mostly sorted (or
+--     reverse-sorted) arrays.
+--
+--   * The algorithm tracks its recursion depth, and if it decides it is
+--     taking too long (depth greater than 2 * lg n), it switches to a heap
+--     sort to maintain O(n lg n) worst case behavior. (This is what makes the
+--     algorithm introsort).
+
+module Data.Vector.Algorithms.Intro
+       ( -- * Sorting
+         sort
+       , sortBy
+       , sortByBounds
+         -- * Selecting
+       , select
+       , selectBy
+       , selectByBounds
+         -- * Partial sorting
+       , partialSort
+       , partialSortBy
+       , partialSortByBounds
+       , Comparison
+       ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import Data.Bits
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison)
+
+import qualified Data.Vector.Algorithms.Insertion as I
+import qualified Data.Vector.Algorithms.Optimal   as O
+import qualified Data.Vector.Algorithms.Heap      as H
+
+-- | Sorts an entire array using the default ordering.
+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
+sort = sortBy compare
+{-# INLINABLE sort #-}
+
+-- | Sorts an entire array using a custom ordering.
+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
+sortBy cmp a = sortByBounds cmp a 0 (length a)
+{-# INLINE sortBy #-}
+
+-- | Sorts a portion of an array [l,u) using a custom ordering
+sortByBounds :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
+sortByBounds cmp a l u
+  | len < 2   = return ()
+  | len == 2  = O.sort2ByOffset cmp a l
+  | len == 3  = O.sort3ByOffset cmp a l
+  | len == 4  = O.sort4ByOffset cmp a l
+  | otherwise = introsort cmp a (ilg len) l u
+ where len = u - l
+{-# INLINE sortByBounds #-}
+
+-- Internal version of the introsort loop which allows partial
+-- sort functions to call with a specified bound on iterations.
+introsort :: (PrimMonad m, MVector v e)
+          => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+introsort cmp a i l u = sort i l u >> I.sortByBounds cmp a l u
+ where
+ sort 0 l u = H.sortByBounds cmp a l u
+ sort d l u
+   | len < threshold = return ()
+   | otherwise = do O.sort3ByIndex cmp a c l (u-1) -- sort the median into the lowest position
+                    p <- unsafeRead a l
+                    mid <- partitionBy cmp a p (l+1) u
+                    unsafeSwap a l (mid - 1)
+                    sort (d-1) mid u
+                    sort (d-1) l   (mid - 1)
+  where
+  len = u - l
+  c   = (u + l) `div` 2
+{-# INLINE introsort #-}
+
+-- | Moves the least k elements to the front of the array in
+-- no particular order.
+select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
+select = selectBy compare
+{-# INLINE select #-}
+
+-- | Moves the least k elements (as defined by the comparison) to
+-- the front of the array in no particular order.
+selectBy :: (PrimMonad m, MVector v e)
+         => Comparison e -> v (PrimState m) e -> Int -> m ()
+selectBy cmp a k = selectByBounds cmp a k 0 (length a)
+{-# INLINE selectBy #-}
+
+-- | Moves the least k elements in the interval [l,u) to the positions
+-- [l,k+l) in no particular order.
+selectByBounds :: (PrimMonad m, MVector v e)
+               => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+selectByBounds cmp a k l u
+  | l >= u    = return ()
+  | otherwise = go (ilg len) l (l + k) u
+ where
+ len = u - l
+ go 0 l m u = H.selectByBounds cmp a (m - l) l u
+ go n l m u = do O.sort3ByIndex cmp a c l (u-1)
+                 p <- unsafeRead a l
+                 mid <- partitionBy cmp a p (l+1) u
+                 unsafeSwap a l (mid - 1)
+                 if m > mid
+                   then go (n-1) mid m u
+                   else if m < mid - 1
+                        then go (n-1) l m (mid - 1)
+                        else return ()
+  where c = (u + l) `div` 2
+{-# INLINE selectByBounds #-}
+
+-- | Moves the least k elements to the front of the array, sorted.
+partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()
+partialSort = partialSortBy compare
+{-# INLINE partialSort #-}
+
+-- | Moves the least k elements (as defined by the comparison) to
+-- the front of the array, sorted.
+partialSortBy :: (PrimMonad m, MVector v e)
+              => Comparison e -> v (PrimState m) e -> Int -> m ()
+partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a)
+{-# INLINE partialSortBy #-}
+
+-- | Moves the least k elements in the interval [l,u) to the positions
+-- [l,k+l), sorted.
+partialSortByBounds :: (PrimMonad m, MVector v e)
+                    => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+partialSortByBounds cmp a k l u
+  | l >= u    = return ()
+  | otherwise = go (ilg len) l (l + k) u
+ where
+ isort = introsort cmp a
+ {-# INLINE [1] isort #-}
+ len = u - l
+ go 0 l m n = H.partialSortByBounds cmp a (m - l) l u
+ go n l m u
+   | l == m    = return ()
+   | otherwise = do O.sort3ByIndex cmp a c l (u-1)
+                    p <- unsafeRead a l
+                    mid <- partitionBy cmp a p (l+1) u
+                    unsafeSwap a l (mid - 1)
+                    case compare m mid of
+                      GT -> do isort (n-1) l (mid - 1)
+                               go (n-1) mid m u
+                      EQ -> isort (n-1) l m
+                      LT -> go n l m (mid - 1)
+  where c = (u + l) `div` 2
+{-# INLINE partialSortByBounds #-}
+
+partitionBy :: forall m v e. (PrimMonad m, MVector v e)
+            => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
+partitionBy cmp a = partUp
+ where
+ -- 6.10 panics without the signatures for partUp and partDown, 6.12 and later
+ -- versions don't need them
+ partUp :: e -> Int -> Int -> m Int
+ partUp p l u
+   | l < u = do e <- unsafeRead a l
+                case cmp e p of
+                  LT -> partUp p (l+1) u
+                  _  -> partDown p l (u-1)
+   | otherwise = return l
+
+ partDown :: e -> Int -> Int -> m Int
+ partDown p l u
+   | l < u = do e <- unsafeRead a u
+                case cmp p e of
+                  LT -> partDown p l (u-1)
+                  _  -> unsafeSwap a l u >> partUp p (l+1) u
+   | otherwise = return l
+{-# INLINE partitionBy #-}
+
+-- computes the number of recursive calls after which heapsort should
+-- be invoked given the lower and upper indices of the array to be sorted
+ilg :: Int -> Int
+ilg m = 2 * loop m 0
+ where
+ loop 0 !k = k - 1
+ loop n !k = loop (n `shiftR` 1) (k+1)
+
+-- the size of array at which the introsort algorithm switches to insertion sort
+threshold :: Int
+threshold = 18
+{-# INLINE threshold #-}
diff --git a/src/Data/Vector/Algorithms/Merge.hs b/src/Data/Vector/Algorithms/Merge.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Merge.hs
@@ -0,0 +1,97 @@
+{-# LANGUAGE BangPatterns #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Merge
+-- Copyright   : (c) 2008-2011 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Portable
+--
+-- This module implements a simple top-down merge sort. The temporary buffer
+-- is preallocated to 1/2 the size of the input array, and shared through
+-- the entire sorting process to ease the amount of allocation performed in
+-- total. This is a stable sort.
+
+module Data.Vector.Algorithms.Merge
+       ( sort
+       , sortBy
+       , Comparison
+       ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Bits
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison, copyOffset)
+
+import qualified Data.Vector.Algorithms.Optimal   as O
+import qualified Data.Vector.Algorithms.Insertion as I
+
+-- | Sorts an array using the default comparison.
+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
+sort = sortBy compare
+{-# INLINABLE sort #-}
+
+-- | Sorts an array using a custom comparison.
+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
+sortBy cmp vec
+  | len <= 1  = return ()
+  | len == 2  = O.sort2ByOffset cmp vec 0
+  | len == 3  = O.sort3ByOffset cmp vec 0
+  | len == 4  = O.sort4ByOffset cmp vec 0
+  | otherwise = do buf <- new len
+                   mergeSortWithBuf cmp vec buf
+ where
+ len = length vec
+{-# INLINE sortBy #-}
+
+mergeSortWithBuf :: (PrimMonad m, MVector v e)
+                 => Comparison e -> v (PrimState m) e -> v (PrimState m) e -> m ()
+mergeSortWithBuf cmp src buf = loop 0 (length src)
+ where
+ loop l u
+   | len < threshold = I.sortByBounds cmp src l u
+   | otherwise       = do loop l mid
+                          loop mid u
+                          merge cmp (unsafeSlice l len src) buf (mid - l)
+  where len = u - l
+        mid = (u + l) `shiftR` 1
+{-# INLINE mergeSortWithBuf #-}
+
+merge :: (PrimMonad m, MVector v e)
+      => Comparison e -> v (PrimState m) e -> v (PrimState m) e
+      -> Int -> m ()
+merge cmp src buf mid = do unsafeCopy tmp lower
+                           eTmp <- unsafeRead tmp 0
+                           eUpp <- unsafeRead upper 0
+                           loop tmp 0 eTmp upper 0 eUpp 0
+ where
+ lower = unsafeSlice 0   mid                src
+ upper = unsafeSlice mid (length src - mid) src
+ tmp   = unsafeSlice 0   mid                buf
+
+ wroteHigh low iLow eLow high iHigh iIns
+   | iHigh >= length high = unsafeCopy (unsafeSlice iIns (length low - iLow) src)
+                                       (unsafeSlice iLow (length low - iLow) low)
+   | otherwise            = do eHigh <- unsafeRead high iHigh
+                               loop low iLow eLow high iHigh eHigh iIns
+
+ wroteLow low iLow high iHigh eHigh iIns
+   | iLow  >= length low  = return ()
+   | otherwise            = do eLow <- unsafeRead low iLow
+                               loop low iLow eLow high iHigh eHigh iIns
+
+ loop !low !iLow !eLow !high !iHigh !eHigh !iIns = case cmp eHigh eLow of
+     LT -> do unsafeWrite src iIns eHigh
+              wroteHigh low iLow eLow high (iHigh + 1) (iIns + 1)
+     _  -> do unsafeWrite src iIns eLow
+              wroteLow low (iLow + 1) high iHigh eHigh (iIns + 1)
+{-# INLINE merge #-}
+
+threshold :: Int
+threshold = 25
+{-# INLINE threshold #-}
diff --git a/src/Data/Vector/Algorithms/Optimal.hs b/src/Data/Vector/Algorithms/Optimal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Optimal.hs
@@ -0,0 +1,245 @@
+{-# LANGUAGE CPP #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Optimal
+-- Copyright   : (c) 2008-2010 Dan Doel
+-- Maintainer  : Dan Doel
+-- Stability   : Experimental
+-- Portability : Portable
+--
+-- Optimal sorts for very small array sizes, or for small numbers of
+-- particular indices in a larger array (to be used, for instance, for
+-- sorting a median of 3 values into the lowest position in an array
+-- for a median-of-3 quicksort).
+
+-- The code herein was adapted from a C algorithm for optimal sorts
+-- of small arrays. The original code was produced for the article
+-- /Sorting Revisited/ by Paul Hsieh, available here:
+--
+--   http://www.azillionmonkeys.com/qed/sort.html
+--
+-- The LICENSE file contains the relevant copyright information for
+-- the reference C code.
+
+module Data.Vector.Algorithms.Optimal
+       ( sort2ByIndex
+       , sort2ByOffset
+       , sort3ByIndex
+       , sort3ByOffset
+       , sort4ByIndex
+       , sort4ByOffset
+       , Comparison
+       ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison)
+
+#include "vector.h"
+
+-- | Sorts the elements at the positions 'off' and 'off + 1' in the given
+-- array using the comparison.
+sort2ByOffset :: (PrimMonad m, MVector v e)
+              => Comparison e -> v (PrimState m) e -> Int -> m ()
+sort2ByOffset cmp a off = sort2ByIndex cmp a off (off + 1)
+{-# INLINABLE sort2ByOffset #-}
+
+-- | Sorts the elements at the two given indices using the comparison. This
+-- is essentially a compare-and-swap, although the first index is assumed to
+-- be the 'lower' of the two.
+sort2ByIndex :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()
+sort2ByIndex cmp a i j = UNSAFE_CHECK(checkIndex) "sort2ByIndex" i (length a)
+                       $ UNSAFE_CHECK(checkIndex) "sort2ByIndex" j (length a) $  do
+  a0 <- unsafeRead a i
+  a1 <- unsafeRead a j
+  case cmp a0 a1 of
+    GT -> unsafeWrite a i a1 >> unsafeWrite a j a0
+    _  -> return ()
+{-# INLINABLE sort2ByIndex #-}
+
+-- | Sorts the three elements starting at the given offset in the array.
+sort3ByOffset :: (PrimMonad m, MVector v e)
+              => Comparison e -> v (PrimState m) e -> Int -> m ()
+sort3ByOffset cmp a off = sort3ByIndex cmp a off (off + 1) (off + 2)
+{-# INLINABLE sort3ByOffset #-}
+
+-- | Sorts the elements at the three given indices. The indices are assumed
+-- to be given from lowest to highest, so if 'l < m < u' then
+-- 'sort3ByIndex cmp a m l u' essentially sorts the median of three into the
+-- lowest position in the array.
+sort3ByIndex :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+sort3ByIndex cmp a i j k = UNSAFE_CHECK(checkIndex) "sort3ByIndex" i (length a)
+                         $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" j (length a)
+                         $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" k (length a) $ do
+  a0 <- unsafeRead a i
+  a1 <- unsafeRead a j
+  a2 <- unsafeRead a k
+  case cmp a0 a1 of
+    GT -> case cmp a0 a2 of
+            GT -> case cmp a2 a1 of
+                    LT -> do unsafeWrite a i a2
+                             unsafeWrite a k a0
+                    _  -> do unsafeWrite a i a1
+                             unsafeWrite a j a2
+                             unsafeWrite a k a0
+            _  -> do unsafeWrite a i a1
+                     unsafeWrite a j a0
+    _  -> case cmp a1 a2 of
+            GT -> case cmp a0 a2 of
+                    GT -> do unsafeWrite a i a2
+                             unsafeWrite a j a0
+                             unsafeWrite a k a1
+                    _  -> do unsafeWrite a j a2
+                             unsafeWrite a k a1
+            _  -> return ()
+{-# INLINABLE sort3ByIndex #-}
+
+-- | Sorts the four elements beginning at the offset.
+sort4ByOffset :: (PrimMonad m, MVector v e)
+              => Comparison e -> v (PrimState m) e -> Int -> m ()
+sort4ByOffset cmp a off = sort4ByIndex cmp a off (off + 1) (off + 2) (off + 3)
+{-# INLINABLE sort4ByOffset #-}
+
+-- The horror...
+
+-- | Sorts the elements at the four given indices. Like the 2 and 3 element
+-- versions, this assumes that the indices are given in increasing order, so
+-- it can be used to sort medians into particular positions and so on.
+sort4ByIndex :: (PrimMonad m, MVector v e)
+             => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> Int -> m ()
+sort4ByIndex cmp a i j k l = UNSAFE_CHECK(checkIndex) "sort4ByIndex" i (length a)
+                           $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" j (length a)
+                           $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" k (length a)
+                           $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" l (length a) $ do
+  a0 <- unsafeRead a i
+  a1 <- unsafeRead a j
+  a2 <- unsafeRead a k
+  a3 <- unsafeRead a l
+  case cmp a0 a1 of
+    GT -> case cmp a0 a2 of
+            GT -> case cmp a1 a2 of
+                    GT -> case cmp a1 a3 of
+                            GT -> case cmp a2 a3 of
+                                    GT -> do unsafeWrite a i a3
+                                             unsafeWrite a j a2
+                                             unsafeWrite a k a1
+                                             unsafeWrite a l a0
+                                    _  -> do unsafeWrite a i a2
+                                             unsafeWrite a j a3
+                                             unsafeWrite a k a1
+                                             unsafeWrite a l a0
+                            _  -> case cmp a0 a3 of
+                                    GT -> do unsafeWrite a i a2
+                                             unsafeWrite a j a1
+                                             unsafeWrite a k a3
+                                             unsafeWrite a l a0
+                                    _  -> do unsafeWrite a i a2
+                                             unsafeWrite a j a1
+                                             unsafeWrite a k a0
+                                             unsafeWrite a l a3
+                    _ -> case cmp a2 a3 of
+                           GT -> case cmp a1 a3 of
+                                   GT -> do unsafeWrite a i a3
+                                            unsafeWrite a j a1
+                                            unsafeWrite a k a2
+                                            unsafeWrite a l a0
+                                   _  -> do unsafeWrite a i a1
+                                            unsafeWrite a j a3
+                                            unsafeWrite a k a2
+                                            unsafeWrite a l a0
+                           _  -> case cmp a0 a3 of
+                                   GT -> do unsafeWrite a i a1
+                                            unsafeWrite a j a2
+                                            unsafeWrite a k a3
+                                            unsafeWrite a l a0
+                                   _  -> do unsafeWrite a i a1
+                                            unsafeWrite a j a2
+                                            unsafeWrite a k a0
+                                            -- unsafeWrite a l a3
+            _  -> case cmp a0 a3 of
+                    GT -> case cmp a1 a3 of
+                            GT -> do unsafeWrite a i a3
+                                     -- unsafeWrite a j a1
+                                     unsafeWrite a k a0
+                                     unsafeWrite a l a2
+                            _  -> do unsafeWrite a i a1
+                                     unsafeWrite a j a3
+                                     unsafeWrite a k a0
+                                     unsafeWrite a l a2
+                    _  -> case cmp a2 a3 of
+                            GT -> do unsafeWrite a i a1
+                                     unsafeWrite a j a0
+                                     unsafeWrite a k a3
+                                     unsafeWrite a l a2
+                            _  -> do unsafeWrite a i a1
+                                     unsafeWrite a j a0
+                                     -- unsafeWrite a k a2
+                                     -- unsafeWrite a l a3
+    _  -> case cmp a1 a2 of
+            GT -> case cmp a0 a2 of
+                    GT -> case cmp a0 a3 of
+                            GT -> case cmp a2 a3 of
+                                    GT -> do unsafeWrite a i a3
+                                             unsafeWrite a j a2
+                                             unsafeWrite a k a0
+                                             unsafeWrite a l a1
+                                    _  -> do unsafeWrite a i a2
+                                             unsafeWrite a j a3
+                                             unsafeWrite a k a0
+                                             unsafeWrite a l a1
+                            _  -> case cmp a1 a3 of
+                                    GT -> do unsafeWrite a i a2
+                                             unsafeWrite a j a0
+                                             unsafeWrite a k a3
+                                             unsafeWrite a l a1
+                                    _  -> do unsafeWrite a i a2
+                                             unsafeWrite a j a0
+                                             unsafeWrite a k a1
+                                             -- unsafeWrite a l a3
+                    _  -> case cmp a2 a3 of
+                            GT -> case cmp a0 a3 of
+                                    GT -> do unsafeWrite a i a3
+                                             unsafeWrite a j a0
+                                             -- unsafeWrite a k a2
+                                             unsafeWrite a l a1
+                                    _  -> do -- unsafeWrite a i a0
+                                             unsafeWrite a j a3
+                                             -- unsafeWrite a k a2
+                                             unsafeWrite a l a1
+                            _  -> case cmp a1 a3 of
+                                    GT -> do -- unsafeWrite a i a0
+                                             unsafeWrite a j a2
+                                             unsafeWrite a k a3
+                                             unsafeWrite a l a1
+                                    _  -> do -- unsafeWrite a i a0
+                                             unsafeWrite a j a2
+                                             unsafeWrite a k a1
+                                             -- unsafeWrite a l a3
+            _  -> case cmp a1 a3 of
+                    GT -> case cmp a0 a3 of
+                            GT -> do unsafeWrite a i a3
+                                     unsafeWrite a j a0
+                                     unsafeWrite a k a1
+                                     unsafeWrite a l a2
+                            _  -> do -- unsafeWrite a i a0
+                                     unsafeWrite a j a3
+                                     unsafeWrite a k a1
+                                     unsafeWrite a l a2
+                    _  -> case cmp a2 a3 of
+                            GT -> do -- unsafeWrite a i a0
+                                     -- unsafeWrite a j a1
+                                     unsafeWrite a k a3
+                                     unsafeWrite a l a2
+                            _  -> do -- unsafeWrite a i a0
+                                     -- unsafeWrite a j a1
+                                     -- unsafeWrite a k a2
+                                     -- unsafeWrite a l a3
+                                     return ()
+{-# INLINABLE sort4ByIndex #-}
diff --git a/src/Data/Vector/Algorithms/Radix.hs b/src/Data/Vector/Algorithms/Radix.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Radix.hs
@@ -0,0 +1,261 @@
+{-# LANGUAGE ScopedTypeVariables, BangPatterns, TypeOperators #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Radix
+-- Copyright   : (c) 2008-2011 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (scoped type variables, bang patterns)
+--
+-- This module provides a radix sort for a subclass of unboxed arrays. The
+-- radix class gives information on
+--   * the number of passes needed for the data type
+--
+--   * the size of the auxiliary arrays
+--
+--   * how to compute the pass-k radix of a value
+--
+-- Radix sort is not a comparison sort, so it is able to achieve O(n) run
+-- time, though it also uses O(n) auxiliary space. In addition, there is a
+-- constant space overhead of 2*size*sizeOf(Int) for the sort, so it is not
+-- advisable to use this sort for large numbers of very small arrays.
+--
+-- A standard example (upon which one could base their own Radix instance)
+-- is Word32:
+--
+--   * We choose to sort on r = 8 bits at a time
+--
+--   * A Word32 has b = 32 bits total
+--
+--   Thus, b/r = 4 passes are required, 2^r = 256 elements are needed in an
+--   auxiliary array, and the radix function is:
+--
+--    > radix k e = (e `shiftR` (k*8)) .&. 255
+
+module Data.Vector.Algorithms.Radix (sort, sortBy, Radix(..)) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import qualified Data.Vector.Primitive.Mutable as PV
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common
+
+import Data.Bits
+import Data.Int
+import Data.Word
+
+
+import Foreign.Storable
+
+class Radix e where
+  -- | The number of passes necessary to sort an array of es
+  passes :: e -> Int
+  -- | The size of an auxiliary array
+  size   :: e -> Int
+  -- | The radix function parameterized by the current pass
+  radix  :: Int -> e -> Int
+
+instance Radix Int where
+  passes _ = sizeOf (undefined :: Int)
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix 0 e = e .&. 255
+  radix i e
+    | i == passes e - 1 = radix' (e `xor` minBound)
+    | otherwise         = radix' e
+   where radix' e = (e `shiftR` (i `shiftL` 3)) .&. 255
+  {-# INLINE radix #-}
+
+instance Radix Int8 where
+  passes _ = 1
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix _ e = 255 .&. fromIntegral e `xor` 128
+  {-# INLINE radix #-}
+
+instance Radix Int16 where
+  passes _ = 2
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix 1 e = fromIntegral (((e `xor` minBound) `shiftR` 8) .&. 255)
+  {-# INLINE radix #-}
+
+instance Radix Int32 where
+  passes _ = 4
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
+  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
+  radix 3 e = fromIntegral (((e `xor` minBound) `shiftR` 24) .&. 255)
+  {-# INLINE radix #-}
+
+instance Radix Int64 where
+  passes _ = 8
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
+  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
+  radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)
+  radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255)
+  radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255)
+  radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255)
+  radix 7 e = fromIntegral (((e `xor` minBound) `shiftR` 56) .&. 255)
+  {-# INLINE radix #-}
+
+instance Radix Word where
+  passes _ = sizeOf (undefined :: Word)
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix i e = fromIntegral ((e `shiftR` (i `shiftL` 3)) .&. 255)
+  {-# INLINE radix #-}
+
+instance Radix Word8 where
+  passes _ = 1
+  {-# INLINE passes #-}
+  size _ = 256
+  {-# INLINE size #-}
+  radix _ = fromIntegral
+  {-# INLINE radix #-}
+
+instance Radix Word16 where
+  passes _ = 2
+  {-# INLINE passes #-}
+  size   _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
+  {-# INLINE radix #-}
+
+instance Radix Word32 where
+  passes _ = 4
+  {-# INLINE passes #-}
+  size   _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
+  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
+  radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)
+  {-# INLINE radix #-}
+
+instance Radix Word64 where
+  passes _ = 8
+  {-# INLINE passes #-}
+  size   _ = 256
+  {-# INLINE size #-}
+  radix 0 e = fromIntegral (e .&. 255)
+  radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)
+  radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)
+  radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)
+  radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255)
+  radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255)
+  radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255)
+  radix 7 e = fromIntegral ((e `shiftR` 56) .&. 255)
+  {-# INLINE radix #-}
+
+instance (Radix i, Radix j) => Radix (i, j) where
+  passes ~(i, j) = passes i + passes j
+  {-# INLINE passes #-}
+  size   ~(i, j) = size i `max` size j
+  {-# INLINE size #-}
+  radix k ~(i, j) | k < passes j = radix k j
+                     | otherwise    = radix (k - passes j) i
+  {-# INLINE radix #-}
+
+-- | Sorts an array based on the Radix instance.
+sort :: forall e m v. (PrimMonad m, MVector v e, Radix e)
+     => v (PrimState m) e -> m ()
+sort arr = sortBy (passes e) (size e) radix arr
+ where
+ e :: e
+ e = undefined
+{-# INLINABLE sort #-}
+
+-- | Radix sorts an array using custom radix information
+-- requires the number of passes to fully sort the array,
+-- the size of of auxiliary arrays necessary (should be
+-- one greater than the maximum value returned by the radix
+-- function), and a radix function, which takes the pass
+-- and an element, and returns the relevant radix.
+sortBy :: (PrimMonad m, MVector v e)
+       => Int               -- ^ the number of passes
+       -> Int               -- ^ the size of auxiliary arrays
+       -> (Int -> e -> Int) -- ^ the radix function
+       -> v (PrimState m) e -- ^ the array to be sorted
+       -> m ()
+sortBy passes size rdx arr = do
+  tmp    <- new (length arr)
+  count  <- new size
+  radixLoop passes rdx arr tmp count
+{-# INLINE sortBy #-}
+
+radixLoop :: (PrimMonad m, MVector v e)
+          => Int                          -- passes
+          -> (Int -> e -> Int)            -- radix function
+          -> v (PrimState m) e            -- array to sort
+          -> v (PrimState m) e            -- temporary array
+          -> PV.MVector (PrimState m) Int -- radix count array
+          -> m ()
+radixLoop passes rdx src dst count = go False 0
+ where
+ len = length src
+ go swap k
+   | k < passes = if swap
+                    then body rdx dst src count k >> go (not swap) (k+1)
+                    else body rdx src dst count k >> go (not swap) (k+1)
+   | otherwise  = when swap (unsafeCopy src dst)
+{-# INLINE radixLoop #-}
+
+body :: (PrimMonad m, MVector v e)
+     => (Int -> e -> Int)            -- radix function
+     -> v (PrimState m) e            -- source array
+     -> v (PrimState m) e            -- destination array
+     -> PV.MVector (PrimState m) Int -- radix count
+     -> Int                          -- current pass
+     -> m ()
+body rdx src dst count k = do
+  countLoop (rdx k) src count
+  accumulate count
+  moveLoop k rdx src dst count
+{-# INLINE body #-}
+
+accumulate :: (PrimMonad m)
+           => PV.MVector (PrimState m) Int -> m ()
+accumulate count = go 0 0
+ where
+ len = length count
+ go i acc
+   | i < len   = do ci <- unsafeRead count i
+                    unsafeWrite count i acc
+                    go (i+1) (acc + ci)
+   | otherwise = return ()
+{-# INLINE accumulate #-}
+
+moveLoop :: (PrimMonad m, MVector v e)
+         => Int -> (Int -> e -> Int) -> v (PrimState m) e
+         -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m ()
+moveLoop k rdx src dst prefix = go 0
+ where
+ len = length src
+ go i
+   | i < len    = do srci <- unsafeRead src i
+                     pf   <- inc prefix (rdx k srci)
+                     unsafeWrite dst pf srci
+                     go (i+1)
+   | otherwise  = return ()
+{-# INLINE moveLoop #-}
+
diff --git a/src/Data/Vector/Algorithms/Search.hs b/src/Data/Vector/Algorithms/Search.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Search.hs
@@ -0,0 +1,135 @@
+{-# LANGUAGE BangPatterns #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Search
+-- Copyright   : (c) 2009-2010 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (bang patterns)
+--
+-- This module implements several methods of searching for indicies to insert
+-- elements into a sorted vector.
+
+module Data.Vector.Algorithms.Search
+       ( binarySearch
+       , binarySearchBy
+       , binarySearchByBounds
+       , binarySearchL
+       , binarySearchLBy
+       , binarySearchLByBounds
+       , binarySearchR
+       , binarySearchRBy
+       , binarySearchRByBounds
+       , binarySearchP
+       , binarySearchPBounds
+       , Comparison
+       ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Bits
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison)
+
+-- | Finds an index in a given sorted vector at which the given element could
+-- be inserted while maintaining the sortedness of the vector.
+binarySearch :: (PrimMonad m, MVector v e, Ord e)
+             => v (PrimState m) e -> e -> m Int
+binarySearch = binarySearchBy compare
+{-# INLINE binarySearch #-}
+
+-- | Finds an index in a given vector, which must be sorted with respect to the
+-- given comparison function, at which the given element could be inserted while
+-- preserving the vector's sortedness.
+binarySearchBy :: (PrimMonad m, MVector v e)
+               => Comparison e -> v (PrimState m) e -> e -> m Int
+binarySearchBy cmp vec e = binarySearchByBounds cmp vec e 0 (length vec)
+{-# INLINE binarySearchBy #-}
+
+-- | Given a vector sorted with respect to a given comparison function in indices
+-- in [l,u), finds an index in [l,u] at which the given element could be inserted
+-- while preserving sortedness.
+binarySearchByBounds :: (PrimMonad m, MVector v e)
+                     => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
+binarySearchByBounds cmp vec e = loop
+ where
+ loop !l !u
+   | u <= l    = return l
+   | otherwise = do e' <- unsafeRead vec k
+                    case cmp e' e of
+                      LT -> loop (k+1) u
+                      EQ -> return k
+                      GT -> loop l     k
+  where k = (u + l) `shiftR` 1
+{-# INLINE binarySearchByBounds #-}
+
+-- | Finds the lowest index in a given sorted vector at which the given element
+-- could be inserted while maintaining the sortedness.
+binarySearchL :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> e -> m Int
+binarySearchL = binarySearchLBy compare
+{-# INLINE binarySearchL #-}
+
+-- | Finds the lowest index in a given vector, which must be sorted with respect to
+-- the given comparison function, at which the given element could be inserted
+-- while preserving the sortedness.
+binarySearchLBy :: (PrimMonad m, MVector v e)
+                => Comparison e -> v (PrimState m) e -> e -> m Int
+binarySearchLBy cmp vec e = binarySearchLByBounds cmp vec e 0 (length vec)
+{-# INLINE binarySearchLBy #-}
+
+-- | Given a vector sorted with respect to a given comparison function on indices
+-- in [l,u), finds the lowest index in [l,u] at which the given element could be
+-- inserted while preserving sortedness.
+binarySearchLByBounds :: (PrimMonad m, MVector v e)
+                      => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
+binarySearchLByBounds cmp vec e = binarySearchPBounds p vec
+ where p e' = case cmp e' e of LT -> False ; _ -> True
+{-# INLINE binarySearchLByBounds #-}
+
+-- | Finds the greatest index in a given sorted vector at which the given element
+-- could be inserted while maintaining sortedness.
+binarySearchR :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> e -> m Int
+binarySearchR = binarySearchRBy compare
+{-# INLINE binarySearchR #-}
+
+-- | Finds the greatest index in a given vector, which must be sorted with respect to
+-- the given comparison function, at which the given element could be inserted
+-- while preserving the sortedness.
+binarySearchRBy :: (PrimMonad m, MVector v e)
+                => Comparison e -> v (PrimState m) e -> e -> m Int
+binarySearchRBy cmp vec e = binarySearchRByBounds cmp vec e 0 (length vec)
+{-# INLINE binarySearchRBy #-}
+
+-- | Given a vector sorted with respect to the given comparison function on indices
+-- in [l,u), finds the greatest index in [l,u] at which the given element could be
+-- inserted while preserving sortedness.
+binarySearchRByBounds :: (PrimMonad m, MVector v e)
+                      => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int
+binarySearchRByBounds cmp vec e = binarySearchPBounds p vec
+ where p e' = case cmp e' e of GT -> True ; _ -> False
+{-# INLINE binarySearchRByBounds #-}
+
+-- | Given a predicate that is guaraneteed to be monotone on the given vector,
+-- finds the first index at which the predicate returns True, or the length of
+-- the array if the predicate is false for the entire array.
+binarySearchP :: (PrimMonad m, MVector v e) => (e -> Bool) -> v (PrimState m) e -> m Int
+binarySearchP p vec = binarySearchPBounds p vec 0 (length vec)
+{-# INLINE binarySearchP #-}
+
+-- | Given a predicate that is guaranteed to be monotone on the indices [l,u) in
+-- a given vector, finds the index in [l,u] at which the predicate turns from
+-- False to True (yielding u if the entire interval is False).
+binarySearchPBounds :: (PrimMonad m, MVector v e)
+                    => (e -> Bool) -> v (PrimState m) e -> Int -> Int -> m Int
+binarySearchPBounds p vec = loop
+ where
+ loop !l !u
+   | u <= l    = return l
+   | otherwise = unsafeRead vec k >>= \e -> if p e then loop l k else loop (k+1) u
+  where k = (u + l) `shiftR` 1
+{-# INLINE binarySearchPBounds #-}
diff --git a/tests/Optimal.hs b/tests/Optimal.hs
deleted file mode 100644
--- a/tests/Optimal.hs
+++ /dev/null
@@ -1,62 +0,0 @@
-{-# LANGUAGE TypeOperators, FlexibleContexts #-}
-
--- Exhaustive test sets for proper sorting and stability of
--- optimal sorts
-
-module Optimal where
-
-import Control.Arrow
-import Control.Monad
-
-import Data.List
-import Data.Function
-
-import Data.Vector.Generic hiding (map, zip, concatMap, (++), replicate, foldM)
-
-interleavings :: [a] -> [a] -> [[a]]
-interleavings [       ] ys        =  [ys]
-interleavings xs        [       ] =  [xs]
-interleavings xs@(x:xt) ys@(y:yt) =  map (x:) (interleavings xt ys)
-                                  ++ map (y:) (interleavings xs yt)
-
-monotones :: Int -> Int -> [[Int]]
-monotones k = atLeastOne 0
- where
- atLeastOne i 0 = [[]]
- atLeastOne i n = map (i:) $ picks i (n-1)
- picks _ 0             = [[]]
- picks i n | i >= k    = [replicate n k]
-           | otherwise = map (i:) (picks i (n-1)) ++ atLeastOne (i+1) n
-
-
-stability :: (Vector v (Int,Int)) => Int -> [v (Int, Int)]
-stability n = concatMap ( map fromList
-                        . foldM interleavings []
-                        . groupBy ((==) `on` fst)
-                        . flip zip [0..])
-              $ monotones (n-2) n
-
-sort2 :: (Vector v Int) => [v Int]
-sort2 = map fromList $ permutations [0,1]
-
-stability2 :: (Vector v (Int,Int)) => [v (Int, Int)]
-stability2 = [fromList [(0, 0), (0, 1)]]
-
-sort3 :: (Vector v Int) => [v Int]
-sort3 = map fromList $ permutations [0..2]
-
-{-
-stability3 :: [UArr (Int :*: Int)]
-stability3 = map toU [ [0:*:0, 0:*:1, 0:*:2]
-                     , [0:*:0, 0:*:1, 1:*:2]
-                     , [0:*:0, 1:*:2, 0:*:1]
-                     , [1:*:2, 0:*:0, 0:*:1]
-                     , [0:*:0, 1:*:1, 1:*:2]
-                     , [1:*:1, 0:*:0, 1:*:2]
-                     , [1:*:1, 1:*:2, 0:*:0]
-                     ]
--}
-
-sort4 :: (Vector v Int) => [v Int]
-sort4 = map fromList $ permutations [0..3]
-
diff --git a/tests/Properties.hs b/tests/Properties.hs
deleted file mode 100644
--- a/tests/Properties.hs
+++ /dev/null
@@ -1,185 +0,0 @@
-{-# LANGUAGE RankNTypes, FlexibleContexts #-}
-
-module Properties where
-
-import Prelude
-
-import Optimal
-
-import Control.Monad
-import Control.Monad.ST
-
-import Data.List
-import Data.Ord
-
-import Data.Vector (Vector)
-import qualified Data.Vector as V
-
-import Data.Vector.Mutable (MVector)
-import qualified Data.Vector.Mutable as MV
-
-import Data.Vector.Generic (modify)
-
-import qualified Data.Vector.Generic.Mutable as G
-
-import Data.Vector.Algorithms.Optimal (Comparison)
-import Data.Vector.Algorithms.Radix (radix, passes, size)
-
-import qualified Data.Map as M
-
-import Test.QuickCheck
-
-import Util
-
-prop_sorted :: (Ord e) => Vector e -> Property
-prop_sorted arr | V.length arr < 2 = property True
-                | otherwise        = check (V.head arr) (V.tail arr)
- where
- check e arr | V.null arr = property True
-             | otherwise  = e <= V.head arr .&. check (V.head arr) (V.tail arr)
-
-prop_empty :: (Ord e) => (forall s. MV.MVector s e -> ST s ()) -> Property
-prop_empty algo = prop_sorted (modify algo $ V.fromList [])
-
-prop_fullsort :: (Ord e)
-              => (forall s mv. G.MVector mv e => mv s e -> ST s ()) -> Vector e -> Property
-prop_fullsort algo arr = prop_sorted $ modify algo arr
-
-{-
-prop_schwartzian :: (UA e, UA k, Ord k)
-                 => (e -> k)
-                 -> (forall e s. (UA e) => (e -> e -> Ordering) -> MUArr e s -> ST s ())
-                 -> UArr e -> Property
-prop_schwartzian f algo arr
-  | lengthU arr < 2 = property True
-  | otherwise       = let srt = modify (algo `usingKeys` f) arr
-                      in check (headU srt) (tailU srt)
- where
- check e arr | nullU arr = property True
-             | otherwise = f e <= f (headU arr) .&. check (headU arr) (tailU arr)
--}
-
-longGen :: (Arbitrary e) => Int -> Gen (Vector e)
-longGen k = liftM2 (\l r -> V.fromList (l ++ r)) (vectorOf k arbitrary) arbitrary
-
-sanity :: Int
-sanity = 100
-
-prop_partialsort :: (Ord e, Arbitrary e, Show e)
-                 => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())
-                 -> Positive Int -> Property
-prop_partialsort = prop_sized $ \algo k ->
-  prop_sorted . V.take k . modify algo
-
-prop_sized_empty :: (Ord e) => (forall s. MV.MVector s e -> Int -> ST s ()) -> Property
-prop_sized_empty algo = prop_empty (flip algo 0) .&&. prop_empty (flip algo 10)
-
-prop_select :: (Ord e, Arbitrary e, Show e)
-            => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())
-            -> Positive Int -> Property
-prop_select = prop_sized $ \algo k arr ->
-  let vec' = modify algo arr
-      l    = V.slice 0 k vec'
-      r    = V.slice k (V.length vec' - k) vec'
-  in V.all (\e -> V.all (e <=) r) l
-
-prop_sized :: (Arbitrary e, Show e, Testable prop)
-           => ((forall s mv. G.MVector mv e => mv s e -> ST s ())
-                 -> Int -> Vector e -> prop)
-           -> (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())
-           -> Positive Int -> Property
-prop_sized prop algo (Positive k) =
-  let k' = k `mod` sanity
-  in forAll (longGen k') $ prop (\marr -> algo marr k') k'
-
-prop_stable :: (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> ST s ())
-            -> Vector Int -> Property
--- prop_stable algo arr = property $ modify algo arr == arr
-prop_stable algo arr = stable $ modify (algo (comparing fst)) $ V.zip arr ix
- where
- ix = V.fromList [1 .. V.length arr]
-
-stable arr | V.null arr = property True
-           | otherwise  = let (e, i) = V.head arr
-                          in V.all (\(e', i') -> e < e' || i < i') (V.tail arr)
-                            .&. stable (V.tail arr)
-
-prop_stable_radix :: (forall e s mv. G.MVector mv e => Int -> Int -> (Int -> e -> Int) 
-                        -> mv s e -> ST s ())
-                  -> Vector Int -> Property
-prop_stable_radix algo arr =
-  stable . modify (algo (passes e) (size e) (\k (e, _) -> radix k e))
-         $ V.zip arr ix
- where
- ix = V.fromList [1 .. V.length arr]
- e = V.head arr
- 
-prop_optimal :: Int
-             -> (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> Int -> ST s ())
-             -> Property
-prop_optimal n algo = label "sorting" sortn .&. label "stability" stabn
- where
- arrn  = V.fromList [0..n-1]
- sortn = all ( (== arrn)
-             . modify (\a -> algo compare a 0)
-             . V.fromList)
-         $ permutations [0..n-1]
- stabn = all ( (== arrn)
-             . snd
-             . V.unzip
-             . modify (\a -> algo (comparing fst) a 0))
-         $ stability n
-
-type Bag e = M.Map e Int
-
-toBag :: (Ord e) => Vector e -> Bag e
-toBag = M.fromListWith (+) . flip zip (repeat 1) . V.toList
-
-prop_permutation :: (Ord e) => (forall s mv. G.MVector mv e => mv s e -> ST s ())
-                 -> Vector e -> Property
-prop_permutation algo arr = property $ 
-                            toBag arr == toBag (modify algo arr)
-
-newtype SortedVec e = Sorted (Vector e)
-
-instance (Show e) => Show (SortedVec e) where
-  show (Sorted a) = show a
-
-instance (Arbitrary e, Ord e) => Arbitrary (SortedVec e) where
-  arbitrary = fmap (Sorted . V.fromList . sort)
-                $ liftM2 (++) (vectorOf 20 arbitrary) arbitrary
-
-ixRanges :: Vector e -> Gen (Int, Int)
-ixRanges vec = do i <- fmap (`mod` len) arbitrary
-                  j <- fmap (`mod` len) arbitrary
-                  return $ if i < j then (i, j) else (j, i)
- where len = V.length vec
-
-prop_search_inrange :: (Ord e)
-                    => (forall s. MVector s e -> e -> Int -> Int -> ST s Int)
-                    -> SortedVec e -> e -> Property
-prop_search_inrange algo (Sorted arr) e = forAll (ixRanges arr) $ \(i, j) ->
-  let k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e i j)
-  in property $ i <= k && k <= j
- where
- len = V.length arr
-
-prop_search_insert :: (e -> e -> Bool) -> (e -> e -> Bool)
-                   -> (forall s. MVector s e -> e -> ST s Int)
-                   -> SortedVec e -> e -> Property
-prop_search_insert lo hi algo (Sorted arr) e =
-  property $ (k == 0   || (arr V.! (k-1)) `lo` e)
-          && (k == len || (arr V.! k) `hi` e)
- where
- len = V.length arr
- k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e)
-
-prop_search_lowbound :: (Ord e)
-                     => (forall s. MVector s e -> e -> ST s Int)
-                     -> SortedVec e -> e -> Property
-prop_search_lowbound = prop_search_insert (<) (>=)
-
-prop_search_upbound :: (Ord e)
-                    => (forall s. MVector s e -> e -> ST s Int)
-                    -> SortedVec e -> e -> Property
-prop_search_upbound = prop_search_insert (<=) (>)
diff --git a/tests/Tests.hs b/tests/Tests.hs
deleted file mode 100644
--- a/tests/Tests.hs
+++ /dev/null
@@ -1,197 +0,0 @@
-{-# LANGUAGE ImpredicativeTypes, RankNTypes, TypeOperators, FlexibleContexts #-}
-
-module Main (main) where
-
-import Properties
-
-import Util
-
-import Test.QuickCheck
-
-import Control.Monad
-import Control.Monad.ST
-
-import Data.Int
-import Data.Word
-
-import qualified Data.ByteString as B
-
-import Data.Vector (Vector)
-import qualified Data.Vector as V
-
-import Data.Vector.Generic.Mutable (MVector)
-import qualified Data.Vector.Generic.Mutable as MV
-
-import qualified Data.Vector.Algorithms.Insertion    as INS
-import qualified Data.Vector.Algorithms.Intro        as INT
-import qualified Data.Vector.Algorithms.Merge        as M
-import qualified Data.Vector.Algorithms.Radix        as R
-import qualified Data.Vector.Algorithms.Heap         as H
-import qualified Data.Vector.Algorithms.Optimal      as O
-import qualified Data.Vector.Algorithms.AmericanFlag as AF
-
-import qualified Data.Vector.Algorithms.Search       as SR
-
-type Algo      e r = forall s mv. MVector mv e => mv s e -> ST s r
-type SizeAlgo  e r = forall s mv. MVector mv e => mv s e -> Int -> ST s r
-type BoundAlgo e r = forall s mv. MVector mv e => mv s e -> Int -> Int -> ST s r
-
-args = stdArgs
-       { maxSuccess = 1000
-       , maxDiscard = 200
-       }
-
-check_Int_sort = forM_ algos $ \(name,algo) ->
-  quickCheckWith args (label name . prop_fullsort algo)
- where
- algos :: [(String, Algo Int ())]
- algos = [ ("introsort", INT.sort)
-         , ("insertion sort", INS.sort)
-         , ("merge sort", M.sort)
-         , ("heapsort", H.sort)
-         ]
-
-check_Int_partialsort = forM_ algos $ \(name,algo) ->
-  quickCheckWith args (label name . prop_partialsort algo)
- where
- algos :: [(String, SizeAlgo Int ())]
- algos = [ ("intro-partialsort", INT.partialSort)
-         , ("heap partialsort", H.partialSort)
-         ]
-
-check_Int_select = forM_ algos $ \(name,algo) ->
-  quickCheckWith args (label name . prop_select algo)
- where
- algos :: [(String, SizeAlgo Int ())]
- algos = [ ("intro-select", INT.select)
-         , ("heap select", H.select)
-         ]
-
-check_radix_sorts = do
-  qc (label "radix Word8"       . prop_fullsort (R.sort :: Algo Word8  ()))
-  qc (label "radix Word16"      . prop_fullsort (R.sort :: Algo Word16 ()))
-  qc (label "radix Word32"      . prop_fullsort (R.sort :: Algo Word32 ()))
-  qc (label "radix Word64"      . prop_fullsort (R.sort :: Algo Word64 ()))
-  qc (label "radix Word"        . prop_fullsort (R.sort :: Algo Word   ()))
-  qc (label "radix Int8"        . prop_fullsort (R.sort :: Algo Int8   ()))
-  qc (label "radix Int16"       . prop_fullsort (R.sort :: Algo Int16  ()))
-  qc (label "radix Int32"       . prop_fullsort (R.sort :: Algo Int32  ()))
-  qc (label "radix Int64"       . prop_fullsort (R.sort :: Algo Int64  ()))
-  qc (label "radix Int"         . prop_fullsort (R.sort :: Algo Int    ()))
-  qc (label "radix (Int, Int)"  . prop_fullsort (R.sort :: Algo (Int, Int) ()))
-
-  qc (label "flag Word8"       . prop_fullsort (AF.sort :: Algo Word8  ()))
-  qc (label "flag Word16"      . prop_fullsort (AF.sort :: Algo Word16 ()))
-  qc (label "flag Word32"      . prop_fullsort (AF.sort :: Algo Word32 ()))
-  qc (label "flag Word64"      . prop_fullsort (AF.sort :: Algo Word64 ()))
-  qc (label "flag Word"        . prop_fullsort (AF.sort :: Algo Word   ()))
-  qc (label "flag Int8"        . prop_fullsort (AF.sort :: Algo Int8   ()))
-  qc (label "flag Int16"       . prop_fullsort (AF.sort :: Algo Int16  ()))
-  qc (label "flag Int32"       . prop_fullsort (AF.sort :: Algo Int32  ()))
-  qc (label "flag Int64"       . prop_fullsort (AF.sort :: Algo Int64  ()))
-  qc (label "flag Int"         . prop_fullsort (AF.sort :: Algo Int    ()))
-  qc (label "flag ByteString"  . prop_fullsort (AF.sort :: Algo B.ByteString ()))
- where
- qc algo = quickCheckWith args algo
-
-{-
-check_schwartzian = do
-  quickCheckWith args (prop_schwartzian i2w INS.sortBy)
- where
- i2w :: Int -> Word
- i2w = fromIntegral
--}
-
-check_stable = do quickCheckWith args (label "merge sort" . prop_stable M.sortBy)
-                  quickCheckWith args (label "radix sort" . prop_stable_radix R.sortBy)
-
-check_optimal = do qc . label "size 2" $ prop_optimal 2 O.sort2ByOffset
-                   qc . label "size 3" $ prop_optimal 3 O.sort3ByOffset
-                   qc . label "size 4" $ prop_optimal 4 O.sort4ByOffset
- where
- qc = quickCheck
-
-check_permutation = do
-  qc $ label "introsort"    . prop_permutation (INT.sort :: Algo Int ())
-  qc $ label "intropartial" . prop_sized (const . prop_permutation)
-                                         (INT.partialSort :: SizeAlgo Int ())
-  qc $ label "introselect"  . prop_sized (const . prop_permutation)
-                                         (INT.select :: SizeAlgo Int ())
-  qc $ label "heapsort"     . prop_permutation (H.sort :: Algo Int ())
-  qc $ label "heappartial"  . prop_sized (const . prop_permutation)
-                                         (H.partialSort :: SizeAlgo Int ())
-  qc $ label "heapselect"   . prop_sized (const . prop_permutation)
-                                         (H.select :: SizeAlgo Int ())
-  qc $ label "mergesort"    . prop_permutation (M.sort :: Algo Int    ())
-  qc $ label "radix I8"     . prop_permutation (R.sort :: Algo Int8   ())
-  qc $ label "radix I16"    . prop_permutation (R.sort :: Algo Int16  ())
-  qc $ label "radix I32"    . prop_permutation (R.sort :: Algo Int32  ())
-  qc $ label "radix I64"    . prop_permutation (R.sort :: Algo Int64  ())
-  qc $ label "radix Int"    . prop_permutation (R.sort :: Algo Int    ())
-  qc $ label "radix W8"     . prop_permutation (R.sort :: Algo Word8  ())
-  qc $ label "radix W16"    . prop_permutation (R.sort :: Algo Word16 ())
-  qc $ label "radix W32"    . prop_permutation (R.sort :: Algo Word32 ())
-  qc $ label "radix W64"    . prop_permutation (R.sort :: Algo Word64 ())
-  qc $ label "radix Word"   . prop_permutation (R.sort :: Algo Word   ())
-  qc $ label "flag I8"      . prop_permutation (AF.sort :: Algo Int8   ())
-  qc $ label "flag I16"     . prop_permutation (AF.sort :: Algo Int16  ())
-  qc $ label "flag I32"     . prop_permutation (AF.sort :: Algo Int32  ())
-  qc $ label "flag I64"     . prop_permutation (AF.sort :: Algo Int64  ())
-  qc $ label "flag Int"     . prop_permutation (AF.sort :: Algo Int    ())
-  qc $ label "flag W8"      . prop_permutation (AF.sort :: Algo Word8  ())
-  qc $ label "flag W16"     . prop_permutation (AF.sort :: Algo Word16 ())
-  qc $ label "flag W32"     . prop_permutation (AF.sort :: Algo Word32 ())
-  qc $ label "flag W64"     . prop_permutation (AF.sort :: Algo Word64 ())
-  qc $ label "flag Word"    . prop_permutation (AF.sort :: Algo Word   ())
-  qc $ label "flag ByteString" . prop_permutation (AF.sort :: Algo B.ByteString ())
- where
- qc prop = quickCheckWith args prop
-
-check_corners = do
-  qc "introsort empty"    $ prop_empty       (INT.sort        :: Algo Int ())
-  qc "intropartial empty" $ prop_sized_empty (INT.partialSort :: SizeAlgo Int ())
-  qc "introselect empty"  $ prop_sized_empty (INT.select      :: SizeAlgo Int ())
-  qc "heapsort empty"     $ prop_empty       (H.sort          :: Algo Int ())
-  qc "heappartial empty"  $ prop_sized_empty (H.partialSort   :: SizeAlgo Int ())
-  qc "heapselect empty"   $ prop_sized_empty (H.select        :: SizeAlgo Int ())
-  qc "mergesort empty"    $ prop_empty       (M.sort          :: Algo Int ())
-  qc "radixsort empty"    $ prop_empty       (R.sort          :: Algo Int ())
-  qc "flagsort empty"     $ prop_empty       (AF.sort         :: Algo Int ())
- where
- qc s prop = quickCheckWith (stdArgs { maxSuccess = 2 }) (label s prop)
-
-type SAlgo e r = forall s mv. MVector mv e => mv s e -> e -> ST s r
-type BoundSAlgo e r = forall s mv. MVector mv e => mv s e -> e -> Int -> Int -> ST s r
-
-check_search_range = do
-  qc $ (label "binarySearchL" .)
-         . prop_search_inrange (SR.binarySearchLByBounds compare :: BoundSAlgo Int Int)
-  qc $ (label "binarySearchL lo-bound" .)
-         . prop_search_lowbound (SR.binarySearchL :: SAlgo Int Int)
-  qc $ (label "binarySearch" .)
-         . prop_search_inrange (SR.binarySearchByBounds compare :: BoundSAlgo Int Int)
-  qc $ (label "binarySearchR" .)
-         . prop_search_inrange (SR.binarySearchRByBounds compare :: BoundSAlgo Int Int)
-  qc $ (label "binarySearchR hi-bound" .)
-         . prop_search_upbound (SR.binarySearchR :: SAlgo Int Int)
- where
- qc prop = quickCheckWith args prop
-
-main = do putStrLn "Int tests:"
-          check_Int_sort
-          check_Int_partialsort
-          check_Int_select
-          putStrLn "Radix sort tests:"
-          check_radix_sorts
---          putStrLn "Schwartzian transform (Int -> Word):"
---          check_schwartzian
-          putStrLn "Stability:"
-          check_stable
-          putStrLn "Optimals:"
-          check_optimal
-          putStrLn "Permutation:"
-          check_permutation
-          putStrLn "Search in range:"
-          check_search_range
-          putStrLn "Corner cases:"
-          check_corners
diff --git a/tests/Util.hs b/tests/Util.hs
deleted file mode 100644
--- a/tests/Util.hs
+++ /dev/null
@@ -1,33 +0,0 @@
-{-# LANGUAGE TypeOperators #-}
-
-module Util where
-
-import Control.Monad
-import Control.Monad.ST
-
-import Data.Word
-import Data.Int
-
-import qualified Data.ByteString as B
-
-import qualified Data.Vector as V
-
-import Data.Vector.Mutable hiding (length)
-
-import Test.QuickCheck
-
-
-mfromList :: [e] -> ST s (MVector s e)
-mfromList l = do v <- new (length l)
-                 fill l 0 v
- where
- fill []     _ v = return v
- fill (x:xs) i v = do write v i x
-                      fill xs (i+1) v
-
-instance (Arbitrary e) => Arbitrary (V.Vector e) where
-  arbitrary = fmap V.fromList arbitrary
-
-instance Arbitrary B.ByteString where
-  arbitrary = B.pack `fmap` arbitrary
-
diff --git a/tests/properties/Optimal.hs b/tests/properties/Optimal.hs
new file mode 100644
--- /dev/null
+++ b/tests/properties/Optimal.hs
@@ -0,0 +1,62 @@
+{-# LANGUAGE TypeOperators, FlexibleContexts #-}
+
+-- Exhaustive test sets for proper sorting and stability of
+-- optimal sorts
+
+module Optimal where
+
+import Control.Arrow
+import Control.Monad
+
+import Data.List
+import Data.Function
+
+import Data.Vector.Generic hiding (map, zip, concatMap, (++), replicate, foldM)
+
+interleavings :: [a] -> [a] -> [[a]]
+interleavings [       ] ys        =  [ys]
+interleavings xs        [       ] =  [xs]
+interleavings xs@(x:xt) ys@(y:yt) =  map (x:) (interleavings xt ys)
+                                  ++ map (y:) (interleavings xs yt)
+
+monotones :: Int -> Int -> [[Int]]
+monotones k = atLeastOne 0
+ where
+ atLeastOne i 0 = [[]]
+ atLeastOne i n = map (i:) $ picks i (n-1)
+ picks _ 0             = [[]]
+ picks i n | i >= k    = [replicate n k]
+           | otherwise = map (i:) (picks i (n-1)) ++ atLeastOne (i+1) n
+
+
+stability :: (Vector v (Int,Int)) => Int -> [v (Int, Int)]
+stability n = concatMap ( map fromList
+                        . foldM interleavings []
+                        . groupBy ((==) `on` fst)
+                        . flip zip [0..])
+              $ monotones (n-2) n
+
+sort2 :: (Vector v Int) => [v Int]
+sort2 = map fromList $ permutations [0,1]
+
+stability2 :: (Vector v (Int,Int)) => [v (Int, Int)]
+stability2 = [fromList [(0, 0), (0, 1)]]
+
+sort3 :: (Vector v Int) => [v Int]
+sort3 = map fromList $ permutations [0..2]
+
+{-
+stability3 :: [UArr (Int :*: Int)]
+stability3 = map toU [ [0:*:0, 0:*:1, 0:*:2]
+                     , [0:*:0, 0:*:1, 1:*:2]
+                     , [0:*:0, 1:*:2, 0:*:1]
+                     , [1:*:2, 0:*:0, 0:*:1]
+                     , [0:*:0, 1:*:1, 1:*:2]
+                     , [1:*:1, 0:*:0, 1:*:2]
+                     , [1:*:1, 1:*:2, 0:*:0]
+                     ]
+-}
+
+sort4 :: (Vector v Int) => [v Int]
+sort4 = map fromList $ permutations [0..3]
+
diff --git a/tests/properties/Properties.hs b/tests/properties/Properties.hs
new file mode 100644
--- /dev/null
+++ b/tests/properties/Properties.hs
@@ -0,0 +1,185 @@
+{-# LANGUAGE RankNTypes, FlexibleContexts #-}
+
+module Properties where
+
+import Prelude
+
+import Optimal
+
+import Control.Monad
+import Control.Monad.ST
+
+import Data.List
+import Data.Ord
+
+import Data.Vector (Vector)
+import qualified Data.Vector as V
+
+import Data.Vector.Mutable (MVector)
+import qualified Data.Vector.Mutable as MV
+
+import Data.Vector.Generic (modify)
+
+import qualified Data.Vector.Generic.Mutable as G
+
+import Data.Vector.Algorithms.Optimal (Comparison)
+import Data.Vector.Algorithms.Radix (radix, passes, size)
+
+import qualified Data.Map as M
+
+import Test.QuickCheck
+
+import Util
+
+prop_sorted :: (Ord e) => Vector e -> Property
+prop_sorted arr | V.length arr < 2 = property True
+                | otherwise        = check (V.head arr) (V.tail arr)
+ where
+ check e arr | V.null arr = property True
+             | otherwise  = e <= V.head arr .&. check (V.head arr) (V.tail arr)
+
+prop_empty :: (Ord e) => (forall s. MV.MVector s e -> ST s ()) -> Property
+prop_empty algo = prop_sorted (modify algo $ V.fromList [])
+
+prop_fullsort :: (Ord e)
+              => (forall s mv. G.MVector mv e => mv s e -> ST s ()) -> Vector e -> Property
+prop_fullsort algo arr = prop_sorted $ modify algo arr
+
+{-
+prop_schwartzian :: (UA e, UA k, Ord k)
+                 => (e -> k)
+                 -> (forall e s. (UA e) => (e -> e -> Ordering) -> MUArr e s -> ST s ())
+                 -> UArr e -> Property
+prop_schwartzian f algo arr
+  | lengthU arr < 2 = property True
+  | otherwise       = let srt = modify (algo `usingKeys` f) arr
+                      in check (headU srt) (tailU srt)
+ where
+ check e arr | nullU arr = property True
+             | otherwise = f e <= f (headU arr) .&. check (headU arr) (tailU arr)
+-}
+
+longGen :: (Arbitrary e) => Int -> Gen (Vector e)
+longGen k = liftM2 (\l r -> V.fromList (l ++ r)) (vectorOf k arbitrary) arbitrary
+
+sanity :: Int
+sanity = 100
+
+prop_partialsort :: (Ord e, Arbitrary e, Show e)
+                 => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())
+                 -> Positive Int -> Property
+prop_partialsort = prop_sized $ \algo k ->
+  prop_sorted . V.take k . modify algo
+
+prop_sized_empty :: (Ord e) => (forall s. MV.MVector s e -> Int -> ST s ()) -> Property
+prop_sized_empty algo = prop_empty (flip algo 0) .&&. prop_empty (flip algo 10)
+
+prop_select :: (Ord e, Arbitrary e, Show e)
+            => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())
+            -> Positive Int -> Property
+prop_select = prop_sized $ \algo k arr ->
+  let vec' = modify algo arr
+      l    = V.slice 0 k vec'
+      r    = V.slice k (V.length vec' - k) vec'
+  in V.all (\e -> V.all (e <=) r) l
+
+prop_sized :: (Arbitrary e, Show e, Testable prop)
+           => ((forall s mv. G.MVector mv e => mv s e -> ST s ())
+                 -> Int -> Vector e -> prop)
+           -> (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())
+           -> Positive Int -> Property
+prop_sized prop algo (Positive k) =
+  let k' = k `mod` sanity
+  in forAll (longGen k') $ prop (\marr -> algo marr k') k'
+
+prop_stable :: (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> ST s ())
+            -> Vector Int -> Property
+-- prop_stable algo arr = property $ modify algo arr == arr
+prop_stable algo arr = stable $ modify (algo (comparing fst)) $ V.zip arr ix
+ where
+ ix = V.fromList [1 .. V.length arr]
+
+stable arr | V.null arr = property True
+           | otherwise  = let (e, i) = V.head arr
+                          in V.all (\(e', i') -> e < e' || i < i') (V.tail arr)
+                            .&. stable (V.tail arr)
+
+prop_stable_radix :: (forall e s mv. G.MVector mv e => Int -> Int -> (Int -> e -> Int) 
+                        -> mv s e -> ST s ())
+                  -> Vector Int -> Property
+prop_stable_radix algo arr =
+  stable . modify (algo (passes e) (size e) (\k (e, _) -> radix k e))
+         $ V.zip arr ix
+ where
+ ix = V.fromList [1 .. V.length arr]
+ e = V.head arr
+ 
+prop_optimal :: Int
+             -> (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> Int -> ST s ())
+             -> Property
+prop_optimal n algo = label "sorting" sortn .&. label "stability" stabn
+ where
+ arrn  = V.fromList [0..n-1]
+ sortn = all ( (== arrn)
+             . modify (\a -> algo compare a 0)
+             . V.fromList)
+         $ permutations [0..n-1]
+ stabn = all ( (== arrn)
+             . snd
+             . V.unzip
+             . modify (\a -> algo (comparing fst) a 0))
+         $ stability n
+
+type Bag e = M.Map e Int
+
+toBag :: (Ord e) => Vector e -> Bag e
+toBag = M.fromListWith (+) . flip zip (repeat 1) . V.toList
+
+prop_permutation :: (Ord e) => (forall s mv. G.MVector mv e => mv s e -> ST s ())
+                 -> Vector e -> Property
+prop_permutation algo arr = property $ 
+                            toBag arr == toBag (modify algo arr)
+
+newtype SortedVec e = Sorted (Vector e)
+
+instance (Show e) => Show (SortedVec e) where
+  show (Sorted a) = show a
+
+instance (Arbitrary e, Ord e) => Arbitrary (SortedVec e) where
+  arbitrary = fmap (Sorted . V.fromList . sort)
+                $ liftM2 (++) (vectorOf 20 arbitrary) arbitrary
+
+ixRanges :: Vector e -> Gen (Int, Int)
+ixRanges vec = do i <- fmap (`mod` len) arbitrary
+                  j <- fmap (`mod` len) arbitrary
+                  return $ if i < j then (i, j) else (j, i)
+ where len = V.length vec
+
+prop_search_inrange :: (Ord e)
+                    => (forall s. MVector s e -> e -> Int -> Int -> ST s Int)
+                    -> SortedVec e -> e -> Property
+prop_search_inrange algo (Sorted arr) e = forAll (ixRanges arr) $ \(i, j) ->
+  let k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e i j)
+  in property $ i <= k && k <= j
+ where
+ len = V.length arr
+
+prop_search_insert :: (e -> e -> Bool) -> (e -> e -> Bool)
+                   -> (forall s. MVector s e -> e -> ST s Int)
+                   -> SortedVec e -> e -> Property
+prop_search_insert lo hi algo (Sorted arr) e =
+  property $ (k == 0   || (arr V.! (k-1)) `lo` e)
+          && (k == len || (arr V.! k) `hi` e)
+ where
+ len = V.length arr
+ k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e)
+
+prop_search_lowbound :: (Ord e)
+                     => (forall s. MVector s e -> e -> ST s Int)
+                     -> SortedVec e -> e -> Property
+prop_search_lowbound = prop_search_insert (<) (>=)
+
+prop_search_upbound :: (Ord e)
+                    => (forall s. MVector s e -> e -> ST s Int)
+                    -> SortedVec e -> e -> Property
+prop_search_upbound = prop_search_insert (<=) (>)
diff --git a/tests/properties/Tests.hs b/tests/properties/Tests.hs
new file mode 100644
--- /dev/null
+++ b/tests/properties/Tests.hs
@@ -0,0 +1,197 @@
+{-# LANGUAGE ImpredicativeTypes, RankNTypes, TypeOperators, FlexibleContexts #-}
+
+module Main (main) where
+
+import Properties
+
+import Util
+
+import Test.QuickCheck
+
+import Control.Monad
+import Control.Monad.ST
+
+import Data.Int
+import Data.Word
+
+import qualified Data.ByteString as B
+
+import Data.Vector (Vector)
+import qualified Data.Vector as V
+
+import Data.Vector.Generic.Mutable (MVector)
+import qualified Data.Vector.Generic.Mutable as MV
+
+import qualified Data.Vector.Algorithms.Insertion    as INS
+import qualified Data.Vector.Algorithms.Intro        as INT
+import qualified Data.Vector.Algorithms.Merge        as M
+import qualified Data.Vector.Algorithms.Radix        as R
+import qualified Data.Vector.Algorithms.Heap         as H
+import qualified Data.Vector.Algorithms.Optimal      as O
+import qualified Data.Vector.Algorithms.AmericanFlag as AF
+
+import qualified Data.Vector.Algorithms.Search       as SR
+
+type Algo      e r = forall s mv. MVector mv e => mv s e -> ST s r
+type SizeAlgo  e r = forall s mv. MVector mv e => mv s e -> Int -> ST s r
+type BoundAlgo e r = forall s mv. MVector mv e => mv s e -> Int -> Int -> ST s r
+
+args = stdArgs
+       { maxSuccess = 1000
+       , maxDiscardRatio = 2
+       }
+
+check_Int_sort = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_fullsort algo)
+ where
+ algos :: [(String, Algo Int ())]
+ algos = [ ("introsort", INT.sort)
+         , ("insertion sort", INS.sort)
+         , ("merge sort", M.sort)
+         , ("heapsort", H.sort)
+         ]
+
+check_Int_partialsort = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_partialsort algo)
+ where
+ algos :: [(String, SizeAlgo Int ())]
+ algos = [ ("intro-partialsort", INT.partialSort)
+         , ("heap partialsort", H.partialSort)
+         ]
+
+check_Int_select = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_select algo)
+ where
+ algos :: [(String, SizeAlgo Int ())]
+ algos = [ ("intro-select", INT.select)
+         , ("heap select", H.select)
+         ]
+
+check_radix_sorts = do
+  qc (label "radix Word8"       . prop_fullsort (R.sort :: Algo Word8  ()))
+  qc (label "radix Word16"      . prop_fullsort (R.sort :: Algo Word16 ()))
+  qc (label "radix Word32"      . prop_fullsort (R.sort :: Algo Word32 ()))
+  qc (label "radix Word64"      . prop_fullsort (R.sort :: Algo Word64 ()))
+  qc (label "radix Word"        . prop_fullsort (R.sort :: Algo Word   ()))
+  qc (label "radix Int8"        . prop_fullsort (R.sort :: Algo Int8   ()))
+  qc (label "radix Int16"       . prop_fullsort (R.sort :: Algo Int16  ()))
+  qc (label "radix Int32"       . prop_fullsort (R.sort :: Algo Int32  ()))
+  qc (label "radix Int64"       . prop_fullsort (R.sort :: Algo Int64  ()))
+  qc (label "radix Int"         . prop_fullsort (R.sort :: Algo Int    ()))
+  qc (label "radix (Int, Int)"  . prop_fullsort (R.sort :: Algo (Int, Int) ()))
+
+  qc (label "flag Word8"       . prop_fullsort (AF.sort :: Algo Word8  ()))
+  qc (label "flag Word16"      . prop_fullsort (AF.sort :: Algo Word16 ()))
+  qc (label "flag Word32"      . prop_fullsort (AF.sort :: Algo Word32 ()))
+  qc (label "flag Word64"      . prop_fullsort (AF.sort :: Algo Word64 ()))
+  qc (label "flag Word"        . prop_fullsort (AF.sort :: Algo Word   ()))
+  qc (label "flag Int8"        . prop_fullsort (AF.sort :: Algo Int8   ()))
+  qc (label "flag Int16"       . prop_fullsort (AF.sort :: Algo Int16  ()))
+  qc (label "flag Int32"       . prop_fullsort (AF.sort :: Algo Int32  ()))
+  qc (label "flag Int64"       . prop_fullsort (AF.sort :: Algo Int64  ()))
+  qc (label "flag Int"         . prop_fullsort (AF.sort :: Algo Int    ()))
+  qc (label "flag ByteString"  . prop_fullsort (AF.sort :: Algo B.ByteString ()))
+ where
+ qc algo = quickCheckWith args algo
+
+{-
+check_schwartzian = do
+  quickCheckWith args (prop_schwartzian i2w INS.sortBy)
+ where
+ i2w :: Int -> Word
+ i2w = fromIntegral
+-}
+
+check_stable = do quickCheckWith args (label "merge sort" . prop_stable M.sortBy)
+                  quickCheckWith args (label "radix sort" . prop_stable_radix R.sortBy)
+
+check_optimal = do qc . label "size 2" $ prop_optimal 2 O.sort2ByOffset
+                   qc . label "size 3" $ prop_optimal 3 O.sort3ByOffset
+                   qc . label "size 4" $ prop_optimal 4 O.sort4ByOffset
+ where
+ qc = quickCheck
+
+check_permutation = do
+  qc $ label "introsort"    . prop_permutation (INT.sort :: Algo Int ())
+  qc $ label "intropartial" . prop_sized (const . prop_permutation)
+                                         (INT.partialSort :: SizeAlgo Int ())
+  qc $ label "introselect"  . prop_sized (const . prop_permutation)
+                                         (INT.select :: SizeAlgo Int ())
+  qc $ label "heapsort"     . prop_permutation (H.sort :: Algo Int ())
+  qc $ label "heappartial"  . prop_sized (const . prop_permutation)
+                                         (H.partialSort :: SizeAlgo Int ())
+  qc $ label "heapselect"   . prop_sized (const . prop_permutation)
+                                         (H.select :: SizeAlgo Int ())
+  qc $ label "mergesort"    . prop_permutation (M.sort :: Algo Int    ())
+  qc $ label "radix I8"     . prop_permutation (R.sort :: Algo Int8   ())
+  qc $ label "radix I16"    . prop_permutation (R.sort :: Algo Int16  ())
+  qc $ label "radix I32"    . prop_permutation (R.sort :: Algo Int32  ())
+  qc $ label "radix I64"    . prop_permutation (R.sort :: Algo Int64  ())
+  qc $ label "radix Int"    . prop_permutation (R.sort :: Algo Int    ())
+  qc $ label "radix W8"     . prop_permutation (R.sort :: Algo Word8  ())
+  qc $ label "radix W16"    . prop_permutation (R.sort :: Algo Word16 ())
+  qc $ label "radix W32"    . prop_permutation (R.sort :: Algo Word32 ())
+  qc $ label "radix W64"    . prop_permutation (R.sort :: Algo Word64 ())
+  qc $ label "radix Word"   . prop_permutation (R.sort :: Algo Word   ())
+  qc $ label "flag I8"      . prop_permutation (AF.sort :: Algo Int8   ())
+  qc $ label "flag I16"     . prop_permutation (AF.sort :: Algo Int16  ())
+  qc $ label "flag I32"     . prop_permutation (AF.sort :: Algo Int32  ())
+  qc $ label "flag I64"     . prop_permutation (AF.sort :: Algo Int64  ())
+  qc $ label "flag Int"     . prop_permutation (AF.sort :: Algo Int    ())
+  qc $ label "flag W8"      . prop_permutation (AF.sort :: Algo Word8  ())
+  qc $ label "flag W16"     . prop_permutation (AF.sort :: Algo Word16 ())
+  qc $ label "flag W32"     . prop_permutation (AF.sort :: Algo Word32 ())
+  qc $ label "flag W64"     . prop_permutation (AF.sort :: Algo Word64 ())
+  qc $ label "flag Word"    . prop_permutation (AF.sort :: Algo Word   ())
+  qc $ label "flag ByteString" . prop_permutation (AF.sort :: Algo B.ByteString ())
+ where
+ qc prop = quickCheckWith args prop
+
+check_corners = do
+  qc "introsort empty"    $ prop_empty       (INT.sort        :: Algo Int ())
+  qc "intropartial empty" $ prop_sized_empty (INT.partialSort :: SizeAlgo Int ())
+  qc "introselect empty"  $ prop_sized_empty (INT.select      :: SizeAlgo Int ())
+  qc "heapsort empty"     $ prop_empty       (H.sort          :: Algo Int ())
+  qc "heappartial empty"  $ prop_sized_empty (H.partialSort   :: SizeAlgo Int ())
+  qc "heapselect empty"   $ prop_sized_empty (H.select        :: SizeAlgo Int ())
+  qc "mergesort empty"    $ prop_empty       (M.sort          :: Algo Int ())
+  qc "radixsort empty"    $ prop_empty       (R.sort          :: Algo Int ())
+  qc "flagsort empty"     $ prop_empty       (AF.sort         :: Algo Int ())
+ where
+ qc s prop = quickCheckWith (stdArgs { maxSuccess = 2 }) (label s prop)
+
+type SAlgo e r = forall s mv. MVector mv e => mv s e -> e -> ST s r
+type BoundSAlgo e r = forall s mv. MVector mv e => mv s e -> e -> Int -> Int -> ST s r
+
+check_search_range = do
+  qc $ (label "binarySearchL" .)
+         . prop_search_inrange (SR.binarySearchLByBounds compare :: BoundSAlgo Int Int)
+  qc $ (label "binarySearchL lo-bound" .)
+         . prop_search_lowbound (SR.binarySearchL :: SAlgo Int Int)
+  qc $ (label "binarySearch" .)
+         . prop_search_inrange (SR.binarySearchByBounds compare :: BoundSAlgo Int Int)
+  qc $ (label "binarySearchR" .)
+         . prop_search_inrange (SR.binarySearchRByBounds compare :: BoundSAlgo Int Int)
+  qc $ (label "binarySearchR hi-bound" .)
+         . prop_search_upbound (SR.binarySearchR :: SAlgo Int Int)
+ where
+ qc prop = quickCheckWith args prop
+
+main = do putStrLn "Int tests:"
+          check_Int_sort
+          check_Int_partialsort
+          check_Int_select
+          putStrLn "Radix sort tests:"
+          check_radix_sorts
+--          putStrLn "Schwartzian transform (Int -> Word):"
+--          check_schwartzian
+          putStrLn "Stability:"
+          check_stable
+          putStrLn "Optimals:"
+          check_optimal
+          putStrLn "Permutation:"
+          check_permutation
+          putStrLn "Search in range:"
+          check_search_range
+          putStrLn "Corner cases:"
+          check_corners
diff --git a/tests/properties/Util.hs b/tests/properties/Util.hs
new file mode 100644
--- /dev/null
+++ b/tests/properties/Util.hs
@@ -0,0 +1,33 @@
+{-# LANGUAGE TypeOperators #-}
+
+module Util where
+
+import Control.Monad
+import Control.Monad.ST
+
+import Data.Word
+import Data.Int
+
+import qualified Data.ByteString as B
+
+import qualified Data.Vector as V
+
+import Data.Vector.Mutable hiding (length)
+
+import Test.QuickCheck
+
+
+mfromList :: [e] -> ST s (MVector s e)
+mfromList l = do v <- new (length l)
+                 fill l 0 v
+ where
+ fill []     _ v = return v
+ fill (x:xs) i v = do write v i x
+                      fill xs (i+1) v
+
+instance (Arbitrary e) => Arbitrary (V.Vector e) where
+  arbitrary = fmap V.fromList arbitrary
+
+instance Arbitrary B.ByteString where
+  arbitrary = B.pack `fmap` arbitrary
+
diff --git a/vector-algorithms.cabal b/vector-algorithms.cabal
--- a/vector-algorithms.cabal
+++ b/vector-algorithms.cabal
@@ -1,73 +1,115 @@
-Name:              vector-algorithms
-Version:           0.5.4.2
-License:           BSD3
-License-File:      LICENSE
-Author:            Dan Doel
-Maintainer:        Dan Doel <dan.doel@gmail.com>
-Homepage:          http://code.haskell.org/~dolio/
-Category:          Data
-Synopsis:          Efficient algorithms for vector arrays
-Description:       Efficient algorithms for vector arrays
-Build-Type:        Simple
-Cabal-Version:     >= 1.2.3
+name:              vector-algorithms
+version:           0.6.0.1
+license:           BSD3
+license-file:      LICENSE
+author:            Dan Doel
+maintainer:        Dan Doel <dan.doel@gmail.com>
+copyright:         (c) 2008,2009,2010,2011,2012,2013 Dan Doel
+homepage:          http://code.haskell.org/~dolio/
+category:          Data
+synopsis:          Efficient algorithms for vector arrays
+description:       Efficient algorithms for vector arrays
+build-type:        Simple
+cabal-version:     >= 1.9.2
 
-Flag BoundsChecks
-  Description: Enable bounds checking
-  Default: True
+flag BoundsChecks
+  description: Enable bounds checking
+  default: True
 
-Flag UnsafeChecks
-  Description: Enable bounds checking in unsafe operations at the cost of a
+flag UnsafeChecks
+  description: Enable bounds checking in unsafe operations at the cost of a
                significant performance penalty.
-  Default: False
+  default: False
 
-Flag InternalChecks
-  Description: Enable internal consistency checks at the cost of a
+flag InternalChecks
+  description: Enable internal consistency checks at the cost of a
                significant performance penalty.
-  Default: False
+  default: False
 
-Library
-    Build-Depends: base >= 3 && < 5,
-                   vector >= 0.6 && < 0.11,
-                   primitive >=0.3 && <0.6,
-                   bytestring >= 0.9 && < 1.0
+flag bench
+  description: Build a benchmarking program to test vector-algorithms
+               performance
+  default: False
 
-    Exposed-Modules:
-        Data.Vector.Algorithms.Optimal
-        Data.Vector.Algorithms.Insertion
-        Data.Vector.Algorithms.Intro
-        Data.Vector.Algorithms.Merge
-        Data.Vector.Algorithms.Radix
-        Data.Vector.Algorithms.Search
-        Data.Vector.Algorithms.Heap
-        Data.Vector.Algorithms.AmericanFlag
+flag properties
+  description: Enable the quickcheck tests
+  default: True
 
-    Other-Modules:
-        Data.Vector.Algorithms.Common
+source-repository head
+  type:     darcs
+  location: http://hub.darcs.net/dolio/vector-algorithms
 
-    Extensions:
-        BangPatterns,
-        TypeOperators,
-        Rank2Types,
-        ScopedTypeVariables,
-        FlexibleContexts,
-        CPP
+library
+  hs-source-dirs: src
 
-    GHC-Options:
-        -Odph
-        -funbox-strict-fields
+  build-depends: base >= 3 && < 5,
+                 vector >= 0.6 && < 0.11,
+                 primitive >=0.3 && <0.6,
+                 bytestring >= 0.9 && < 1.0
 
-    Include-Dirs:
-        include
+  exposed-modules:
+    Data.Vector.Algorithms.Optimal
+    Data.Vector.Algorithms.Insertion
+    Data.Vector.Algorithms.Intro
+    Data.Vector.Algorithms.Merge
+    Data.Vector.Algorithms.Radix
+    Data.Vector.Algorithms.Search
+    Data.Vector.Algorithms.Heap
+    Data.Vector.Algorithms.AmericanFlag
 
-    Install-Includes:
-        vector.h
+  other-modules:
+    Data.Vector.Algorithms.Common
 
-    if flag(BoundsChecks)
-        cpp-options: -DVECTOR_BOUNDS_CHECKS
+  ghc-options:
+    -Odph
+    -funbox-strict-fields
 
-    if flag(UnsafeChecks)
-        cpp-options: -DVECTOR_UNSAFE_CHECKS
+  include-dirs:
+    include
 
-    if flag(InternalChecks)
-        cpp-options: -DVECTOR_INTERNAL_CHECKS
+  install-includes:
+    vector.h
 
+  if flag(BoundsChecks)
+    cpp-options: -DVECTOR_BOUNDS_CHECKS
+
+  if flag(UnsafeChecks)
+    cpp-options: -DVECTOR_UNSAFE_CHECKS
+
+  if flag(InternalChecks)
+    cpp-options: -DVECTOR_INTERNAL_CHECKS
+
+executable vector-algorithms-bench
+  hs-source-dirs: bench
+
+  if !flag(bench)
+    buildable: False
+
+  main-is: Main.hs
+
+  other-modules:
+    Blocks
+
+  build-depends: base, mwc-random, vector, vector-algorithms, mtl
+  ghc-options: -Wall -Odph
+
+test-suite properties
+  hs-source-dirs: tests/properties
+  type: exitcode-stdio-1.0
+  main-is: Tests.hs
+
+  other-modules:
+    Optimal
+    Properties
+    Util
+
+  if !flag(properties)
+    buildable: False
+  else
+    build-depends:
+      base,
+      bytestring,
+      containers,
+      QuickCheck >= 2,
+      vector,
+      vector-algorithms