diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Andrew Martin (c) 2017
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Andrew Martin nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,1 @@
+# mergesort
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/bench/Main.hs b/bench/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/Main.hs
@@ -0,0 +1,149 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE MagicHash #-}
+{-# LANGUAGE TypeApplications #-}
+
+import Gauge.Main 
+import Type.Reflection (typeRep,TypeRep)
+import Data.Primitive (ByteArray(..),PrimArray(..),Prim)
+import Control.Monad.ST (ST,runST)
+import Data.Int
+import Data.Word
+import System.Random (mkStdGen,randoms,Random)
+import GHC.Prim (proxy#, Proxy#)
+import qualified GHC.OldList as L
+import qualified Data.Primitive as P
+import qualified Data.Primitive.Sort
+import qualified GHC.Exts as E
+
+main :: IO ()
+main = defaultMain
+  [ bgroup "contiguous"
+    [ benchType (typeRep :: TypeRep Int8) (primArrayToByteArray . Data.Primitive.Sort.sort @PrimArray @Int8 . byteArrayToPrimArray)
+    , benchType (typeRep :: TypeRep Word) (primArrayToByteArray . Data.Primitive.Sort.sort @PrimArray @Word . byteArrayToPrimArray)
+    ]
+  , bgroup "tagged-unique"
+    [ bench "mini" (whnf (\(k,v) -> evalPair (Data.Primitive.Sort.sortUniqueTagged k v)) (sizedInts Mini, sizedInts Mini))
+    , bench "tiny" (whnf (\(k,v) -> evalPair (Data.Primitive.Sort.sortUniqueTagged k v)) (sizedInts Tiny, sizedInts Tiny))
+    , bench "small" (whnf (\(k,v) -> evalPair (Data.Primitive.Sort.sortUniqueTagged k v)) (sizedInts Small, sizedInts Small))
+    ]
+  ]
+
+evalPair :: (PrimArray a, PrimArray b) -> ()
+evalPair (!_,!_) = ()
+
+primArrayToByteArray :: PrimArray a -> ByteArray
+primArrayToByteArray (PrimArray x) = ByteArray x
+
+byteArrayToPrimArray :: ByteArray -> PrimArray a
+byteArrayToPrimArray (ByteArray x) = PrimArray x
+
+data Size = Mini | Tiny | Small | Medium | Large | Gigantic
+  deriving (Enum,Bounded)
+
+data Arrangement = Unsorted | Presorted | Reversed
+  deriving (Enum,Bounded)
+
+allSizes :: [Size]
+allSizes = [minBound..maxBound]
+
+allArrangements :: [Arrangement]
+allArrangements = [minBound..maxBound]
+
+showSize :: Size -> String
+showSize x = case x of
+  Mini -> "mini"
+  Tiny -> "tiny"
+  Small -> "small"
+  Medium -> "medium"
+  Large -> "large"
+  Gigantic -> "gigantic"
+
+numSize :: Size -> Int
+numSize x = case x of
+  Mini -> 10
+  Tiny -> 100
+  Small -> 1000
+  Medium -> 10000
+  Large -> 100000
+  Gigantic -> 1000000
+
+sizedInts :: Size -> PrimArray Int
+sizedInts x = case x of
+  Mini -> intsMini
+  Tiny -> intsTiny
+  Small -> intsSmall
+  Medium -> intsMedium
+  Large -> intsLarge
+  Gigantic -> intsGigantic
+
+intsMini,intsTiny,intsSmall,intsMedium,intsLarge,intsGigantic :: PrimArray Int
+intsMini = E.fromList (L.take 10 (randoms (mkStdGen 23) :: [Int]))
+intsTiny = E.fromList (L.take 100 (randoms (mkStdGen 87) :: [Int]))
+intsSmall = E.fromList (L.take 1000 (randoms (mkStdGen 19) :: [Int]))
+intsMedium = E.fromList (L.take 10000 (randoms (mkStdGen 47) :: [Int]))
+intsLarge = E.fromList (L.take 100000 (randoms (mkStdGen 53) :: [Int]))
+intsGigantic = E.fromList (L.take 1000000 (randoms (mkStdGen 12) :: [Int]))
+
+showArrangement :: Arrangement -> String
+showArrangement x = case x of
+  Unsorted -> "unsorted"
+  Presorted -> "presorted"
+  Reversed -> "reversed"
+
+buildArrangement :: (Prim a, Num a, Random a, Enum a, Bounded a)
+  => Arrangement -> TypeRep a -> Int -> ByteArray
+buildArrangement x = case x of
+  Unsorted -> unsorted
+  Presorted -> presorted
+  Reversed -> reversed
+
+benchType :: (Prim a, Num a, Random a, Enum a, Bounded a)
+  => TypeRep a -> (ByteArray -> ByteArray) -> Benchmark
+benchType rep sort = bgroup
+  (show rep)
+  (map (\arrange -> benchArrangement rep arrange sort) allArrangements)
+
+benchArrangement :: (Prim a, Num a, Random a, Enum a, Bounded a)
+  => TypeRep a -> Arrangement -> (ByteArray -> ByteArray) -> Benchmark
+benchArrangement rep arrange sort = bgroup
+  (showArrangement arrange)
+  (map (\sz -> let arr = buildArrangement arrange rep (numSize sz) in benchSize arr sz sort) allSizes)
+
+benchSize :: ByteArray -> Size -> (ByteArray -> ByteArray) -> Benchmark
+benchSize arr sz sort =
+  bench (showSize sz) (whnf sort arr)
+
+unsorted :: forall a. (Prim a, Random a) => TypeRep a -> Int -> ByteArray
+unsorted typ n = byteArrayFromList
+  (L.take n (randoms (mkStdGen 42) :: [a]))
+
+presorted :: forall a. (Prim a, Num a, Enum a, Bounded a) => TypeRep a -> Int -> ByteArray
+presorted typ n = byteArrayFromList
+  (L.take n (iterate (+1) (minBound :: a)))
+
+reversed :: forall a. (Prim a, Num a, Enum a, Bounded a)
+  => TypeRep a -> Int -> ByteArray
+reversed typ n = byteArrayFromList
+  (L.take n (iterate (subtract 1) (maxBound :: a)))
+
+
+
+byteArrayFromList :: Prim a => [a] -> ByteArray
+byteArrayFromList xs = byteArrayFromListN (L.length xs) xs
+
+byteArrayFromListN :: forall a. Prim a => Int -> [a] -> ByteArray
+byteArrayFromListN len vs = runST run where
+  run :: forall s. ST s ByteArray
+  run = do
+    arr <- P.newByteArray (len * P.sizeOf (undefined :: a))
+    let go :: [a] -> Int -> ST s ()
+        go !xs !ix = case xs of
+          [] -> return ()
+          a : as -> do
+            P.writeByteArray arr ix a
+            go as (ix + 1)
+    go vs 0
+    P.unsafeFreezeByteArray arr
+
diff --git a/primitive-sort.cabal b/primitive-sort.cabal
new file mode 100644
--- /dev/null
+++ b/primitive-sort.cabal
@@ -0,0 +1,81 @@
+cabal-version: 2.0
+name: primitive-sort
+version: 0.1.0.0
+synopsis: Sort primitive arrays
+description:
+  This library provides a stable sorting algorithm for primitive arrays.
+  When extra capabilities are available, the sort is parallelized.
+  .
+  The algorithm currently uses mergesort on large chunks and switches
+  to insertion sort on small chunks. The are also novel improvements
+  to increase the performance if the input array is already mostly sorted.
+homepage: https://github.com/andrewthad/primitive-sort
+license: BSD3
+license-file: LICENSE
+author: Andrew Martin
+maintainer: andrew.thaddeus@gmail.com
+copyright: 2018 Andrew Martin
+category: software
+build-type: Simple
+extra-source-files: README.md
+
+library
+  hs-source-dirs: src
+  exposed-modules:
+    Data.Primitive.Sort
+  build-depends:
+      base >= 0.4.9 && < 5
+    , primitive >= 0.6.4.0
+    , ghc-prim
+    , contiguous >= 0.1 && < 0.2
+  ghc-options: -O2
+  default-language: Haskell2010
+
+test-suite test
+  type: exitcode-stdio-1.0
+  hs-source-dirs: test
+  main-is: Main.hs
+  build-depends:
+      base
+    , primitive-sort
+    , tasty
+    , tasty-hunit
+    , tasty-smallcheck
+    , tasty-quickcheck
+    , primitive
+    , containers
+    , smallcheck
+    , QuickCheck
+    , HUnit
+  ghc-options: -threaded -rtsopts -O2 -with-rtsopts=-N
+  default-language: Haskell2010
+
+test-suite doctest
+  type: exitcode-stdio-1.0
+  hs-source-dirs: test
+  main-is: Doctest.hs
+  build-depends:
+      base
+    , primitive-sort
+    , doctest >= 0.10
+    , QuickCheck
+  default-language: Haskell2010
+
+benchmark bench
+  type: exitcode-stdio-1.0
+  build-depends:
+      base
+    , primitive-sort
+    , gauge
+    , random
+    , primitive
+    , ghc-prim
+  ghc-options: -threaded -rtsopts -O2 -with-rtsopts=-N
+  default-language: Haskell2010
+  hs-source-dirs: bench
+  main-is: Main.hs
+
+source-repository head
+  type: git
+  location: https://github.com/andrewthad/primitive-sort
+
diff --git a/src/Data/Primitive/Sort.hs b/src/Data/Primitive/Sort.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Primitive/Sort.hs
@@ -0,0 +1,796 @@
+{-# LANGUAGE MultiWayIf #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE MagicHash #-}
+{-# LANGUAGE UnboxedTuples #-}
+
+{-# OPTIONS_GHC -Wall #-}
+
+-- | Sort primitive arrays with a stable sorting algorithm. All functions
+-- in this module are marked as @INLINABLE@, so they will specialize
+-- when used in a monomorphic setting.
+module Data.Primitive.Sort
+  ( -- * Immutable
+    sort
+  , sortUnique
+  , sortTagged
+  , sortUniqueTagged
+    -- * Mutable
+  , sortMutable
+  , sortUniqueMutable
+  , sortTaggedMutable
+  , sortUniqueTaggedMutable
+  ) where
+
+import Control.Monad.ST
+import Control.Applicative
+import GHC.ST (ST(..))
+import GHC.IO (IO(..))
+import GHC.Int (Int(..))
+import Control.Monad
+import GHC.Prim
+import Control.Concurrent (getNumCapabilities)
+import Data.Primitive.Contiguous (Contiguous,Mutable,Element)
+import qualified Data.Primitive.Contiguous as C
+
+-- | Sort an immutable array. Duplicate elements are preserved.
+--
+-- >>> sort ([5,6,7,9,5,4,5,7] :: Array Int)
+-- fromListN 8 [4,5,5,5,6,7,7,9]
+sort :: (Contiguous arr, Element arr a, Ord a)
+  => arr a
+  -> arr a
+{-# INLINABLE sort #-}
+sort !src = runST $ do
+  let len = C.size src
+  dst <- C.new (C.size src)
+  C.copy dst 0 src 0 len
+  res <- sortMutable dst
+  C.unsafeFreeze res
+
+-- | Sort a tagged immutable array. Each element from the @keys@ array is
+-- paired up with an element from the @values@ array at the matching
+-- index. The sort permutes the @values@ array so that a value end up
+-- in the same position as its corresponding key. The two argument array
+-- should be of the same length, but if one is shorter than the other,
+-- the longer one will be truncated so that the lengths match.
+--
+-- >>> sortTagged ([5,6,7,5,5,7] :: Array Int) ([1,2,3,4,5,6] :: Array Int)
+-- (fromListN 6 [5,5,5,6,7,7],fromListN 6 [1,4,5,2,3,6])
+--
+-- Since the sort is stable, the values corresponding to a key that
+-- appears multiple times have their original order preserved.
+sortTagged :: forall k v karr varr. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => karr k -- ^ keys
+  -> varr v -- ^ values
+  -> (karr k,varr v)
+{-# INLINABLE sortTagged #-}
+sortTagged !src !srcTags = runST $ do
+  let len = min (C.size src) (C.size srcTags)
+  dst <- C.new len
+  C.copy dst 0 src 0 len
+  dstTags <- C.new len
+  C.copy dstTags 0 srcTags 0 len
+  (res,resTags) <- sortTaggedMutableN len dst dstTags
+  res' <- C.unsafeFreeze res
+  resTags' <- C.unsafeFreeze resTags
+  return (res',resTags')
+
+-- | Sort a tagged immutable array. Only a single copy of each
+-- duplicate key is preserved, along with the last value from @values@
+-- that corresponded to it. The two argument arrays
+-- should be of the same length, but if one is shorter than the other,
+-- the longer one will be truncated so that the lengths match.
+--
+-- >>> sortUniqueTagged ([5,6,7,5,5,7] :: Array Int) ([1,2,3,4,5,6] :: Array Int)
+-- (fromListN 3 [5,6,7],fromListN 3 [5,2,6])
+sortUniqueTagged :: forall k v karr varr. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => karr k -- ^ keys
+  -> varr v -- ^ values
+  -> (karr k,varr v)
+{-# INLINABLE sortUniqueTagged #-}
+sortUniqueTagged !src !srcTags = runST $ do
+  let len = min (C.size src) (C.size srcTags)
+  dst <- C.new len
+  C.copy dst 0 src 0 len
+  dstTags <- C.new len
+  C.copy dstTags 0 srcTags 0 len
+  (res0,resTags0) <- sortTaggedMutableN len dst dstTags
+  (res1,resTags1) <- uniqueTaggedMutableN len res0 resTags0
+  res' <- C.unsafeFreeze res1
+  resTags' <- C.unsafeFreeze resTags1
+  return (res',resTags')
+
+-- | Sort the mutable array. This operation preserves duplicate
+-- elements. The argument may either be modified in-place, or another
+-- array may be allocated and returned. The argument
+-- may not be reused after being passed to this function.
+sortMutable :: (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a
+  -> ST s (Mutable arr s a)
+{-# INLINABLE sortMutable #-}
+sortMutable !dst = do
+  len <- C.sizeMutable dst
+  if len < threshold
+    then insertionSortRange dst 0 len
+    else do
+      work <- C.new len
+      C.copyMutable work 0 dst 0 len 
+      caps <- unsafeEmbedIO getNumCapabilities
+      let minElemsPerThread = 20000
+          maxThreads = unsafeQuot len minElemsPerThread
+          preThreads = min caps maxThreads
+          threads = if preThreads == 1 then 1 else preThreads * 8
+      -- I cannot understand why, but GHC's runtime does better
+      -- when we let this schedule 8 times as many threads as
+      -- we have capabilities. However, we only get this benefit
+      -- when we actually have more than one capability.
+      splitMergeParallel dst work threads 0 len
+  return dst
+
+-- | Sort an array of a key type @k@, rearranging the values of
+-- type @v@ according to the element they correspond to in the
+-- key array. The argument arrays may not be reused after they
+-- are passed to the function.
+sortTaggedMutable :: (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k
+  -> Mutable varr s v
+  -> ST s (Mutable karr s k, Mutable varr s v)
+{-# INLINABLE sortTaggedMutable #-}
+sortTaggedMutable !dst0 !dstTags0 = do
+  (!dst,!dstTags,!len) <- alignArrays dst0 dstTags0
+  sortTaggedMutableN len dst dstTags
+
+alignArrays :: (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k
+  -> Mutable varr s v
+  -> ST s (Mutable karr s k, Mutable varr s v,Int)
+{-# INLINABLE alignArrays #-}
+alignArrays dst0 dstTags0 = do
+  lenDst <- C.sizeMutable dst0
+  lenDstTags <- C.sizeMutable dstTags0
+  -- This cleans up mismatched lengths.
+  if lenDst == lenDstTags
+    then return (dst0,dstTags0,lenDst)
+    else if lenDst < lenDstTags
+      then do
+        dstTags <- C.resize dstTags0 lenDst
+        return (dst0,dstTags,lenDst)
+      else do
+        dst <- C.resize dst0 lenDstTags
+        return (dst,dstTags0,lenDstTags)
+
+sortUniqueTaggedMutable :: (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k -- ^ keys
+  -> Mutable varr s v -- ^ values
+  -> ST s (Mutable karr s k, Mutable varr s v)
+{-# INLINABLE sortUniqueTaggedMutable #-}
+sortUniqueTaggedMutable dst0 dstTags0 = do
+  (!dst1,!dstTags1,!len) <- alignArrays dst0 dstTags0
+  (!dst2,!dstTags2) <- sortTaggedMutableN len dst1 dstTags1
+  uniqueTaggedMutableN len dst2 dstTags2
+
+sortTaggedMutableN :: (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Int
+  -> Mutable karr s k
+  -> Mutable varr s v
+  -> ST s (Mutable karr s k, Mutable varr s v)
+{-# INLINABLE sortTaggedMutableN #-}
+sortTaggedMutableN !len !dst !dstTags = if len < thresholdTagged
+  then do
+    insertionSortTaggedRange dst dstTags 0 len
+    return (dst,dstTags)
+  else do
+    work <- C.cloneMutable dst 0 len 
+    workTags <- C.cloneMutable dstTags 0 len 
+    caps <- unsafeEmbedIO getNumCapabilities
+    let minElemsPerThread = 20000
+        maxThreads = unsafeQuot len minElemsPerThread
+        preThreads = min caps maxThreads
+        threads = if preThreads == 1 then 1 else preThreads * 8
+    splitMergeParallelTagged dst work dstTags workTags threads 0 len
+    return (dst,dstTags)
+
+-- | Sort an immutable array. Only a single copy of each duplicated
+-- element is preserved.
+--
+-- >>> sortUnique ([5,6,7,9,5,4,5,7] :: Array Int)
+-- fromListN 5 [4,5,6,7,9]
+sortUnique :: (Contiguous arr, Element arr a, Ord a)
+  => arr a -> arr a
+{-# INLINABLE sortUnique #-}
+sortUnique src = runST $ do
+  let len = C.size src
+  dst <- C.new len
+  C.copy dst 0 src 0 len
+  res <- sortUniqueMutable dst
+  C.unsafeFreeze res
+
+-- | Sort an immutable array. Only a single copy of each duplicated
+-- element is preserved. This operation may run in-place, or it may
+-- need to allocate a new array, so the argument may not be reused
+-- after this function is applied to it. 
+sortUniqueMutable :: (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a
+  -> ST s (Mutable arr s a)
+{-# INLINABLE sortUniqueMutable #-}
+sortUniqueMutable marr = do
+  res <- sortMutable marr
+  uniqueMutable res
+
+-- | Discards adjacent equal elements from an array. This operation
+-- may run in-place, or it may need to allocate a new array, so the
+-- argument may not be reused after this function is applied to it.
+uniqueMutable :: forall arr s a. (Contiguous arr, Element arr a, Eq a)
+  => Mutable arr s a -> ST s (Mutable arr s a)
+{-# INLINABLE uniqueMutable #-}
+uniqueMutable !marr = do
+  !len <- C.sizeMutable marr
+  if len > 1
+    then do
+      !a0 <- C.read marr 0
+      let findFirstDuplicate :: a -> Int -> ST s Int
+          findFirstDuplicate !prev !ix = if ix < len
+            then do
+              a <- C.read marr ix
+              if a == prev
+                then return ix
+                else findFirstDuplicate a (ix + 1)
+            else return ix
+      dupIx <- findFirstDuplicate a0 1
+      if dupIx == len
+        then return marr
+        else do
+          let deduplicate :: a -> Int -> Int -> ST s Int
+              deduplicate !prev !srcIx !dstIx = if srcIx < len
+                then do
+                  a <- C.read marr srcIx
+                  if a == prev
+                    then deduplicate a (srcIx + 1) dstIx
+                    else do
+                      C.write marr dstIx a
+                      deduplicate a (srcIx + 1) (dstIx + 1)
+                else return dstIx
+          !a <- C.read marr dupIx
+          !reducedLen <- deduplicate a (dupIx + 1) dupIx
+          C.resize marr reducedLen
+    else return marr
+
+uniqueTaggedMutableN :: forall karr varr s k v. (Contiguous karr, Element karr k, Eq k, Contiguous varr, Element varr v)
+  => Int
+  -> Mutable karr s k
+  -> Mutable varr s v
+  -> ST s (Mutable karr s k, Mutable varr s v)
+{-# INLINABLE uniqueTaggedMutableN #-}
+uniqueTaggedMutableN !len !marr !marrTags = if len > 1
+  then do
+    !a0 <- C.read marr 0
+    let findFirstDuplicate :: k -> Int -> ST s Int
+        findFirstDuplicate !prev !ix = if ix < len
+          then do
+            a <- C.read marr ix
+            if a == prev
+              then return ix
+              else findFirstDuplicate a (ix + 1)
+          else return ix
+    dupIx <- findFirstDuplicate a0 1
+    if dupIx == len
+      then return (marr,marrTags)
+      else do
+        C.read marrTags dupIx >>= C.write marrTags (dupIx - 1)
+        let deduplicate :: k -> Int -> Int -> ST s Int
+            deduplicate !prev !srcIx !dstIx = if srcIx < len
+              then do
+                a <- C.read marr srcIx
+                if a == prev
+                  then do
+                    C.read marrTags srcIx >>= C.write marrTags (dstIx - 1)
+                    deduplicate a (srcIx + 1) dstIx
+                  else do
+                    C.read marrTags srcIx >>= C.write marrTags dstIx
+                    C.write marr dstIx a
+                    deduplicate a (srcIx + 1) (dstIx + 1)
+              else return dstIx
+        !a <- C.read marr dupIx
+        !reducedLen <- deduplicate a (dupIx + 1) dupIx
+        liftA2 (,) (C.resize marr reducedLen) (C.resize marrTags reducedLen)
+  else return (marr,marrTags)
+
+unsafeEmbedIO :: IO a -> ST s a
+unsafeEmbedIO (IO f) = ST (unsafeCoerce# f)
+
+half :: Int -> Int
+half x = unsafeQuot x 2
+
+splitMergeParallel :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a -- source and destination
+  -> Mutable arr s a -- work array
+  -> Int -- spark limit, should be power of two
+  -> Int -- start
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE splitMergeParallel #-}
+splitMergeParallel !arr !work !level !start !end = if level > 1
+  then if end - start < threshold
+    then insertionSortRange arr start end
+    else do
+      let !mid = unsafeQuot (end + start) 2
+          !levelDown = half level
+      tandem 
+        (splitMergeParallel work arr levelDown start mid)
+        (splitMergeParallel work arr levelDown mid end)
+      mergeParallel work arr level start mid end
+  else splitMerge arr work start end
+
+splitMergeParallelTagged :: forall karr varr s k v. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k -- source and destination
+  -> Mutable karr s k -- work array
+  -> Mutable varr s v -- source and destination tags
+  -> Mutable varr s v -- work tags
+  -> Int -- spark limit, should be power of two
+  -> Int -- start
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE splitMergeParallelTagged #-}
+splitMergeParallelTagged !arr !work !arrTags !workTags !level !start !end = if level > 1
+  then do
+    let !mid = unsafeQuot (end + start) 2
+        !levelDown = half level
+    tandem 
+      (splitMergeParallelTagged work arr workTags arrTags levelDown start mid)
+      (splitMergeParallelTagged work arr workTags arrTags levelDown mid end)
+    mergeParallelTagged work arr workTags arrTags level start mid end
+  else splitMergeTagged arr work arrTags workTags start end
+
+splitMerge :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a -- source and destination
+  -> Mutable arr s a -- work array
+  -> Int -- start
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE splitMerge #-}
+splitMerge !arr !work !start !end = if end - start < 2
+  then return ()
+  else if end - start > threshold
+    then do
+      let !mid = unsafeQuot (end + start) 2
+      splitMerge work arr start mid
+      splitMerge work arr mid end
+      mergeNonContiguous work arr start mid mid end start
+    else insertionSortRange arr start end
+
+splitMergeTagged :: (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k -- source and destination
+  -> Mutable karr s k -- work array
+  -> Mutable varr s v
+  -> Mutable varr s v
+  -> Int -- start
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE splitMergeTagged #-}
+splitMergeTagged !arr !work !arrTags !workTags !start !end = if end - start < 2
+  then return ()
+  else if end - start > thresholdTagged
+    then do
+      let !mid = unsafeQuot (end + start) 2
+      splitMergeTagged work arr workTags arrTags start mid
+      splitMergeTagged work arr workTags arrTags mid end
+      mergeNonContiguousTagged work arr workTags arrTags start mid mid end start
+    else insertionSortTaggedRange arr arrTags start end
+
+-- Precondition: threads is greater than 0
+mergeParallel :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a -- source
+  -> Mutable arr s a -- dest
+  -> Int -- threads
+  -> Int -- start
+  -> Int -- middle
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE mergeParallel #-}
+mergeParallel !src !dst !threads !start !mid !end = do
+  !lock <- newLock
+  let go :: Int -- previous A end
+         -> Int -- previous B end
+         -> Int -- how many chunk have we already iterated over
+         -> ST s Int
+      go !prevEndA !prevEndB !ix = 
+        if | prevEndA == mid && prevEndB == end -> return ix
+           | prevEndA == mid -> do
+               forkST_ $ do
+                 let !startA = mid
+                     !endA = mid
+                     !startB = prevEndB
+                     !endB = end
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguous src dst startA endA startB endB startDst
+                 putLock lock
+               go mid end (ix + 1)
+           | prevEndB == end -> do
+               forkST_ $ do
+                 let !startA = prevEndA
+                     !endA = mid
+                     !startB = end
+                     !endB = end
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguous src dst startA endA startB endB startDst
+                 putLock lock
+               go mid end (ix + 1)
+           | ix == threads - 1 -> do
+               forkST_ $ do
+                 let !startA = prevEndA
+                     !endA = mid
+                     !startB = prevEndB
+                     !endB = end
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguous src dst startA endA startB endB startDst
+                 putLock lock
+               return (ix + 1)
+           | otherwise -> do
+               -- We use the left half for this lookup. We could instead
+               -- use both halves and take the median.
+               !endElem <- C.read src (start + chunk * (ix + 1))
+               !endA <- findIndexOfGtElem src (endElem :: a) prevEndA mid
+               !endB <- findIndexOfGtElem src endElem prevEndB end
+               forkST_ $ do
+                 let !startA = prevEndA
+                     !startB = prevEndB
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguous src dst startA endA startB endB startDst
+                 putLock lock
+               go endA endB (ix + 1)
+  !endElem <- C.read src (start + chunk) 
+  !endA <- findIndexOfGtElem src (endElem :: a) start mid
+  !endB <- findIndexOfGtElem src endElem mid end
+  forkST_ $ do
+    let !startA = start
+        !startB = mid
+        !startDst = (startA - start) + (startB - mid) + start
+    mergeNonContiguous src dst startA endA startB endB startDst
+    putLock lock
+  total <- go endA endB 1
+  replicateM_ total (takeLock lock)
+  where
+  !chunk = unsafeQuot (end - start) threads
+
+-- Precondition: threads is greater than 0
+-- This function is just a copy of mergeParallel but with
+-- the tags arrays passed to mergeNonContiguousTagged
+mergeParallelTagged :: forall karr varr s k v. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k -- source
+  -> Mutable karr s k -- dest
+  -> Mutable varr s v -- source tags
+  -> Mutable varr s v -- dest tags
+  -> Int -- threads
+  -> Int -- start
+  -> Int -- middle
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE mergeParallelTagged #-}
+mergeParallelTagged !src !dst !srcTags !dstTags !threads !start !mid !end = do
+  !lock <- newLock
+  let go :: Int -- previous A end
+         -> Int -- previous B end
+         -> Int -- how many chunk have we already iterated over
+         -> ST s Int
+      go !prevEndA !prevEndB !ix = 
+        if | prevEndA == mid && prevEndB == end -> return ix
+           | prevEndA == mid -> do
+               forkST_ $ do
+                 let !startA = mid
+                     !endA = mid
+                     !startB = prevEndB
+                     !endB = end
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguousTagged src dst srcTags dstTags startA endA startB endB startDst
+                 putLock lock
+               go mid end (ix + 1)
+           | prevEndB == end -> do
+               forkST_ $ do
+                 let !startA = prevEndA
+                     !endA = mid
+                     !startB = end
+                     !endB = end
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguousTagged src dst srcTags dstTags startA endA startB endB startDst
+                 putLock lock
+               go mid end (ix + 1)
+           | ix == threads - 1 -> do
+               forkST_ $ do
+                 let !startA = prevEndA
+                     !endA = mid
+                     !startB = prevEndB
+                     !endB = end
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguousTagged src dst srcTags dstTags startA endA startB endB startDst
+                 putLock lock
+               return (ix + 1)
+           | otherwise -> do
+               -- We use the left half for this lookup. We could instead
+               -- use both halves and take the median.
+               !endElem <- C.read src (start + chunk * (ix + 1))
+               !endA <- findIndexOfGtElem src (endElem :: k) prevEndA mid
+               !endB <- findIndexOfGtElem src endElem prevEndB end
+               forkST_ $ do
+                 let !startA = prevEndA
+                     !startB = prevEndB
+                     !startDst = (startA - start) + (startB - mid) + start
+                 mergeNonContiguousTagged src dst srcTags dstTags startA endA startB endB startDst
+                 putLock lock
+               go endA endB (ix + 1)
+  !endElem <- C.read src (start + chunk) 
+  !endA <- findIndexOfGtElem src (endElem :: k) start mid
+  !endB <- findIndexOfGtElem src endElem mid end
+  forkST_ $ do
+    let !startA = start
+        !startB = mid
+        !startDst = (startA - start) + (startB - mid) + start
+    mergeNonContiguousTagged src dst srcTags dstTags startA endA startB endB startDst
+    putLock lock
+  total <- go endA endB 1
+  replicateM_ total (takeLock lock)
+  where
+  !chunk = unsafeQuot (end - start) threads
+
+unsafeQuot :: Int -> Int -> Int
+unsafeQuot (I# a) (I# b) = I# (quotInt# a b)
+
+-- If the needle is bigger than everything in the slice
+-- of the array, this returns the end index (which is out
+-- of bounds). Callers of this function should be able
+-- to handle that.
+findIndexOfGtElem :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a -> a -> Int -> Int -> ST s Int
+{-# INLINABLE findIndexOfGtElem #-}
+findIndexOfGtElem !v !needle !start !end = go start end
+  where
+  go :: Int -> Int -> ST s Int
+  go !lo !hi = if lo < hi
+    then do
+      let !mid = lo + half (hi - lo)
+      !val <- C.read v mid
+      if | val == needle -> gallopToGtIndex v needle (mid + 1) hi
+         | val < needle -> go (mid + 1) hi
+         | otherwise -> go lo mid
+    else return lo
+
+-- | TODO: should probably turn this into a real galloping search
+gallopToGtIndex :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a -> a -> Int -> Int -> ST s Int
+{-# INLINABLE gallopToGtIndex #-}
+gallopToGtIndex !v !val !start !end = go start
+  where
+  go :: Int -> ST s Int
+  go !ix = if ix < end
+    then do
+      !a <- C.read v ix
+      if a > val
+        then return ix
+        else go (ix + 1)
+    else return end
+
+-- stepA assumes that we previously incremented ixA.
+-- Consequently, we do not need to check that ixB
+-- is still in bounds. As a precondition, both
+-- indices are guarenteed to start in bounds.
+mergeNonContiguous :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a -- source
+  -> Mutable arr s a -- dest
+  -> Int -- start A
+  -> Int -- end A
+  -> Int -- start B
+  -> Int -- end B
+  -> Int -- start destination
+  -> ST s ()
+{-# INLINABLE mergeNonContiguous #-}
+mergeNonContiguous !src !dst !startA !endA !startB !endB !startDst =
+  if startB < endB
+    then stepA startA startB startDst
+    else if startA < endA
+      then stepB startA startB startDst
+      else return ()
+  where
+  continue :: Int -> Int -> Int -> ST s ()
+  continue ixA ixB ixDst = do
+    !a <- C.read src ixA
+    !b <- C.read src ixB
+    if (a :: a) <= b
+      then do
+        C.write dst ixDst a
+        stepA (ixA + 1) ixB (ixDst + 1)
+      else do
+        C.write dst ixDst b
+        stepB ixA (ixB + 1) (ixDst + 1)
+  stepB :: Int -> Int -> Int -> ST s ()
+  stepB !ixA !ixB !ixDst = if ixB < endB
+    then continue ixA ixB ixDst
+    else finishA ixA ixDst
+  stepA :: Int -> Int -> Int -> ST s ()
+  stepA !ixA !ixB !ixDst = if ixA < endA
+    then continue ixA ixB ixDst
+    else finishB ixB ixDst
+  finishB :: Int -> Int -> ST s ()
+  finishB !ixB !ixDst = C.copyMutable dst ixDst src ixB (endB - ixB)
+  finishA :: Int -> Int -> ST s ()
+  finishA !ixA !ixDst = C.copyMutable dst ixDst src ixA (endA - ixA)
+
+mergeNonContiguousTagged :: forall karr varr k v s. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k -- source
+  -> Mutable karr s k -- dest
+  -> Mutable varr s v -- source tags
+  -> Mutable varr s v -- dest tags
+  -> Int -- start A
+  -> Int -- end A
+  -> Int -- start B
+  -> Int -- end B
+  -> Int -- start destination
+  -> ST s ()
+{-# INLINABLE mergeNonContiguousTagged #-}
+mergeNonContiguousTagged !src !dst !srcTags !dstTags !startA !endA !startB !endB !startDst =
+  if startB < endB
+    then stepA startA startB startDst
+    else if startA < endA
+      then stepB startA startB startDst
+      else return ()
+  where
+  continue :: Int -> Int -> Int -> ST s ()
+  continue ixA ixB ixDst = do
+    !a <- C.read src ixA
+    !b <- C.read src ixB
+    if a <= b
+      then do
+        C.write dst ixDst a
+        (C.read srcTags ixA :: ST s v) >>= C.write dstTags ixDst
+        stepA (ixA + 1) ixB (ixDst + 1)
+      else do
+        C.write dst ixDst b
+        (C.read srcTags ixB :: ST s v) >>= C.write dstTags ixDst
+        stepB ixA (ixB + 1) (ixDst + 1)
+  stepB :: Int -> Int -> Int -> ST s ()
+  stepB !ixA !ixB !ixDst = if ixB < endB
+    then continue ixA ixB ixDst
+    else finishA ixA ixDst
+  stepA :: Int -> Int -> Int -> ST s ()
+  stepA !ixA !ixB !ixDst = if ixA < endA
+    then continue ixA ixB ixDst
+    else finishB ixB ixDst
+  finishB :: Int -> Int -> ST s ()
+  finishB !ixB !ixDst = do
+    C.copyMutable dst ixDst src ixB (endB - ixB)
+    C.copyMutable dstTags ixDst srcTags ixB (endB - ixB)
+  finishA :: Int -> Int -> ST s ()
+  finishA !ixA !ixDst = do
+    C.copyMutable dst ixDst src ixA (endA - ixA)
+    C.copyMutable dstTags ixDst srcTags ixA (endA - ixA)
+
+threshold :: Int
+threshold = 16
+
+thresholdTagged :: Int
+thresholdTagged = 16
+
+insertionSortRange :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a
+  -> Int -- start
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE insertionSortRange #-}
+insertionSortRange !arr !start !end = go start
+  where
+  go :: Int -> ST s ()
+  go !ix = if ix < end
+    then do
+      !a <- C.read arr ix
+      insertElement arr (a :: a) start ix
+      go (ix + 1)
+    else return ()
+    
+insertElement :: forall arr s a. (Contiguous arr, Element arr a, Ord a)
+  => Mutable arr s a
+  -> a
+  -> Int
+  -> Int
+  -> ST s ()
+{-# INLINABLE insertElement #-}
+insertElement !arr !a !start !end = go end
+  where
+  go :: Int -> ST s ()
+  go !ix = if ix > start
+    then do
+      !b <- C.read arr (ix - 1)
+      if b <= a
+        then do
+          C.copyMutable arr (ix + 1) arr ix (end - ix)
+          C.write arr ix a
+        else go (ix - 1)
+    else do
+      C.copyMutable arr (ix + 1) arr ix (end - ix)
+      C.write arr ix a
+
+insertionSortTaggedRange :: forall karr varr s k v. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k
+  -> Mutable varr s v
+  -> Int -- start
+  -> Int -- end
+  -> ST s ()
+{-# INLINABLE insertionSortTaggedRange #-}
+insertionSortTaggedRange !karr !varr !start !end = go start
+  where
+  go :: Int -> ST s ()
+  go !ix = if ix < end
+    then do
+      !a <- C.read karr ix
+      !v <- C.read varr ix
+      insertElementTagged karr varr a v start ix
+      go (ix + 1)
+    else return ()
+    
+insertElementTagged :: forall karr varr s k v. (Contiguous karr, Element karr k, Ord k, Contiguous varr, Element varr v)
+  => Mutable karr s k
+  -> Mutable varr s v
+  -> k
+  -> v
+  -> Int
+  -> Int
+  -> ST s ()
+{-# INLINABLE insertElementTagged #-}
+insertElementTagged !karr !varr !a !v !start !end = go end
+  where
+  go :: Int -> ST s ()
+  go !ix = if ix > start
+    then do
+      !b <- C.read karr (ix - 1)
+      if b <= a
+        then do
+          C.copyMutable karr (ix + 1) karr ix (end - ix)
+          C.write karr ix a
+          C.copyMutable varr (ix + 1) varr ix (end - ix)
+          C.write varr ix v
+        else go (ix - 1)
+    else do
+      C.copyMutable karr (ix + 1) karr ix (end - ix)
+      C.write karr ix a
+      C.copyMutable varr (ix + 1) varr ix (end - ix)
+      C.write varr ix v
+
+
+forkST_ :: ST s a -> ST s ()
+forkST_ action = ST $ \s1 -> case forkST# action s1 of
+  (# s2, _ #) -> (# s2, () #)
+
+forkST# :: a -> State# s -> (# State# s, ThreadId# #)
+forkST# = unsafeCoerce# fork#
+
+data Lock s = Lock (MVar# s ())
+
+newLock :: ST s (Lock s)
+newLock = ST $ \s1 -> case newMVar# s1 of
+  (# s2, v #) -> (# s2, Lock v #)
+
+takeLock :: Lock s -> ST s ()
+takeLock (Lock mvar#) = ST $ \ s# -> takeMVar# mvar# s#
+
+putLock  :: Lock s -> ST s ()
+putLock (Lock mvar#) = ST $ \ s# ->
+  case putMVar# mvar# () s# of
+    s2# -> (# s2#, () #)
+
+-- | Execute the first computation on the main thread and
+--   the second one on another thread in parallel. Blocks
+--   until both are finished.
+tandem :: ST s () -> ST s () -> ST s ()
+tandem a b = do
+  lock <- newLock
+  forkST_ (b >> putLock lock)
+  a
+  takeLock lock
+
+-- $setup
+--
+-- These are to make doctest work correctly.
+--
+-- >>> :set -XOverloadedLists
+-- >>> import Data.Primitive.Array (Array)
+--
+
diff --git a/test/Doctest.hs b/test/Doctest.hs
new file mode 100644
--- /dev/null
+++ b/test/Doctest.hs
@@ -0,0 +1,6 @@
+import Test.DocTest
+
+main :: IO ()
+main = doctest
+  [ "src/Data/Primitive/Sort.hs"
+  ]
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,299 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeApplications #-}
+
+{-# OPTIONS_GHC -Wall -fno-warn-orphans #-}
+
+import Test.Tasty
+import Test.Tasty.HUnit (testCase)
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC
+import Test.QuickCheck as Q
+import qualified Test.QuickCheck.Property as QP
+import Type.Reflection (TypeRep,typeRep)
+import qualified Test.SmallCheck.Series as SCS
+import Test.HUnit.Base ((@?=))
+
+import Data.List
+import Data.Word
+import Data.Int
+import Data.Primitive (ByteArray,Prim)
+import Data.Proxy (Proxy(..))
+import Control.Monad.ST (ST,runST)
+import Test.SmallCheck.Series (Serial(..),Series)
+import Control.Exception (Exception,toException)
+import Control.Applicative (liftA2)
+import Data.Primitive (ByteArray(..),PrimArray(..),Prim,Array)
+
+import qualified GHC.Exts as E
+import qualified GHC.OldList as L
+import qualified Data.Set as S
+import qualified Data.Map as M
+import qualified Data.Primitive as P
+import qualified Data.Primitive.Sort
+-- import qualified Sort.Merge.Int8
+-- import qualified Sort.Merge.Word16
+-- import qualified Sort.Merge.Word
+
+main :: IO ()
+main = defaultMain $ testGroup "Sort"
+  [ testGroup "Contiguous"
+    [ tests (typeRep :: TypeRep Int8) (primArrayToByteArray . Data.Primitive.Sort.sort @PrimArray @Int8 . byteArrayToPrimArray)
+    , tests (typeRep :: TypeRep Word) (primArrayToByteArray . Data.Primitive.Sort.sort @PrimArray @Word . byteArrayToPrimArray)
+    , SC.testProperty "sortUnique == Set.toList . Set.fromList" $ \(list :: [Int]) ->
+        let actual = E.toList (Data.Primitive.Sort.sortUnique (E.fromList list :: Array Int))
+            expected = S.toList (S.fromList list)
+         in if actual == expected
+              then Right "unused"
+              else Left ("expected " ++ show expected ++ " but got " ++ show actual)
+    , testCase "sortTagged" $
+        Data.Primitive.Sort.sortTagged
+          (E.fromList [2, 1, 0] :: Array Int)
+          (E.fromList [True, True, False] :: Array Bool)
+        @?=
+        (E.fromList [0,1,2], E.fromList [False,True,True] :: Array Bool)
+    , testCase "sortUniqueTagged" $
+        Data.Primitive.Sort.sortUniqueTagged
+          (E.fromList [2, 1, 0] :: Array Int)
+          (E.fromList [True, True, False] :: Array Bool)
+        @?=
+        (E.fromList [0,1,2], E.fromList [False,True,True] :: Array Bool)
+    , SC.testProperty "sortUniqueTagged == Map.toList . Map.fromList" $ \(list :: [(Int,Bool)]) ->
+        let keys = E.fromList (map fst list) :: Array Int
+            vals = E.fromList (map snd list) :: Array Bool
+            (actualKeys,actualVals) = Data.Primitive.Sort.sortUniqueTagged keys vals
+            actual = zip (E.toList actualKeys) (E.toList actualVals)
+            expected = M.toList (M.fromList list)
+         in if actual == expected
+              then Right "unused"
+              else Left ("expected " ++ show expected ++ " but got " ++ show actual)
+    ]
+  , testGroup "Tagged"
+    [ testsTagged (typeRep :: TypeRep Word16) (typeRep :: TypeRep Word32)
+        (\k v -> pairPrimArrayToByteArray $ uncurry (Data.Primitive.Sort.sortTagged @Word16 @Word32) $ pairByteArrayToPrimArray (k,v))
+    ]
+  ]
+
+primArrayToByteArray :: PrimArray a -> ByteArray
+primArrayToByteArray (PrimArray x) = ByteArray x
+
+byteArrayToPrimArray :: ByteArray -> PrimArray a
+byteArrayToPrimArray (ByteArray x) = PrimArray x
+
+pairPrimArrayToByteArray :: (PrimArray a, PrimArray b) -> (ByteArray,ByteArray)
+pairPrimArrayToByteArray (PrimArray x,PrimArray y) = (ByteArray x,ByteArray y)
+
+pairByteArrayToPrimArray :: (ByteArray,ByteArray) -> (PrimArray a, PrimArray b) 
+pairByteArrayToPrimArray (ByteArray x,ByteArray y) = (PrimArray x,PrimArray y)
+
+tests :: forall n. (Prim n, Ord n, Show n, Arbitrary n, Serial IO n) => TypeRep n -> (ByteArray -> ByteArray) -> TestTree
+tests p sortArray = testGroup (show p) [properties (Proxy :: Proxy n) sortArray, unitTests (Proxy :: Proxy n) sortArray]
+
+testsTagged :: forall n a. (Prim a, Ord a, Show a, Arbitrary a, Serial IO a, Prim n, Ord n, Show n, Arbitrary n, Serial IO n, Num n, Enum n)
+  => TypeRep a -> TypeRep n -> (ByteArray -> ByteArray -> (ByteArray, ByteArray)) -> TestTree
+testsTagged p n sortArray = testGroup (show p ++ " tagged with " ++ show n) 
+  [ propertiesTagged (Proxy :: Proxy a) (Proxy :: Proxy n) sortArray
+  ]
+
+properties :: (Prim n, Ord n, Show n, Arbitrary n, Serial IO n) => Proxy n -> (ByteArray -> ByteArray) -> TestTree
+properties p sortArray = testGroup "Properties"
+  [ scProps p sortArray
+  , qcProps p sortArray
+  ]
+
+propertiesTagged :: (Prim a, Ord a, Show a, Arbitrary a, Serial IO a, Prim n, Ord n, Show n, Arbitrary n, Serial IO n, Num n, Enum n)
+  => Proxy a -> Proxy n -> (ByteArray -> ByteArray -> (ByteArray, ByteArray)) -> TestTree
+propertiesTagged p n sortArray = testGroup "Properties"
+  [ scPropsTagged p n sortArray
+  , qcPropsTagged p n sortArray
+  ]
+
+scProps :: forall n. (Prim n, Ord n, Show n, Serial IO n) => Proxy n -> (ByteArray -> ByteArray) -> TestTree
+scProps _ sortArray = testGroup "(checked by SmallCheck)"
+  [ SC.testProperty "sort == sort . reverse" $ \list ->
+      eqByteArray (sortArray (byteArrayFromList (list :: [n]))) (sortArray (byteArrayFromList (reverse list)))
+  , SC.testProperty "sort == Data.List.sort" $ \list ->
+      (==) (byteArrayToList (sortArray (byteArrayFromList (list :: [n])))) (Data.List.sort list)
+  ]
+
+scPropsTagged :: forall n a. (Prim a, Ord a, Show a, Serial IO a, Prim n, Ord n, Show n, Serial IO n, Num n, Enum n)
+  => Proxy a -> Proxy n -> (ByteArray -> ByteArray -> (ByteArray,ByteArray)) -> TestTree
+scPropsTagged _ _ sortArray = testGroup "(checked by SmallCheck)"
+  [ SC.testProperty "sort == Data.List.sort" $ \list ->
+      let taggedList = tagWithIndices list :: [Tag a n]
+          actual = taggedByteArrayToList (uncurry sortArray (taggedByteArrayFromList (taggedList :: [Tag a n])))
+          expected = Data.List.sort taggedList
+       in if actual == expected
+            then Right "unused"
+            else Left ("expected " ++ show expected ++ " but got " ++ show actual)
+  ]
+
+
+qcProps :: forall n. (Prim n, Arbitrary n, Show n, Ord n) => Proxy n -> (ByteArray -> ByteArray) -> TestTree
+qcProps p sortArray = testGroup "(checked by QuickCheck)"
+  [ testGroup "sort == sort . reverse"
+    [ sizedQuickCheckReverse p sortArray "small" 20 10 100
+    , sizedQuickCheckReverse p sortArray "medium" 5 10000 100000
+    , sizedQuickCheckReverse p sortArray "large" 2 100000 200000
+    ]
+  , testGroup "sort == Data.List.sort"
+    [ sizedQuickCheckCorrect p sortArray "small" 20 10 100
+    , sizedQuickCheckCorrect p sortArray "medium" 5 10000 100000
+    , sizedQuickCheckCorrect p sortArray "large" 2 100000 200000
+    ]
+  ]
+
+qcPropsTagged :: forall n a. (Prim a, Arbitrary a, Show a, Ord a, Prim n, Arbitrary n, Show n, Ord n)
+  => Proxy a -> Proxy n -> (ByteArray -> ByteArray -> (ByteArray,ByteArray)) -> TestTree
+qcPropsTagged p n sortArray = testGroup "(checked by QuickCheck)"
+  [ testGroup "sort == Data.List.sort"
+    [ sizedQuickCheckCorrectTagged p n sortArray "small" 20 10 100
+    , sizedQuickCheckCorrectTagged p n sortArray "medium" 5 10000 100000
+    , sizedQuickCheckCorrectTagged p n sortArray "large" 2 100000 200000
+    ]
+  ]
+
+sizedQuickCheckReverse :: forall n. (Arbitrary n, Prim n)
+  => Proxy n -> (ByteArray -> ByteArray) -> String -> Int -> Int -> Int -> TestTree
+sizedQuickCheckReverse _ sortArray szName countTests szMin szMax = 
+  adjustOption (\_ -> QC.QuickCheckTests countTests) $
+    QC.testProperty szName $ do
+      sz <- Q.choose (szMin,szMax)
+      list <- Q.vector sz
+      return (eqByteArray (sortArray (byteArrayFromList (list :: [n]))) (sortArray (byteArrayFromList (reverse list))))
+
+sizedQuickCheckCorrect :: forall n. (Arbitrary n, Prim n, Ord n, Show n)
+  => Proxy n -> (ByteArray -> ByteArray) -> String -> Int -> Int -> Int -> TestTree
+sizedQuickCheckCorrect _ sortArray szName countTests szMin szMax = 
+  adjustOption (\_ -> QC.QuickCheckTests countTests) $
+    QC.testProperty szName $ do
+      sz <- Q.choose (szMin,szMax)
+      list <- Q.vector sz
+      let actual = byteArrayToList (sortArray (byteArrayFromList (list :: [n])))
+          expected = Data.List.sort list
+      return $ if actual == expected
+        then property QP.succeeded
+        else if sz < 100
+          then property (QP.exception ("expected " ++ show expected ++ " but got " ++ show actual) (toException MyException))
+          else property QP.failed
+
+sizedQuickCheckCorrectTagged :: forall n a. (Arbitrary a, Prim a, Ord a, Show a, Arbitrary n, Prim n, Ord n, Show n)
+  => Proxy a -> Proxy n -> (ByteArray -> ByteArray -> (ByteArray,ByteArray)) -> String -> Int -> Int -> Int -> TestTree
+sizedQuickCheckCorrectTagged _ _ sortArray szName countTests szMin szMax = 
+  adjustOption (\_ -> QC.QuickCheckTests countTests) $
+    QC.testProperty szName $ do
+      sz <- Q.choose (szMin,szMax)
+      list <- Q.vector sz
+      let actual = taggedByteArrayToList (uncurry sortArray (taggedByteArrayFromList (list :: [Tag a n])))
+          expected = Data.List.sort list
+      return $ if actual == expected
+        then property QP.succeeded
+        else if sz < 100
+          then property (QP.exception ("expected " ++ show expected ++ " but got " ++ show actual) (toException MyException))
+          else property QP.failed
+
+data MyException = MyException
+  deriving (Show,Eq)
+instance Exception MyException
+
+unitTests :: forall n. Prim n => Proxy n -> (ByteArray -> ByteArray) -> TestTree
+unitTests _ _ = testGroup "Unit Tests"
+  [ -- testCase "List comparison (different length)" $
+    --   [1, 2, 3] `compare` [1,2] @?= GT
+  ]
+
+
+byteArrayToList :: forall a. Prim a => ByteArray -> [a]
+byteArrayToList arr = go 0 where
+  !len = div (P.sizeofByteArray arr) (P.sizeOf (undefined :: a))
+  go :: Int -> [a]
+  go !ix = if ix < len
+    then P.indexByteArray arr ix : go (ix + 1)
+    else []
+
+taggedByteArrayToList :: forall n a. (Prim a, Prim n) => (ByteArray, ByteArray) -> [Tag a n]
+taggedByteArrayToList (arr,tags) = go 0 where
+  !len = div (P.sizeofByteArray arr) (P.sizeOf (undefined :: a))
+  go :: Int -> [Tag a n]
+  go !ix = if ix < len
+    then Tag (P.indexByteArray arr ix) (P.indexByteArray tags ix) : go (ix + 1)
+    else []
+
+byteArrayFromList :: Prim a => [a] -> ByteArray
+byteArrayFromList xs = byteArrayFromListN (L.length xs) xs
+
+taggedByteArrayFromList :: (Prim a, Prim n) => [Tag a n] -> (ByteArray,ByteArray)
+taggedByteArrayFromList xs = taggedByteArrayFromListN (L.length xs) xs
+
+byteArrayFromListN :: forall a. Prim a => Int -> [a] -> ByteArray
+byteArrayFromListN len vs = runST run where
+  run :: forall s. ST s ByteArray
+  run = do
+    arr <- P.newByteArray (len * P.sizeOf (undefined :: a))
+    let go :: [a] -> Int -> ST s ()
+        go !xs !ix = case xs of
+          [] -> return ()
+          a : as -> do
+            P.writeByteArray arr ix a
+            go as (ix + 1)
+    go vs 0
+    P.unsafeFreezeByteArray arr
+
+taggedByteArrayFromListN :: forall n a. (Prim a, Prim n)
+  => Int -> [Tag a n] -> (ByteArray,ByteArray)
+taggedByteArrayFromListN len vs = runST run where
+  run :: forall s. ST s (ByteArray,ByteArray)
+  run = do
+    arr <- P.newByteArray (len * P.sizeOf (undefined :: a))
+    tags <- P.newByteArray (len * P.sizeOf (undefined :: n))
+    let go :: [Tag a n] -> Int -> ST s ()
+        go !xs !ix = case xs of
+          [] -> return ()
+          Tag a n : as -> do
+            P.writeByteArray arr ix a
+            P.writeByteArray tags ix n
+            go as (ix + 1)
+    go vs 0
+    liftA2 (,) (P.unsafeFreezeByteArray arr) (P.unsafeFreezeByteArray tags)
+
+
+eqByteArray :: ByteArray -> ByteArray -> Bool
+eqByteArray paA paB =
+  let !sizA = P.sizeofByteArray paA
+      !sizB = P.sizeofByteArray paB
+      go !ix = if ix < sizA
+        then if P.indexByteArray paA ix == (P.indexByteArray paB ix :: Word8)
+          then go (ix + 1)
+          else False
+        else True
+  in if sizA == sizB
+       then go 0
+       else False
+
+data Tag a b = Tag a b
+  deriving (Show)
+
+instance Eq a => Eq (Tag a b) where
+  Tag a1 _ == Tag a2 _ = a1 == a2
+
+instance Ord a => Ord (Tag a b) where
+  compare (Tag a1 _) (Tag a2 _) = compare a1 a2
+
+instance (Serial m a, Serial m b) => Serial m (Tag a b) where
+  series = fmap tagFromTuple series
+
+instance (Arbitrary a, Arbitrary b) => Arbitrary (Tag a b) where
+  arbitrary = liftA2 Tag arbitrary arbitrary
+
+tagFromTuple :: (a,b) -> Tag a b
+tagFromTuple (a,b) = Tag a b
+
+tagWithIndices :: (Num n, Enum n) => [a] -> [Tag a n]
+tagWithIndices xs = map tagFromTuple (zip xs [0,1..])
+
+
+