diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,37 @@
+## Version 0.9.1.0 (2025-02-05)
+
+- More inlining for `sort` and `nib` functions.
+
+## Version 0.9.0.3 (2024-11-25)
+
+- Fix an off-by-one error Heap.partialSort functions.
+- Support latest ghcs.
+
+## Version 0.9.0.2 (2024-05-23)
+
+- Add `TypeOperators` pragma where needed.
+
+## Version 0.9.0.1 (2022-07-28)
+
+- Allow building with vector-0.13.*.
+
+## Version 0.9.0.0 (2022-05-19)
+
+- Add nub related functions.
+- Add sortUniq related functions (sorts, then removes duplicates).
+
+## Version 0.8.0.4 (2020-12-06)
+
+- Fix out of range access in Intro.partialSort.
+- Update QuickCheck dependency bounds.
+
+## Version 0.8.0.3 (2019-12-02)
+
+- Fix out-of-bounds access in Timsort.
+
+## Version 0.8.0.2 (2019-11-28)
+
+- Bump upper bounds on primitive and QuickCheck.
+- Expose 'terminate' function from 'AmericanFlag' module.
+- Fix an off-by-one error in Data.Vector.Algorithms.Heaps.heapInsert.
+
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/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,5 @@
-Copyright (c) 2008-2010 Dan Doel
+Copyright (c) 2015 Dan Doel
+Copyright (c) 2015 Tim Baumann
 
 All rights reserved.
 
@@ -32,7 +33,7 @@
 ------------------------------------------------------------------------------
 
 The code in Data.Array.Vector.Algorithms.Mutable.Optimal is adapted from a C
-algorithm for the same purpose. The folowing is the copyright notice for said
+algorithm for the same purpose. The following is the copyright notice for said
 C code:
 
 Copyright (c) 2004 Paul Hsieh
diff --git a/bench/Blocks.hs b/bench/Blocks.hs
deleted file mode 100644
--- a/bench/Blocks.hs
+++ /dev/null
@@ -1,62 +0,0 @@
-{-# LANGUAGE Rank2Types #-}
-
-module Blocks where
-
-import Control.Monad
-import Control.Monad.ST
-
-import Data.Vector.Unboxed.Mutable
-
-import System.CPUTime
-
-import System.Random.Mersenne
-
--- Some conveniences for doing evil stuff in the ST monad.
--- All the tests get run in IO, but uvector stuff happens
--- in ST, so we temporarily coerce.
-clock :: IO Integer
-clock = getCPUTime
-
--- Strategies for filling the initial arrays
-rand :: (MTRandom e) => MTGen -> Int -> IO e
-rand g _ = random g
-
-ascend :: Num e => Int -> IO e
-ascend = return . fromIntegral
-
-descend :: Num e => e -> Int -> IO e
-descend m n = return $ m - fromIntegral n
-
-modulo :: Integral e => e -> Int -> IO e
-modulo m n = return $ fromIntegral n `mod` m
-
--- This is the worst case for the median-of-three quicksort
--- used in the introsort implementation.
-medianKiller :: Integral e => e -> Int -> IO e
-medianKiller m n'
-  | n < k     = return $ if even n then n + 1 else n + k
-  | otherwise = return $ (n - k + 1) * 2
- where
- n = fromIntegral n'
- k = m `div` 2
-{-# INLINE medianKiller #-}
-
-initialize :: (Unbox e) => MVector RealWorld e -> Int -> (Int -> IO e) -> IO ()
-initialize arr len fill = init $ len - 1
- where init n = fill n >>= unsafeWrite arr n >> when (n > 0) (init $ n - 1)
-{-# INLINE initialize #-}
-
-speedTest :: (Unbox e) => Int
-                       -> (Int -> IO e)
-                       -> (MVector RealWorld e -> IO ())
-                       -> IO Integer
-speedTest n fill algo = do
-  arr <- new n
-  initialize arr n fill
-  t0 <- clock
-  algo arr
-  t1 <- clock
-  return $ t1 - t0
-{-# INLINE speedTest #-}
-
-
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
deleted file mode 100644
--- a/bench/Main.hs
+++ /dev/null
@@ -1,195 +0,0 @@
-{-# LANGUAGE Rank2Types #-}
-
-module Main (main) where
-
-import Prelude hiding (read, length)
-import qualified Prelude as P
-
-import Control.Monad.ST
-import Control.Monad.Error
-
-import Data.Char
-import Data.Ord  (comparing)
-import Data.List (maximumBy)
-
-import Data.Vector.Unboxed.Mutable
-
-import qualified Data.Vector.Algorithms.Insertion    as INS
-import qualified Data.Vector.Algorithms.Intro        as INT
-import qualified Data.Vector.Algorithms.Heap         as H
-import qualified Data.Vector.Algorithms.Merge        as M
-import qualified Data.Vector.Algorithms.Radix        as R
-import qualified Data.Vector.Algorithms.AmericanFlag as AF
-
-import System.Environment
-import System.Console.GetOpt
-import System.Random.Mersenne
-
-import Blocks
-
--- Does nothing. For testing the speed/heap allocation of the building blocks.
-noalgo :: (Unbox e) => MVector RealWorld e -> IO ()
-noalgo _ = return ()
-
--- Allocates a temporary buffer, like mergesort for similar purposes as noalgo.
-alloc :: (Unbox e) => MVector RealWorld e -> IO ()
-alloc arr | len <= 4  = arr `seq` return ()
-          | otherwise = (new (len `div` 2) :: IO (MVector RealWorld Int)) >> return ()
- where len = length arr
-
-displayTime :: String -> Integer -> IO ()
-displayTime s elapsed = putStrLn $
-    s ++ " : " ++ show (fromIntegral elapsed / 1e12) ++ " seconds"
-
-run :: String -> IO Integer -> IO ()
-run s t = t >>= displayTime s
-
-sortSuite :: String -> MTGen -> Int -> (MVector RealWorld Int -> IO ()) -> IO ()
-sortSuite str g n sort = do
-  putStrLn $ "Testing: " ++ str
-  run "Random            " $ speedTest n (rand g >=> modulo n) sort
-  run "Sorted            " $ speedTest n ascend sort
-  run "Reverse-sorted    " $ speedTest n (descend n) sort
-  run "Random duplicates " $ speedTest n (rand g >=> modulo 1000) sort
-  let m = 4 * (n `div` 4)
-  run "Median killer     " $ speedTest m (medianKiller m) sort
-
-partialSortSuite :: String -> MTGen -> Int -> Int
-                 -> (MVector RealWorld Int -> Int -> IO ()) -> IO ()
-partialSortSuite str g n k sort = sortSuite str g n (\a -> sort a k)
-
--- -----------------
--- Argument handling
--- -----------------
-
-data Algorithm = DoNothing
-               | Allocate
-               | InsertionSort
-               | IntroSort
-               | IntroPartialSort
-               | IntroSelect
-               | HeapSort
-               | HeapPartialSort
-               | HeapSelect
-               | MergeSort
-               | RadixSort
-               | AmericanFlagSort
-               deriving (Show, Read, Enum, Bounded)
-
-data Options = O { algos :: [Algorithm], elems :: Int, portion :: Int, usage :: Bool } deriving (Show)
-
-defaultOptions :: Options
-defaultOptions = O [] 10000 1000 False
-
-type OptionsT = Options -> Either String Options
-
-options :: [OptDescr OptionsT]
-options = [ Option ['A']     ["algorithm"] (ReqArg parseAlgo "ALGO")
-               ("Specify an algorithm to be run. Options:\n" ++ algoOpts)
-          , Option ['n']     ["num-elems"] (ReqArg parseN    "INT")
-               "Specify the size of arrays in algorithms."
-          , Option ['k']     ["portion"]   (ReqArg parseK    "INT")
-               "Specify the number of elements to partial sort/select in\nrelevant algorithms."
-          , Option ['?','v'] ["help"]      (NoArg $ \o -> Right $ o { usage = True })
-               "Show options."
-          ]
- where
- allAlgos :: [Algorithm]
- allAlgos = [minBound .. maxBound]
- algoOpts = fmt allAlgos
- fmt (x:y:zs) = '\t' : pad (show x) ++ show y ++ "\n" ++ fmt zs
- fmt [x]      = '\t' : show x ++ "\n"
- fmt []       = ""
- size         = ("    " ++) . maximumBy (comparing P.length) . map show $ allAlgos
- pad str      = zipWith const (str ++ repeat ' ') size
-
-parseAlgo :: String -> Options -> Either String Options
-parseAlgo "None" o = Right $ o { algos = [] }
-parseAlgo "All"  o = Right $ o { algos = [DoNothing .. AmericanFlagSort] }
-parseAlgo s      o = leftMap (\e -> "Unrecognized algorithm `" ++ e ++ "'")
-                     . fmap (\v -> o { algos = v : algos o }) $ readEither s
-
-leftMap :: (a -> b) -> Either a c -> Either b c
-leftMap f (Left a)  = Left (f a)
-leftMap _ (Right c) = Right c
-
-parseNum :: (Int -> Options) -> String -> Either String Options
-parseNum f = leftMap (\e -> "Invalid numeric argument `" ++ e ++ "'") . fmap f . readEither
-
-parseN, parseK :: String -> Options -> Either String Options
-parseN s o = parseNum (\n -> o { elems   = n }) s
-parseK s o = parseNum (\k -> o { portion = k }) s
-
-readEither :: Read a => String -> Either String a
-readEither s = case reads s of
-  [(x,t)] | all isSpace t -> Right x
-  _                       -> Left s
-
-runTest :: MTGen -> 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
-  InsertionSort      -> sortSuite        "insertion sort"        g n   insertionSort
-  IntroSort          -> sortSuite        "introsort"             g n   introSort
-  IntroPartialSort   -> partialSortSuite "partial introsort"     g n k introPSort
-  IntroSelect        -> partialSortSuite "introselect"           g n k introSelect
-  HeapSort           -> sortSuite        "heap sort"             g n   heapSort
-  HeapPartialSort    -> partialSortSuite "partial heap sort"     g n k heapPSort
-  HeapSelect         -> partialSortSuite "heap select"           g n k heapSelect
-  MergeSort          -> sortSuite        "merge sort"            g n   mergeSort
-  RadixSort          -> sortSuite        "radix sort"            g n   radixSort
-  AmericanFlagSort   -> sortSuite        "flag sort"             g n   flagSort
-  _                  -> putStrLn $ "Currently unsupported algorithm: " ++ show alg
-
-mergeSort :: MVector RealWorld Int -> IO ()
-mergeSort v = M.sort v
-{-# NOINLINE mergeSort #-}
-
-introSort :: MVector RealWorld Int -> IO ()
-introSort v = INT.sort v
-{-# NOINLINE introSort #-}
-
-introPSort :: MVector RealWorld Int -> Int -> IO ()
-introPSort v k = INT.partialSort v k
-{-# NOINLINE introPSort #-}
-
-introSelect :: MVector RealWorld Int -> Int -> IO ()
-introSelect v k = INT.select v k
-{-# NOINLINE introSelect #-}
-
-heapSort :: MVector RealWorld Int -> IO ()
-heapSort v = H.sort v
-{-# NOINLINE heapSort #-}
-
-heapPSort :: MVector RealWorld Int -> Int -> IO ()
-heapPSort v k = H.partialSort v k
-{-# NOINLINE heapPSort #-}
-
-heapSelect :: MVector RealWorld Int -> Int -> IO ()
-heapSelect v k = H.select v k
-{-# NOINLINE heapSelect #-}
-
-insertionSort :: MVector RealWorld Int -> IO ()
-insertionSort v = INS.sort v
-{-# NOINLINE insertionSort #-}
-
-radixSort :: MVector RealWorld Int -> IO ()
-radixSort v = R.sort v
-{-# NOINLINE radixSort #-}
-
-flagSort :: MVector RealWorld Int -> IO ()
-flagSort v = AF.sort v
-{-# 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
-
-
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/simple/Blocks.hs b/bench/simple/Blocks.hs
new file mode 100644
--- /dev/null
+++ b/bench/simple/Blocks.hs
@@ -0,0 +1,62 @@
+{-# LANGUAGE Rank2Types #-}
+
+module Blocks where
+
+import Control.Monad
+import Control.Monad.ST
+
+import Data.Vector.Unboxed.Mutable
+
+import System.CPUTime
+
+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
+-- in ST, so we temporarily coerce.
+clock :: IO Integer
+clock = getCPUTime
+
+-- Strategies for filling the initial arrays
+rand :: Variate e => GenIO -> Int -> IO e
+rand g _ = uniform g
+
+ascend :: Num e => Int -> IO e
+ascend = return . fromIntegral
+
+descend :: Num e => e -> Int -> IO e
+descend m n = return $ m - fromIntegral n
+
+modulo :: Integral e => e -> Int -> IO e
+modulo m n = return $ fromIntegral n `mod` m
+
+-- This is the worst case for the median-of-three quicksort
+-- used in the introsort implementation.
+medianKiller :: Integral e => e -> Int -> IO e
+medianKiller m n'
+  | n < k     = return $ if even n then n + 1 else n + k
+  | otherwise = return $ (n - k + 1) * 2
+ where
+ n = fromIntegral n'
+ k = m `div` 2
+{-# INLINE medianKiller #-}
+
+initialize :: (Unbox e) => MVector RealWorld e -> Int -> (Int -> IO e) -> IO ()
+initialize arr len fill = initial $ len - 1
+ where initial n = fill n >>= unsafeWrite arr n >> when (n > 0) (initial $ n - 1)
+{-# INLINE initialize #-}
+
+speedTest :: (Unbox e) => MVector RealWorld e
+                       -> Int
+                       -> (Int -> IO e)
+                       -> (MVector RealWorld e -> IO ())
+                       -> IO Integer
+speedTest arr n fill algo = do
+  initialize arr n fill
+  t0 <- clock
+  algo arr
+  t1 <- clock
+  return $ t1 - t0
+{-# INLINE speedTest #-}
+
+
diff --git a/bench/simple/Main.hs b/bench/simple/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/simple/Main.hs
@@ -0,0 +1,202 @@
+{-# LANGUAGE Rank2Types #-}
+
+module Main (main) where
+
+import Prelude hiding (read, length)
+import qualified Prelude as P
+
+import Control.Monad
+import Control.Monad.ST
+
+import Data.Char
+import Data.Ord  (comparing)
+import Data.List (maximumBy)
+
+import qualified Data.Vector.Unboxed.Mutable as UVector
+import Data.Vector.Unboxed.Mutable (MVector, Unbox)
+
+import qualified Data.Vector.Algorithms.Insertion    as INS
+import qualified Data.Vector.Algorithms.Intro        as INT
+import qualified Data.Vector.Algorithms.Heap         as H
+import qualified Data.Vector.Algorithms.Merge        as M
+import qualified Data.Vector.Algorithms.Radix        as R
+import qualified Data.Vector.Algorithms.AmericanFlag as AF
+import qualified Data.Vector.Algorithms.Tim          as T
+
+import System.Environment
+import System.Console.GetOpt
+import System.Random.MWC
+
+import Blocks
+
+-- Does nothing. For testing the speed/heap allocation of the building blocks.
+noalgo :: (Unbox e) => MVector RealWorld e -> IO ()
+noalgo _ = return ()
+
+-- Allocates a temporary buffer, like mergesort for similar purposes as noalgo.
+alloc :: (Unbox e) => MVector RealWorld e -> IO ()
+alloc arr | len <= 4  = arr `seq` return ()
+          | otherwise = (UVector.new (len `div` 2) :: IO (MVector RealWorld Int)) >> return ()
+ where len = UVector.length arr
+
+displayTime :: String -> Integer -> IO ()
+displayTime s elapsed = putStrLn $
+    s ++ " : " ++ show (fromIntegral elapsed / (1e12 :: Double)) ++ " seconds"
+
+run :: String -> IO Integer -> IO ()
+run s t = t >>= displayTime s
+
+sortSuite :: String -> GenIO -> Int -> (MVector RealWorld Int -> IO ()) -> IO ()
+sortSuite str g n sort = do
+  arr <- UVector.new n
+  putStrLn $ "Testing: " ++ str
+  run "Random            " $ speedTest arr n (rand g >=> modulo n) sort
+  run "Sorted            " $ speedTest arr n ascend sort
+  run "Reverse-sorted    " $ speedTest arr n (descend n) sort
+  run "Random duplicates " $ speedTest arr n (rand g >=> modulo 1000) sort
+  let m = 4 * (n `div` 4)
+  run "Median killer     " $ speedTest arr m (medianKiller m) sort
+
+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)
+
+-- -----------------
+-- Argument handling
+-- -----------------
+
+data Algorithm = DoNothing
+               | Allocate
+               | InsertionSort
+               | IntroSort
+               | IntroPartialSort
+               | IntroSelect
+               | HeapSort
+               | HeapPartialSort
+               | HeapSelect
+               | MergeSort
+               | RadixSort
+               | AmericanFlagSort
+               | TimSort
+               deriving (Show, Read, Enum, Bounded)
+
+data Options = O { algos :: [Algorithm], elems :: Int, portion :: Int, usage :: Bool } deriving (Show)
+
+defaultOptions :: Options
+defaultOptions = O [] 10000 1000 False
+
+type OptionsT = Options -> Either String Options
+
+options :: [OptDescr OptionsT]
+options = [ Option ['A']     ["algorithm"] (ReqArg parseAlgo "ALGO")
+               ("Specify an algorithm to be run. Options:\n" ++ algoOpts)
+          , Option ['n']     ["num-elems"] (ReqArg parseN    "INT")
+               "Specify the size of arrays in algorithms."
+          , Option ['k']     ["portion"]   (ReqArg parseK    "INT")
+               "Specify the number of elements to partial sort/select in\nrelevant algorithms."
+          , Option ['?','v'] ["help"]      (NoArg $ \o -> Right $ o { usage = True })
+               "Show options."
+          ]
+ where
+ allAlgos :: [Algorithm]
+ allAlgos = [minBound .. maxBound]
+ algoOpts = fmt allAlgos
+ fmt (x:y:zs) = '\t' : pad (show x) ++ show y ++ "\n" ++ fmt zs
+ fmt [x]      = '\t' : show x ++ "\n"
+ fmt []       = ""
+ size         = ("    " ++) . maximumBy (comparing P.length) . map show $ allAlgos
+ pad str      = zipWith const (str ++ repeat ' ') size
+
+parseAlgo :: String -> Options -> Either String Options
+parseAlgo "None" o = Right $ o { algos = [] }
+parseAlgo "All"  o = Right $ o { algos = [DoNothing .. AmericanFlagSort] }
+parseAlgo s      o = leftMap (\e -> "Unrecognized algorithm `" ++ e ++ "'")
+                     . fmap (\v -> o { algos = v : algos o }) $ readEither s
+
+leftMap :: (a -> b) -> Either a c -> Either b c
+leftMap f (Left a)  = Left (f a)
+leftMap _ (Right c) = Right c
+
+parseNum :: (Int -> Options) -> String -> Either String Options
+parseNum f = leftMap (\e -> "Invalid numeric argument `" ++ e ++ "'") . fmap f . readEither
+
+parseN, parseK :: String -> Options -> Either String Options
+parseN s o = parseNum (\n -> o { elems   = n }) s
+parseK s o = parseNum (\k -> o { portion = k }) s
+
+readEither :: Read a => String -> Either String a
+readEither s = case reads s of
+  [(x,t)] | all isSpace t -> Right x
+  _                       -> Left s
+
+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
+  InsertionSort      -> sortSuite        "insertion sort"        g n   insertionSort
+  IntroSort          -> sortSuite        "introsort"             g n   introSort
+  IntroPartialSort   -> partialSortSuite "partial introsort"     g n k introPSort
+  IntroSelect        -> partialSortSuite "introselect"           g n k introSelect
+  HeapSort           -> sortSuite        "heap sort"             g n   heapSort
+  HeapPartialSort    -> partialSortSuite "partial heap sort"     g n k heapPSort
+  HeapSelect         -> partialSortSuite "heap select"           g n k heapSelect
+  MergeSort          -> sortSuite        "merge sort"            g n   mergeSort
+  RadixSort          -> sortSuite        "radix sort"            g n   radixSort
+  AmericanFlagSort   -> sortSuite        "flag sort"             g n   flagSort
+  TimSort            -> sortSuite        "tim sort"              g n   timSort
+
+mergeSort :: MVector RealWorld Int -> IO ()
+mergeSort v = M.sort v
+{-# NOINLINE mergeSort #-}
+
+introSort :: MVector RealWorld Int -> IO ()
+introSort v = INT.sort v
+{-# NOINLINE introSort #-}
+
+introPSort :: MVector RealWorld Int -> Int -> IO ()
+introPSort v k = INT.partialSort v k
+{-# NOINLINE introPSort #-}
+
+introSelect :: MVector RealWorld Int -> Int -> IO ()
+introSelect v k = INT.select v k
+{-# NOINLINE introSelect #-}
+
+heapSort :: MVector RealWorld Int -> IO ()
+heapSort v = H.sort v
+{-# NOINLINE heapSort #-}
+
+heapPSort :: MVector RealWorld Int -> Int -> IO ()
+heapPSort v k = H.partialSort v k
+{-# NOINLINE heapPSort #-}
+
+heapSelect :: MVector RealWorld Int -> Int -> IO ()
+heapSelect v k = H.select v k
+{-# NOINLINE heapSelect #-}
+
+insertionSort :: MVector RealWorld Int -> IO ()
+insertionSort v = INS.sort v
+{-# NOINLINE insertionSort #-}
+
+radixSort :: MVector RealWorld Int -> IO ()
+radixSort v = R.sort v
+{-# NOINLINE radixSort #-}
+
+flagSort :: MVector RealWorld Int -> IO ()
+flagSort v = AF.sort v
+{-# NOINLINE flagSort #-}
+
+timSort :: MVector RealWorld Int -> IO ()
+timSort v = T.sort v
+{-# NOINLINE timSort #-}
+
+main :: IO ()
+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 "vector-algorithms-bench" options
+    (_, _, errs) -> putStrLn $ usageInfo (concat errs) options
+
+
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.hs b/src/Data/Vector/Algorithms.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms.hs
@@ -0,0 +1,77 @@
+{-# language BangPatterns, RankNTypes, ScopedTypeVariables #-}
+module Data.Vector.Algorithms where
+
+import Prelude hiding (length)
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.ST (runST)
+
+import Data.Vector.Generic.Mutable
+import qualified Data.Vector.Generic as V
+import qualified Data.Vector.Unboxed.Mutable as UMV
+import qualified Data.Bit as Bit
+
+import Data.Vector.Algorithms.Common (Comparison)
+import Data.Vector.Algorithms.Intro (sortUniqBy)
+import qualified Data.Vector.Algorithms.Search  as S
+
+-- | The `nub` function which removes duplicate elements from a vector.
+nub :: forall v e . (V.Vector v e, Ord e) => v e -> v e
+nub = nubBy compare
+{-# INLINE nub #-}
+
+-- | A version of `nub` with a custom comparison predicate.
+--
+-- /Note:/ This function makes use of `sortByUniq` using the intro
+-- sort algorithm.
+nubBy ::
+  forall v e . (V.Vector v e) =>
+  Comparison e -> v e -> v e
+nubBy cmp vec = runST $ do
+  mv <- V.unsafeThaw vec -- safe as the nubByMut algorithm copies the input
+  destMV <- nubByMut sortUniqBy cmp mv
+  v <- V.unsafeFreeze destMV
+  pure (V.force v)
+{-# INLINE nubBy #-}
+
+-- | The `nubByMut` function takes in an in-place sort algorithm
+-- and uses it to do a de-deduplicated sort. It then uses this to
+-- remove duplicate elements from the input.
+--
+-- /Note:/ Since this algorithm needs the original input and so
+-- copies before sorting in-place. As such, it is safe to use on
+-- immutable inputs.
+nubByMut ::
+  forall m v e . (PrimMonad m, MVector v e) =>
+  (Comparison e -> v (PrimState m) e -> m (v (PrimState m) e))
+  -> Comparison e -> v (PrimState m) e -> m (v (PrimState m) e)
+nubByMut alg cmp inp = do
+  let len = length inp
+  inp' <- clone inp
+  sortUniqs <- alg cmp inp'
+  let uniqLen = length sortUniqs
+  bitmask <- UMV.replicate uniqLen (Bit.Bit False) -- bitmask to track which elements have
+                                                   -- already been seen.
+  dest ::  v (PrimState m) e <- unsafeNew uniqLen  -- return vector
+  let
+    go :: Int -> Int -> m ()
+    go !srcInd !destInd
+      | srcInd == len = pure ()
+      | destInd == uniqLen = pure ()
+      | otherwise = do
+          curr    <- unsafeRead inp srcInd                -- read current element
+          sortInd <- S.binarySearchBy cmp sortUniqs curr  -- find sorted index
+          bit <- UMV.unsafeRead bitmask sortInd           -- check if we have already seen
+                                                          -- this element in bitvector
+          case bit of
+            -- if we have seen it then iterate
+            Bit.Bit True -> go (srcInd + 1) destInd
+            -- if we haven't then write it into output
+            -- and mark that it has been seen
+            Bit.Bit False -> do
+              UMV.unsafeWrite bitmask sortInd (Bit.Bit True)
+              unsafeWrite dest destInd curr
+              go (srcInd + 1) (destInd + 1)
+  go 0 0
+  pure dest
+{-# INLINABLE nubByMut #-}
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,402 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# 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
+                                           , sortUniq
+                                           , sortBy
+                                           , sortUniqBy
+                                           , terminate
+                                           , Lexicographic(..)
+                                           ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import Data.Proxy
+
+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
+
+import Foreign.Storable
+
+-- | 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
+  -- | Computes the length of a representative of a stripe. It should take 'n'
+  -- passes to sort values of extent 'n'. The extent may not be uniform across
+  -- all values of the type.
+  extent    :: e -> Int
+
+  -- | The size of the bucket array necessary for sorting es
+  size      :: Proxy e -> Int
+  -- | Determines which bucket a given element should inhabit for a
+  -- particular iteration.
+  index     :: Int -> e -> Int
+
+instance Lexicographic Word8 where
+  extent _ = 1
+  {-# INLINE extent #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index _ n = fromIntegral n
+  {-# INLINE index #-}
+
+instance Lexicographic Word16 where
+  extent _ = 2
+  {-# INLINE extent #-}
+  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
+  extent _ = 4
+  {-# INLINE extent #-}
+  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
+  extent _ = 8
+  {-# INLINE extent #-}
+  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
+  extent _ = sizeOf (0 :: Word)
+  {-# INLINE extent #-}
+  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
+  extent _ = 1
+  {-# INLINE extent #-}
+  size _ = 256
+  {-# INLINE size #-}
+  index _ n = 255 .&. fromIntegral n `xor` 128
+  {-# INLINE index #-}
+
+instance Lexicographic Int16 where
+  extent _ = 2
+  {-# INLINE extent #-}
+  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
+  extent _ = 4
+  {-# INLINE extent #-}
+  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
+  extent _ = 8
+  {-# INLINE extent #-}
+  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
+  extent _ = sizeOf (0 :: Int)
+  {-# INLINE extent #-}
+  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
+  extent = B.length
+  {-# INLINE extent #-}
+  size _ = 257
+  {-# INLINE size #-}
+  index i b
+    | i >= B.length b = 0
+    | otherwise       = fromIntegral (B.index b i) + 1
+  {-# INLINE index #-}
+
+instance (Lexicographic a, Lexicographic b) => Lexicographic (a, b) where
+  extent (a,b) = extent a + extent b
+  {-# INLINE extent #-}
+  size _ = size (Proxy :: Proxy a) `max` size (Proxy :: Proxy b)
+  {-# INLINE size #-}
+  index i (a,b)
+    | i >= extent a = index i b
+    | otherwise     = index i a
+  {-# INLINE index #-}
+
+instance (Lexicographic a, Lexicographic b) => Lexicographic (Either a b) where
+  extent (Left  a) = 1 + extent a
+  extent (Right b) = 1 + extent b
+  {-# INLINE extent #-}
+  size _ = size (Proxy :: Proxy a) `max` size (Proxy :: Proxy b)
+  {-# INLINE size #-}
+  index 0 (Left  _) = 0
+  index 0 (Right _) = 1
+  index n (Left  a) = index (n-1) a
+  index n (Right b) = index (n-1) b
+  {-# INLINE index #-}
+
+-- | Given a representative of a stripe and an index number, this
+-- function determines whether to stop sorting.
+terminate :: Lexicographic e => e -> Int -> Bool
+terminate e i = i >= extent e
+{-# INLINE terminate #-}
+
+-- | 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 p) index v
+ where p :: Proxy e
+       p = Proxy
+{-# INLINE sort #-}
+
+-- | A variant on `sort` that returns a vector of unique elements.
+sortUniq :: forall e m v. (PrimMonad m, MVector v e, Lexicographic e, Ord e)
+     => v (PrimState m) e -> m (v (PrimState m) e)
+sortUniq v = sortUniqBy compare terminate (size p) index v
+ where p :: Proxy e
+       p = Proxy
+{-# INLINE sortUniq #-}
+
+-- | 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 #-}
+
+-- | A variant on `sortBy` which returns a vector of unique elements.
+sortUniqBy :: (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 (v (PrimState m) e)
+sortUniqBy cmp stop buckets radix v
+  | length v == 0 = return v
+  | otherwise     = do count <- new buckets
+                       pile <- new buckets
+                       countLoop (radix 0) v count
+                       flagLoop cmp stop radix count pile v
+                       uniqueMutableBy cmp v
+{-# INLINE sortUniqBy #-}
+
+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 already 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,132 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- 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
+  ( type Comparison
+  , copyOffset
+  , inc
+  , countLoop
+  , midPoint
+  , uniqueMutableBy
+  )
+  where
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Vector.Generic.Mutable
+import Data.Word (Word)
+
+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 #-}
+
+midPoint :: Int -> Int -> Int
+midPoint a b =
+  toInt $ (toWord a + toWord b) `div` 2
+  where
+    toWord :: Int -> Word
+    toWord = fromIntegral
+
+    toInt :: Word -> Int
+    toInt = fromIntegral
+{-# INLINE midPoint #-}
+
+-- Adapted from Andrew Martin's uniquqMutable in the primitive-sort package
+uniqueMutableBy :: forall m v a . (PrimMonad m, MVector v a)
+  => Comparison a -> v (PrimState m) a -> m (v (PrimState m) a)
+uniqueMutableBy cmp mv = do
+  let !len = basicLength mv
+  if len > 1
+    then do
+      !a0 <- unsafeRead mv 0
+      let findFirstDuplicate :: a -> Int -> m Int
+          findFirstDuplicate !prev !ix = if ix < len
+            then do
+              a <- unsafeRead mv ix
+              if cmp a prev == EQ
+                then return ix
+                else findFirstDuplicate a (ix + 1)
+            else return ix
+      dupIx <- findFirstDuplicate a0 1
+      if dupIx == len
+        then return mv
+        else do
+          let deduplicate :: a -> Int -> Int -> m Int
+              deduplicate !prev !srcIx !dstIx = if srcIx < len
+                then do
+                  a <- unsafeRead mv srcIx
+                  if cmp a prev == EQ
+                    then deduplicate a (srcIx + 1) dstIx
+                    else do
+                      unsafeWrite mv dstIx a
+                      deduplicate a (srcIx + 1) (dstIx + 1)
+                else return dstIx
+          !a <- unsafeRead mv dupIx
+          !reducedLen <- deduplicate a (dupIx + 1) dupIx
+          resizeVector mv reducedLen
+    else return mv
+{-# INLINABLE uniqueMutableBy #-}
+
+-- Used internally in uniqueMutableBy: copies the elements of a vector to one
+-- of a smaller size.
+resizeVector
+  :: (MVector v a, PrimMonad m)
+  =>  v (PrimState m) a -> Int -> m (v (PrimState m) a)
+resizeVector !src !sz = do
+  dst <- unsafeNew sz
+  copyToSmaller dst src
+  pure dst
+{-# inline resizeVector #-}
+
+-- Used internally in resizeVector: copy a vector from a larger to
+-- smaller vector. Should not be used if the source vector
+-- is smaller than the target vector.
+copyToSmaller
+  :: (MVector v a, PrimMonad m)
+  => v (PrimState m) a -> v (PrimState m) a -> m ()
+copyToSmaller !dst !src = stToPrim $ do_copy 0
+    where
+      !n = basicLength dst
+
+      do_copy i | i < n = do
+                            x <- basicUnsafeRead src i
+                            basicUnsafeWrite dst i x
+                            do_copy (i+1)
+                | otherwise = return ()
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,348 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Heap
+-- Copyright   : (c) 2008-2015 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
+       , sortUniq
+       , sortBy
+       , sortUniqBy
+       , sortByBounds
+         -- * Selection
+       , select
+       , selectBy
+       , selectByBounds
+         -- * Partial sorts
+       , partialSort
+       , partialSortBy
+       , partialSortByBounds
+         -- * Heap operations
+       , heapify
+       , pop
+       , popTo
+       , sortHeap
+       , heapInsert
+       , 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, uniqueMutableBy)
+
+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
+{-# INLINE sort #-}
+
+-- | A variant on `sort` that returns a vector of unique elements.
+sortUniq :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m (v (PrimState m) e)
+sortUniq = sortUniqBy compare
+{-# INLINE sortUniq #-}
+
+-- | 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 #-}
+
+-- | A variant on `sortBy` which returns a vector of unique elements.
+sortUniqBy :: (PrimMonad m, MVector v e)
+  => Comparison e -> v (PrimState m) e -> m (v (PrimState m) e)
+sortUniqBy cmp a = do
+  sortByBounds cmp a 0 (length a)
+  uniqueMutableBy cmp a
+{-# INLINE sortUniqBy #-}
+
+-- | 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 -- ^ lower index, l
+  -> Int -- ^ upper index, u
+  -> 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 -- ^ number of elements to select, k
+  -> 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 -- ^ number of elements to select, k
+  -> 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 -- ^ number of elements to select, k
+  -> Int -- ^ lower index, l
+  -> Int -- ^ upper index, u
+  -> 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.
+--
+-- The remaining values of the array will be in no particular order.
+partialSort
+  :: (PrimMonad m, MVector v e, Ord e)
+  => v (PrimState m) e
+  -> Int -- ^ number of elements to sort, k
+  -> m ()
+partialSort = partialSortBy compare
+{-# INLINE partialSort #-}
+
+-- | Moves the lowest k elements (as defined by the comparison) to
+-- the front of the array, sorted.
+--
+-- The remaining values of the array will be in no particular order.
+partialSortBy
+  :: (PrimMonad m, MVector v e)
+  => Comparison e
+  -> v (PrimState m) e
+  -> Int -- ^ number of elements to sort, k
+  -> 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.
+--
+-- The remaining values in [l,u) will be in no particular order. Values outside
+-- the range [l,u) will be unaffected.
+partialSortByBounds
+  :: (PrimMonad m, MVector v e)
+  => Comparison e
+  -> v (PrimState m) e
+  -> Int -- ^ number of elements to sort, k
+  -> Int -- ^ lower index, l
+  -> Int -- ^ upper index, u
+  -> 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 + 1) l u
+                    sortHeap cmp a l (l + 4) (l + k + 1)
+                    O.sort4ByOffset cmp a l
+ where
+ len = u - l
+{-# INLINE partialSortByBounds #-}
+
+-- | Constructs a heap in a portion of an array [l, u), using the values therein.
+--
+-- Note: 'heapify' is more efficient than constructing a heap by repeated
+-- insertion. Repeated insertion has complexity O(n*log n) while 'heapify' is able
+-- to construct a heap in O(n), where n is the number of elements in the heap.
+heapify
+  :: (PrimMonad m, MVector v e)
+  => Comparison e
+  -> v (PrimState m) e
+  -> Int -- ^ lower index, l
+  -> Int -- ^ upper index, u
+  -> 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 -- ^ lower heap index, l
+  -> Int -- ^ upper heap index, u
+  -> 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 -- ^ lower heap index, l
+  -> Int -- ^ upper heap index, u
+  -> Int -- ^ index to pop to, t
+  -> 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 -- ^ lower heap index, l
+  -> Int -- ^ lower bound of final sorted portion, m
+  -> Int -- ^ upper heap index, u
+  -> 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 #-}
+
+-- | Given a heap stored in a portion of an array [l,u) and an element e,
+-- inserts the element into the heap, resulting in a heap in [l,u].
+--
+-- Note: it is best to only use this operation when incremental construction of
+-- a heap is required. 'heapify' is capable of building a heap in O(n) time,
+-- while repeated insertion takes O(n*log n) time.
+heapInsert
+  :: (PrimMonad m, MVector v e)
+  => Comparison e
+  -> v (PrimState m) e
+  -> Int -- ^ lower heap index, l
+  -> Int -- ^ upper heap index, u
+  -> e -- ^ element to be inserted, e
+  -> m ()
+heapInsert cmp v l u e = sift (u - l)
+ where
+ sift k
+   | k <= 0    = unsafeWrite v l e
+   | otherwise = let pi = shiftR (k-1) 2
+                  in unsafeRead v (l + pi) >>= \p -> case cmp p e of
+                       LT -> unsafeWrite v (l + k) p >> sift pi
+                       _  -> unsafeWrite v (l + k) e
+{-# INLINE heapInsert #-}
+
+-- 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,96 @@
+{-# LANGUAGE TypeFamilies #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- 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
+       , sortUniq
+       , sortBy
+       , sortUniqBy
+       , sortByBounds
+       , sortByBounds'
+       , Comparison
+       ) where
+
+
+import Prelude hiding (read, length)
+
+import Control.Monad.Primitive
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison, uniqueMutableBy)
+
+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
+{-# INLINE sort #-}
+
+-- | A variant on `sort` that returns a vector of unique elements.
+sortUniq :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m (v (PrimState m) e)
+sortUniq = sortUniqBy compare
+{-# INLINE sortUniq #-}
+
+-- | 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 #-}
+
+-- | A variant on `sortBy` which returns a vector of unique elements.
+sortUniqBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m (v (PrimState m) e)
+sortUniqBy cmp a = do
+  sortByBounds cmp a 0 (length a)
+  uniqueMutableBy cmp a
+{-# INLINE sortUniqBy #-}
+
+-- | 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,263 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Intro
+-- Copyright   : (c) 2008-2015 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
+       , sortUniq
+       , sortBy
+       , sortUniqBy
+       , 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, midPoint, uniqueMutableBy)
+
+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
+{-# INLINE sort #-}
+
+-- | A variant on `sort` that returns a vector of unique elements.
+sortUniq :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m (v (PrimState m) e)
+sortUniq = sortUniqBy compare
+{-# INLINE sortUniq #-}
+
+-- | A variant on `sortBy` which returns a vector of unique elements.
+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 an entire array using a custom ordering returning a vector of
+-- the unique elements.
+sortUniqBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m (v (PrimState m) e)
+sortUniqBy cmp a = do
+  sortByBounds cmp a 0 (length a)
+  uniqueMutableBy cmp a
+{-# INLINE sortUniqBy #-}
+
+-- | 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 -- ^ lower index, l
+  -> Int -- ^ upper index, u
+  -> 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   = midPoint u l
+{-# 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 -- ^ number of elements to select, k
+  -> 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 -- ^ number of elements to select, k
+  -> 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 -- ^ number of elements to select, k
+  -> Int -- ^ lower bound, l
+  -> Int -- ^ upper bound, u
+  -> 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 = midPoint u l
+{-# 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 -- ^ number of elements to sort, k
+  -> 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 -- ^ number of elements to sort, k
+  -> 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 -- ^ number of elements to sort, k
+  -> Int -- ^ lower index, l
+  -> Int -- ^ upper index, u
+  -> m ()
+partialSortByBounds cmp a k l u
+  | l >= u    = return ()
+  | otherwise = let k' = min (u-l) k
+                      -- N.B. Clamp k to the length of the range
+                      -- being sorted.
+                in 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 = midPoint u l
+{-# 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
+ 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,119 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- 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
+       , sortUniq
+       , sortBy
+       , sortUniqBy
+       , 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, midPoint, uniqueMutableBy)
+
+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
+{-# INLINE sort #-}
+
+-- | A variant on `sort` that returns a vector of unique elements.
+sortUniq :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m (v (PrimState m) e)
+sortUniq = sortUniqBy compare
+{-# INLINE sortUniq #-}
+
+-- | Sorts an array using a custom comparison.
+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()
+sortBy cmp vec = if len <= 4
+                    then if len <= 2
+                            then if len /= 2
+                                    then return ()
+                                    else O.sort2ByOffset cmp vec 0
+                            else if len == 3
+                                    then O.sort3ByOffset cmp vec 0
+                                    else O.sort4ByOffset cmp vec 0
+                    else if len < threshold
+                            then I.sortByBounds cmp vec 0 len
+                            else do buf <- new halfLen
+                                    mergeSortWithBuf cmp vec buf
+ where
+ len     = length vec
+ -- odd lengths have a larger half that needs to fit, so use ceiling, not floor
+ halfLen = (len + 1) `div` 2
+{-# INLINE sortBy #-}
+
+-- | A variant on `sortBy` which returns a vector of unique elements.
+sortUniqBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m (v (PrimState m) e)
+sortUniqBy cmp vec = do
+  sortBy cmp vec
+  uniqueMutableBy cmp vec
+{-# INLINE sortUniqBy #-}
+
+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 = midPoint u l
+{-# 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,252 @@
+{-# 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)
+
+#if MIN_VERSION_vector(0,13,0)
+import qualified Data.Vector.Internal.Check as Ck
+# define CHECK_INDEX(name, i, n) Ck.checkIndex Ck.Unsafe (i) (n)
+#else
+# define CHECK_INDEX(name, i, n) UNSAFE_CHECK(checkIndex) name (i) (n)
+#endif
+
+#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 = CHECK_INDEX("sort2ByIndex", i, length a)
+                       $ CHECK_INDEX("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 = CHECK_INDEX("sort3ByIndex", i, length a)
+                         $ CHECK_INDEX("sort3ByIndex", j, length a)
+                         $ CHECK_INDEX("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 = CHECK_INDEX("sort4ByIndex", i, length a)
+                           $ CHECK_INDEX("sort4ByIndex", j, length a)
+                           $ CHECK_INDEX("sort4ByIndex", k, length a)
+                           $ CHECK_INDEX("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,264 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- 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
+{-# INLINE 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,209 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Search
+-- Copyright   : (c) 2009-2015 Dan Doel, 2015 Tim Baumann
+-- 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
+       , gallopingSearchLeftP
+       , gallopingSearchLeftPBounds
+       , gallopingSearchRightP
+       , gallopingSearchRightPBounds
+       , 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, midPoint)
+
+-- | 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 = midPoint u l
+{-# 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 guaranteed 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 = midPoint u l
+{-# INLINE binarySearchPBounds #-}
+
+-- | Given a predicate that is guaranteed to be monotone on the vector elements
+-- in order, finds the index at which the predicate turns from False to True.
+-- The length of the vector is returned if the predicate is False for the entire
+-- vector.
+--
+-- Begins searching at the start of the vector, in increasing steps of size 2^n.
+gallopingSearchLeftP
+  :: (PrimMonad m, MVector v e) => (e -> Bool) -> v (PrimState m) e -> m Int
+gallopingSearchLeftP p vec = gallopingSearchLeftPBounds p vec 0 (length vec)
+{-# INLINE gallopingSearchLeftP #-}
+
+-- | Given a predicate that is guaranteed to be monotone on the vector elements
+-- in order, finds the index at which the predicate turns from False to True.
+-- The length of the vector is returned if the predicate is False for the entire
+-- vector.
+--
+-- Begins searching at the end of the vector, in increasing steps of size 2^n.
+gallopingSearchRightP
+  :: (PrimMonad m, MVector v e) => (e -> Bool) -> v (PrimState m) e -> m Int
+gallopingSearchRightP p vec = gallopingSearchRightPBounds p vec 0 (length vec)
+{-# INLINE gallopingSearchRightP #-}
+
+-- | 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).
+-- Begins searching at l, going right in increasing (2^n)-steps.
+gallopingSearchLeftPBounds :: (PrimMonad m, MVector v e)
+                           => (e -> Bool)
+                           -> v (PrimState m) e
+                           -> Int -- ^ l
+                           -> Int -- ^ u
+                           -> m Int
+gallopingSearchLeftPBounds p vec l u
+  | u <= l    = return l
+  | otherwise = do x <- unsafeRead vec l
+                   if p x then return l else iter (l+1) l 2
+ where
+ binSearch = binarySearchPBounds p vec
+ iter !i !j !_stepSize | i >= u - 1 = do
+   x <- unsafeRead vec (u-1)
+   if p x then binSearch (j+1) (u-1) else return u
+ iter !i !j !stepSize = do
+   x <- unsafeRead vec i
+   if p x then binSearch (j+1) i else iter (i+stepSize) i (2*stepSize)
+{-# INLINE gallopingSearchLeftPBounds #-}
+
+-- | 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).
+-- Begins searching at u, going left in increasing (2^n)-steps.
+gallopingSearchRightPBounds :: (PrimMonad m, MVector v e)
+                            => (e -> Bool)
+                            -> v (PrimState m) e
+                            -> Int -- ^ l
+                            -> Int -- ^ u
+                            -> m Int
+gallopingSearchRightPBounds p vec l u
+  | u <= l    = return l
+  | otherwise = iter (u-1) (u-1) (-1)
+ where
+ binSearch = binarySearchPBounds p vec
+ iter !i !j !_stepSize | i <= l = do
+   x <- unsafeRead vec l
+   if p x then return l else binSearch (l+1) j
+ iter !i !j !stepSize = do
+   x <- unsafeRead vec i
+   if p x then iter (i+stepSize) i (2*stepSize) else binSearch (i+1) j
+{-# INLINE gallopingSearchRightPBounds #-}
diff --git a/src/Data/Vector/Algorithms/Tim.hs b/src/Data/Vector/Algorithms/Tim.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Algorithms/Tim.hs
@@ -0,0 +1,382 @@
+{-# LANGUAGE BangPatterns #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Tim
+-- Copyright   : (c) 2013-2015 Dan Doel, 2015 Tim Baumann
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (bang patterns)
+--
+-- Timsort is a complex, adaptive, bottom-up merge sort. It is designed to
+-- minimize comparisons as much as possible, even at some cost in overhead.
+-- Thus, it may not be ideal for sorting simple primitive types, for which
+-- comparison is cheap. It may, however, be significantly faster for sorting
+-- arrays of complex values (strings would be an example, though an algorithm
+-- not based on comparison would probably be superior in that particular
+-- case).
+--
+-- For more information on the details of the algorithm, read on.
+--
+-- The first step of the algorithm is to identify runs of elements. These can
+-- either be non-decreasing or strictly decreasing sequences of elements in
+-- the input. Strictly decreasing sequences are used rather than
+-- non-increasing so that they can be easily reversed in place without the
+-- sort becoming unstable.
+--
+-- If the natural runs are too short, they are padded to a minimum value. The
+-- minimum is chosen based on the length of the array, and padded runs are put
+-- in order using insertion sort. The length of the minimum run size is
+-- determined as follows:
+--
+--   * If the length of the array is less than 64, the minimum size is the
+--     length of the array, and insertion sort is used for the entirety
+--
+--   * Otherwise, a value between 32 and 64 is chosen such that N/min is
+--     either equal to or just below a power of two. This avoids having a
+--     small chunk left over to merge into much larger chunks at the end.
+--
+-- This is accomplished by taking the the mininum to be the lowest six bits
+-- containing the highest set bit, and adding one if any other bits are set.
+-- For instance:
+--
+--     length: 00000000 00000000 00000000 00011011 = 25
+--     min:    00000000 00000000 00000000 00011011 = 25
+--
+--     length: 00000000 11111100 00000000 00000000 = 63 * 2^18
+--     min:    00000000 00000000 00000000 00111111 = 63
+--
+--     length: 00000000 11111100 00000000 00000001 = 63 * 2^18 + 1
+--     min:    00000000 00000000 00000000 01000000 = 64
+--
+-- Once chunks can be produced, the next step is merging them. The indices of
+-- all runs are stored in a stack. When we identify a new run, we push it onto
+-- the stack. However, certain invariants are maintained of the stack entries.
+-- Namely:
+--
+--   if stk = _ :> z :> y :> x
+--     length x + length y < length z
+--
+--   if stk = _ :> y :> x
+--     length x < length y
+--
+-- This ensures that the chunks stored are decreasing, and that the chunk
+-- sizes follow something like the fibonacci sequence, ensuring there at most
+-- log-many chunks at any time. If pushing a new chunk on the stack would
+-- violate either of the invariants, we first perform a merge.
+--
+-- If length x + length y >= length z, then y is merged with the smaller of x
+-- and z (if they are tied, x is chosen, because it is more likely to be
+-- cached). If, further,  length x >= length y then they are merged. These steps
+-- are repeated until the invariants are established.
+--
+-- The last important piece of the algorithm is the merging. At first, two
+-- chunks are merged element-wise. However, while doing so, counts are kept of
+-- the number of elements taken from one chunk without any from its partner. If
+-- this count exceeds a threshold, the merge switches to searching for elements
+-- from one chunk in the other, and copying chunks at a time. If these chunks
+-- start falling below the threshold, the merge switches back to element-wise.
+--
+-- The search used in the merge is also special. It uses a galloping strategy,
+-- where exponentially increasing indices are tested, and once two such indices
+-- are determined to bracket the desired value, binary search is used to find
+-- the exact index within that range. This is asymptotically the same as simply
+-- using binary search, but is likely to do fewer comparisons than binary search
+-- would.
+--
+-- One aspect that is not yet implemented from the original Tim sort is the
+-- adjustment of the above threshold. When galloping saves time, the threshold
+-- is lowered, and when it doesn't, it is raised. This may be implemented in the
+-- future.
+
+module Data.Vector.Algorithms.Tim
+       ( sort
+       , sortUniq
+       , sortBy
+       , sortUniqBy
+       ) where
+
+import Prelude hiding (length, reverse)
+
+import Control.Monad.Primitive
+import Control.Monad (when)
+import Data.Bits
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Search ( gallopingSearchRightPBounds
+                                     , gallopingSearchLeftPBounds
+                                     )
+import Data.Vector.Algorithms.Insertion (sortByBounds', Comparison)
+import Data.Vector.Algorithms.Common (uniqueMutableBy)
+
+-- | Sorts an array using the default comparison.
+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
+sort = sortBy compare
+{-# INLINE sort #-}
+
+-- | A variant on `sort` that returns a vector of unique elements.
+sortUniq :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m (v (PrimState m) e)
+sortUniq = sortUniqBy compare
+{-# INLINE sortUniq #-}
+
+-- | Sorts an array using a custom comparison.
+sortBy :: (PrimMonad m, MVector v e)
+       => Comparison e -> v (PrimState m) e -> m ()
+sortBy cmp vec
+  | mr == len = iter [0] 0 (error "no merge buffer needed!")
+  | otherwise = new 256 >>= iter [] 0
+ where
+ len = length vec
+ mr = minrun len
+ iter s i tmpBuf
+   | i >= len  = performRemainingMerges s tmpBuf
+   | otherwise = do (order, runLen) <- nextRun cmp vec i len
+                    when (order == Descending) $
+                      reverse $ unsafeSlice i runLen vec
+                    let runEnd = min len (i + max runLen mr)
+                    sortByBounds' cmp vec i (i+runLen) runEnd
+                    (s', tmpBuf') <- performMerges (i : s) runEnd tmpBuf
+                    iter s' runEnd tmpBuf'
+ runLengthInvariantBroken a b c i = (b - a <= i - b) || (c - b <= i - c)
+ performMerges [b,a] i tmpBuf
+   | i - b >= b - a = merge cmp vec a b i tmpBuf >>= performMerges [a] i
+ performMerges (c:b:a:ss) i tmpBuf
+   | runLengthInvariantBroken a b c i =
+     if i - c <= b - a
+       then merge cmp vec b c i tmpBuf >>= performMerges (b:a:ss) i
+       else do tmpBuf' <- merge cmp vec a b c tmpBuf
+               (ass', tmpBuf'') <- performMerges (a:ss) c tmpBuf'
+               performMerges (c:ass') i tmpBuf''
+ performMerges s _ tmpBuf = return (s, tmpBuf)
+ performRemainingMerges (b:a:ss) tmpBuf =
+   merge cmp vec a b len tmpBuf >>= performRemainingMerges (a:ss)
+ performRemainingMerges _ _ = return ()
+{-# INLINE sortBy #-}
+
+-- | A variant on `sortBy` which returns a vector of unique elements.
+sortUniqBy :: (PrimMonad m, MVector v e)
+       => Comparison e -> v (PrimState m) e -> m (v (PrimState m) e)
+sortUniqBy cmp vec = do
+  sortBy cmp vec
+  uniqueMutableBy cmp vec
+{-# INLINE sortUniqBy #-}
+
+-- | Computes the minimum run size for the sort. The goal is to choose a size
+-- such that there are almost if not exactly 2^n chunks of that size in the
+-- array.
+minrun :: Int -> Int
+minrun n0 = (n0 `unsafeShiftR` extra) + if (lowMask .&. n0) > 0 then 1 else 0
+ where
+ -- smear the bits down from the most significant bit
+ !n1 = n0 .|. unsafeShiftR n0 1
+ !n2 = n1 .|. unsafeShiftR n1 2
+ !n3 = n2 .|. unsafeShiftR n2 4
+ !n4 = n3 .|. unsafeShiftR n3 8
+ !n5 = n4 .|. unsafeShiftR n4 16
+ !n6 = n5 .|. unsafeShiftR n5 32
+
+ -- mask for the bits lower than the 6 highest bits
+ !lowMask = n6 `unsafeShiftR` 6
+
+ !extra = popCount lowMask
+{-# INLINE minrun #-}
+
+data Order = Ascending | Descending deriving (Eq, Show)
+
+-- | Identify the next run (that is a monotonically increasing or strictly
+-- decreasing sequence) in the slice [l,u) in vec. Returns the order and length
+-- of the run.
+nextRun :: (PrimMonad m, MVector v e)
+        => Comparison e
+        -> v (PrimState m) e
+        -> Int -- ^ l
+        -> Int -- ^ u
+        -> m (Order, Int)
+nextRun _ _ i len | i+1 >= len = return (Ascending, 1)
+nextRun cmp vec i len = do x <- unsafeRead vec i
+                           y <- unsafeRead vec (i+1)
+                           if x `gt` y then desc y 2 else asc  y 2
+ where
+ gt a b = cmp a b == GT
+ desc _ !k | i + k >= len = return (Descending, k)
+ desc x !k = do y <- unsafeRead vec (i+k)
+                if x `gt` y then desc y (k+1) else return (Descending, k)
+ asc _ !k | i + k >= len = return (Ascending, k)
+ asc x !k = do y <- unsafeRead vec (i+k)
+               if x `gt` y then return (Ascending, k) else asc y (k+1)
+{-# INLINE nextRun #-}
+
+-- | Tests if a temporary buffer has a given size. If not, allocates a new
+-- buffer and returns it instead of the old temporary buffer.
+ensureCapacity :: (PrimMonad m, MVector v e)
+               => Int -> v (PrimState m) e -> m (v (PrimState m) e)
+ensureCapacity l tmpBuf
+  | l <= length tmpBuf = return tmpBuf
+  | otherwise          = new (2*l)
+{-# INLINE ensureCapacity #-}
+
+-- | Copy the slice [i,i+len) from vec to tmpBuf. If tmpBuf is not large enough,
+-- a new buffer is allocated and used. Returns the buffer.
+cloneSlice :: (PrimMonad m, MVector v e)
+           => Int -- ^ i
+           -> Int -- ^ len
+           -> v (PrimState m) e -- ^ vec
+           -> v (PrimState m) e -- ^ tmpBuf
+           -> m (v (PrimState m) e)
+cloneSlice i len vec tmpBuf = do
+  tmpBuf' <- ensureCapacity len tmpBuf
+  unsafeCopy (unsafeSlice 0 len tmpBuf') (unsafeSlice i len vec)
+  return tmpBuf'
+{-# INLINE cloneSlice #-}
+
+-- | Number of consecutive times merge chooses the element from the same run
+-- before galloping mode is activated.
+minGallop :: Int
+minGallop = 7
+{-# INLINE minGallop #-}
+
+-- | Merge the adjacent sorted slices [l,m) and [m,u) in vec. This is done by
+-- copying the slice [l,m) to a temporary buffer. Returns the (enlarged)
+-- temporary buffer.
+mergeLo :: (PrimMonad m, MVector v e)
+        => Comparison e
+        -> v (PrimState m) e -- ^ vec
+        -> Int -- ^ l
+        -> Int -- ^ m
+        -> Int -- ^ u
+        -> v (PrimState m) e -- ^ tmpBuf
+        -> m (v (PrimState m) e)
+mergeLo cmp vec l m u tempBuf' = do
+  tmpBuf <- cloneSlice l tmpBufLen vec tempBuf'
+  vi <- unsafeRead tmpBuf 0
+  vj <- unsafeRead vec m
+  iter tmpBuf 0 m l vi vj minGallop minGallop
+  return tmpBuf
+ where
+ gt  a b = cmp a b == GT
+ gte a b = cmp a b /= LT
+ tmpBufLen = m - l
+
+ finalize tmpBuf i k = do
+   let from = unsafeSlice i (tmpBufLen-i) tmpBuf
+       to   = unsafeSlice k (tmpBufLen-i) vec
+   unsafeCopy to from
+
+ iter _ i _ _ _ _ _ _ | i >= tmpBufLen = return ()
+ iter tmpBuf i j k _ _ _ _ | j >= u = finalize tmpBuf i k
+ iter tmpBuf i j k _ vj 0 _ = do
+   i' <- gallopingSearchLeftPBounds (`gt` vj) tmpBuf i tmpBufLen
+   let gallopLen = i' - i
+       from = unsafeSlice i gallopLen tmpBuf
+       to   = unsafeSlice k gallopLen vec
+   unsafeCopy to from
+   when (i' < tmpBufLen) $ do
+     vi' <- unsafeRead tmpBuf i'
+     iter tmpBuf i' j (k+gallopLen) vi' vj minGallop minGallop
+ iter tmpBuf i j k vi _ _ 0 = do
+   j' <- gallopingSearchLeftPBounds (`gte` vi) vec j u
+   let gallopLen = j' - j
+       from = slice j gallopLen vec
+       to   = slice k gallopLen vec
+   unsafeMove to from
+   if j' >= u then finalize tmpBuf i (k + gallopLen) else do
+     vj' <- unsafeRead vec j'
+     iter tmpBuf i j' (k+gallopLen) vi vj' minGallop minGallop
+ iter tmpBuf i j k vi vj ga gb
+   | vj `gte` vi = do unsafeWrite vec k vi
+                      when (i + 1 < tmpBufLen) $ do
+                        vi' <- unsafeRead tmpBuf (i+1)
+                        iter tmpBuf (i+1) j (k+1) vi' vj (ga-1) minGallop
+   | otherwise   = do unsafeWrite vec k vj
+                      if j + 1 >= u then finalize tmpBuf i (k + 1) else do
+                        vj' <- unsafeRead vec (j+1)
+                        iter tmpBuf i (j+1) (k+1) vi vj' minGallop (gb-1)
+{-# INLINE mergeLo #-}
+
+-- | Merge the adjacent sorted slices [l,m) and [m,u) in vec. This is done by
+-- copying the slice [j,k) to a temporary buffer. Returns the (enlarged)
+-- temporary buffer.
+mergeHi :: (PrimMonad m, MVector v e)
+        => Comparison e
+        -> v (PrimState m) e -- ^ vec
+        -> Int -- ^ l
+        -> Int -- ^ m
+        -> Int -- ^ u
+        -> v (PrimState m) e -- ^ tmpBuf
+        -> m (v (PrimState m) e)
+mergeHi cmp vec l m u tmpBuf' = do
+  tmpBuf <- cloneSlice m tmpBufLen vec tmpBuf'
+  vi <- unsafeRead vec (m-1)
+  vj <- unsafeRead tmpBuf (tmpBufLen-1)
+  iter tmpBuf (m-1) (tmpBufLen-1) (u-1) vi vj minGallop minGallop
+  return tmpBuf
+ where
+ gt  a b = cmp a b == GT
+ gte a b = cmp a b /= LT
+ tmpBufLen = u - m
+
+ finalize tmpBuf j = do
+   let from = unsafeSlice 0 (j+1) tmpBuf
+       to   = unsafeSlice l (j+1) vec
+   unsafeCopy to from
+
+ iter _ _ j _ _ _ _ _ | j < 0 = return ()
+ iter tmpBuf i j _ _ _ _ _ | i < l = finalize tmpBuf j
+ iter tmpBuf i j k _ vj 0 _ = do
+   i' <- gallopingSearchRightPBounds (`gt` vj) vec l i
+   let gallopLen = i - i'
+       from = slice (i'+1) gallopLen vec
+       to   = slice (k-gallopLen+1) gallopLen vec
+   unsafeMove to from
+   if i' < l then finalize tmpBuf j else do
+     vi' <- unsafeRead vec i'
+     iter tmpBuf i' j (k-gallopLen) vi' vj minGallop minGallop
+ iter tmpBuf i j k vi _ _ 0 = do
+   j' <- gallopingSearchRightPBounds (`gte` vi) tmpBuf 0 j
+   let gallopLen = j - j'
+       from = slice (j'+1) gallopLen tmpBuf
+       to   = slice (k-gallopLen+1) gallopLen vec
+   unsafeCopy to from
+   when (j' >= 0) $ do
+     vj' <- unsafeRead tmpBuf j'
+     iter tmpBuf i j' (k-gallopLen) vi vj' minGallop minGallop
+ iter tmpBuf i j k vi vj ga gb
+   | vi `gt` vj = do unsafeWrite vec k vi
+                     if i - 1 < l then finalize tmpBuf j else do
+                       vi' <- unsafeRead vec (i-1)
+                       iter tmpBuf (i-1) j (k-1) vi' vj (ga-1) minGallop
+   | otherwise  = do unsafeWrite vec k vj
+                     when (j - 1 >= 0) $ do
+                       vj' <- unsafeRead tmpBuf (j-1)
+                       iter tmpBuf i (j-1) (k-1) vi vj' minGallop (gb-1)
+{-# INLINE mergeHi #-}
+
+-- | Merge the adjacent sorted slices A=[l,m) and B=[m,u) in vec. This begins
+-- with galloping searches to find the index of vec[m] in A and the index of
+-- vec[m-1] in B to reduce the sizes of A and B. Then it uses `mergeHi` or
+-- `mergeLo` depending on whether A or B is larger. Returns the (enlarged)
+-- temporary buffer.
+merge :: (PrimMonad m, MVector v e)
+      => Comparison e
+      -> v (PrimState m) e -- ^ vec
+      -> Int -- ^ l
+      -> Int -- ^ m
+      -> Int -- ^ u
+      -> v (PrimState m) e -- ^ tmpBuf
+      -> m (v (PrimState m) e)
+merge cmp vec l m u tmpBuf = do
+  vm <- unsafeRead vec m
+  l' <- gallopingSearchLeftPBounds (`gt` vm) vec l m
+  if l' >= m
+    then return tmpBuf
+    else do
+      vn <- unsafeRead vec (m-1)
+      u' <- gallopingSearchRightPBounds (`gte` vn) vec m u
+      if u' <= m
+        then return tmpBuf
+        else (if (m-l') <= (u'-m) then mergeLo else mergeHi) cmp vec l' m u' tmpBuf
+ where
+ gt  a b = cmp a b == GT
+ gte a b = cmp a b /= LT
+{-# INLINE merge #-}
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 qualified Data.List as 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 []
+                        . List.groupBy ((==) `on` fst)
+                        . flip zip [0..])
+              $ monotones (n-2) n
+
+sort2 :: (Vector v Int) => [v Int]
+sort2 = map fromList $ List.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 $ List.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 $ List.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,224 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE TypeOperators #-}
+
+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 qualified Data.Vector.Generic as GV
+
+import Data.Vector.Algorithms.Optimal (Comparison)
+import Data.Vector.Algorithms.Radix (radix, passes, size)
+import qualified Data.Vector.Algorithms as Alg
+
+import qualified Data.Map as M
+
+import Test.QuickCheck hiding (Sorted)
+
+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_sorted_uniq :: (Ord e) => Vector e -> Property
+prop_sorted_uniq 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
+
+runFreeze
+  :: forall e . (Ord e)
+  => (forall s mv . G.MVector mv e => mv s e -> ST s (mv s e))
+  -> (forall s v mv. (GV.Vector v e, mv ~ GV.Mutable v) => mv s e -> ST s (v e))
+runFreeze alg mv = do
+  mv <- alg mv
+  GV.unsafeFreeze mv
+
+prop_full_sortUniq
+  :: (Ord e, Show e)
+  => (forall s . MV.MVector s e -> ST s (Vector e))
+  -> Vector e -> Property
+prop_full_sortUniq algo arr = runST $ do
+  mv <- V.unsafeThaw arr
+  arr' <- algo mv
+  pure (prop_sorted_uniq 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 v -> do
+  let newVec = modify algo v
+      vhead = V.take k newVec
+      vtail = V.drop k newVec
+  prop_sorted vhead
+    .&&.
+      -- Every element in the head should be < every element in the tail.
+      if V.null vtail then 1 == 1 else V.maximum vhead <= V.minimum vtail
+
+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 (<=) (>)
+
+prop_nub :: (Ord e, Show e) => Vector e -> Property
+prop_nub v =
+  V.fromList (nub (V.toList v)) === Alg.nub v
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,230 @@
+{-# LANGUAGE RankNTypes, TypeOperators, FlexibleContexts, TypeApplications #-}
+
+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 qualified Data.Vector.Mutable as BoxedMV
+
+import qualified Data.Vector.Generic as G
+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.Tim          as T
+
+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
+type MonoAlgo  e r = forall s . BoxedMV.MVector s e -> ST s r
+
+newtype WrappedAlgo      e r = WrapAlgo      { unWrapAlgo      :: Algo      e r }
+newtype WrappedSizeAlgo  e r = WrapSizeAlgo  { unWrapSizeAlgo  :: SizeAlgo  e r }
+newtype WrappedBoundAlgo e r = WrapBoundAlgo { unWrapBoundAlgo :: BoundAlgo e r }
+newtype WrappedMonoAlgo  e r = MonoAlgo      { unWrapMonoAlgo  :: MonoAlgo  e r }
+
+args = stdArgs
+       { maxSuccess = 1000
+       , maxDiscardRatio = 2
+       }
+
+check_Int_sort = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_fullsort (unWrapAlgo algo))
+ where
+ algos :: [(String, WrappedAlgo Int ())]
+ algos = [ ("introsort", WrapAlgo INT.sort)
+         , ("insertion sort", WrapAlgo INS.sort)
+         , ("merge sort", WrapAlgo M.sort)
+         , ("heapsort", WrapAlgo H.sort)
+         , ("timsort", WrapAlgo T.sort)
+         ]
+
+check_Int_sortUniq = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_full_sortUniq (unWrapMonoAlgo algo))
+ where
+ algos :: [(String, WrappedMonoAlgo Int (Vector Int))]
+ algos = [ ("intro_sortUniq", MonoAlgo (runFreeze INT.sortUniq))
+         , ("insertion sortUniq", MonoAlgo (runFreeze INS.sortUniq))
+         , ("merge sortUniq", MonoAlgo (runFreeze M.sortUniq))
+         , ("heap_sortUniq", MonoAlgo (runFreeze H.sortUniq))
+         , ("tim_sortUniq", MonoAlgo (runFreeze T.sortUniq))
+         ]
+
+check_Int_partialsort = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_partialsort (unWrapSizeAlgo algo))
+ where
+ algos :: [(String, WrappedSizeAlgo Int ())]
+ algos = [ ("intro-partialsort", WrapSizeAlgo INT.partialSort)
+         , ("heap partialsort", WrapSizeAlgo H.partialSort)
+         ]
+
+check_Int_select = forM_ algos $ \(name,algo) ->
+  quickCheckWith args (label name . prop_select (unWrapSizeAlgo algo))
+ where
+ algos :: [(String, WrappedSizeAlgo Int ())]
+ algos = [ ("intro-select", WrapSizeAlgo INT.select)
+         , ("heap select", WrapSizeAlgo H.select)
+         ]
+
+check_nub = quickCheckWith args (label "nub Int" . (prop_nub @Int))
+
+
+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)
+                  quickCheckWith args (label "tim sort" . prop_stable T.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 "heapsort"     . prop_permutation (H.sort :: Algo Int ())
+
+  qc $ label "mergesort"    . prop_permutation (M.sort :: Algo Int    ())
+  qc $ label "timsort"      . prop_permutation (T.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 ())
+  qc $ label "intropartial" . prop_sized (\x -> const (prop_permutation x))
+                                         (INT.partialSort :: SizeAlgo Int ())
+  qc $ label "introselect"  . prop_sized (\x -> const (prop_permutation x))
+                                         (INT.select :: SizeAlgo Int ())
+  qc $ label "heappartial"  . prop_sized (\x -> const (prop_permutation x))
+                                         (H.partialSort :: SizeAlgo Int ())
+  qc $ label "heapselect"   . prop_sized (\x -> const (prop_permutation x))
+                                         (H.select :: SizeAlgo Int ())
+
+ 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 "timsort empty"      $ prop_empty       (T.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_sortUniq
+          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
+          putStrLn "Algorithms:"
+          check_nub
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,161 @@
-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
+cabal-version:     >= 1.10
+name:              vector-algorithms
+version:           0.9.1.0
+license:           BSD3
+license-file:      LICENSE
+author:            Dan Doel
+maintainer:        Dan Doel <dan.doel@gmail.com>
+                   Erik de Castro Lopo <erikd@mega-nerd.com>
+copyright:         (c) 2008,2009,2010,2011,2012,2013,2014,2015 Dan Doel
+                   (c) 2015 Tim Baumann
+homepage:          https://github.com/erikd/vector-algorithms/
+category:          Data
+synopsis:          Efficient algorithms for vector arrays
+description:       Efficient algorithms for sorting vector arrays. At some stage
+                   other vector algorithms may be added.
+build-type:        Simple
 
-Flag BoundsChecks
-  Description: Enable bounds checking
-  Default: True
+extra-source-files: CHANGELOG.md
 
-Flag UnsafeChecks
-  Description: Enable bounds checking in unsafe operations at the cost of a
+tested-with:
+  GHC == 9.12.1
+  GHC == 9.10.1
+  GHC == 9.8.2
+  GHC == 9.6.3
+  GHC == 9.4.7
+  GHC == 9.2.8
+  GHC == 9.0.2
+  GHC == 8.10.7
+  GHC == 8.8.4
+  GHC == 8.6.5
+  GHC == 8.4.4
+  GHC == 8.2.2
+
+flag BoundsChecks
+  description: Enable bounds checking
+  default: True
+
+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: True
 
-    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 dump-simpl
+--   description: Dumps the simplified core during compilation
+--   default: False
 
-    Other-Modules:
-        Data.Vector.Algorithms.Common
+flag llvm
+  description: Build using llvm
+  default: False
 
-    Extensions:
-        BangPatterns,
-        TypeOperators,
-        Rank2Types,
-        ScopedTypeVariables,
-        FlexibleContexts,
-        CPP
+source-repository head
+  type:     git
+  location: https://github.com/erikd/vector-algorithms/
 
-    GHC-Options:
-        -Odph
-        -funbox-strict-fields
+library
+  hs-source-dirs: src
+  default-language: Haskell2010
 
-    Include-Dirs:
-        include
+  build-depends: base >= 4.8 && < 5,
+                 bitvec >= 1.0 && < 1.2,
+                 vector >= 0.6 && < 0.14,
+                 primitive >= 0.6.2.0 && < 0.10,
+                 bytestring >= 0.9 && < 1
 
-    Install-Includes:
-        vector.h
+  if ! impl (ghc >= 7.8)
+    build-depends: tagged >= 0.4 && < 0.9
 
-    if flag(BoundsChecks)
-        cpp-options: -DVECTOR_BOUNDS_CHECKS
+  exposed-modules:
+    Data.Vector.Algorithms
+    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
+    Data.Vector.Algorithms.Tim
 
-    if flag(UnsafeChecks)
-        cpp-options: -DVECTOR_UNSAFE_CHECKS
+  other-modules:
+    Data.Vector.Algorithms.Common
 
-    if flag(InternalChecks)
-        cpp-options: -DVECTOR_INTERNAL_CHECKS
+  ghc-options:
+    -funbox-strict-fields
 
+  -- Cabal/Hackage complains about these
+  -- if flag(dump-simpl)
+  --   ghc-options: -ddump-simpl -ddump-to-file
+
+  if flag(llvm)
+    ghc-options: -fllvm
+
+  include-dirs:
+    include
+
+  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
+
+benchmark simple-bench
+  hs-source-dirs: bench/simple
+  type: exitcode-stdio-1.0
+  default-language: Haskell2010
+
+  if !flag(bench)
+    buildable: False
+
+  main-is: Main.hs
+
+  other-modules:
+    Blocks
+
+  build-depends: base, mwc-random, vector, vector-algorithms
+  ghc-options: -Wall
+
+  -- Cabal/Hackage complains about these
+  -- if flag(dump-simpl)
+  --   ghc-options: -ddump-simpl -ddump-to-file
+
+  if flag(llvm)
+    ghc-options: -fllvm
+
+test-suite properties
+  hs-source-dirs: tests/properties
+  type: exitcode-stdio-1.0
+  main-is: Tests.hs
+  default-language: Haskell2010
+
+  other-modules:
+    Optimal
+    Properties
+    Util
+
+  build-depends:
+    base >= 4.9,
+    bytestring,
+    containers,
+    QuickCheck > 2.9 && < 2.16,
+    vector,
+    vector-algorithms
+
+  if flag(llvm)
+    ghc-options: -fllvm