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.
 
diff --git a/bench/Main.hs b/bench/Main.hs
--- a/bench/Main.hs
+++ b/bench/Main.hs
@@ -20,6 +20,7 @@
 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
@@ -74,6 +75,7 @@
                | MergeSort
                | RadixSort
                | AmericanFlagSort
+               | TimSort
                deriving (Show, Read, Enum, Bounded)
 
 data Options = O { algos :: [Algorithm], elems :: Int, portion :: Int, usage :: Bool } deriving (Show)
@@ -139,6 +141,7 @@
   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
   _                  -> putStrLn $ "Currently unsupported algorithm: " ++ show alg
 
 mergeSort :: MVector RealWorld Int -> IO ()
@@ -180,6 +183,10 @@
 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 ->
diff --git a/src/Data/Vector/Algorithms/AmericanFlag.hs b/src/Data/Vector/Algorithms/AmericanFlag.hs
--- a/src/Data/Vector/Algorithms/AmericanFlag.hs
+++ b/src/Data/Vector/Algorithms/AmericanFlag.hs
@@ -292,7 +292,7 @@
                             then unsafeRead count (r-1)
                             else return 0
                     case () of
-                      -- if the current element is alunsafeReady in the right pile,
+                      -- 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
diff --git a/src/Data/Vector/Algorithms/Heap.hs b/src/Data/Vector/Algorithms/Heap.hs
--- a/src/Data/Vector/Algorithms/Heap.hs
+++ b/src/Data/Vector/Algorithms/Heap.hs
@@ -4,7 +4,7 @@
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Heap
--- Copyright   : (c) 2008-2011 Dan Doel
+-- Copyright   : (c) 2008-2015 Dan Doel
 -- Maintainer  : Dan Doel <dan.doel@gmail.com>
 -- Stability   : Experimental
 -- Portability : Non-portable (type operators)
@@ -34,6 +34,7 @@
        , pop
        , popTo
        , sortHeap
+       , heapInsert
        , Comparison
        ) where
 
@@ -61,8 +62,13 @@
 {-# 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
+  :: (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
@@ -74,22 +80,37 @@
 
 -- | 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
+  :: (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
+-- | 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
+  :: (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 -> Int -> Int -> m ()
+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 ()
@@ -105,21 +126,42 @@
 {-# 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 ()
+--
+-- 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.
-partialSortBy :: (PrimMonad m, MVector v e)
-              => Comparison e -> v (PrimState m) e -> Int -> m ()
+--
+-- 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.
-partialSortByBounds :: (PrimMonad m, MVector v e)
-                    => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()
+--
+-- 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
@@ -138,9 +180,18 @@
  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 ()
+-- | 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
@@ -152,15 +203,26 @@
 
 -- | 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
+  :: (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 -> Int -> Int -> m ()
+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
@@ -170,14 +232,44 @@
 -- | 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
+  :: (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 = l + shiftR (k-1) 2
+                  in unsafeRead v 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
diff --git a/src/Data/Vector/Algorithms/Intro.hs b/src/Data/Vector/Algorithms/Intro.hs
--- a/src/Data/Vector/Algorithms/Intro.hs
+++ b/src/Data/Vector/Algorithms/Intro.hs
@@ -6,7 +6,7 @@
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Intro
--- Copyright   : (c) 2008-2011 Dan Doel
+-- Copyright   : (c) 2008-2015 Dan Doel
 -- Maintainer  : Dan Doel <dan.doel@gmail.com>
 -- Stability   : Experimental
 -- Portability : Non-portable (type operators, bang patterns)
@@ -73,8 +73,13 @@
 {-# 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
+  :: (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
@@ -106,21 +111,35 @@
 
 -- | 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
+  :: (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 -> m ()
+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 -> Int -> Int -> m ()
+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
@@ -140,21 +159,35 @@
 {-# 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
+  :: (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 -> m ()
+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 -> Int -> Int -> m ()
+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 = go (ilg len) l (l + k) u
@@ -181,8 +214,6 @@
             => 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
diff --git a/src/Data/Vector/Algorithms/Search.hs b/src/Data/Vector/Algorithms/Search.hs
--- a/src/Data/Vector/Algorithms/Search.hs
+++ b/src/Data/Vector/Algorithms/Search.hs
@@ -4,7 +4,7 @@
 -- ---------------------------------------------------------------------------
 -- |
 -- Module      : Data.Vector.Algorithms.Search
--- Copyright   : (c) 2009-2010 Dan Doel
+-- Copyright   : (c) 2009-2015 Dan Doel, 2015 Tim Baumann
 -- Maintainer  : Dan Doel <dan.doel@gmail.com>
 -- Stability   : Experimental
 -- Portability : Non-portable (bang patterns)
@@ -24,6 +24,10 @@
        , binarySearchRByBounds
        , binarySearchP
        , binarySearchPBounds
+       , gallopingSearchLeftP
+       , gallopingSearchLeftPBounds
+       , gallopingSearchRightP
+       , gallopingSearchRightPBounds
        , Comparison
        ) where
 
@@ -134,3 +138,72 @@
    | otherwise = unsafeRead vec k >>= \e -> if p e then loop l k else loop (k+1) u
   where k = (u + l) `shiftR` 1
 {-# INLINE binarySearchPBounds #-}
+
+-- | 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,352 @@
+{-# 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
+       , sortBy
+       ) 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)
+
+-- | 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
+  | 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 #-}
+
+-- | 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
+ iter _ i _ _ _ _ _ _ | i >= tmpBufLen = return ()
+ iter tmpBuf i j k _ _ _ _ | j >= u = do
+   let from = unsafeSlice i (tmpBufLen-i) tmpBuf
+       to   = unsafeSlice k (tmpBufLen-i) vec
+   unsafeCopy to from
+ 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
+   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
+   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
+                      vi' <- unsafeRead tmpBuf (i+1)
+                      iter tmpBuf (i+1) j (k+1) vi' vj (ga-1) minGallop
+   | otherwise   = do unsafeWrite vec k vj
+                      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
+ iter _ _ j _ _ _ _ _ | j < 0 = return ()
+ iter tmpBuf i j _ _ _ _ _ | i < l = do
+   let from = unsafeSlice 0 (j+1) tmpBuf
+       to   = unsafeSlice l (j+1) vec
+   unsafeCopy to from
+ 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
+   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
+   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
+                     vi' <- unsafeRead vec (i-1)
+                     iter tmpBuf (i-1) j (k-1) vi' vj (ga-1) minGallop
+   | otherwise  = do unsafeWrite vec k vj
+                     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/properties/Tests.hs b/tests/properties/Tests.hs
--- a/tests/properties/Tests.hs
+++ b/tests/properties/Tests.hs
@@ -29,6 +29,7 @@
 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
 
@@ -49,6 +50,7 @@
          , ("insertion sort", INS.sort)
          , ("merge sort", M.sort)
          , ("heapsort", H.sort)
+         , ("timsort", T.sort)
          ]
 
 check_Int_partialsort = forM_ algos $ \(name,algo) ->
@@ -104,7 +106,9 @@
 
 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
@@ -123,6 +127,7 @@
   qc $ label "heapselect"   . prop_sized (const . prop_permutation)
                                          (H.select :: SizeAlgo 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  ())
@@ -155,6 +160,7 @@
   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
diff --git a/vector-algorithms.cabal b/vector-algorithms.cabal
--- a/vector-algorithms.cabal
+++ b/vector-algorithms.cabal
@@ -1,10 +1,11 @@
 name:              vector-algorithms
-version:           0.6.0.4
+version:           0.7
 license:           BSD3
 license-file:      LICENSE
 author:            Dan Doel
 maintainer:        Dan Doel <dan.doel@gmail.com>
-copyright:         (c) 2008,2009,2010,2011,2012,2013,2014 Dan Doel
+copyright:         (c) 2008,2009,2010,2011,2012,2013,2014,2015 Dan Doel
+                   (c) 2015 Tim Baumann
 homepage:          http://code.haskell.org/~dolio/
 category:          Data
 synopsis:          Efficient algorithms for vector arrays
@@ -56,6 +57,7 @@
     Data.Vector.Algorithms.Search
     Data.Vector.Algorithms.Heap
     Data.Vector.Algorithms.AmericanFlag
+    Data.Vector.Algorithms.Tim
 
   other-modules:
     Data.Vector.Algorithms.Common