diff --git a/Data/Vector/Algorithms/AmericanFlag.hs b/Data/Vector/Algorithms/AmericanFlag.hs
new file mode 100644
--- /dev/null
+++ b/Data/Vector/Algorithms/AmericanFlag.hs
@@ -0,0 +1,300 @@
+{-# 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 ()
+--
+-- 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).
+
+module Data.Vector.Algorithms.AmericanFlag ( sort
+                                           , sortBy
+                                           ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import Data.Word
+import Data.Int
+import Data.Bits
+
+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
+
+class Lexicographic e where
+  terminate :: e -> Int -> Bool
+  size      :: e -> Int
+  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 #-}
+
+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
+{-# INLINE sort #-}
+
+sortBy :: (PrimMonad m, MVector v e)
+       => Comparison e
+       -> (e -> Int -> Bool) -- a stopping predicate
+       -> Int                -- size of auxiliary arrays
+       -> (Int -> e -> Int)  -- big-endian radix
+       -> 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/Common.hs b/Data/Vector/Algorithms/Common.hs
--- a/Data/Vector/Algorithms/Common.hs
+++ b/Data/Vector/Algorithms/Common.hs
@@ -1,8 +1,9 @@
+{-# LANGUAGE FlexibleContexts #-}
 
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Common
--- Copyright   : (c) 2008-2010 Dan Doel
+-- Copyright   : (c) 2008-2011 Dan Doel
 -- Maintainer  : Dan Doel
 -- Stability   : Experimental
 -- Portability : Portable
@@ -17,6 +18,8 @@
 
 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
 
@@ -25,3 +28,20 @@
 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
new file mode 100644
--- /dev/null
+++ b/Data/Vector/Algorithms/Heap.hs
@@ -0,0 +1,240 @@
+{-# LANGUAGE TypeOperators #-}
+
+-- ---------------------------------------------------------------------------
+-- |
+-- Module      : Data.Vector.Algorithms.Heap
+-- Copyright   : (c) 2008-2011 Dan Doel
+-- Maintainer  : Dan Doel <dan.doel@gmail.com>
+-- Stability   : Experimental
+-- Portability : Non-portable (type operators)
+--
+-- This module implements operations for working with a quaternary heap stored
+-- in an unboxed array. Most heapsorts are defined in terms of a binary heap,
+-- in which each internal node has at most two children. By contrast, a
+-- quaternary heap has internal nodes with up to four children. This reduces
+-- the number of comparisons in a heapsort slightly, and improves locality
+-- (again, slightly) by flattening out the heap.
+
+module Data.Vector.Algorithms.Heap
+       ( -- * Sorting
+         sort
+       , sortBy
+       , sortByBounds
+         -- * Selection
+       , select
+       , selectBy
+       , selectByBounds
+         -- * Partial sorts
+       , partialSort
+       , partialSortBy
+       , partialSortByBounds
+         -- * Heap operations
+       , heapify
+       , pop
+       , popTo
+       , sortHeap
+       , Comparison
+       ) where
+
+import Prelude hiding (read, length)
+
+import Control.Monad
+import Control.Monad.Primitive
+
+import Data.Bits
+
+import Data.Vector.Generic.Mutable
+
+import Data.Vector.Algorithms.Common (Comparison)
+
+import qualified Data.Vector.Algorithms.Optimal as O
+
+-- | Sorts an entire array using the default ordering.
+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()
+sort = sortBy compare
+{-# INLINE 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/Intro.hs b/Data/Vector/Algorithms/Intro.hs
--- a/Data/Vector/Algorithms/Intro.hs
+++ b/Data/Vector/Algorithms/Intro.hs
@@ -3,7 +3,7 @@
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Intro
--- Copyright   : (c) 2008-2010 Dan Doel
+-- Copyright   : (c) 2008-2011 Dan Doel
 -- Maintainer  : Dan Doel <dan.doel@gmail.com>
 -- Stability   : Experimental
 -- Portability : Non-portable (type operators, bang patterns)
@@ -57,7 +57,7 @@
 
 import qualified Data.Vector.Algorithms.Insertion as I
 import qualified Data.Vector.Algorithms.Optimal   as O
-import qualified Data.Vector.Algorithms.TriHeap   as H
+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 ()
@@ -118,7 +118,9 @@
 -- [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 = go (ilg len) l (l + k) u
+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
@@ -150,7 +152,9 @@
 -- [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 = go (ilg len) l (l + k) u
+partialSortByBounds cmp a k l u
+  | l >= u    = return ()
+  | otherwise = go (ilg len) l (l + k) u
  where
  len = u - l
  go 0 l m n = H.partialSortByBounds cmp a (m - l) l u
diff --git a/Data/Vector/Algorithms/Merge.hs b/Data/Vector/Algorithms/Merge.hs
--- a/Data/Vector/Algorithms/Merge.hs
+++ b/Data/Vector/Algorithms/Merge.hs
@@ -1,7 +1,7 @@
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Merge
--- Copyright   : (c) 2008-2010 Dan Doel
+-- Copyright   : (c) 2008-2011 Dan Doel
 -- Maintainer  : Dan Doel <dan.doel@gmail.com>
 -- Stability   : Experimental
 -- Portability : Portable
@@ -72,17 +72,22 @@
  upper = unsafeSlice mid (length src - mid) src
  tmp   = unsafeSlice 0   mid                buf
 
- loop low iLow eLow high iHigh eHigh iIns
-   | iLow  >= length low  = return ()
+ wroteHigh low iLow eLow high iHigh iIns
    | iHigh >= length high = unsafeCopy (unsafeSlice iIns (length low - iLow) src)
                                        (unsafeSlice iLow (length low - iLow) low)
-   | otherwise            = case cmp eHigh eLow of
-                             LT -> do unsafeWrite src iIns eHigh
-                                      eHigh <- unsafeRead high (iHigh + 1)
-                                      loop low iLow eLow high (iHigh + 1) eHigh (iIns + 1)
-                             _  -> do unsafeWrite src iIns eLow
-                                      eLow <- unsafeRead low (iLow + 1)
-                                      loop low (iLow + 1) eLow high iHigh eHigh (iIns + 1)
+   | 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
diff --git a/Data/Vector/Algorithms/Radix.hs b/Data/Vector/Algorithms/Radix.hs
--- a/Data/Vector/Algorithms/Radix.hs
+++ b/Data/Vector/Algorithms/Radix.hs
@@ -3,7 +3,7 @@
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Radix
--- Copyright   : (c) 2008-2010 Dan Doel
+-- Copyright   : (c) 2008-2011 Dan Doel
 -- Maintainer  : Dan Doel <dan.doel@gmail.com>
 -- Stability   : Experimental
 -- Portability : Non-portable (scoped type variables, bang patterns)
@@ -43,7 +43,7 @@
 import qualified Data.Vector.Primitive.Mutable as PV
 import Data.Vector.Generic.Mutable
 
-import Data.Vector.Algorithms.Common (Comparison)
+import Data.Vector.Algorithms.Common
 
 import Data.Bits
 import Data.Int
@@ -200,8 +200,7 @@
 sortBy passes size rdx arr = do
   tmp    <- new (length arr)
   count  <- new size
-  prefix <- new size
-  radixLoop passes rdx arr tmp count prefix
+  radixLoop passes rdx arr tmp count
 {-# INLINE sortBy #-}
 
 radixLoop :: (PrimMonad m, MVector v e)
@@ -210,15 +209,14 @@
           -> v (PrimState m) e            -- array to sort
           -> v (PrimState m) e            -- temporary array
           -> PV.MVector (PrimState m) Int -- radix count array
-          -> PV.MVector (PrimState m) Int -- placement array
           -> m ()
-radixLoop passes rdx src dst count prefix = go False 0
+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 prefix k >> go (not swap) (k+1)
-                    else body rdx src dst count prefix k >> go (not swap) (k+1)
+                    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 #-}
 
@@ -227,41 +225,25 @@
      -> v (PrimState m) e            -- source array
      -> v (PrimState m) e            -- destination array
      -> PV.MVector (PrimState m) Int -- radix count
-     -> PV.MVector (PrimState m) Int -- placement
      -> Int                          -- current pass
      -> m ()
-body rdx src dst count prefix k = do
-  set count 0
-  countLoop k rdx src count
-  unsafeWrite prefix 0 0
-  prefixLoop count prefix
-  moveLoop k rdx src dst prefix
+body rdx src dst count k = do
+  countLoop (rdx k) src count
+  accumulate count
+  moveLoop k rdx src dst count
 {-# INLINE body #-}
 
-countLoop :: (PrimMonad m, MVector v e)
-          => Int -> (Int -> e -> Int)
-          -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m ()
-countLoop k rdx src count = go 0
- where
- len = length src
- go i
-   | i < len    = unsafeRead src i >>= inc count . rdx k >> go (i+1)
-   | otherwise  = return ()
-{-# INLINE countLoop #-}
-
-prefixLoop :: (PrimMonad m)
-           => PV.MVector (PrimState m) Int -> PV.MVector (PrimState m) Int
-           -> m ()
-prefixLoop count prefix = go 1 0
+accumulate :: (PrimMonad m)
+           => PV.MVector (PrimState m) Int -> m ()
+accumulate count = go 0 0
  where
  len = length count
- go i pi
-   | i < len   = do ci <- unsafeRead count (i-1)
-                    let pi' = pi + ci
-                    unsafeWrite prefix i pi'
-                    go (i+1) pi'
+ go i acc
+   | i < len   = do ci <- unsafeRead count i
+                    unsafeWrite count i acc
+                    go (i+1) (acc + ci)
    | otherwise = return ()
-{-# INLINE prefixLoop #-}
+{-# INLINE accumulate #-}
 
 moveLoop :: (PrimMonad m, MVector v e)
          => Int -> (Int -> e -> Int) -> v (PrimState m) e
@@ -277,6 +259,3 @@
    | otherwise  = return ()
 {-# INLINE moveLoop #-}
 
-inc :: (PrimMonad m) => PV.MVector (PrimState m) Int -> Int -> m Int
-inc arr i = unsafeRead arr i >>= \e -> unsafeWrite arr i (e+1) >> return e
-{-# INLINE inc #-}
diff --git a/Data/Vector/Algorithms/TriHeap.hs b/Data/Vector/Algorithms/TriHeap.hs
deleted file mode 100644
--- a/Data/Vector/Algorithms/TriHeap.hs
+++ /dev/null
@@ -1,218 +0,0 @@
-{-# LANGUAGE TypeOperators #-}
-
--- ---------------------------------------------------------------------------
--- |
--- Module      : Data.Vector.Algorithms.TriHeap
--- Copyright   : (c) 2008-2010 Dan Doel
--- Maintainer  : Dan Doel <dan.doel@gmail.com>
--- Stability   : Experimental
--- Portability : Non-portable (type operators)
---
--- This module implements operations for working with a trinary 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
--- trinary heap has internal nodes with up to three 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.TriHeap
-       ( -- * 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.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
-{-# INLINE 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) `div` 3
-  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 * 3 + 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
-  | 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
-{-# INLINE maximumChild #-}
diff --git a/bench/Main.hs b/bench/Main.hs
--- a/bench/Main.hs
+++ b/bench/Main.hs
@@ -14,11 +14,12 @@
 
 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.TriHeap   as TH
-import qualified Data.Vector.Algorithms.Merge     as M
-import qualified Data.Vector.Algorithms.Radix     as R
+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
@@ -67,11 +68,12 @@
                | IntroSort
                | IntroPartialSort
                | IntroSelect
-               | TriHeapSort
-               | TriHeapPartialSort
-               | TriHeapSelect
+               | HeapSort
+               | HeapPartialSort
+               | HeapSelect
                | MergeSort
                | RadixSort
+               | AmericanFlagSort
                deriving (Show, Read, Enum, Bounded)
 
 data Options = O { algos :: [Algorithm], elems :: Int, portion :: Int, usage :: Bool } deriving (Show)
@@ -103,7 +105,7 @@
 
 parseAlgo :: String -> Options -> Either String Options
 parseAlgo "None" o = Right $ o { algos = [] }
-parseAlgo "All"  o = Right $ o { algos = [DoNothing .. RadixSort] }
+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
 
@@ -131,11 +133,12 @@
   IntroSort          -> sortSuite        "introsort"             g n   introSort
   IntroPartialSort   -> partialSortSuite "partial introsort"     g n k introPSort
   IntroSelect        -> partialSortSuite "introselect"           g n k introSelect
-  TriHeapSort        -> sortSuite        "tri-heap sort"         g n   triHeapSort
-  TriHeapPartialSort -> partialSortSuite "partial tri-heap sort" g n k triHeapPSort
-  TriHeapSelect      -> partialSortSuite "tri-heap select"       g n k triHeapSelect
+  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 ()
@@ -154,17 +157,17 @@
 introSelect v k = INT.select v k
 {-# NOINLINE introSelect #-}
 
-triHeapSort :: MVector RealWorld Int -> IO ()
-triHeapSort v = TH.sort v
-{-# NOINLINE triHeapSort #-}
+heapSort :: MVector RealWorld Int -> IO ()
+heapSort v = H.sort v
+{-# NOINLINE heapSort #-}
 
-triHeapPSort :: MVector RealWorld Int -> Int -> IO ()
-triHeapPSort v k = TH.partialSort v k
-{-# NOINLINE triHeapPSort #-}
+heapPSort :: MVector RealWorld Int -> Int -> IO ()
+heapPSort v k = H.partialSort v k
+{-# NOINLINE heapPSort #-}
 
-triHeapSelect :: MVector RealWorld Int -> Int -> IO ()
-triHeapSelect v k = TH.select v k
-{-# NOINLINE triHeapSelect #-}
+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
@@ -173,6 +176,10 @@
 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
diff --git a/tests/Optimal.hs b/tests/Optimal.hs
--- a/tests/Optimal.hs
+++ b/tests/Optimal.hs
@@ -11,7 +11,7 @@
 import Data.List
 import Data.Function
 
-import Data.Vector.Generic hiding (map, zip, concatMap, (++), replicate)
+import Data.Vector.Generic hiding (map, zip, concatMap, (++), replicate, foldM)
 
 interleavings :: [a] -> [a] -> [[a]]
 interleavings [       ] ys        =  [ys]
diff --git a/tests/Properties.hs b/tests/Properties.hs
--- a/tests/Properties.hs
+++ b/tests/Properties.hs
@@ -18,6 +18,8 @@
 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)
@@ -36,9 +38,12 @@
  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 $ apply algo arr
+prop_fullsort algo arr = prop_sorted $ modify algo arr
 
 {-
 prop_schwartzian :: (UA e, UA k, Ord k)
@@ -47,7 +52,7 @@
                  -> UArr e -> Property
 prop_schwartzian f algo arr
   | lengthU arr < 2 = property True
-  | otherwise       = let srt = apply (algo `usingKeys` f) arr
+  | otherwise       = let srt = modify (algo `usingKeys` f) arr
                       in check (headU srt) (tailU srt)
  where
  check e arr | nullU arr = property True
@@ -64,13 +69,16 @@
                  => (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 . apply algo
+  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' = apply algo 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
@@ -86,8 +94,8 @@
 
 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 $ apply algo arr == arr
-prop_stable algo arr = stable $ apply (algo (comparing fst)) $ V.zip arr ix
+-- 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]
 
@@ -100,7 +108,7 @@
                         -> mv s e -> ST s ())
                   -> Vector Int -> Property
 prop_stable_radix algo arr =
-  stable . apply (algo (passes e) (size e) (\k (e, _) -> radix k e))
+  stable . modify (algo (passes e) (size e) (\k (e, _) -> radix k e))
          $ V.zip arr ix
  where
  ix = V.fromList [1 .. V.length arr]
@@ -113,13 +121,13 @@
  where
  arrn  = V.fromList [0..n-1]
  sortn = all ( (== arrn)
-             . apply (\a -> algo compare a 0)
+             . modify (\a -> algo compare a 0)
              . V.fromList)
          $ permutations [0..n-1]
  stabn = all ( (== arrn)
              . snd
              . V.unzip
-             . apply (\a -> algo (comparing fst) a 0))
+             . modify (\a -> algo (comparing fst) a 0))
          $ stability n
 
 type Bag e = M.Map e Int
@@ -130,7 +138,7 @@
 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 (apply algo arr)
+                            toBag arr == toBag (modify algo arr)
 
 newtype SortedVec e = Sorted (Vector e)
 
diff --git a/tests/Tests.hs b/tests/Tests.hs
--- a/tests/Tests.hs
+++ b/tests/Tests.hs
@@ -20,21 +20,22 @@
 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.TriHeap   as TH
-import qualified Data.Vector.Algorithms.Optimal   as O
+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
+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 = 300
+       { maxSuccess = 1000
        , maxDiscard = 200
        }
 
@@ -45,7 +46,7 @@
  algos = [ ("introsort", INT.sort)
          , ("insertion sort", INS.sort)
          , ("merge sort", M.sort)
-         , ("tri-heapsort", TH.sort)
+         , ("heapsort", H.sort)
          ]
 
 check_Int_partialsort = forM_ algos $ \(name,algo) ->
@@ -53,7 +54,7 @@
  where
  algos :: [(String, SizeAlgo Int ())]
  algos = [ ("intro-partialsort", INT.partialSort)
-         , ("tri-heap partialsort", TH.partialSort)
+         , ("heap partialsort", H.partialSort)
          ]
 
 check_Int_select = forM_ algos $ \(name,algo) ->
@@ -61,21 +62,33 @@
  where
  algos :: [(String, SizeAlgo Int ())]
  algos = [ ("intro-select", INT.select)
-         , ("tri-heap select", TH.select)
+         , ("heap select", H.select)
          ]
 
 check_radix_sorts = do
-  qc (label "Word8"       . prop_fullsort (R.sort :: Algo Word8  ()))
-  qc (label "Word16"      . prop_fullsort (R.sort :: Algo Word16 ()))
-  qc (label "Word32"      . prop_fullsort (R.sort :: Algo Word32 ()))
-  qc (label "Word64"      . prop_fullsort (R.sort :: Algo Word64 ()))
-  qc (label "Word"        . prop_fullsort (R.sort :: Algo Word   ()))
-  qc (label "Int8"        . prop_fullsort (R.sort :: Algo Int8   ()))
-  qc (label "Int16"       . prop_fullsort (R.sort :: Algo Int16  ()))
-  qc (label "Int32"       . prop_fullsort (R.sort :: Algo Int32  ()))
-  qc (label "Int64"       . prop_fullsort (R.sort :: Algo Int64  ()))
-  qc (label "Int"         . prop_fullsort (R.sort :: Algo Int    ()))
-  qc (label "(Int, Int)"  . prop_fullsort (R.sort :: Algo (Int, Int) ()))
+  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 (Int, Int)"  . prop_fullsort (R.sort :: Algo (Int, Int) ()))
  where
  qc algo = quickCheckWith args algo
 
@@ -102,11 +115,11 @@
                                          (INT.partialSort :: SizeAlgo Int ())
   qc $ label "introselect"  . prop_sized (const . prop_permutation)
                                          (INT.select :: SizeAlgo Int ())
-  qc $ label "heapsort"     . prop_permutation (TH.sort :: Algo Int ())
+  qc $ label "heapsort"     . prop_permutation (H.sort :: Algo Int ())
   qc $ label "heappartial"  . prop_sized (const . prop_permutation)
-                                         (TH.partialSort :: SizeAlgo Int ())
+                                         (H.partialSort :: SizeAlgo Int ())
   qc $ label "heapselect"   . prop_sized (const . prop_permutation)
-                                         (TH.select :: SizeAlgo Int ())
+                                         (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  ())
@@ -118,9 +131,32 @@
   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   ())
  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 BoundSAlgo e r = forall s mv. MVector mv e => mv s e -> e -> Int -> Int -> ST s r
 
 check_search_range = do
@@ -147,3 +183,5 @@
           check_permutation
           putStrLn "Search in range:"
           check_search_range
+          putStrLn "Corner cases:"
+          check_corners
diff --git a/tests/Util.hs b/tests/Util.hs
--- a/tests/Util.hs
+++ b/tests/Util.hs
@@ -26,29 +26,3 @@
 instance (Arbitrary e) => Arbitrary (V.Vector e) where
   arbitrary = fmap V.fromList arbitrary
 
-instance Arbitrary Int8 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Int16 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Int32 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Int64 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Word8 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Word16 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Word32 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Word64 where
-  arbitrary = fromInteger `fmap` arbitrary
-
-instance Arbitrary Word where
-  arbitrary = fromInteger `fmap` arbitrary
diff --git a/vector-algorithms.cabal b/vector-algorithms.cabal
--- a/vector-algorithms.cabal
+++ b/vector-algorithms.cabal
@@ -1,5 +1,5 @@
 Name:              vector-algorithms
-Version:           0.4
+Version:           0.5.0
 License:           BSD3
 License-File:      LICENSE
 Author:            Dan Doel
@@ -35,7 +35,8 @@
         Data.Vector.Algorithms.Merge
         Data.Vector.Algorithms.Radix
         Data.Vector.Algorithms.Search
-        Data.Vector.Algorithms.TriHeap
+        Data.Vector.Algorithms.Heap
+        Data.Vector.Algorithms.AmericanFlag
 
     Other-Modules:
         Data.Vector.Algorithms.Common
@@ -45,6 +46,7 @@
         TypeOperators,
         Rank2Types,
         ScopedTypeVariables,
+        FlexibleContexts,
         CPP
 
     GHC-Options: