uvector-algorithms (empty) → 0.1
raw patch · 11 files changed
+1096/−0 lines, 11 filesdep +basedep +uvectorsetup-changed
Dependencies added: base, uvector
Files
- Data/Array/Vector/Algorithms/Common.hs +32/−0
- Data/Array/Vector/Algorithms/Immutable.hs +24/−0
- Data/Array/Vector/Algorithms/Insertion.hs +74/−0
- Data/Array/Vector/Algorithms/Intro.hs +187/−0
- Data/Array/Vector/Algorithms/Merge.hs +84/−0
- Data/Array/Vector/Algorithms/Optimal.hs +191/−0
- Data/Array/Vector/Algorithms/Radix.hs +207/−0
- Data/Array/Vector/Algorithms/TriHeap.hs +189/−0
- LICENSE +65/−0
- Setup.lhs +3/−0
- uvector-algorithms.cabal +40/−0
+ Data/Array/Vector/Algorithms/Common.hs view
@@ -0,0 +1,32 @@+-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Common+-- Copyright : (c) 2008 Dan Doel+-- Maintainer : Dan Doel+-- Stability : Experimental+-- Portability : Portable+--+-- Common operations and utility functions for all sorts++module Data.Array.Vector.Algorithms.Common where++import Control.Monad.ST++import Data.Array.Vector++-- | A type of comparisons between two values of a given type.+type Comparison e = e -> e -> Ordering++-- | Swaps the elements at two positions in an array.+swap :: (UA e) => MUArr e s -> Int -> Int -> ST s ()+swap arr i j = do ei <- readMU arr i+ readMU arr j >>= writeMU arr i+ writeMU arr j ei+{-# INLINE swap #-}++mcopyMU :: (UA e) => MUArr e s -> MUArr e s -> Int -> Int -> Int -> ST s ()+mcopyMU from to iFrom iTo len = go 0+ where+ go n | n < len = readMU from (iFrom + n) >>= writeMU to (iTo + n) >> go (n+1)+ | otherwise = return ()+{-# INLINE mcopyMU #-}
+ Data/Array/Vector/Algorithms/Immutable.hs view
@@ -0,0 +1,24 @@+{-# LANGUAGE Rank2Types #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Immutable+-- Copyright : (c) 2008 Dan Doel+-- Maintainer : Dan Doel <dan.doel@gmail.com>+-- Stability : Experimental+-- Portability : Non-portable (rank-2 types)+--+-- The purpose of this module is to apply the algorithms on mutable arrays+-- in other packages to immutable arrays. The idea is to copy the immutable+-- array into a mutable intermediate, perform the algorithm on the mutable+-- array, and freeze it, yielding a new immutable array.++module Data.Array.Vector.Algorithms.Immutable ( apply ) where++import Control.Monad.ST++import Data.Array.Vector++-- | Safely applies a mutable array algorithm to an immutable array.+apply :: (UA e) => (forall s. MUArr e s -> ST s ()) -> UArr e -> UArr e+apply algo v = newU (lengthU v) (\arr -> copyMU arr 0 v >> algo arr)
+ Data/Array/Vector/Algorithms/Insertion.hs view
@@ -0,0 +1,74 @@++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Insertion+-- Copyright : (c) 2008 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.Array.Vector.Algorithms.Insertion+ ( sort+ , sortBy+ , sortByBounds+ , sortByBounds'+ ) where+++import Control.Monad.ST++import Data.Array.Vector+import Data.Array.Vector.Algorithms.Common++import qualified Data.Array.Vector.Algorithms.Optimal as O++-- | Sorts an entire array using the default comparison for the type+sort :: (UA e, Ord e) => MUArr e s -> ST s ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an entire array using a given comparison+sortBy :: (UA e) => Comparison e -> MUArr e s -> ST s ()+sortBy cmp a = sortByBounds cmp a 0 (lengthMU a)+{-# INLINE sortBy #-}++-- | Sorts the portion of an array delimited by [l,u)+sortByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+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' :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+sortByBounds' cmp a l m u = sort m+ where+ sort i+ | i < u = do v <- readMU 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 :: (UA e) => Comparison e -> MUArr e s -> Int -> e -> Int -> ST s ()+insert cmp a l = loop+ where+ loop val j+ | j <= l = writeMU a l val+ | otherwise = do e <- readMU a (j - 1)+ case cmp val e of+ LT -> writeMU a j e >> loop val (j - 1)+ _ -> writeMU a j val+{-# INLINE insert #-}
+ Data/Array/Vector/Algorithms/Intro.hs view
@@ -0,0 +1,187 @@+{-# LANGUAGE TypeOperators, BangPatterns #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Intro+-- Copyright : (c) 2008 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.Array.Vector.Algorithms.Intro+ ( -- * Sorting+ sort+ , sortBy+ , sortByBounds + -- * Selecting+ , select+ , selectBy+ , selectByBounds+ -- * Partial sorting+ , partialSort+ , partialSortBy+ , partialSortByBounds+ ) where++import Control.Monad+import Control.Monad.ST++import Data.Array.Vector+import Data.Array.Vector.Algorithms.Common+import Data.Bits++import qualified Data.Array.Vector.Algorithms.Insertion as I+import qualified Data.Array.Vector.Algorithms.Optimal as O+import qualified Data.Array.Vector.Algorithms.TriHeap as H++-- | Sorts an entire array using the default ordering.+sort :: (UA e, Ord e) => MUArr e s -> ST s ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an entire array using a custom ordering.+sortBy :: (UA e) => Comparison e -> MUArr e s -> ST s ()+sortBy cmp a = sortByBounds cmp a 0 (lengthMU a)+{-# INLINE sortBy #-}++-- | Sorts a portion of an array [l,u) using a custom ordering+sortByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+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 <- readMU a l+ mid <- partitionBy cmp a p (l+1) u+ swap 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 :: (UA e, Ord e) => MUArr e s -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+selectBy cmp a k = selectByBounds cmp a k 0 (lengthMU a)+{-# INLINE selectBy #-}++-- | Moves the least k elements in the interval [l,u) to the positions+-- [l,k+l) in no particular order.+selectByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+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 <- readMU a l+ mid <- partitionBy cmp a p (l+1) u+ swap 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 :: (UA e, Ord e) => MUArr e s -> Int -> ST s ()+partialSort = partialSortBy compare+{-# INLINE partialSort #-}++-- | Moves the least k elements (as defined by the comparison) to+-- the front of the array, sorted.+partialSortBy :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+partialSortBy cmp a k = partialSortByBounds cmp a k 0 (lengthMU a)+{-# INLINE partialSortBy #-}++-- | Moves the least k elements in the interval [l,u) to the positions+-- [l,k+l), sorted.+partialSortByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+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 = do O.sort3ByIndex cmp a c l (u-1)+ p <- readMU a l+ mid <- partitionBy cmp a p (l+1) u+ swap 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 :: (UA e) => Comparison e -> MUArr e s -> e -> Int -> Int -> ST s Int+partitionBy cmp a = partUp+ where+ partUp p l u+ | l < u = do e <- readMU 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 <- readMU a u+ case cmp p e of+ LT -> partDown p l (u-1)+ _ -> swap 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/Array/Vector/Algorithms/Merge.hs view
@@ -0,0 +1,84 @@+-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Merge+-- Copyright : (c) 2008 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.Array.Vector.Algorithms.Merge (sort, sortBy, sortByBounds) where++import Control.Monad.ST++import Data.Bits+import Data.Array.Vector+import Data.Array.Vector.Algorithms.Common++import qualified Data.Array.Vector.Algorithms.Optimal as O+import qualified Data.Array.Vector.Algorithms.Insertion as I++-- | Sorts an array using the default comparison.+sort :: (Ord e, UA e) => MUArr e s -> ST s ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an array using a custom comparison.+sortBy :: (UA e) => Comparison e -> MUArr e s -> ST s ()+sortBy cmp arr = sortByBounds cmp arr 0 (lengthMU arr)+{-# INLINE sortBy #-}++-- | Sorts a portion of an array [l,u) using a custom comparison.+sortByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+sortByBounds cmp arr l u+ | len < 1 = return ()+ | len == 2 = O.sort2ByOffset cmp arr l+ | len == 3 = O.sort3ByOffset cmp arr l+ | len == 4 = O.sort4ByOffset cmp arr l+ | otherwise = do tmp <- newMU size+ mergeSortWithBuf cmp arr tmp l u+ where+ len = u - l+ size = (u + l) `div` 2 - l+{-# INLINE sortByBounds #-}++mergeSortWithBuf :: (UA e) => Comparison e -> MUArr e s -> MUArr e s -> Int -> Int -> ST s ()+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 #-}++merge :: (UA e) => Comparison e -> MUArr e s -> MUArr e s -> Int -> Int -> Int -> ST s ()+merge cmp arr tmp l m u = do mcopyMU arr tmp l 0 uTmp+ eTmp <- readMU tmp 0+ eArr <- readMU 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 = mcopyMU tmp arr iTmp iIns (uTmp - iTmp)+ | otherwise = case cmp eArr eTmp of+ LT -> do writeMU arr iIns eArr+ eArr <- readMU arr (iArr+1)+ loop iTmp eTmp (iArr+1) eArr (iIns+1)+ _ -> do writeMU arr iIns eTmp+ eTmp <- readMU tmp (iTmp+1)+ loop (iTmp+1) eTmp iArr eArr (iIns+1)+{-# INLINE merge #-}++threshold :: Int+threshold = 25+{-# INLINE threshold #-}
+ Data/Array/Vector/Algorithms/Optimal.hs view
@@ -0,0 +1,191 @@++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Optimal+-- Copyright : (c) 2008 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.Array.Vector.Algorithms.Optimal+ ( sort2ByIndex+ , sort2ByOffset+ , sort3ByIndex+ , sort3ByOffset+ , sort4ByIndex+ , sort4ByOffset+ ) where++import Control.Monad.ST++import Data.Array.Vector++import Data.Array.Vector.Algorithms.Common+-- | Sorts the elements at the positions 'off' and 'off + 1' in the given+-- array using the comparison.+sort2ByOffset :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+sort2ByIndex cmp a i j = do+ a0 <- readMU a i+ a1 <- readMU a j+ case cmp a0 a1 of+ GT -> writeMU a i a1 >> writeMU a j a0+ _ -> return ()+{-# INLINE sort2ByIndex #-}++-- | Sorts the three elements starting at the given offset in the array.+sort3ByOffset :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+sort3ByIndex cmp a i j k = do+ a0 <- readMU a i+ a1 <- readMU a j+ a2 <- readMU a k+ case cmp a0 a1 of+ GT -> case cmp a0 a2 of+ GT -> case cmp a2 a1 of+ GT -> do writeMU a i a1+ writeMU a j a2+ writeMU a k a0+ _ -> do writeMU a i a2+ writeMU a k a0+ _ -> do writeMU a i a1+ writeMU a j a0+ _ -> case cmp a1 a2 of+ GT -> case cmp a2 a0 of+ GT -> do writeMU a j a2+ writeMU a k a1+ _ -> do writeMU a i a2+ writeMU a k a1+ writeMU a j a0+ _ -> return ()+{-# INLINE sort3ByIndex #-}++-- | Sorts the four elements beginning at the offset.+sort4ByOffset :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> Int -> ST s ()+sort4ByIndex cmp a i j k l = do+ a0 <- readMU a i+ a1 <- readMU a j+ a2 <- readMU a k+ a3 <- readMU a l+ case cmp a0 a1 of+ LT -> case cmp a1 a2 of+ LT -> case cmp a1 a3 of+ LT -> case cmp a2 a3 of+ GT -> do writeMU a k a3+ writeMU a l a2+ _ -> return ()+ _ -> do case cmp a0 a3 of+ LT -> writeMU a j a3+ _ -> do writeMU a j a0+ writeMU a i a3+ writeMU a l a2+ writeMU a k a1+ _ -> case cmp a0 a2 of+ LT -> case cmp a2 a3 of+ LT -> case cmp a1 a3 of+ LT -> do writeMU a j a2+ writeMU a k a1+ _ -> do writeMU a l a1+ writeMU a j a2+ writeMU a k a3+ _ -> case cmp a0 a3 of+ LT -> do writeMU a l a1+ writeMU a j a3+ _ -> do writeMU a i a3+ writeMU a l a1+ writeMU a j a0+ _ -> case cmp a0 a3 of+ LT -> do writeMU a i a2+ case cmp a1 a3 of+ LT -> writeMU a k a1+ _ -> do writeMU a k a3+ writeMU a l a1+ writeMU a j a0+ _ -> case cmp a2 a3 of+ LT -> do writeMU a i a2+ writeMU a k a0+ writeMU a j a3+ writeMU a l a1+ _ -> do writeMU a j a2+ writeMU a k a0+ writeMU a i a3+ writeMU a l a1+ _ -> case cmp a0 a2 of+ LT -> case cmp a0 a3 of+ LT -> do writeMU a i a1+ writeMU a j a0+ case cmp a2 a3 of+ GT -> do writeMU a k a3+ writeMU a l a2+ _ -> return ()+ _ -> do case cmp a1 a3 of+ LT -> do writeMU a i a1+ writeMU a j a3+ _ -> writeMU a i a3+ writeMU a l a2+ writeMU a k a0+ _ -> case cmp a1 a2 of+ LT -> case cmp a2 a3 of+ LT -> do writeMU a i a1+ writeMU a j a2+ case cmp a0 a3 of+ LT -> writeMU a k a0+ _ -> do writeMU a k a3+ writeMU a l a0+ _ -> do case cmp a1 a3 of+ LT -> do writeMU a i a1+ writeMU a j a3+ _ -> writeMU a i a3+ writeMU a l a0+ _ -> case cmp a1 a3 of+ LT -> do writeMU a i a2+ case cmp a0 a3 of+ LT -> writeMU a k a0+ _ -> do writeMU a k a3+ writeMU a l a0+ _ -> case cmp a2 a3 of+ LT -> do writeMU a i a2+ writeMU a k a1+ writeMU a j a3+ writeMU a l a0+ _ -> do writeMU a i a3+ writeMU a l a0+ writeMU a j a2+ writeMU a k a1+{-# INLINE sort4ByIndex #-}
+ Data/Array/Vector/Algorithms/Radix.hs view
@@ -0,0 +1,207 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.Radix+-- Copyright : (c) 2008 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.Array.Vector.Algorithms.Radix (sort, Radix(..)) where++import Control.Monad+import Control.Monad.ST++import Data.Array.Vector+import Data.Array.Vector.Algorithms.Common++import Data.Bits+import Data.Int+import Data.Word+++import Foreign.Storable++class UA e => 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 + 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 = fromIntegral e + 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 + 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 + 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 + 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 #-}++-- | Sorts an array based on the Radix instance.+sort :: forall e s. Radix e => MUArr e s -> ST s ()+sort arr = do+ tmp <- newMU len+ count <- newMU (size e)+ prefix <- newMU (size e)+ go False arr tmp count prefix 0+ where+ len = lengthMU arr+ e :: e+ e = undefined+ go !swap src dst count prefix k+ | k < passes e = do zero 0 count+ countLoop 0 k src count+ writeMU prefix 0 0+ prefixLoop 1 0 count prefix+ moveLoop 0 k src dst prefix+ go (not swap) dst src count prefix (k+1)+ | otherwise = when swap (mcopyMU src dst 0 0 len)+ zero i a+ | i < size e = writeMU a i 0 >> zero (i+1) a+ | otherwise = return ()+ countLoop i k src count+ | i < len = readMU src i >>= inc count . radix k >> countLoop (i+1) k src count+ | otherwise = return ()+ prefixLoop i pi count prefix+ | i < size e = do ci <- readMU count (i-1)+ let pi' = pi + ci+ writeMU prefix i pi'+ prefixLoop (i+1) pi' count prefix+ | otherwise = return ()+ moveLoop i k src dst prefix+ | i < len = do srci <- readMU src i+ pf <- inc prefix (radix k srci)+ writeMU dst pf srci+ moveLoop (i+1) k src dst prefix+ | otherwise = return ()+{-# INLINE sort #-}++inc :: MUArr Int s -> Int -> ST s Int+inc arr i = readMU arr i >>= \e -> writeMU arr i (e+1) >> return e+{-# INLINE inc #-}
+ Data/Array/Vector/Algorithms/TriHeap.hs view
@@ -0,0 +1,189 @@+{-# LANGUAGE TypeOperators #-}++-- ---------------------------------------------------------------------------+-- |+-- Module : Data.Array.Vector.Algorithms.TriHeap+-- Copyright : (c) 2008 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.Array.Vector.Algorithms.TriHeap+ ( -- * Sorting+ sort+ , sortBy+ , sortByBounds+ -- * Selection+ , select+ , selectBy+ , selectByBounds+ -- * Partial sorts+ , partialSort+ , partialSortBy+ , partialSortByBounds+ -- * Heap operations+ , heapify+ , pop+ , popTo+ , sortHeap ) where++import Control.Monad+import Control.Monad.ST++import Data.Array.Vector+import Data.Array.Vector.Algorithms.Common++import qualified Data.Array.Vector.Algorithms.Optimal as O++-- | Sorts an entire array using the default ordering.+sort :: (UA e, Ord e) => MUArr e s -> ST s ()+sort = sortBy compare+{-# INLINE sort #-}++-- | Sorts an entire array using a custom ordering.+sortBy :: (UA e) => Comparison e -> MUArr e s -> ST s ()+sortBy cmp a = sortByBounds cmp a 0 (lengthMU a)+{-# INLINE sortBy #-}++-- | Sorts a portion of an array [l,u) using a custom ordering+sortByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+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 :: (UA e, Ord e) => MUArr e s -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+selectBy cmp a k = selectByBounds cmp a k 0 (lengthMU 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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+selectByBounds cmp a k l u = heapify cmp a l (l + k) >> go l (l + k) u+ where+ go l m u+ | u < m = return ()+ | otherwise = do el <- readMU a l+ eu <- readMU 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 :: (UA e, Ord e) => MUArr e s -> Int -> ST s ()+partialSort = partialSortBy compare+{-# INLINE partialSort #-}++-- | Moves the lowest k elements (as defined by the comparison) to+-- the front of the array, sorted.+partialSortBy :: (UA e) => Comparison e -> MUArr e s -> Int -> ST s ()+partialSortBy cmp a k = partialSortByBounds cmp a k 0 (lengthMU a)+{-# INLINE partialSortBy #-}++-- | Moves the lowest k elements in the portion [l,u) of the array+-- into positions [l,k+l), sorted.+partialSortByBounds :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+partialSortByBounds cmp a k l u = do selectByBounds cmp a k l u+ sortHeap cmp a l (l + 4) (l + k)+ O.sort4ByOffset cmp a l+ -- this technically does extra work for k < 4, but+ -- I'm not sure that's a significant concern.+{-# INLINE partialSortByBounds #-}++-- | Constructs a heap in a portion of an array [l, u)+heapify :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+heapify cmp a l u = loop $ (len - 1) `div` 3+ where+ len = u - l+ loop k+ | k < 0 = return ()+ | otherwise = readMU 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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> ST s ()+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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+popTo cmp a l u t = do al <- readMU a l+ at <- readMU a t+ writeMU 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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s ()+sortHeap cmp a l m u = loop (u-1) >> swap 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 :: (UA e) => Comparison e -> MUArr e s -> e -> Int -> Int -> Int -> ST s ()+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 -> writeMU a (root + off) ac >> sift val child' len+ _ -> writeMU a (root + off) val+ | otherwise = writeMU 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 :: (UA e) => Comparison e -> MUArr e s -> Int -> Int -> Int -> ST s (Int :*: e)+maximumChild cmp a off child1 len+ | child3 < len = do ac1 <- readMU a (child1 + off)+ ac2 <- readMU a (child2 + off)+ ac3 <- readMU 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 <- readMU a (child1 + off)+ ac2 <- readMU a (child2 + off)+ return $ case cmp ac1 ac2 of+ LT -> child2 :*: ac2+ _ -> child1 :*: ac1+ | otherwise = do ac1 <- readMU a (child1 + off) ; return (child1 :*: ac1)+ where+ child2 = child1 + 1+ child3 = child1 + 2+{-# INLINE maximumChild #-}
+ LICENSE view
@@ -0,0 +1,65 @@+Copyright (c) 2008 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
+ uvector-algorithms.cabal view
@@ -0,0 +1,40 @@+Name: uvector-algorithms+Version: 0.1+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 uvector unboxed arrays+Description: Efficient algorithms for uvector unboxed arrays+ be sure to compile with -O2, and -fvia-C -optc-O3 is+ recommended.+Build-Type: Simple+Cabal-Version: >= 1.2++Library+ Build-Depends: base, uvector++ Exposed-Modules:+ Data.Array.Vector.Algorithms.Immutable+ Data.Array.Vector.Algorithms.Optimal+ Data.Array.Vector.Algorithms.Insertion+ Data.Array.Vector.Algorithms.Intro+ Data.Array.Vector.Algorithms.Merge+ Data.Array.Vector.Algorithms.Radix+ Data.Array.Vector.Algorithms.TriHeap++ Other-Modules:+ Data.Array.Vector.Algorithms.Common++ Extensions:+ BangPatterns,+ TypeOperators,+ Rank2Types,+ ScopedTypeVariables++ GHC-Options:+ -O2+ -fvia-C -optc-O3+ -funbox-strict-fields