vector-algorithms (empty) → 0.3
raw patch · 22 files changed
+2374/−0 lines, 22 filesdep +basedep +primitivedep +vectorsetup-changed
Dependencies added: base, primitive, vector
Files
- Data/Vector/Algorithms/Combinators.hs +76/−0
- Data/Vector/Algorithms/Common.hs +27/−0
- Data/Vector/Algorithms/Insertion.hs +81/−0
- Data/Vector/Algorithms/Intro.hs +200/−0
- Data/Vector/Algorithms/Merge.hs +155/−0
- Data/Vector/Algorithms/Optimal.hs +244/−0
- Data/Vector/Algorithms/Radix.hs +282/−0
- Data/Vector/Algorithms/Search.hs +88/−0
- Data/Vector/Algorithms/TriHeap.hs +218/−0
- LICENSE +65/−0
- Setup.lhs +3/−0
- bench/Blocks.hs +62/−0
- bench/LICENSE +30/−0
- bench/Main.hs +188/−0
- bench/RadSieve.hs +97/−0
- bench/vector-algorithms-bench.cabal +22/−0
- include/vector.h +31/−0
- tests/Optimal.hs +62/−0
- tests/Properties.hs +167/−0
- tests/Tests.hs +151/−0
- tests/Util.hs +54/−0
- vector-algorithms.cabal +71/−0
+ Data/Vector/Algorithms/Combinators.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE Rank2Types, TypeOperators #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Combinators+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel <dan.doel@gmail.com>+-- Stability : Experimental+-- Portability : Non-portable (rank-2 types)+--+-- The purpose of this module is to supply various combinators for commonly+-- used idioms for the algorithms in this package. Examples at the time of+-- this writing include running an algorithm keyed on some function of the+-- elements (but only computing said function once per element), and safely+-- applying the algorithms on mutable arrays to immutable arrays.++module Data.Vector.Algorithms.Combinators+ ( apply+-- , usingKeys+-- , usingIxKeys+ ) where++import Prelude hiding (length)++import Control.Monad.ST++import Data.Ord++import Data.Vector.Generic++import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic.New as N++-- | Safely applies a mutable array algorithm to an immutable array.+apply :: (Vector v e) => (forall s mv. M.MVector mv e => mv s e -> ST s ()) -> v e -> v e+apply algo v = new . N.New $+ M.unstream (stream v) >>= \arr -> algo arr >> return arr++{-+-- | Uses a function to compute a key for each element which the+-- algorithm should use in lieu of the actual element. For instance:+--+-- > usingKeys sortBy f arr+--+-- should produce the same results as:+--+-- > sortBy (comparing f) arr+--+-- the difference being that usingKeys computes each key only once+-- which can be more efficient for expensive key functions.+usingKeys :: (UA e, UA k, Ord k)+ => (forall e'. (UA e') => Comparison e' -> MUArr e' s -> ST s ())+ -> (e -> k)+ -> MUArr e s+ -> ST s ()+usingKeys algo f arr = usingIxKeys algo (const f) arr+{-# INLINE usingKeys #-}++-- | As usingKeys, only the key function has access to the array index+-- at which each element is stored.+usingIxKeys :: (UA e, UA k, Ord k)+ => (forall e'. (UA e') => Comparison e' -> MUArr e' s -> ST s ())+ -> (Int -> e -> k)+ -> MUArr e s+ -> ST s ()+usingIxKeys algo f arr = do+ keys <- newMU (lengthMU arr)+ fill len keys+ algo (comparing fstS) (unsafeZipMU keys arr)+ where+ len = lengthMU arr+ fill k keys+ | k < 0 = return ()+ | otherwise = readMU arr k >>= writeMU keys k . f k >> fill (k-1) keys+{-# INLINE usingIxKeys #-}+-}
+ Data/Vector/Algorithms/Common.hs view
@@ -0,0 +1,27 @@++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Common+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel+-- Stability : Experimental+-- Portability : Portable+--+-- Common operations and utility functions for all sorts++module Data.Vector.Algorithms.Common where++import Prelude hiding (read, length)++import Control.Monad.Primitive++import Data.Vector.Generic.Mutable++-- | A type of comparisons between two values of a given type.+type Comparison e = e -> e -> Ordering++copyOffset :: (PrimMonad m, MVector v e)+ => v (PrimState m) e -> v (PrimState m) e -> Int -> Int -> Int -> m ()+copyOffset from to iFrom iTo len =+ unsafeCopy (unsafeSlice iTo len to) (unsafeSlice iFrom len from)+{-# INLINE copyOffset #-}
+ Data/Vector/Algorithms/Insertion.hs view
@@ -0,0 +1,81 @@++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Insertion+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel+-- Stability : Experimental+-- Portability : Portable+--+-- A simple insertion sort. Though it's O(n^2), its iterative nature can be+-- beneficial for small arrays. It is used to sort small segments of an array+-- by some of the more heavy-duty, recursive algorithms.++module Data.Vector.Algorithms.Insertion+ ( sort+ , sortBy+ , sortByBounds+ , sortByBounds'+ , Comparison+ ) where+++import Prelude hiding (read, length)++import Control.Monad.Primitive++import Data.Vector.Generic.Mutable++import Data.Vector.Algorithms.Common (Comparison)++import qualified Data.Vector.Algorithms.Optimal as O++-- | Sorts an entire array using the default comparison for the type+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an entire array using a given comparison+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()+sortBy cmp a = sortByBounds cmp a 0 (length a)+{-# INLINE sortBy #-}++-- | Sorts the portion of an array delimited by [l,u)+sortByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()+sortByBounds cmp a l u+ | len < 2 = return ()+ | len == 2 = O.sort2ByOffset cmp a l+ | len == 3 = O.sort3ByOffset cmp a l+ | len == 4 = O.sort4ByOffset cmp a l+ | otherwise = O.sort4ByOffset cmp a l >> sortByBounds' cmp a l (l + 4) u+ where+ len = u - l+{-# INLINE sortByBounds #-}++-- | Sorts the portion of the array delimited by [l,u) under the assumption+-- that [l,m) is already sorted.+sortByBounds' :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()+sortByBounds' cmp a l m u = sort m+ where+ sort i+ | i < u = do v <- unsafeRead a i+ insert cmp a l v i+ sort (i+1)+ | otherwise = return ()+{-# INLINE sortByBounds' #-}++-- Given a sorted array in [l,u), inserts val into its proper position,+-- yielding a sorted [l,u]+insert :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> e -> Int -> m ()+insert cmp a l = loop+ where+ loop val j+ | j <= l = unsafeWrite a l val+ | otherwise = do e <- unsafeRead a (j - 1)+ case cmp val e of+ LT -> unsafeWrite a j e >> loop val (j - 1)+ _ -> unsafeWrite a j val+{-# INLINE insert #-}
+ Data/Vector/Algorithms/Intro.hs view
@@ -0,0 +1,200 @@+{-# LANGUAGE TypeOperators, BangPatterns #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Intro+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel <dan.doel@gmail.com>+-- Stability : Experimental+-- Portability : Non-portable (type operators, bang patterns)+--+-- This module implements various algorithms based on the introsort algorithm,+-- originally described by David R. Musser in the paper /Introspective Sorting+-- and Selection Algorithms/. It is also in widespread practical use, as the+-- standard unstable sort used in the C++ Standard Template Library.+--+-- Introsort is at its core a quicksort. The version implemented here has the+-- following optimizations that make it perform better in practice:+--+-- * Small segments of the array are left unsorted until a final insertion+-- sort pass. This is faster than recursing all the way down to+-- one-element arrays.+--+-- * The pivot for segment [l,u) is chosen as the median of the elements at+-- l, u-1 and (u+l)/2. This yields good behavior on mostly sorted (or+-- reverse-sorted) arrays.+--+-- * The algorithm tracks its recursion depth, and if it decides it is+-- taking too long (depth greater than 2 * lg n), it switches to a heap+-- sort to maintain O(n lg n) worst case behavior. (This is what makes the+-- algorithm introsort).++module Data.Vector.Algorithms.Intro+ ( -- * Sorting+ sort+ , sortBy+ , sortByBounds + -- * Selecting+ , select+ , selectBy+ , selectByBounds+ -- * Partial sorting+ , partialSort+ , partialSortBy+ , partialSortByBounds+ , Comparison+ ) where++import Prelude hiding (read, length)++import Control.Monad+import Control.Monad.Primitive++import Data.Bits+import Data.Vector.Generic.Mutable++import Data.Vector.Algorithms.Common (Comparison)++import qualified Data.Vector.Algorithms.Insertion as I+import qualified Data.Vector.Algorithms.Optimal as O+import qualified Data.Vector.Algorithms.TriHeap as H++-- | Sorts an entire array using the default ordering.+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an entire array using a custom ordering.+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()+sortBy cmp a = sortByBounds cmp a 0 (length a)+{-# INLINE sortBy #-}++-- | Sorts a portion of an array [l,u) using a custom ordering+sortByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()+sortByBounds cmp a l u+ | len < 2 = return ()+ | len == 2 = O.sort2ByOffset cmp a l+ | len == 3 = O.sort3ByOffset cmp a l+ | len == 4 = O.sort4ByOffset cmp a l+ | otherwise = introsort cmp a (ilg len) l u+ where len = u - l+{-# INLINE sortByBounds #-}++-- Internal version of the introsort loop which allows partial+-- sort functions to call with a specified bound on iterations.+introsort :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()+introsort cmp a i l u = sort i l u >> I.sortByBounds cmp a l u+ where+ sort 0 l u = H.sortByBounds cmp a l u+ sort d l u+ | len < threshold = return ()+ | otherwise = do O.sort3ByIndex cmp a c l (u-1) -- sort the median into the lowest position+ p <- unsafeRead a l+ mid <- partitionBy cmp a p (l+1) u+ unsafeSwap a l (mid - 1)+ sort (d-1) mid u+ sort (d-1) l (mid - 1)+ where+ len = u - l+ c = (u + l) `div` 2+{-# INLINE introsort #-}++-- | Moves the least k elements to the front of the array in+-- no particular order.+select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()+select = selectBy compare+{-# INLINE select #-}++-- | Moves the least k elements (as defined by the comparison) to+-- the front of the array in no particular order.+selectBy :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> m ()+selectBy cmp a k = selectByBounds cmp a k 0 (length a)+{-# INLINE selectBy #-}++-- | Moves the least k elements in the interval [l,u) to the positions+-- [l,k+l) in no particular order.+selectByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()+selectByBounds cmp a k l u = go (ilg len) l (l + k) u+ where+ len = u - l+ go 0 l m u = H.selectByBounds cmp a (m - l) l u+ go n l m u = do O.sort3ByIndex cmp a c l (u-1)+ p <- unsafeRead a l+ mid <- partitionBy cmp a p (l+1) u+ unsafeSwap a l (mid - 1)+ if m > mid+ then go (n-1) mid m u+ else if m < mid - 1+ then go (n-1) l m (mid - 1)+ else return ()+ where c = (u + l) `div` 2+{-# INLINE selectByBounds #-}++-- | Moves the least k elements to the front of the array, sorted.+partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m ()+partialSort = partialSortBy compare+{-# INLINE partialSort #-}++-- | Moves the least k elements (as defined by the comparison) to+-- the front of the array, sorted.+partialSortBy :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> m ()+partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a)+{-# INLINE partialSortBy #-}++-- | Moves the least k elements in the interval [l,u) to the positions+-- [l,k+l), sorted.+partialSortByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()+partialSortByBounds cmp a k l u = go (ilg len) l (l + k) u+ where+ len = u - l+ go 0 l m n = H.partialSortByBounds cmp a (m - l) l u+ go n l m u+ | l == m = return ()+ | otherwise = do O.sort3ByIndex cmp a c l (u-1)+ p <- unsafeRead a l+ mid <- partitionBy cmp a p (l+1) u+ unsafeSwap a l (mid - 1)+ case compare m mid of+ GT -> do introsort cmp a (n-1) l (mid - 1)+ go (n-1) mid m u+ EQ -> introsort cmp a (n-1) l m+ LT -> go n l m (mid - 1)+ where c = (u + l) `div` 2+{-# INLINE partialSortByBounds #-}++partitionBy :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int+partitionBy cmp a = partUp+ where+ partUp p l u+ | l < u = do e <- unsafeRead a l+ case cmp e p of+ LT -> partUp p (l+1) u+ _ -> partDown p l (u-1)+ | otherwise = return l+ partDown p l u+ | l < u = do e <- unsafeRead a u+ case cmp p e of+ LT -> partDown p l (u-1)+ _ -> unsafeSwap a l u >> partUp p (l+1) u+ | otherwise = return l+{-# INLINE partitionBy #-}++-- computes the number of recursive calls after which heapsort should+-- be invoked given the lower and upper indices of the array to be sorted+ilg :: Int -> Int+ilg m = 2 * loop m 0+ where+ loop 0 !k = k - 1+ loop n !k = loop (n `shiftR` 1) (k+1)++-- the size of array at which the introsort algorithm switches to insertion sort+threshold :: Int+threshold = 18+{-# INLINE threshold #-}
+ Data/Vector/Algorithms/Merge.hs view
@@ -0,0 +1,155 @@+-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Merge+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel <dan.doel@gmail.com>+-- Stability : Experimental+-- Portability : Portable+--+-- This module implements a simple top-down merge sort. The temporary buffer+-- is preallocated to 1/2 the size of the input array, and shared through+-- the entire sorting process to ease the amount of allocation performed in+-- total. This is a stable sort.++module Data.Vector.Algorithms.Merge+ ( sort+ , sortBy+ , sortByBounds+ , Comparison+ ) where++import Prelude hiding (read, length)++import Control.Monad.Primitive++import Data.Bits+import Data.Vector.Generic.Mutable++import Data.Vector.Algorithms.Common (Comparison, copyOffset)++import qualified Data.Vector.Algorithms.Optimal as O+import qualified Data.Vector.Algorithms.Insertion as I++-- | Sorts an array using the default comparison.+sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an array using a custom comparison.+sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m ()+sortBy cmp vec = sortByBounds cmp vec 0 (length vec)+{-+sortBy cmp vec+ | len <= 1 = return ()+ | len == 2 = O.sort2ByOffset cmp vec 0+ | len == 3 = O.sort3ByOffset cmp vec 0+ | len == 4 = O.sort4ByOffset cmp vec 0+ | otherwise = do buf <- new len+ mergeSortWithBuf cmp vec buf+ where+ len = length vec+-}+{-# INLINE sortBy #-}+++-- | Sorts a portion of an array [l,u) using a custom comparison.+sortByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()+sortByBounds cmp vec l u+ | len <= 1 = return ()+ | len == 2 = O.sort2ByOffset cmp vec l+ | len == 3 = O.sort3ByOffset cmp vec l+ | len == 4 = O.sort4ByOffset cmp vec l+ | otherwise = do tmp <- new size+ mergeSortWithBuf cmp vec tmp l u+ where+ len = u - l+ size = (u + l) `div` 2 - l+{-# INLINE sortByBounds #-}++mergeSortWithBuf :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> v (PrimState m) e+ -> Int -> Int -> m ()+mergeSortWithBuf cmp arr tmp = loop+ where+ loop l u+ | len < threshold = I.sortByBounds cmp arr l u+ | otherwise = do loop l mid+ loop mid u+ merge cmp arr tmp l mid u+ where+ len = u - l+ mid = (u + l) `shiftR` 1+{-# INLINE mergeSortWithBuf #-}++{-+mergeSortWithBuf :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> v (PrimState m) e -> m ()+mergeSortWithBuf cmp src buf+ | length src < threshold = I.sortByBounds cmp src 0 (length src)+ | otherwise = do mergeSortWithBuf cmp srcL bufL+ mergeSortWithBuf cmp srcU bufU+ merge cmp src buf mid+ where+ len = length src+ mid = len `shiftR` 1++ srcL = unsafeSlice 0 mid src+ srcU = unsafeSlice mid (len - mid) src+ + bufL = unsafeSlice 0 mid buf+ bufU = unsafeSlice mid (len - mid) buf+-}++merge :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> v (PrimState m) e+ -> Int -> Int -> Int -> m ()+merge cmp arr tmp l m u = do copyOffset arr tmp l 0 uTmp+ eTmp <- unsafeRead tmp 0+ eArr <- unsafeRead arr m+ loop 0 eTmp m eArr l+ where+ uTmp = m - l+ uArr = u+ loop iTmp eTmp iArr eArr iIns+ | iTmp >= uTmp = return ()+ | iArr >= uArr = copyOffset tmp arr iTmp iIns (uTmp - iTmp)+ | otherwise = case cmp eArr eTmp of+ LT -> do unsafeWrite arr iIns eArr+ eArr <- unsafeRead arr (iArr+1)+ loop iTmp eTmp (iArr+1) eArr (iIns+1)+ _ -> do unsafeWrite arr iIns eTmp+ eTmp <- unsafeRead tmp (iTmp+1)+ loop (iTmp+1) eTmp iArr eArr (iIns+1)+{-# INLINE merge #-}++{-+merge :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> v (PrimState m) e+ -> Int -> m ()+merge cmp src buf mid = do unsafeCopy tmp lower+ eTmp <- unsafeRead tmp 0+ eUpp <- unsafeRead upper 0+ loop tmp 0 eTmp upper 0 eUpp 0+ where+ lower = unsafeSlice 0 mid src+ upper = unsafeSlice mid (length src - mid) src+ tmp = unsafeSlice 0 mid buf++ loop low iLow eLow high iHigh eHigh iIns+ | iLow >= length low = return ()+ | 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)+{-# INLINE merge #-}+-}++threshold :: Int+threshold = 25+{-# INLINE threshold #-}
+ Data/Vector/Algorithms/Optimal.hs view
@@ -0,0 +1,244 @@++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Optimal+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel+-- Stability : Experimental+-- Portability : Portable+--+-- Optimal sorts for very small array sizes, or for small numbers of+-- particular indices in a larger array (to be used, for instance, for+-- sorting a median of 3 values into the lowest position in an array+-- for a median-of-3 quicksort).++-- The code herein was adapted from a C algorithm for optimal sorts+-- of small arrays. The original code was produced for the article+-- /Sorting Revisited/ by Paul Hsieh, available here:+--+-- http://www.azillionmonkeys.com/qed/sort.html+--+-- The LICENSE file contains the relevant copyright information for+-- the reference C code.++module Data.Vector.Algorithms.Optimal+ ( sort2ByIndex+ , sort2ByOffset+ , sort3ByIndex+ , sort3ByOffset+ , sort4ByIndex+ , sort4ByOffset+ , Comparison+ ) where++import Prelude hiding (read, length)++import Control.Monad.Primitive++import Data.Vector.Generic.Mutable++import Data.Vector.Algorithms.Common (Comparison)++#include "vector.h"++-- | Sorts the elements at the positions 'off' and 'off + 1' in the given+-- array using the comparison.+sort2ByOffset :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> m ()+sort2ByOffset cmp a off = sort2ByIndex cmp a off (off + 1)+{-# INLINE sort2ByOffset #-}++-- | Sorts the elements at the two given indices using the comparison. This+-- is essentially a compare-and-swap, although the first index is assumed to+-- be the 'lower' of the two.+sort2ByIndex :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> m ()+sort2ByIndex cmp a i j = UNSAFE_CHECK(checkIndex) "sort2ByIndex" i (length a)+ $ UNSAFE_CHECK(checkIndex) "sort2ByIndex" j (length a) $ do+ a0 <- unsafeRead a i+ a1 <- unsafeRead a j+ case cmp a0 a1 of+ GT -> unsafeWrite a i a1 >> unsafeWrite a j a0+ _ -> return ()+{-# INLINE sort2ByIndex #-}++-- | Sorts the three elements starting at the given offset in the array.+sort3ByOffset :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> m ()+sort3ByOffset cmp a off = sort3ByIndex cmp a off (off + 1) (off + 2)+{-# INLINE sort3ByOffset #-}++-- | Sorts the elements at the three given indices. The indices are assumed+-- to be given from lowest to highest, so if 'l < m < u' then+-- 'sort3ByIndex cmp a m l u' essentially sorts the median of three into the+-- lowest position in the array.+sort3ByIndex :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m ()+sort3ByIndex cmp a i j k = UNSAFE_CHECK(checkIndex) "sort3ByIndex" i (length a)+ $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" j (length a) + $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" k (length a) $ do+ a0 <- unsafeRead a i+ a1 <- unsafeRead a j+ a2 <- unsafeRead a k+ case cmp a0 a1 of+ GT -> case cmp a0 a2 of+ GT -> case cmp a2 a1 of+ LT -> do unsafeWrite a i a2+ unsafeWrite a k a0+ _ -> do unsafeWrite a i a1+ unsafeWrite a j a2+ unsafeWrite a k a0+ _ -> do unsafeWrite a i a1+ unsafeWrite a j a0+ _ -> case cmp a1 a2 of+ GT -> case cmp a0 a2 of+ GT -> do unsafeWrite a i a2+ unsafeWrite a j a0+ unsafeWrite a k a1+ _ -> do unsafeWrite a j a2+ unsafeWrite a k a1+ _ -> return ()+{-# INLINE sort3ByIndex #-}++-- | Sorts the four elements beginning at the offset.+sort4ByOffset :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> m ()+sort4ByOffset cmp a off = sort4ByIndex cmp a off (off + 1) (off + 2) (off + 3)+{-# INLINE sort4ByOffset #-}++-- The horror...++-- | Sorts the elements at the four given indices. Like the 2 and 3 element+-- versions, this assumes that the indices are given in increasing order, so+-- it can be used to sort medians into particular positions and so on.+sort4ByIndex :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> Int -> m ()+sort4ByIndex cmp a i j k l = UNSAFE_CHECK(checkIndex) "sort4ByIndex" i (length a)+ $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" j (length a)+ $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" k (length a)+ $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" l (length a) $ do+ a0 <- unsafeRead a i+ a1 <- unsafeRead a j+ a2 <- unsafeRead a k+ a3 <- unsafeRead a l+ case cmp a0 a1 of+ GT -> case cmp a0 a2 of+ GT -> case cmp a1 a2 of+ GT -> case cmp a1 a3 of+ GT -> case cmp a2 a3 of+ GT -> do unsafeWrite a i a3+ unsafeWrite a j a2+ unsafeWrite a k a1+ unsafeWrite a l a0+ _ -> do unsafeWrite a i a2+ unsafeWrite a j a3+ unsafeWrite a k a1+ unsafeWrite a l a0+ _ -> case cmp a0 a3 of+ GT -> do unsafeWrite a i a2+ unsafeWrite a j a1+ unsafeWrite a k a3+ unsafeWrite a l a0+ _ -> do unsafeWrite a i a2+ unsafeWrite a j a1+ unsafeWrite a k a0+ unsafeWrite a l a3+ _ -> case cmp a2 a3 of+ GT -> case cmp a1 a3 of+ GT -> do unsafeWrite a i a3+ unsafeWrite a j a1+ unsafeWrite a k a2+ unsafeWrite a l a0+ _ -> do unsafeWrite a i a1+ unsafeWrite a j a3+ unsafeWrite a k a2+ unsafeWrite a l a0+ _ -> case cmp a0 a3 of+ GT -> do unsafeWrite a i a1+ unsafeWrite a j a2+ unsafeWrite a k a3+ unsafeWrite a l a0+ _ -> do unsafeWrite a i a1+ unsafeWrite a j a2+ unsafeWrite a k a0+ -- unsafeWrite a l a3+ _ -> case cmp a0 a3 of+ GT -> case cmp a1 a3 of+ GT -> do unsafeWrite a i a3+ -- unsafeWrite a j a1+ unsafeWrite a k a0+ unsafeWrite a l a2+ _ -> do unsafeWrite a i a1+ unsafeWrite a j a3+ unsafeWrite a k a0+ unsafeWrite a l a2+ _ -> case cmp a2 a3 of+ GT -> do unsafeWrite a i a1+ unsafeWrite a j a0+ unsafeWrite a k a3+ unsafeWrite a l a2+ _ -> do unsafeWrite a i a1+ unsafeWrite a j a0+ -- unsafeWrite a k a2+ -- unsafeWrite a l a3+ _ -> case cmp a1 a2 of+ GT -> case cmp a0 a2 of+ GT -> case cmp a0 a3 of+ GT -> case cmp a2 a3 of+ GT -> do unsafeWrite a i a3+ unsafeWrite a j a2+ unsafeWrite a k a0+ unsafeWrite a l a1+ _ -> do unsafeWrite a i a2+ unsafeWrite a j a3+ unsafeWrite a k a0+ unsafeWrite a l a1+ _ -> case cmp a1 a3 of+ GT -> do unsafeWrite a i a2+ unsafeWrite a j a0+ unsafeWrite a k a3+ unsafeWrite a l a1+ _ -> do unsafeWrite a i a2+ unsafeWrite a j a0+ unsafeWrite a k a1+ -- unsafeWrite a l a3+ _ -> case cmp a2 a3 of+ GT -> case cmp a0 a3 of+ GT -> do unsafeWrite a i a3+ unsafeWrite a j a0+ -- unsafeWrite a k a2+ unsafeWrite a l a1+ _ -> do -- unsafeWrite a i a0+ unsafeWrite a j a3+ -- unsafeWrite a k a2+ unsafeWrite a l a1+ _ -> case cmp a1 a3 of+ GT -> do -- unsafeWrite a i a0+ unsafeWrite a j a2+ unsafeWrite a k a3+ unsafeWrite a l a1+ _ -> do -- unsafeWrite a i a0+ unsafeWrite a j a2+ unsafeWrite a k a1+ -- unsafeWrite a l a3+ _ -> case cmp a1 a3 of+ GT -> case cmp a0 a3 of+ GT -> do unsafeWrite a i a3+ unsafeWrite a j a0+ unsafeWrite a k a1+ unsafeWrite a l a2+ _ -> do -- unsafeWrite a i a0+ unsafeWrite a j a3+ unsafeWrite a k a1+ unsafeWrite a l a2+ _ -> case cmp a2 a3 of+ GT -> do -- unsafeWrite a i a0+ -- unsafeWrite a j a1+ unsafeWrite a k a3+ unsafeWrite a l a2+ _ -> do -- unsafeWrite a i a0+ -- unsafeWrite a j a1+ -- unsafeWrite a k a2+ -- unsafeWrite a l a3+ return ()+{-# INLINE sort4ByIndex #-}
+ Data/Vector/Algorithms/Radix.hs view
@@ -0,0 +1,282 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns, TypeOperators #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Radix+-- Copyright : (c) 2008-2010 Dan Doel+-- Maintainer : Dan Doel <dan.doel@gmail.com>+-- Stability : Experimental+-- Portability : Non-portable (scoped type variables, bang patterns)+--+-- This module provides a radix sort for a subclass of unboxed arrays. The +-- radix class gives information on+-- * the number of passes needed for the data type+--+-- * the size of the auxiliary arrays+--+-- * how to compute the pass-k radix of a value+--+-- Radix sort is not a comparison sort, so it is able to achieve O(n) run+-- time, though it also uses O(n) auxiliary space. In addition, there is a+-- constant space overhead of 2*size*sizeOf(Int) for the sort, so it is not+-- advisable to use this sort for large numbers of very small arrays.+--+-- A standard example (upon which one could base their own Radix instance)+-- is Word32:+--+-- * We choose to sort on r = 8 bits at a time+--+-- * A Word32 has b = 32 bits total+--+-- Thus, b/r = 4 passes are required, 2^r = 256 elements are needed in an+-- auxiliary array, and the radix function is:+--+-- > radix k e = (e `shiftR` (k*8)) .&. 256++module Data.Vector.Algorithms.Radix (sort, sortBy, Radix(..)) where++import Prelude hiding (read, length)++import Control.Monad+import Control.Monad.Primitive++import qualified Data.Vector.Primitive.Mutable as PV+import Data.Vector.Generic.Mutable++import Data.Vector.Algorithms.Common (Comparison)++import Data.Bits+import Data.Int+import Data.Word+++import Foreign.Storable++class Radix e where+ -- | The number of passes necessary to sort an array of es+ passes :: e -> Int+ -- | The size of an auxiliary array+ size :: e -> Int+ -- | The radix function parameterized by the current pass+ radix :: Int -> e -> Int++instance Radix Int where+ passes _ = sizeOf (undefined :: Int)+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = e .&. 255+ radix i e+ | i == passes e - 1 = radix' (e `xor` minBound)+ | otherwise = radix' e+ where radix' e = (e `shiftR` (i `shiftL` 3)) .&. 255+ {-# INLINE radix #-}++instance Radix Int8 where+ passes _ = 1+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix _ e = 255 .&. fromIntegral e `xor` 128 + {-# INLINE radix #-}++instance Radix Int16 where+ passes _ = 2+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix 1 e = fromIntegral (((e `xor` minBound) `shiftR` 8) .&. 255)+ {-# INLINE radix #-}++instance Radix Int32 where+ passes _ = 4+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)+ radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)+ radix 3 e = fromIntegral (((e `xor` minBound) `shiftR` 24) .&. 255)+ {-# INLINE radix #-}++instance Radix Int64 where+ passes _ = 8+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)+ radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)+ radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)+ radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255)+ radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255)+ radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255)+ radix 7 e = fromIntegral (((e `xor` minBound) `shiftR` 56) .&. 255)+ {-# INLINE radix #-}++instance Radix Word where+ passes _ = sizeOf (undefined :: Word)+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix i e = fromIntegral ((e `shiftR` (i `shiftL` 3)) .&. 255)+ {-# INLINE radix #-}++instance Radix Word8 where+ passes _ = 1+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix _ = fromIntegral+ {-# INLINE radix #-}++instance Radix Word16 where+ passes _ = 2+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)+ {-# INLINE radix #-}++instance Radix Word32 where+ passes _ = 4+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)+ radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)+ radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)+ {-# INLINE radix #-}++instance Radix Word64 where+ passes _ = 8+ {-# INLINE passes #-}+ size _ = 256+ {-# INLINE size #-}+ radix 0 e = fromIntegral (e .&. 255)+ radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255)+ radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255)+ radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255)+ radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255)+ radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255)+ radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255)+ radix 7 e = fromIntegral ((e `shiftR` 56) .&. 255)+ {-# INLINE radix #-}++instance (Radix i, Radix j) => Radix (i, j) where+ passes ~(i, j) = passes i + passes j+ {-# INLINE passes #-}+ size ~(i, j) = size i `max` size j+ {-# INLINE size #-}+ radix k ~(i, j) | k < passes j = radix k j+ | otherwise = radix (k - passes j) i+ {-# INLINE radix #-}++-- | Sorts an array based on the Radix instance.+sort :: forall e m v. (PrimMonad m, MVector v e, Radix e)+ => v (PrimState m) e -> m ()+sort arr = sortBy (passes e) (size e) radix arr+ where+ e :: e+ e = undefined+{-# INLINE sort #-}++-- | Radix sorts an array using custom radix information+-- requires the number of passes to fully sort the array,+-- the size of of auxiliary arrays necessary (should be+-- one greater than the maximum value returned by the radix+-- function), and a radix function, which takes the pass+-- and an element, and returns the relevant radix.+sortBy :: (PrimMonad m, MVector v e)+ => Int -- ^ the number of passes+ -> Int -- ^ the size of auxiliary arrays+ -> (Int -> e -> Int) -- ^ the radix function+ -> v (PrimState m) e -- ^ the array to be sorted+ -> m ()+sortBy passes size rdx arr = do+ tmp <- new (length arr)+ count <- new size+ prefix <- new size+ radixLoop passes rdx arr tmp count prefix+{-# INLINE sortBy #-}++radixLoop :: (PrimMonad m, MVector v e)+ => Int -- passes+ -> (Int -> e -> Int) -- radix function+ -> v (PrimState m) e -- array to sort+ -> v (PrimState m) e -- temporary array+ -> PV.MVector (PrimState m) Int -- radix count array+ -> PV.MVector (PrimState m) Int -- placement array+ -> m ()+radixLoop passes rdx src dst count prefix = 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)+ | otherwise = when swap (unsafeCopy src dst)+{-# INLINE radixLoop #-}++body :: (PrimMonad m, MVector v e)+ => (Int -> e -> Int) -- radix function+ -> v (PrimState m) e -- source array+ -> v (PrimState m) e -- destination array+ -> PV.MVector (PrimState m) Int -- radix count+ -> 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+{-# 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+ 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'+ | otherwise = return ()+{-# INLINE prefixLoop #-}++moveLoop :: (PrimMonad m, MVector v e)+ => Int -> (Int -> e -> Int) -> v (PrimState m) e+ -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m ()+moveLoop k rdx src dst prefix = go 0+ where+ len = length src+ go i+ | i < len = do srci <- unsafeRead src i+ pf <- inc prefix (rdx k srci)+ unsafeWrite dst pf srci+ go (i+1)+ | otherwise = return ()+{-# INLINE moveLoop #-}++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 #-}
+ Data/Vector/Algorithms/Search.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE BangPatterns #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Algorithms.Search+-- Copyright : (c) 2009-2010 Dan Doel+-- Maintainer : Dan Doel <dan.doel@gmail.com>+-- Stability : Experimental+-- Portability : Non-portable (bang patterns)+--+-- This module implements several methods of searching for indicies to insert+-- elements into a sorted array.++module Data.Vector.Algorithms.Search+ ( binarySearch+ , binarySearchBy+ , binarySearchByBounds+ , binarySearchL+ , binarySearchLBy+ , binarySearchLByBounds+ , Comparison+ ) where++import Prelude hiding (read, length)++import Control.Monad.Primitive++import Data.Bits++import Data.Vector.Generic.Mutable++import Data.Vector.Algorithms.Common (Comparison)++-- | Finds an index in a givesn sorted array at which the given element could+-- be inserted while maintaining the sortedness of the array.+binarySearch :: (PrimMonad m, MVector v e, Ord e)+ => v (PrimState m) e -> e -> m Int+binarySearch = binarySearchBy compare++-- | Finds an index in a given array, which must be sorted with respect to the+-- given comparison function, at which the given element could be inserted while+-- preserving the array's sortedness.+binarySearchBy :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> e -> m Int+binarySearchBy cmp arr e = binarySearchByBounds cmp arr e 0 (length arr)++-- | Given an array sorted with respect to a given comparison function in indices+-- in [l,u), finds an index in [l,u] at which the given element could be inserted+-- while preserving sortedness.+binarySearchByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int+binarySearchByBounds cmp arr e l u+ | u <= l = return l+ | otherwise = do e' <- unsafeRead arr k+ case cmp e' e of+ LT -> binarySearchByBounds cmp arr e (k+1) u+ EQ -> return k+ GT -> binarySearchByBounds cmp arr e l k+ where k = (u + l) `shiftR` 1+{-# INLINE binarySearchByBounds #-}++-- | Finds the lowest index in a given sorted array at which the given element+-- could be inserted while maintaining the sortedness.+binarySearchL :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> e -> m Int+binarySearchL = binarySearchLBy compare+{-# INLINE binarySearchL #-}++-- | Finds the lowest index in a given array, which must be sorted with respect to +-- the given comparison function, at which the given element could be inserted+-- while preserving the sortedness.+binarySearchLBy :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> e -> m Int+binarySearchLBy cmp arr e = binarySearchLByBounds cmp arr e 0 (length arr)+{-# INLINE binarySearchLBy #-}++-- | Given an array sorted with respect to a given comparison function on indices+-- in [l,u), finds the lowest index in [l,u] at which the given element could be+-- inserted while preserving sortedness.+binarySearchLByBounds :: (PrimMonad m, MVector v e)+ => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int+binarySearchLByBounds cmp arr e !l !u+ | u <= l = return l+ | otherwise = do e' <- unsafeRead arr k+ case cmp e' e of+ LT -> binarySearchLByBounds cmp arr e (k+1) u+ _ -> binarySearchLByBounds cmp arr e l k+ where k = (u + l) `shiftR` 1+{-# INLINE binarySearchLByBounds #-}
+ Data/Vector/Algorithms/TriHeap.hs view
@@ -0,0 +1,218 @@+{-# 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 #-}
+ LICENSE view
@@ -0,0 +1,65 @@+Copyright (c) 2008-2010 Dan Doel++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.++------------------------------------------------------------------------------++The code in Data.Array.Vector.Algorithms.Mutable.Optimal is adapted from a C+algorithm for the same purpose. The folowing is the copyright notice for said+C code:++Copyright (c) 2004 Paul Hsieh+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 sorttest nor the names of its 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.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ bench/Blocks.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE Rank2Types #-}++module Blocks where++import Control.Monad+import Control.Monad.ST++import Data.Vector.Unboxed.Mutable++import System.CPUTime++import System.Random.Mersenne++-- Some conveniences for doing evil stuff in the ST monad.+-- All the tests get run in IO, but uvector stuff happens+-- in ST, so we temporarily coerce.+clock :: IO Integer+clock = getCPUTime++-- Strategies for filling the initial arrays+rand :: (MTRandom e) => MTGen -> Int -> IO e+rand g _ = random g++ascend :: Num e => Int -> IO e+ascend = return . fromIntegral++descend :: Num e => e -> Int -> IO e+descend m n = return $ m - fromIntegral n++modulo :: Integral e => e -> Int -> IO e+modulo m n = return $ fromIntegral n `mod` m++-- This is the worst case for the median-of-three quicksort+-- used in the introsort implementation.+medianKiller :: Integral e => e -> Int -> IO e+medianKiller m n'+ | n < k = return $ if even n then n + 1 else n + k+ | otherwise = return $ (n - k + 1) * 2+ where+ n = fromIntegral n'+ k = m `div` 2+{-# INLINE medianKiller #-}++initialize :: (Unbox e) => MVector RealWorld e -> Int -> (Int -> IO e) -> IO ()+initialize arr len fill = init $ len - 1+ where init n = fill n >>= unsafeWrite arr n >> when (n > 0) (init $ n - 1)+{-# INLINE initialize #-}++speedTest :: (Unbox e) => Int+ -> (Int -> IO e)+ -> (MVector RealWorld e -> IO ())+ -> IO Integer+speedTest n fill algo = do+ arr <- new n+ initialize arr n fill+ t0 <- clock+ algo arr+ t1 <- clock+ return $ t1 - t0+{-# INLINE speedTest #-}++
+ bench/LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2009 Dan Doel++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.
+ bench/Main.hs view
@@ -0,0 +1,188 @@+{-# LANGUAGE Rank2Types #-}++module Main (main) where++import Prelude hiding (read, length)+import qualified Prelude as P++import Control.Monad.ST+import Control.Monad.Error++import Data.Char+import Data.Ord (comparing)+import Data.List (maximumBy)++import Data.Vector.Unboxed.Mutable++import qualified Data.Vector.Algorithms.Insertion as INS+import qualified Data.Vector.Algorithms.Intro as INT+import qualified Data.Vector.Algorithms.TriHeap as TH+import qualified Data.Vector.Algorithms.Merge as M+import qualified Data.Vector.Algorithms.Radix as R++import System.Environment+import System.Console.GetOpt+import System.Random.Mersenne++import Blocks++-- Does nothing. For testing the speed/heap allocation of the building blocks.+noalgo :: (Unbox e) => MVector RealWorld e -> IO ()+noalgo _ = return ()++-- Allocates a temporary buffer, like mergesort for similar purposes as noalgo.+alloc :: (Unbox e) => MVector RealWorld e -> IO ()+alloc arr | len <= 4 = arr `seq` return ()+ | otherwise = (new (len `div` 2) :: IO (MVector RealWorld Int)) >> return ()+ where len = length arr++displayTime :: String -> Integer -> IO ()+displayTime s elapsed = putStrLn $+ s ++ " : " ++ show (fromIntegral elapsed / 1e12) ++ " seconds"++run :: String -> IO Integer -> IO ()+run s t = t >>= displayTime s++sortSuite :: String -> MTGen -> Int -> (MVector RealWorld Int -> IO ()) -> IO ()+sortSuite str g n sort = do+ putStrLn $ "Testing: " ++ str+ run "Random " $ speedTest n (rand g >=> modulo n) sort+ run "Sorted " $ speedTest n ascend sort+ run "Reverse-sorted " $ speedTest n (descend n) sort+ run "Random duplicates " $ speedTest n (rand g >=> modulo 1000) sort+ let m = 4 * (n `div` 4)+ run "Median killer " $ speedTest m (medianKiller m) sort++partialSortSuite :: String -> MTGen -> Int -> Int+ -> (MVector RealWorld Int -> Int -> IO ()) -> IO ()+partialSortSuite str g n k sort = sortSuite str g n (\a -> sort a k)++-- -----------------+-- Argument handling+-- -----------------++data Algorithm = DoNothing+ | Allocate+ | InsertionSort+ | IntroSort+ | IntroPartialSort+ | IntroSelect+ | TriHeapSort+ | TriHeapPartialSort+ | TriHeapSelect+ | MergeSort+ | RadixSort+ deriving (Show, Read, Enum, Bounded)++data Options = O { algos :: [Algorithm], elems :: Int, portion :: Int, usage :: Bool } deriving (Show)++defaultOptions :: Options+defaultOptions = O [] 10000 1000 False++type OptionsT = Options -> Either String Options++options :: [OptDescr OptionsT]+options = [ Option ['A'] ["algorithm"] (ReqArg parseAlgo "ALGO")+ ("Specify an algorithm to be run. Options:\n" ++ algoOpts)+ , Option ['n'] ["num-elems"] (ReqArg parseN "INT")+ "Specify the size of arrays in algorithms."+ , Option ['k'] ["portion"] (ReqArg parseK "INT")+ "Specify the number of elements to partial sort/select in\nrelevant algorithms."+ , Option ['?','v'] ["help"] (NoArg $ \o -> Right $ o { usage = True })+ "Show options."+ ]+ where+ allAlgos :: [Algorithm]+ allAlgos = [minBound .. maxBound]+ algoOpts = fmt allAlgos+ fmt (x:y:zs) = '\t' : pad (show x) ++ show y ++ "\n" ++ fmt zs+ fmt [x] = '\t' : show x ++ "\n"+ fmt [] = ""+ size = (" " ++) . maximumBy (comparing P.length) . map show $ allAlgos+ pad str = zipWith const (str ++ repeat ' ') size++parseAlgo :: String -> Options -> Either String Options+parseAlgo "None" o = Right $ o { algos = [] }+parseAlgo "All" o = Right $ o { algos = [DoNothing .. RadixSort] }+parseAlgo s o = leftMap (\e -> "Unrecognized algorithm `" ++ e ++ "'")+ . fmap (\v -> o { algos = v : algos o }) $ readEither s++leftMap :: (a -> b) -> Either a c -> Either b c+leftMap f (Left a) = Left (f a)+leftMap _ (Right c) = Right c++parseNum :: (Int -> Options) -> String -> Either String Options+parseNum f = leftMap (\e -> "Invalid numeric argument `" ++ e ++ "'") . fmap f . readEither++parseN, parseK :: String -> Options -> Either String Options+parseN s o = parseNum (\n -> o { elems = n }) s+parseK s o = parseNum (\k -> o { portion = k }) s++readEither :: Read a => String -> Either String a+readEither s = case reads s of+ [(x,t)] | all isSpace t -> Right x+ _ -> Left s++runTest :: MTGen -> Int -> Int -> Algorithm -> IO ()+runTest g n k alg = case alg of+ DoNothing -> sortSuite "no algorithm" g n noalgo+ Allocate -> sortSuite "allocate" g n alloc+ InsertionSort -> sortSuite "insertion sort" g n insertionSort+ IntroSort -> sortSuite "introsort" g n introSort+ IntroPartialSort -> partialSortSuite "partial introsort" g n k introPSort+ IntroSelect -> partialSortSuite "introselect" g n k introSelect+ 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+ MergeSort -> sortSuite "merge sort" g n mergeSort+ RadixSort -> sortSuite "radix sort" g n radixSort+ _ -> putStrLn $ "Currently unsupported algorithm: " ++ show alg++mergeSort :: MVector RealWorld Int -> IO ()+mergeSort v = M.sort v+{-# NOINLINE mergeSort #-}++introSort :: MVector RealWorld Int -> IO ()+introSort v = INT.sort v+{-# NOINLINE introSort #-}++introPSort :: MVector RealWorld Int -> Int -> IO ()+introPSort v k = INT.partialSort v k+{-# NOINLINE introPSort #-}++introSelect :: MVector RealWorld Int -> Int -> IO ()+introSelect v k = INT.select v k+{-# NOINLINE introSelect #-}++triHeapSort :: MVector RealWorld Int -> IO ()+triHeapSort v = TH.sort v+{-# NOINLINE triHeapSort #-}++triHeapPSort :: MVector RealWorld Int -> Int -> IO ()+triHeapPSort v k = TH.partialSort v k+{-# NOINLINE triHeapPSort #-}++triHeapSelect :: MVector RealWorld Int -> Int -> IO ()+triHeapSelect v k = TH.select v k+{-# NOINLINE triHeapSelect #-}++insertionSort :: MVector RealWorld Int -> IO ()+insertionSort v = INS.sort v+{-# NOINLINE insertionSort #-}++radixSort :: MVector RealWorld Int -> IO ()+radixSort v = R.sort v+{-# NOINLINE radixSort #-}++main :: IO ()+main = do args <- getArgs+ gen <- getStdGen+ case getOpt Permute options args of+ (fs, _, []) -> case foldl (>>=) (Right defaultOptions) fs of+ Left err -> putStrLn $ usageInfo err options+ Right opts | not (usage opts) ->+ mapM_ (runTest gen (elems opts) (portion opts)) (algos opts)+ | otherwise -> putStrLn $ usageInfo "uvector-algorithms-bench" options+ (_, _, errs) -> putStrLn $ usageInfo (concat errs) options++
+ bench/RadSieve.hs view
@@ -0,0 +1,97 @@+-- ------------------------------------------------------------------+--+-- Module : RadSieve+-- Copyright : (c) 2009 Dan Doel+--+-- ------------------------------------------------------------------+-- An implementation of a radical sieve, inspired by solving Project+-- Euler problem #124.+--+-- Reproduction fo the problem text:+--+-- The radical of n, rad(n), is the product of distinct prime factors+-- of n. For example, 504 = 23 × 32 × 7, so rad(504) = 2 × 3 × 7 = 42.+--+-- If we calculate rad(n) for 1 ≤ n ≤ 10, then sort them on rad(n),+-- and sorting on n if the radical values are equal, we get:+--+-- Unsorted Sorted+-- n rad(n) n rad(n) k+-- 1 1 1 1 1+-- 2 2 2 2 2+-- 3 3 4 2 3+-- 4 2 8 2 4+-- 5 5 3 3 5+-- 6 6 9 3 6+-- 7 7 5 5 7+-- 8 2 6 6 8+-- 9 3 7 7 9+-- 10 10 10 10 10+--+-- Let E(k) be the kth element in the sorted n column; for example,+-- E(4) = 8 and E(6) = 9.+--+-- If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000).++module RadSieve where++import Control.Monad+import Control.Monad.ST++import Data.Array.Vector++-- Radicals can be sieved as follows:+-- set a[1,n] = 1+-- for i from 2 to n+-- if a[i] == 1 -- i must be prime+-- then a[j*i] *= i for positive integers j, j*i <= n+-- else do nothing -- i is composite, so its prime factors+-- -- have been accounted for+--+-- This sieves for radicals up to the given integer.+radSieve :: Int -> ST s (MUArr Int s)+radSieve n = do arr <- newMU (n + 1)+ fill arr n+ sieve arr 1+ return arr+ where+ fill arr i | i < 0 = return ()+ | otherwise = writeMU arr i 1 >> fill arr (i-1)+ sieve arr i | n < i = return ()+ | otherwise = do e <- readMU arr i+ when (e == 1) $ mark arr i i+ sieve arr (i+1)+ mark arr p j | n < j = return ()+ | otherwise = readMU arr j >>= writeMU arr j . (*p)+ >> mark arr p (j+p)++-- Computes the answer to the above Project Euler problem. The correct+-- answer is only generated for a stable sorting function.+stableSortedRad :: Int -> Int+ -> (forall s e. UA e => Comparison e -> MUArr e s -> ST s ()) + -> Int+stableSortedRad n k sortBy = runST (do rads <- radSieve n+ index <- newMU (n + 1)+ fillUp index n+ sortBy (comparing fstS)+ (unsafeZipMU rads index)+ readMU k index)+ where+ fillUp arr k | k < 0 = return ()+ | otherwise = writeMU arr k k >> fillUp arr (k-1)++-- Computes the answer to the above Project Euler problem. This version+-- will generate the correct answer even for unstable sorts, but may be+-- marginally slower.+unstableSortedRad :: Int -> Int+ -> (forall s e. UA e => Comparison e -> MUArr e s -> ST s ()) + -> Int+unstableSortedRad n k sortBy = runST (do rads <- radSieve n+ index <- newMU (n + 1)+ fillUp index n+ sortBy compare (unsafeZipMU rads index)+ readMU k index)+ where+ fillUp arr k | k < 0 = return ()+ | otherwise = writeMU arr k k >> fillUp arr (k-1)+
+ bench/vector-algorithms-bench.cabal view
@@ -0,0 +1,22 @@+name: vector-algorithms-bench+version: 0.3+license: BSD3+license-file: LICENSE+author: Dan Doel+maintainer: Dan Doel <dan.doel@gmail.com>+homepage: http://code.haskell.org/~doio/+category: Benchmark+synopsis: Benchmarks for vector-algorithms+description: A suite of various benchmarks for verifying the+ performance of the algorithms in vector-algorithms.+build-type: Simple+cabal-version: >= 1.2++executable vec-bench+ build-depends: base, mersenne-random, vector, vector-algorithms, mtl++ ghc-options: -Wall -Odph+ main-is: Main.hs++ extensions:+ Rank2Types
+ include/vector.h view
@@ -0,0 +1,31 @@+#define PHASE_STREAM [1]+#define PHASE_INNER [0]++#define INLINE_STREAM INLINE PHASE_STREAM+#define INLINE_INNER INLINE PHASE_INNER++#ifndef NOT_VECTOR_MODULE+import qualified Data.Vector.Internal.Check as Ck+#endif++#define ERROR(f) (Ck.f __FILE__ __LINE__)+#define ASSERT (Ck.assert __FILE__ __LINE__)+#define ENSURE (Ck.f __FILE__ __LINE__)+#define CHECK(f) (Ck.f __FILE__ __LINE__)++#define BOUNDS_ERROR(f) (ERROR(f) Ck.Bounds)+#define BOUNDS_ASSERT (ASSERT Ck.Bounds)+#define BOUNDS_ENSURE (ENSURE Ck.Bounds)+#define BOUNDS_CHECK(f) (CHECK(f) Ck.Bounds)++#define UNSAFE_ERROR(f) (ERROR(f) Ck.Unsafe)+#define UNSAFE_ASSERT (ASSERT Ck.Unsafe)+#define UNSAFE_ENSURE (ENSURE Ck.Unsafe)+#define UNSAFE_CHECK(f) (CHECK(f) Ck.Unsafe)++#define INTERNAL_ERROR(f) (ERROR(f) Ck.Internal)+#define INTERNAL_ASSERT (ASSERT Ck.Internal)+#define INTERNAL_ENSURE (ENSURE Ck.Internal)+#define INTERNAL_CHECK(f) (CHECK(f) Ck.Internal)++
+ tests/Optimal.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE TypeOperators, FlexibleContexts #-}++-- Exhaustive test sets for proper sorting and stability of+-- optimal sorts++module Optimal where++import Control.Arrow+import Control.Monad++import Data.List+import Data.Function++import Data.Vector.Generic hiding (map, zip, concatMap, (++), replicate)++interleavings :: [a] -> [a] -> [[a]]+interleavings [ ] ys = [ys]+interleavings xs [ ] = [xs]+interleavings xs@(x:xt) ys@(y:yt) = map (x:) (interleavings xt ys)+ ++ map (y:) (interleavings xs yt)++monotones :: Int -> Int -> [[Int]]+monotones k = atLeastOne 0+ where+ atLeastOne i 0 = [[]]+ atLeastOne i n = map (i:) $ picks i (n-1)+ picks _ 0 = [[]]+ picks i n | i >= k = [replicate n k]+ | otherwise = map (i:) (picks i (n-1)) ++ atLeastOne (i+1) n+++stability :: (Vector v (Int,Int)) => Int -> [v (Int, Int)]+stability n = concatMap ( map fromList+ . foldM interleavings []+ . groupBy ((==) `on` fst)+ . flip zip [0..])+ $ monotones (n-2) n++sort2 :: (Vector v Int) => [v Int]+sort2 = map fromList $ permutations [0,1]++stability2 :: (Vector v (Int,Int)) => [v (Int, Int)]+stability2 = [fromList [(0, 0), (0, 1)]]++sort3 :: (Vector v Int) => [v Int]+sort3 = map fromList $ permutations [0..2]++{-+stability3 :: [UArr (Int :*: Int)]+stability3 = map toU [ [0:*:0, 0:*:1, 0:*:2]+ , [0:*:0, 0:*:1, 1:*:2]+ , [0:*:0, 1:*:2, 0:*:1]+ , [1:*:2, 0:*:0, 0:*:1]+ , [0:*:0, 1:*:1, 1:*:2]+ , [1:*:1, 0:*:0, 1:*:2]+ , [1:*:1, 1:*:2, 0:*:0]+ ]+-}++sort4 :: (Vector v Int) => [v Int]+sort4 = map fromList $ permutations [0..3]+
+ tests/Properties.hs view
@@ -0,0 +1,167 @@+{-# LANGUAGE RankNTypes, FlexibleContexts #-}++module Properties where++import Prelude++import Optimal++import Control.Monad+import Control.Monad.ST++import Data.List+import Data.Ord++import Data.Vector (Vector)+import qualified Data.Vector as V++import Data.Vector.Mutable (MVector)+import qualified Data.Vector.Mutable as MV++import qualified Data.Vector.Generic.Mutable as G++import Data.Vector.Algorithms.Optimal (Comparison)+import Data.Vector.Algorithms.Radix (radix, passes, size)+import Data.Vector.Algorithms.Combinators++import qualified Data.Map as M++import Test.QuickCheck++import Util++prop_sorted :: (Ord e) => Vector e -> Property+prop_sorted arr | V.length arr < 2 = property True+ | otherwise = check (V.head arr) (V.tail arr)+ where+ check e arr | V.null arr = property True+ | otherwise = e <= V.head arr .&. check (V.head arr) (V.tail arr)++prop_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_schwartzian :: (UA e, UA k, Ord k)+ => (e -> k)+ -> (forall e s. (UA e) => (e -> e -> Ordering) -> MUArr e s -> ST s ())+ -> UArr e -> Property+prop_schwartzian f algo arr+ | lengthU arr < 2 = property True+ | otherwise = let srt = apply (algo `usingKeys` f) arr+ in check (headU srt) (tailU srt)+ where+ check e arr | nullU arr = property True+ | otherwise = f e <= f (headU arr) .&. check (headU arr) (tailU arr)+-}++longGen :: (Arbitrary e) => Int -> Gen (Vector e)+longGen k = liftM2 (\l r -> V.fromList (l ++ r)) (vectorOf k arbitrary) arbitrary++sanity :: Int+sanity = 100++prop_partialsort :: (Ord e, Arbitrary e, Show e)+ => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())+ -> Positive Int -> Property+prop_partialsort = prop_sized $ \algo k ->+ prop_sorted . V.take k . apply algo++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+ l = V.slice 0 k vec'+ r = V.slice k (V.length vec' - k) vec'+ in V.all (\e -> V.all (e <=) r) l++prop_sized :: (Arbitrary e, Show e, Testable prop)+ => ((forall s mv. G.MVector mv e => mv s e -> ST s ())+ -> Int -> Vector e -> prop)+ -> (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ())+ -> Positive Int -> Property+prop_sized prop algo (Positive k) =+ let k' = k `mod` sanity+ in forAll (longGen k') $ prop (\marr -> algo marr k') k'++prop_stable :: (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> ST s ())+ -> Vector Int -> Property+-- prop_stable algo arr = property $ apply algo arr == arr+prop_stable algo arr = stable $ apply (algo (comparing fst)) $ V.zip arr ix+ where+ ix = V.fromList [1 .. V.length arr]++stable arr | V.null arr = property True+ | otherwise = let (e, i) = V.head arr+ in V.all (\(e', i') -> e < e' || i < i') (V.tail arr)+ .&. stable (V.tail arr)++prop_stable_radix :: (forall e s mv. G.MVector mv e => Int -> Int -> (Int -> e -> Int) + -> mv s e -> ST s ())+ -> Vector Int -> Property+prop_stable_radix algo arr =+ stable . apply (algo (passes e) (size e) (\k (e, _) -> radix k e))+ $ V.zip arr ix+ where+ ix = V.fromList [1 .. V.length arr]+ e = V.head arr+ +prop_optimal :: Int+ -> (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> Int -> ST s ())+ -> Property+prop_optimal n algo = label "sorting" sortn .&. label "stability" stabn+ where+ arrn = V.fromList [0..n-1]+ sortn = all ( (== arrn)+ . apply (\a -> algo compare a 0)+ . V.fromList)+ $ permutations [0..n-1]+ stabn = all ( (== arrn)+ . snd+ . V.unzip+ . apply (\a -> algo (comparing fst) a 0))+ $ stability n++type Bag e = M.Map e Int++toBag :: (Ord e) => Vector e -> Bag e+toBag = M.fromListWith (+) . flip zip (repeat 1) . V.toList++prop_permutation :: (Ord e) => (forall s mv. G.MVector mv e => mv s e -> ST s ())+ -> Vector e -> Property+prop_permutation algo arr = property $ + toBag arr == toBag (apply algo arr)++newtype SortedVec e = Sorted (Vector e)++instance (Show e) => Show (SortedVec e) where+ show (Sorted a) = show a++instance (Arbitrary e, Ord e) => Arbitrary (SortedVec e) where+ arbitrary = fmap (Sorted . V.fromList . sort)+ $ liftM2 (++) (vectorOf 20 arbitrary) arbitrary++ixRanges :: Vector e -> Gen (Int, Int)+ixRanges vec = do i <- fmap (`mod` len) arbitrary+ j <- fmap (`mod` len) arbitrary+ return $ if i < j then (i, j) else (j, i)+ where len = V.length vec++prop_search_inrange :: (Ord e)+ => (forall s. MVector s e -> e -> Int -> Int -> ST s Int)+ -> SortedVec e -> e -> Property+prop_search_inrange algo (Sorted arr) e = forAll (ixRanges arr) $ \(i, j) ->+ let k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e i j)+ in property $ i <= k && k <= j+ where+ len = V.length arr++prop_search_lowbound :: (Ord e)+ => (forall s. MVector s e -> e -> ST s Int)+ -> SortedVec e -> e -> Property+prop_search_lowbound algo (Sorted arr) e = property $ (k == 0 || arr V.! (k-1) < e)+ && (k == len || arr V.! k >= e)+ where+ len = V.length arr+ k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e)
+ tests/Tests.hs view
@@ -0,0 +1,151 @@+{-# LANGUAGE ImpredicativeTypes, RankNTypes, TypeOperators, FlexibleContexts #-}++module Main (main) where++import Properties++import Util++import Test.QuickCheck++import Control.Monad+import Control.Monad.ST++import Data.Int+import Data.Word++import Data.Vector (Vector)+import qualified Data.Vector as V++import Data.Vector.Generic.Mutable (MVector)+import qualified Data.Vector.Generic.Mutable as MV++import Data.Vector.Algorithms.Combinators++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.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+ , maxDiscard = 200+ }++check_Int_sort = forM_ algos $ \(name,algo) ->+ quickCheckWith args (label name . prop_fullsort algo)+ where+ algos :: [(String, Algo Int ())]+ algos = [ ("introsort", INT.sort)+ , ("insertion sort", INS.sort)+ , ("merge sort", M.sort)+ , ("tri-heapsort", TH.sort)+ ]++check_Int_partialsort = forM_ algos $ \(name,algo) ->+ quickCheckWith args (label name . prop_partialsort algo)+ where+ algos :: [(String, SizeAlgo Int ())]+ algos = [ ("intro-partialsort", INT.partialSort)+ , ("tri-heap partialsort", TH.partialSort)+ ]++check_Int_select = forM_ algos $ \(name,algo) ->+ quickCheckWith args (label name . prop_select algo)+ where+ algos :: [(String, SizeAlgo Int ())]+ algos = [ ("intro-select", INT.select)+ , ("tri-heap select", TH.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) ()))+ where+ qc algo = quickCheckWith args algo++{-+check_schwartzian = do+ quickCheckWith args (prop_schwartzian i2w INS.sortBy)+ where+ i2w :: Int -> Word+ i2w = fromIntegral+-}++check_stable = do quickCheckWith args (label "merge sort" . prop_stable M.sortBy)+ quickCheckWith args (label "radix sort" . prop_stable_radix R.sortBy)++check_optimal = do qc . label "size 2" $ prop_optimal 2 O.sort2ByOffset+ qc . label "size 3" $ prop_optimal 3 O.sort3ByOffset+ qc . label "size 4" $ prop_optimal 4 O.sort4ByOffset+ where+ qc = quickCheck++check_permutation = do+ qc $ label "introsort" . prop_permutation (INT.sort :: Algo Int ())+ qc $ label "intropartial" . prop_sized (const . prop_permutation)+ (INT.partialSort :: SizeAlgo Int ())+ qc $ label "introselect" . prop_sized (const . prop_permutation)+ (INT.select :: SizeAlgo Int ())+ qc $ label "heapsort" . prop_permutation (TH.sort :: Algo Int ())+ qc $ label "heappartial" . prop_sized (const . prop_permutation)+ (TH.partialSort :: SizeAlgo Int ())+ qc $ label "heapselect" . prop_sized (const . prop_permutation)+ (TH.select :: SizeAlgo Int ())+ qc $ label "mergesort" . prop_permutation (M.sort :: Algo Int ())+ qc $ label "radix I8" . prop_permutation (R.sort :: Algo Int8 ())+ qc $ label "radix I16" . prop_permutation (R.sort :: Algo Int16 ())+ qc $ label "radix I32" . prop_permutation (R.sort :: Algo Int32 ())+ qc $ label "radix I64" . prop_permutation (R.sort :: Algo Int64 ())+ qc $ label "radix Int" . prop_permutation (R.sort :: Algo Int ())+ qc $ label "radix W8" . prop_permutation (R.sort :: Algo Word8 ())+ qc $ label "radix W16" . prop_permutation (R.sort :: Algo Word16 ())+ qc $ label "radix W32" . prop_permutation (R.sort :: Algo Word32 ())+ qc $ label "radix W64" . prop_permutation (R.sort :: Algo Word64 ())+ qc $ label "radix Word" . prop_permutation (R.sort :: Algo Word ())+ where+ qc prop = quickCheckWith args prop++type BoundSAlgo e r = forall s mv. MVector mv e => mv s e -> e -> Int -> Int -> ST s r++check_search_range = do+ qc $ (label "binarySearchL" .) + . prop_search_inrange (SR.binarySearchLByBounds compare :: BoundSAlgo Int Int)+ qc $ (label "binarySearch" .)+ . prop_search_inrange (SR.binarySearchByBounds compare :: BoundSAlgo Int Int)+ where+ qc prop = quickCheckWith args prop++main = do putStrLn "Int tests:"+ check_Int_sort+ check_Int_partialsort+ check_Int_select+ putStrLn "Radix sort tests:"+ check_radix_sorts+-- putStrLn "Schwartzian transform (Int -> Word):"+-- check_schwartzian+ putStrLn "Stability:"+ check_stable+ putStrLn "Optimals:"+ check_optimal+ putStrLn "Permutation:"+ check_permutation+ putStrLn "Search in range:"+ check_search_range
+ tests/Util.hs view
@@ -0,0 +1,54 @@+{-# LANGUAGE TypeOperators #-}++module Util where++import Control.Monad+import Control.Monad.ST++import Data.Word+import Data.Int++import qualified Data.Vector as V++import Data.Vector.Mutable hiding (length)++import Test.QuickCheck+++mfromList :: [e] -> ST s (MVector s e)+mfromList l = do v <- new (length l)+ fill l 0 v+ where+ fill [] _ v = return v+ fill (x:xs) i v = do write v i x+ fill xs (i+1) v++instance (Arbitrary e) => Arbitrary (V.Vector e) where+ arbitrary = fmap V.fromList arbitrary++instance Arbitrary 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
+ vector-algorithms.cabal view
@@ -0,0 +1,71 @@+Name: vector-algorithms+Version: 0.3+License: BSD3+License-File: LICENSE+Author: Dan Doel+Maintainer: Dan Doel <dan.doel@gmail.com>+Homepage: http://code.haskell.org/~dolio/+Category: Data+Synopsis: Efficient algorithms for vector arrays+Description: Efficient algorithms for vector arrays+ be sure to compile with -O2, and -fvia-C -optc-O3 is+ recommended.+Build-Type: Simple+Cabal-Version: >= 1.2.3++Flag BoundsChecks+ Description: Enable bounds checking+ Default: True++Flag UnsafeChecks+ Description: Enable bounds checking in unsafe operations at the cost of a+ significant performance penalty.+ Default: False++Flag InternalChecks+ Description: Enable internal consistency checks at the cost of a+ significant performance penalty.+ Default: False++Library+ Build-Depends: base >= 3 && < 5, vector >= 0.5 && < 0.6, primitive >=0.2 && <0.3++ Exposed-Modules:+ Data.Vector.Algorithms.Combinators+ Data.Vector.Algorithms.Optimal+ Data.Vector.Algorithms.Insertion+ Data.Vector.Algorithms.Intro+ Data.Vector.Algorithms.Merge+ Data.Vector.Algorithms.Radix+ Data.Vector.Algorithms.Search+ Data.Vector.Algorithms.TriHeap++ Other-Modules:+ Data.Vector.Algorithms.Common++ Extensions:+ BangPatterns,+ TypeOperators,+ Rank2Types,+ ScopedTypeVariables,+ CPP++ GHC-Options:+ -Odph+ -funbox-strict-fields++ Include-Dirs:+ include++ Install-Includes:+ vector.h++ if flag(BoundsChecks)+ cpp-options: -DVECTOR_BOUNDS_CHECKS++ if flag(UnsafeChecks)+ cpp-options: -DVECTOR_UNSAFE_CHECKS++ if flag(InternalChecks)+ cpp-options: -DVECTOR_INTERNAL_CHECKS+