diff --git a/Data/Array/Vector/Algorithms/Common.hs b/Data/Array/Vector/Algorithms/Common.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Common.hs @@ -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 #-} diff --git a/Data/Array/Vector/Algorithms/Immutable.hs b/Data/Array/Vector/Algorithms/Immutable.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Immutable.hs @@ -0,0 +1,24 @@ +{-# LANGUAGE Rank2Types #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Array.Vector.Algorithms.Immutable +-- Copyright : (c) 2008 Dan Doel +-- Maintainer : Dan Doel +-- 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) diff --git a/Data/Array/Vector/Algorithms/Insertion.hs b/Data/Array/Vector/Algorithms/Insertion.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Insertion.hs @@ -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 #-} diff --git a/Data/Array/Vector/Algorithms/Intro.hs b/Data/Array/Vector/Algorithms/Intro.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Intro.hs @@ -0,0 +1,187 @@ +{-# LANGUAGE TypeOperators, BangPatterns #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Array.Vector.Algorithms.Intro +-- Copyright : (c) 2008 Dan Doel +-- Maintainer : Dan Doel +-- 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 #-} diff --git a/Data/Array/Vector/Algorithms/Merge.hs b/Data/Array/Vector/Algorithms/Merge.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Merge.hs @@ -0,0 +1,84 @@ +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Array.Vector.Algorithms.Merge +-- Copyright : (c) 2008 Dan Doel +-- Maintainer : Dan Doel +-- 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 #-} diff --git a/Data/Array/Vector/Algorithms/Optimal.hs b/Data/Array/Vector/Algorithms/Optimal.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Optimal.hs @@ -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 #-} diff --git a/Data/Array/Vector/Algorithms/Radix.hs b/Data/Array/Vector/Algorithms/Radix.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/Radix.hs @@ -0,0 +1,207 @@ +{-# LANGUAGE ScopedTypeVariables, BangPatterns #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Array.Vector.Algorithms.Radix +-- Copyright : (c) 2008 Dan Doel +-- Maintainer : Dan Doel +-- 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 #-} diff --git a/Data/Array/Vector/Algorithms/TriHeap.hs b/Data/Array/Vector/Algorithms/TriHeap.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Vector/Algorithms/TriHeap.hs @@ -0,0 +1,189 @@ +{-# LANGUAGE TypeOperators #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Array.Vector.Algorithms.TriHeap +-- Copyright : (c) 2008 Dan Doel +-- Maintainer : Dan Doel +-- 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 #-} diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/uvector-algorithms.cabal b/uvector-algorithms.cabal new file mode 100644 --- /dev/null +++ b/uvector-algorithms.cabal @@ -0,0 +1,40 @@ +Name: uvector-algorithms +Version: 0.1 +License: BSD3 +License-File: LICENSE +Author: Dan Doel +Maintainer: Dan Doel +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