diff --git a/Data/Array/CArray.hs b/Data/Array/CArray.hs new file mode 100644 --- /dev/null +++ b/Data/Array/CArray.hs @@ -0,0 +1,85 @@ +-- | +-- Module : Data.Array.CArray +-- Copyright : (c) 2008 Jed Brown +-- License : BSD-style +-- +-- Maintainer : jed@59A2.org +-- Stability : experimental +-- Portability : non-portable +-- +-- This module provides the immutable 'CArray' which uses pinned memory on the +-- GC'd heap. Elements are stored according to the class 'Storable'. You can +-- obtain a pointer to the array contents to manipulate elements from +-- languages like C. +-- +-- 'CArray' is 16-byte aligned by default. If you create a 'CArray' with +-- 'unsafeForeignPtrToCArray' then it may not be aligned. This will be an issue +-- if you intend to use SIMD instructions. +-- +-- 'CArray' is similar to 'Data.Array.Unboxed.UArray' but slower if you stay +-- within Haskell. 'CArray' can handle more types and can be used by external +-- libraries. +----------------------------------------------------------------------------- + +module Data.Array.CArray ( + -- * CArray type + CArray, + + -- * Foreign support + withCArray, + unsafeForeignPtrToCArray, + createCArray, + + -- * Multi-dimensional + + -- ** Fast reshaping + reshape, + flatten, + + -- ** Query + rank, + shape, + size, + + -- * Mapping + + -- ** General + ixmapWithIndP, + ixmapWithInd, + ixmapWithP, + ixmapWith, + ixmapP, + + -- ** Slicing + sliceStrideWithP, + sliceStrideWith, + sliceStrideP, + sliceStride, + sliceWithP, + sliceWith, + sliceP, + slice, + + -- * Lifting + liftArrayP, + liftArray, + liftArray2P, + liftArray2, + liftArray3P, + liftArray3, + + -- * Norms + normp, + norm2, + normSup, + + -- * Types + Shapable, + Abs, + + -- * The overloaded immutable array interface + module Data.Array.IArray +) where + +import Data.Array.IArray +import Data.Array.CArray.Base diff --git a/Data/Array/CArray/Base.hs b/Data/Array/CArray/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Array/CArray/Base.hs @@ -0,0 +1,618 @@ +{-# OPTIONS_GHC -frewrite-rules #-} +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, MagicHash, + FlexibleInstances, FlexibleContexts, UnboxedTuples, ScopedTypeVariables, + DeriveDataTypeable, CPP #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Array.CArray.Base +-- Copyright : (c) 2001 The University of Glasgow +-- (c) 2008 Jed Brown +-- License : BSD-style +-- +-- Maintainer : jed@59A2.org +-- Stability : experimental +-- Portability : non-portable +-- +-- This module provides both the immutable 'CArray' and mutable 'IOCArray'. The +-- underlying storage is exactly the same - pinned memory on the GC'd heap. +-- Elements are stored according to the class 'Storable'. You can obtain a +-- pointer to the array contents to manipulate elements from languages like C. +-- +-- 'CArray' is 16-byte aligned by default. If you create a 'CArray' with +-- 'unsafeForeignPtrToCArray' then it may not be aligned. This will be an issue +-- if you intend to use SIMD instructions. +-- +-- 'CArray' is similar to 'Data.Array.Unboxed.UArray' but slower if you stay +-- within Haskell. 'CArray' can handle more types and can be used by external +-- libraries. +-- +-- 'IOCArray' is equivalent to 'Data.Array.Storable.StorableArray' and similar +-- to 'Data.Array.IO.IOUArray' but slower. 'IOCArray' has O(1) versions of +-- 'unsafeFreeze' and 'unsafeThaw' when converting to/from 'CArray'. +----------------------------------------------------------------------------- + +module Data.Array.CArray.Base where + +import Control.Applicative +import Control.Monad +import Data.Array.Base +import Data.Array.MArray +import Data.Array.IArray +import Data.Complex +import Data.List +import System.IO.Unsafe (unsafePerformIO) +import Foreign.Storable +import Foreign.ForeignPtr +import Foreign.Ptr +import Foreign.Marshal.Alloc +import Foreign.Marshal.Array (copyArray,peekArray,pokeArray) +import Data.Word (Word8,Word) + +import Data.Generics (Data(..), Typeable(..)) +import GHC.Base (realWorld#, Addr#) +import GHC.IOBase (IO(..)) +import GHC.Ptr (Ptr(..)) +import GHC.ForeignPtr (ForeignPtr(..), mallocPlainForeignPtrBytes) + +-- | The immutable array type. +data CArray i e = CArray !i !i Int !(ForeignPtr e) + deriving (Data, Typeable) + +-- | Absolutely equivalent representation, but used for the mutable interface. +data IOCArray i e = IOCArray !i !i Int !(ForeignPtr e) + deriving (Data, Typeable) + +instance Storable e => MArray IOCArray e IO where + getBounds (IOCArray l u _ _) = return (l,u) + + getNumElements (IOCArray l u n _) = return n + + newArray (l,u) init = do + fp <- mallocForeignPtrArrayAligned size + withForeignPtr fp $ \a -> + sequence_ [pokeElemOff a i init | i <- [0..size-1]] + return (IOCArray l u size fp) + where size = rangeSize (l,u) + + unsafeNewArray_ (l,u) = do + let n = rangeSize (l,u) + fp <- mallocForeignPtrArrayAligned n + return (IOCArray l u n fp) + + newArray_ = unsafeNewArray_ + + unsafeRead (IOCArray _ _ _ fp) i = + withForeignPtr fp $ \a -> peekElemOff a i + + unsafeWrite (IOCArray _ _ _ fp) i e = + withForeignPtr fp $ \a -> pokeElemOff a i e + +-- |The pointer to the array contents is obtained by 'withCArray'. +-- The idea is similar to 'ForeignPtr' (used internally here). +-- The pointer should be used only during execution of the 'IO' action +-- retured by the function passed as argument to 'withCArray'. +withCArray :: CArray i e -> (Ptr e -> IO a) -> IO a +withCArray (CArray _ _ _ fp) f = withForeignPtr fp f + +withIOCArray :: IOCArray i e -> (Ptr e -> IO a) -> IO a +withIOCArray (IOCArray _ _ _ fp) f = withForeignPtr fp f + +-- |If you want to use it afterwards, ensure that you +-- 'touchCArray' after the last use of the pointer, +-- so the array is not freed too early. +touchIOCArray :: IOCArray i e -> IO () +touchIOCArray (IOCArray _ _ _ fp) = touchForeignPtr fp + +-- |Construct a 'CArray' from an arbitrary 'ForeignPtr'. It is +-- the caller's responsibility to ensure that the 'ForeignPtr' points to +-- an area of memory sufficient for the specified bounds. +unsafeForeignPtrToCArray + :: Ix i => ForeignPtr e -> (i,i) -> IO (CArray i e) +unsafeForeignPtrToCArray p (l,u) = + return (CArray l u (rangeSize (l,u)) p) + +unsafeForeignPtrToIOCArray + :: Ix i => ForeignPtr e -> (i,i) -> IO (IOCArray i e) +unsafeForeignPtrToIOCArray p (l,u) = + return (IOCArray l u (rangeSize (l,u)) p) + +copy :: (Ix i, Storable e) => CArray i e -> IO (CArray i e) +copy ain@(CArray l u n fp) = + createCArray (l,u) $ \op -> + withCArray ain $ \ip -> + copyArray op ip n + +freezeIOCArray :: (Ix i, Storable e) => IOCArray i e -> IO (CArray i e) +freezeIOCArray = unsafeFreezeIOCArray >=> copy + +unsafeFreezeIOCArray :: (Ix i) => IOCArray i e -> IO (CArray i e) +unsafeFreezeIOCArray (IOCArray l u n fp) = return (CArray l u n fp) + +thawIOCArray :: (Ix i, Storable e) => CArray i e -> IO (IOCArray i e) +thawIOCArray = copy >=> unsafeThawIOCArray + +unsafeThawIOCArray :: (Ix i) => CArray i e -> IO (IOCArray i e) +unsafeThawIOCArray (CArray l u n fp) = return (IOCArray l u n fp) + +-- Since we can remove the (Storable e) restriction for these, the rules are +-- compact and general. +{-# RULES +"unsafeFreeze/IOCArray" unsafeFreeze = unsafeFreezeIOCArray +"unsafeThaw/IOCArray" unsafeThaw = unsafeThawIOCArray + #-} + +-- Since we can't parameterize the rules with the (Storable e) constraint, we +-- have to specialize manually. This is unfortunate since it is less general. +{-# RULES +"freeze/IOCArray/Int" freeze = freezeIOCArray :: (Ix i) => IOCArray i Int -> IO (CArray i Int) +"freeze/IOCArray/Float" freeze = freezeIOCArray :: (Ix i) => IOCArray i Float -> IO (CArray i Float) +"freeze/IOCArray/Double" freeze = freezeIOCArray :: (Ix i) => IOCArray i Double -> IO (CArray i Double) +"thaw/IOCArray/Int" thaw = thawIOCArray :: (Ix i) => CArray i Int -> IO (IOCArray i Int) +"thaw/IOCArray/Float" thaw = thawIOCArray :: (Ix i) => CArray i Float -> IO (IOCArray i Float) +"thaw/IOCArray/Double" thaw = thawIOCArray :: (Ix i) => CArray i Double -> IO (IOCArray i Double) + #-} + + +instance Storable e => IArray CArray e where + {-# INLINE bounds #-} + bounds (CArray l u _ _) = (l,u) + {-# INLINE numElements #-} + numElements (CArray _ _ n _) = n + {-# NOINLINE unsafeArray #-} + unsafeArray lu ies = unsafePerformIO $ unsafeArrayCArray lu ies (zeroElem (undefined :: e)) + {-# INLINE unsafeAt #-} + unsafeAt (CArray _ _ _ fp) i = unsafeInlinePerformIO $ + withForeignPtr fp $ \a -> peekElemOff a i + {-# NOINLINE unsafeReplace #-} + unsafeReplace arr ies = unsafePerformIO $ unsafeReplaceCArray arr ies + {-# NOINLINE unsafeAccum #-} + unsafeAccum f arr ies = unsafePerformIO $ unsafeAccumCArray f arr ies + {-# NOINLINE unsafeAccumArray #-} + unsafeAccumArray f init lu ies = unsafePerformIO $ unsafeAccumArrayCArray f init lu ies + + +-- | Hackish way to get the zero element for a Storable type. +{-# NOINLINE zeroElem #-} +zeroElem :: Storable a => a -> a +zeroElem u = unsafePerformIO $ do + allocaBytes size $ \p -> do + sequence_ [pokeByteOff p off (0 :: Word8) | off <- [0 .. (size - 1)] ] + peek (castPtr p) + where size = sizeOf u + +{-# INLINE unsafeArrayCArray #-} +unsafeArrayCArray :: (MArray IOCArray e IO, Storable e, Ix i) + => (i,i) -> [(Int, e)] -> e -> IO (CArray i e) +unsafeArrayCArray lu ies default_elem = do + marr <- newArray lu default_elem + sequence_ [unsafeWrite marr i e | (i, e) <- ies] + unsafeFreezeIOCArray marr + +{-# INLINE unsafeReplaceCArray #-} +unsafeReplaceCArray :: (MArray IOCArray e IO, Storable e, Ix i) + => CArray i e -> [(Int, e)] -> IO (CArray i e) +unsafeReplaceCArray arr ies = do + marr <- thawIOCArray arr + sequence_ [unsafeWrite marr i e | (i, e) <- ies] + unsafeFreezeIOCArray marr + +{-# INLINE unsafeAccumCArray #-} +unsafeAccumCArray :: (MArray IOCArray e IO, Storable e, Ix i) + => (e -> e' -> e) -> CArray i e -> [(Int, e')] + -> IO (CArray i e) +unsafeAccumCArray f arr ies = do + marr <- thawIOCArray arr + sequence_ [do + old <- unsafeRead marr i + unsafeWrite marr i (f old new) + | (i, new) <- ies] + unsafeFreezeIOCArray marr + +{-# INLINE unsafeAccumArrayCArray #-} +unsafeAccumArrayCArray :: (MArray IOCArray e IO, Storable e, Ix i) + => (e -> e' -> e) -> e -> (i,i) -> [(Int, e')] + -> IO (CArray i e) +unsafeAccumArrayCArray f init lu ies = do + marr <- newArray lu init + sequence_ [do + old <- unsafeRead marr i + unsafeWrite marr i (f old new) + | (i, new) <- ies] + unsafeFreezeIOCArray marr + +{-# INLINE eqCArray #-} +eqCArray :: (IArray CArray e, Ix i, Eq e) + => CArray i e -> CArray i e -> Bool +eqCArray arr1@(CArray l1 u1 n1 _) arr2@(CArray l2 u2 n2 _) = + if n1 == 0 then n2 == 0 else + l1 == l2 && u1 == u2 && + and [unsafeAt arr1 i == unsafeAt arr2 i | i <- [0 .. n1 - 1]] + +{-# INLINE cmpCArray #-} +cmpCArray :: (IArray CArray e, Ix i, Ord e) + => CArray i e -> CArray i e -> Ordering +cmpCArray arr1 arr2 = compare (assocs arr1) (assocs arr2) + +{-# INLINE cmpIntCArray #-} +cmpIntCArray :: (IArray CArray e, Ord e) + => CArray Int e -> CArray Int e -> Ordering +cmpIntCArray arr1@(CArray l1 u1 n1 _) arr2@(CArray l2 u2 n2 _) = + if n1 == 0 then if n2 == 0 then EQ else LT else + if n2 == 0 then GT else + case compare l1 l2 of + EQ -> foldr cmp (compare u1 u2) [0 .. (n1 `min` n2) - 1] + other -> other + where + cmp i rest = case compare (unsafeAt arr1 i) (unsafeAt arr2 i) of + EQ -> rest + other -> other + +{-# RULES "cmpCArray/Int" cmpCArray = cmpIntCArray #-} + +instance (Ix ix, Eq e, IArray CArray e) => Eq (CArray ix e) where + (==) = eqCArray + +instance (Ix ix, Ord e, IArray CArray e) => Ord (CArray ix e) where + compare = cmpCArray + +instance (Ix ix, Show ix, Show e, IArray CArray e) => Show (CArray ix e) where + showsPrec = showsIArray + +-- +-- General purpose array operations which happen to be very fast for CArray. +-- + +-- | O(1) reshape an array. The number of elements in the new shape must not +-- exceed the number in the old shape. The elements are in C-style ordering. +reshape :: (Ix i, Ix j) => (j,j) -> CArray i e -> CArray j e +reshape (l',u') (CArray l u n fp) | n' > n = error "reshape: new size too large" + |otherwise = CArray l' u' n' fp + where n' = rangeSize (l', u') + +-- | O(1) make a rank 1 array from an arbitrary shape. +-- It has the property that 'reshape (0, size a - 1) a == flatten a'. +flatten :: Ix i => CArray i e -> CArray Int e +flatten (CArray l u n fp) = CArray 0 (n - 1) n fp + +-- | Determine the rank of an array. +rank :: (Shapable i, Ix i, IArray a e) => a i e -> Int +rank = sRank . fst . bounds + +-- +-- Useful for multi-dimensional arrays. Not specific to CArray. + +-- | Canonical representation of the shape. +-- The following properties hold: +-- 'length . shape = rank' +-- 'product . shape = size' +shape :: (Shapable i, Ix i, IArray a e) => a i e -> [Int] +shape = uncurry sShape . bounds + +-- | How much the offset changes when you move one element in the given +-- direction. Since arrays are in row-major order, 'last . shapeToStride = const 1' +shapeToStride :: [Int] -> [Int] +shapeToStride = scanr (*) 1 . tail + +-- | Number of elements in the Array. +size :: (Ix i, IArray a e) => a i e -> Int +size = numElements + +-- +-- None of the following are specific to CArray. Some could have slightly +-- faster versions specialized to CArray. In general, slicing is expensive +-- because the slice is not contiguous in memory, so must be copied. There are +-- many specialized versions. +-- + +-- | Generic slice and map. This takes the new range, the inverse map on +-- indices, and function to produce the next element. It is the most general +-- operation in its class. +ixmapWithIndP :: (Ix i, Ix i', IArray a e, IArray a' e') + => (i',i') -> (i' -> i) -> (i -> e -> i' -> e') -> a i e -> a' i' e' +ixmapWithIndP lu f g arr = listArray lu + [ let i = f i' in g i (arr ! i) i' | i' <- range lu ] + +-- | Less polymorphic version. +ixmapWithInd :: (Ix i, Ix i', IArray a e, IArray a e') + => (i',i') -> (i' -> i) -> (i -> e -> i' -> e') -> a i e -> a i' e' +ixmapWithInd = ixmapWithIndP + +-- | Perform an operation on the elements, independent of their location. +ixmapWithP :: (Ix i, Ix i', IArray a e, IArray a' e') + => (i',i') -> (i' -> i) -> (e -> e') -> a i e -> a' i' e' +ixmapWithP lu f g arr = listArray lu + [ g (arr ! f i') | i' <- range lu ] + +-- | Less polymorphic version. +ixmapWith :: (Ix i, Ix i', IArray a e, IArray a e') + => (i',i') -> (i' -> i) -> (e -> e') -> a i e -> a i' e' +ixmapWith = ixmapWithP + +-- | More polymorphic version of 'ixmap'. +ixmapP :: (Ix i, Ix i', IArray a e, IArray a' e) + => (i',i') -> (i' -> i) -> a i e -> a' i' e +ixmapP lu f arr = ixmapWithP lu f id arr + +-- | More friendly sub-arrays with element mapping. +sliceStrideWithP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e') + => (i',i') -> (i,i,i) -> (e -> e') -> a i e -> a' i' e' +sliceStrideWithP lu (start,next,end) f arr + | all (inRange (bounds arr)) [start,next,end] = listArray lu elems + | otherwise = error "sliceStrideWith: out of bounds" + where is = offsetShapeFromThenTo (shape arr) (index' start) (index' next) (index' end) + elems = map (f . (unsafeAt arr)) is + index' = indexes arr + +-- | Less polymorphic version. +sliceStrideWith :: (Ix i, Shapable i, Ix i', IArray a e, IArray a e') + => (i',i') -> (i,i,i) -> (e -> e') -> a i e -> a i' e' +sliceStrideWith = sliceStrideWithP + +-- | Strided sub-array without element mapping. +sliceStrideP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e) + => (i',i') -> (i,i,i) -> a i e -> a' i' e +sliceStrideP lu sne = sliceStrideWithP lu sne id + +-- | Less polymorphic version. +sliceStride :: (Ix i, Shapable i, Ix i', IArray a e) + => (i',i') -> (i,i,i) -> a i e -> a i' e +sliceStride = sliceStrideP + +-- | Contiguous sub-array with element mapping. +sliceWithP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e') + => (i',i') -> (i,i) -> (e -> e') -> a i e -> a' i' e' +sliceWithP lu (start,end) f arr + | all (inRange (bounds arr)) [start,end] = listArray lu elems + | otherwise = error "sliceWith: out of bounds" + where is = offsetShapeFromTo (shape arr) (index' start) (index' end) + elems = map (f . (unsafeAt arr)) is + index' = indexes arr + +-- | Less polymorphic version. +sliceWith :: (Ix i, Shapable i, Ix i', IArray a e, IArray a e') + => (i',i') -> (i,i) -> (e -> e') -> a i e -> a i' e' +sliceWith = sliceWithP + +-- | Contiguous sub-array without element mapping. +sliceP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e) + => (i',i') -> (i,i) -> a i e -> a' i' e +sliceP lu se = sliceWithP lu se id + +-- | Less polymorphic version. +slice :: (Ix i, Shapable i, Ix i', IArray a e) + => (i',i') -> (i,i) -> a i e -> a i' e +slice = sliceP + +-- | In-place map on CArray. Note that this is /IN PLACE/ so you should not +-- retain any reference to the original. It flagrantly breaks referential +-- transparency! +{-# INLINE mapCArrayInPlace #-} +mapCArrayInPlace :: (Ix i, IArray CArray e, Storable e) => (e -> e) -> CArray i e -> CArray i e +mapCArrayInPlace f a = unsafeInlinePerformIO $ do + withCArray a $ \p -> + forM_ [0 .. size a - 1] $ \i -> + peekElemOff p i >>= pokeElemOff p i . f + return a + +----------------------------------------- +-- These are meant to be internal only +indexes :: (Ix i, Shapable i, IArray a e) => a i e -> i -> [Int] +indexes a i = map pred $ (sShape . fst . bounds) a i + +offsetShapeFromThenTo s a b c = foldr (liftA2 (+)) [0] (ilists stride a b c) + where ilists = zipWith4 (\s a b c -> map (*s) $ enumFromThenTo a b c) + stride = shapeToStride s + +offsetShapeFromTo = offsetShapeFromTo' id + +offsetShapeFromTo' f s a b = foldr (liftA2 (+)) [0] (f $ ilists stride a b) + where ilists = zipWith3 (\s a b -> map (*s) $ enumFromTo a b) + stride = shapeToStride s + +offsets :: (Ix a, Shapable a) => (a, a) -> a -> [Int] +offsets lu i = reverse . osets (index lu i) . reverse . scanl1 (*) . uncurry sShape $ lu + where osets 0 [] = [] + osets i (b:bs) = r : osets d bs + where (d,r) = i `divMod` b + osets _ _ = error "osets" +----------------------------------------- + +-- | p-norm on the array taken as a vector +normp :: (Ix i, RealFloat e', Abs e e', IArray a e) => e' -> a i e -> e' +normp p a | 1 <= p && not (isInfinite p) = (** (1/p)) $ foldl' (\z e -> z + (abs_ e) ** p) 0 (elems a) + | otherwise = error "normp: p < 1" + +-- | 2-norm on the array taken as a vector (Frobenius norm for matrices) +norm2 :: (Ix i, Floating e', Abs e e', IArray a e) => a i e -> e' +norm2 a = sqrt $ foldl' (\z e -> z + abs_ e ^ 2) 0 (elems a) + +-- | Sup norm on the array taken as a vector +normSup :: (Ix i, Num e', Ord e', Abs e e', IArray a e) => a i e -> e' +normSup a = foldl' (\z e -> z `max` abs_ e) 0 (elems a) + +-- | Polymorphic version of amap. +liftArrayP :: (Ix i, IArray a e, IArray a1 e1) + => (e -> e1) -> a i e -> a1 i e1 +liftArrayP f a = listArray (bounds a) (map f (elems a)) + +-- | Equivalent to amap. Here for consistency only. +liftArray :: (Ix i, IArray a e, IArray a e1) + => (e -> e1) -> a i e -> a i e1 +liftArray = liftArrayP + +-- | Polymorphic 2-array lift. +liftArray2P :: (Ix i, IArray a e, IArray a1 e1, IArray a2 e2) + => (e -> e1 -> e2) -> a i e -> a1 i e1 -> a2 i e2 +liftArray2P f a b | aBounds == bounds b = + listArray aBounds (zipWith f (elems a) (elems b)) + | otherwise = error "liftArray2: array bounds must match" + where aBounds = bounds a + +-- | Less polymorphic version. +liftArray2 :: (Ix i, IArray a e, IArray a e1, IArray a e2) + => (e -> e1 -> e2) -> a i e -> a i e1 -> a i e2 +liftArray2 = liftArray2P + +-- | Polymorphic 3-array lift. +liftArray3P :: (Ix i, IArray a e, IArray a1 e1, IArray a2 e2, IArray a3 e3) + => (e -> e1 -> e2 -> e3) -> a i e -> a1 i e1 -> a2 i e2 -> a3 i e3 +liftArray3P f a b c | aBounds == bounds b && aBounds == bounds c = + listArray aBounds (zipWith3 f (elems a) (elems b) (elems c)) + | otherwise = error "liftArray2: array bounds must match" + where aBounds = bounds a + +-- | Less polymorphic version. +liftArray3 :: (Ix i, IArray a e, IArray a e1, IArray a e2, IArray a e3) + => (e -> e1 -> e2 -> e3) -> a i e -> a i e1 -> a i e2 -> a i e3 +liftArray3 = liftArray3P + + +-- | We need this type class to distinguish between different tuples of Ix. +-- There are Shapable instances for homogenous Int tuples, but may Haddock +-- doesn't see them. +class Shapable i where + sRank :: i -> Int + sShape :: i -> i -> [Int] + sBounds :: [Int] -> (i,i) + +instance Shapable Int where + sRank _ = 1 + sShape a a' = [rangeSize (a,a')] + sBounds [a] = (0,a-1) + +instance Shapable (Int,Int) where + sRank _ = 2 + sShape (a,b) (a',b') = [rangeSize (a,a'), rangeSize (b,b')] + sBounds [a,b] = ((0,0),(a-1,b-1)) + +instance Shapable (Int,Int,Int) where + sRank _ = 3 + sShape (a,b,c) (a',b',c') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c')] + sBounds [a,b,c] = ((0,0,0),(a-1,b-1,c-1)) + +instance Shapable (Int,Int,Int,Int) where + sRank _ = 4 + sShape (a,b,c,d) (a',b',c',d') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c'), rangeSize (d,d')] + sBounds [a,b,c,d] = ((0,0,0,0),(a-1,b-1,c-1,d-1)) + +instance Shapable (Int,Int,Int,Int,Int) where + sRank _ = 5 + sShape (a,b,c,d,e) (a',b',c',d',e') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c'), rangeSize (d,d') + , rangeSize (e,e')] + sBounds [a,b,c,d,e] = ((0,0,0,0,0),(a-1,b-1,c-1,d-1,e-1)) + +instance Shapable (Int,Int,Int,Int,Int,Int) where + sRank _ = 6 + sShape (a,b,c,d,e,f) (a',b',c',d',e',f') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c'), rangeSize (d,d') + , rangeSize (e,e'), rangeSize (f,f')] + sBounds [a,b,c,d,e,f] = ((0,0,0,0,0,0),(a-1,b-1,c-1,d-1,e-1,f-1)) + +instance Shapable (Int,Int,Int,Int,Int,Int,Int) where + sRank _ = 7 + sShape (a,b,c,d,e,f,g) (a',b',c',d',e',f',g') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c'), rangeSize (d,d') + , rangeSize (e,e'), rangeSize (f,f'), rangeSize (g,g')] + sBounds [a,b,c,d,e,f,g] = ((0,0,0,0,0,0,0),(a-1,b-1,c-1,d-1,e-1,f-1,g-1)) + +instance Shapable (Int,Int,Int,Int,Int,Int,Int,Int) where + sRank _ = 8 + sShape (a,b,c,d,e,f,g,h) (a',b',c',d',e',f',g',h') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c'), rangeSize (d,d') + , rangeSize (e,e'), rangeSize (f,f'), rangeSize (g,g'), rangeSize (h,h')] + sBounds [a,b,c,d,e,f,g,h] = ((0,0,0,0,0,0,0,0),(a-1,b-1,c-1,d-1,e-1,f-1,g-1,h-1)) + +instance Shapable (Int,Int,Int,Int,Int,Int,Int,Int,Int) where + sRank _ = 9 + sShape (a,b,c,d,e,f,g,h,i) (a',b',c',d',e',f',g',h',i') = + [rangeSize (a,a'), rangeSize (b,b'), rangeSize (c,c'), rangeSize (d,d') + , rangeSize (e,e'), rangeSize (f,f'), rangeSize (g,g'), rangeSize (h,h') + , rangeSize (i,i')] + sBounds [a,b,c,d,e,f,g,h,i] = ((0,0,0,0,0,0,0,0,0) + ,(a-1,b-1,c-1,d-1,e-1,f-1,g-1,h-1,i-1)) + + +-- Represent the Complex value so that it conforms to the C99, C++, and FFTW +-- Complex format. This instance should probably be in the base libs. +instance (RealFloat a, Storable a) => Storable (Complex a) where + sizeOf _ = 2 * sizeOf (undefined :: a) + alignment _ = sizeOf (undefined :: a) + peek p = do + [r,i] <- peekArray 2 (castPtr p) + return (r :+ i) + poke p (r :+ i) = pokeArray (castPtr p) [r,i] + + +-- | Hack so that norms have a sensible type. +class Abs a b | a -> b where + abs_ :: a -> b +instance Abs (Complex Double) Double where + abs_ = magnitude +instance Abs (Complex Float) Float where + abs_ = magnitude +instance Abs Double Double where + abs_ = abs +instance Abs Float Float where + abs_ = abs + + +-- | This variant of 'unsafePerformIO' is quite /mind-bogglingly unsafe/. It +-- unstitches the dependency chain that holds the IO monad together and breaks +-- all your ordinary intuitions about IO, sequencing and side effects. Avoid +-- it unless you really know what you are doing. +-- +-- It is only safe for operations which are genuinely pure (not just +-- externally pure) for example reading from an immutable foreign data +-- structure. In particular, you should do no memory allocation inside an +-- 'unsafeInlinePerformIO' block. This is because an allocation is a constant +-- and is likely to be floated out and shared. More generally, any part of any +-- IO action that does not depend on a function argument is likely to be +-- floated to the top level and have its result shared. +-- +-- It is more efficient because in addition to the checks that +-- 'unsafeDupablePerformIO' omits, we also inline. Additionally we do not +-- pretend that the body is lazy which allows the strictness analyser to see +-- the strictness in the body. In turn this allows some re-ordering of +-- operations and any corresponding side-effects. +-- +-- With GHC it compiles to essentially no code and it exposes the body to +-- further inlining. +-- +{-# INLINE unsafeInlinePerformIO #-} +unsafeInlinePerformIO :: IO a -> a +#ifdef __GLASGOW_HASKELL__ +unsafeInlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r +#else +unsafeInlinePerformIO = unsafePerformIO +#endif + +-- | Allocate an array which is 16-byte aligned. Essential for SIMD instructions. +mallocForeignPtrArrayAligned :: Storable a => Int -> IO (ForeignPtr a) +mallocForeignPtrArrayAligned n = doMalloc undefined n + where + doMalloc :: Storable b => b -> Int -> IO (ForeignPtr b) + doMalloc dummy size = mallocForeignPtrBytesAligned (size * sizeOf dummy) + +-- | Allocate memory which is 16-byte aligned. This is essential for SIMD +-- instructions. We know that mallocPlainForeignPtrBytes will give word-aligned +-- memory, so we pad enough to be able to return the desired amount of memory +-- after aligning our pointer. +mallocForeignPtrBytesAligned :: Int -> IO (ForeignPtr a) +mallocForeignPtrBytesAligned n = do + (ForeignPtr addr contents) <- mallocPlainForeignPtrBytes (n + pad) + let (Ptr addr') = alignPtr (Ptr addr) 16 + return (ForeignPtr addr' contents) + where pad = 16 - sizeOf (undefined :: Word) + +-- | Make a new CArray with an IO action. +createCArray :: (Ix i, Storable e) => (i,i) -> (Ptr e -> IO ()) -> IO (CArray i e) +createCArray lu f = do + fp <- mallocForeignPtrArrayAligned (rangeSize lu) + withForeignPtr fp f + unsafeForeignPtrToCArray fp lu + +unsafeCreateCArray :: (Ix i, Storable e) => (i,i) -> (Ptr e -> IO ()) -> CArray i e +unsafeCreateCArray lu = unsafePerformIO . createCArray lu diff --git a/Data/Array/IOCArray.hs b/Data/Array/IOCArray.hs new file mode 100644 --- /dev/null +++ b/Data/Array/IOCArray.hs @@ -0,0 +1,40 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Array.IOCArray +-- Copyright : (c) 2008 Jed Brown +-- License : BSD-style +-- +-- Maintainer : jed@59A2.org +-- Stability : experimental +-- Portability : non-portable +-- +-- This module provides both the mutable 'IOCArray' which uses pinned memory on +-- the GC'd heap. Elements are stored according to the class 'Storable'. You +-- can obtain a pointer to the array contents to manipulate elements from +-- languages like C. +-- +-- 'IOCArray' is 16-byte aligned by default. If you create a 'IOCArray' with +-- 'unsafeForeignPtrToIOCArray' then it may not be aligned. This will be an +-- issue if you intend to use SIMD instructions. +-- +-- 'IOCArray' is equivalent to 'Data.Array.Storable.StorableArray' and similar +-- to 'Data.Array.IO.IOUArray' but slower. 'IOCArray' has O(1) versions of +-- 'unsafeFreeze' and 'unsafeThaw' when converting to/from 'CArray'. +----------------------------------------------------------------------------- + + +module Data.Array.IOCArray ( + -- * IOCArray type + IOCArray, + + -- * Foreign support + withIOCArray, + touchIOCArray, + unsafeForeignPtrToIOCArray, + + -- * The overloaded mutable array interface + module Data.Array.MArray +) where + +import Data.Array.CArray.Base +import Data.Array.MArray diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2008 Jed Brown + +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. diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,29 @@ +This package provides immutable and mutable arrays that can be used in foreign +calls. They are 16-byte aligned by default to facilitate use of SIMD +instructions. To build this package, use: + + runhaskell Setup.lhs configure + runhaskell Setup.lhs build + runhaskell Setup.lhs haddock (optional) + runhaskell Setup.lhs install + +Then run the tests: + + cd tests + ghc -O2 --make tests.hs -o tests && ./tests # checks QC properties + +In addition, there are versions two of shootout entries which use arrays. +Modified versions of these are in the tests directory, using various array +implementations. To build, benchmark, and check that results match, run: + + ./runtests.sh + + +Exposed Modules: + +Data.Array.CArray Immutable interface, enhanced for foreign calls, + multiple dimensions, mapping, and norms. + +Data.Array.IOCArray Mutable interface, enhanced for foreign calls + +Data.Array.CArray.Base Internals 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/carray.cabal b/carray.cabal new file mode 100644 --- /dev/null +++ b/carray.cabal @@ -0,0 +1,23 @@ +name: carray +version: 0.1.0.0 +synopsis: A C-compatible array library. +description: + A C-compatible array library. + . + Provides both an immutable and mutable (in the IO monad) interface. + Includes utilities for multi-dimensional arrays, slicing and norms. + Memory is 16-byte aligned by default to enable use of SIMD instructions. + . +category: Data +license: BSD3 +license-file: LICENSE +author: Jed Brown +maintainer: Jed Brown +stability: experimental +build-depends: base, array +build-type: Simple +exposed-modules: Data.Array.CArray + Data.Array.IOCArray + Data.Array.CArray.Base +extensions: +ghc-options: diff --git a/tests/meteor-contest-c.hs b/tests/meteor-contest-c.hs new file mode 100644 --- /dev/null +++ b/tests/meteor-contest-c.hs @@ -0,0 +1,268 @@ +{-# OPTIONS -O2 -fbang-patterns -optc-O3 #-} + +-- The Computer Language Benchmarks Game +-- http://shootout.alioth.debian.org/ +-- +-- Sterling Clover's translation of Tim Hochberg's Clean implementation + +module Main where +import System.Environment +import Data.Bits +import Data.List +import Data.Array.CArray +import Control.Arrow + +--- The Board --- +n_elem = 5 +n_col = 5 +n_row = 10 + +m_top :: Mask +m_top = 0x1F + +cells :: [Cell] +cells = [0..49] + +colors :: [Color] +colors = [0..9] + +cellAt x y = x + n_col * y +coordOf i = snd &&& fst $ i `quotRem` n_col +isValid x y = 0 <= x && x < n_col && 0 <= y && y < n_row + +--- Piece Operations --- +data Direction = E | SE | SW | W | NW | NE deriving (Enum, Eq, Ord) +type Piece = [Direction] +type CellCoord = (Int, Int) +type Mask = Int; type Color = Int; type Row = Int; +type Col = Int; type Tag = Int; type Cell = Int +type Solution = [Mask] + +pieces :: Array Int Piece +pieces = array (0,9) $ zip [0..9] $ + [[E, E, E, SE], + [SE, SW, W, SW], + [W, W, SW, SE], + [E, E, SW, SE], + [NW, W, NW, SE, SW], + [E, E, NE, W], + [NW, NE, NE, W], + [NE, SE, E, NE], + [SE, SE, E, SE], + [E, NW, NW, NW]] + +permutations :: Piece -> [Piece] +permutations p = take 12 (perms p) + where + perms p = p:(flip p) : perms (rotate p) + rotate piece = map r piece + where r E = NE + r NE = NW + r NW = W + r W = SW + r SW = SE + r SE = E + flip piece = map f piece + where f E = W + f NE = NW + f NW = NE + f W = E + f SW = SE + f SE = SW + +--- Mask Operations ---- +untag :: Mask -> Mask +untag mask = mask .&. 0x1ffffff + +retag :: Mask -> Tag -> Mask +retag mask n = untag mask .|. n `shiftL` 25 + +tagof :: Mask -> Tag +tagof mask = mask `shiftR` 25 + +tag :: Mask -> Tag -> Mask +tag mask n = mask .|. n `shiftL` 25 + +count1s :: Mask -> Int +count1s i + | i == 0 = 0 + | i .&. 1 == 1 = 1 + count1s (i `shiftR` 1) + | otherwise = count1s (i `shiftR` 1) + +first0 :: Mask -> Int +first0 i + | i .&. 1 == 0 = 0 + | otherwise = 1 + first0 (i `shiftR` 1) + +--- Making the Bitmasks --- +mod2 x = x .&. 1 +packSize a b = a*5+b +unpackSize n = quotRem n 5 + +move :: Direction -> CellCoord -> CellCoord +move E (x, y) = (x+1, y) +move W (x, y) = (x-1, y) +move NE (x, y) = (x+(mod2 y), y-1) +move NW (x, y) = (x+(mod2 y)-1, y-1) +move SE (x, y) = (x+(mod2 y), y+1) +move SW (x, y) = (x+(mod2 y)-1, y+1) + +pieceBounds :: Piece -> Bool -> (Int, Int, Int, Int) +pieceBounds piece isodd = bnds piece 0 y0 0 y0 0 y0 + where + y0 | isodd = 1 | otherwise = 0 + bnds [] _ _ xmin ymin xmax ymax = (xmin, ymin, xmax, ymax) + bnds (d:rest) x y xmin ymin xmax ymax = + bnds rest x' y' (min x' xmin) (min y' ymin) (max x' xmax) (max y' ymax) + where (x', y') = move d (x, y) + +pieceMask :: Piece -> (Mask, Mask) +pieceMask piece + | odd y1 = (tag (msk piece x2 y2 0) (packSize w2 h2), + tag (msk piece x1 (y1+1) 0 `shiftR` n_col) (packSize w1 h1)) + | otherwise = (tag (msk piece x1 y1 0) (packSize w1 h1), + tag (msk piece x2 (y2+1) 0 `shiftR` n_col) (packSize w2 h2)) + where + (xmin, ymin, xmax, ymax) = pieceBounds piece False + (x1, y1) = (-xmin, -ymin) + w1 = xmax - xmin + h1 = ymax - ymin + (xmin', ymin', xmax', ymax') = pieceBounds piece True + (x2, y2) = (-xmin', (-ymin')+1) + w2 = xmax' - xmin' + h2 = ymax' - ymin' + msk :: Piece -> Col -> Row -> Mask -> Mask + msk [] x y m = m `setBit` cellAt x y + msk (d:rest) x y m = msk rest x' y' (m `setBit` cellAt x y) + where (x', y') = move d (x, y) + +templatesForColor :: Color -> ([Mask], [Mask]) +templatesForColor c = (unzip . map pieceMask) perms + where perms | c == 5 = take 6 ps | otherwise = ps + ps = permutations $ pieces ! c + +--- Looking for Islands --- +noLineIslands :: Mask -> Cell -> Cell -> Int -> Bool +noLineIslands mask start stop step + | (fnd testBit . fnd ((not .) . testBit) . fnd testBit) start > stop = True + | otherwise = False + where + fnd test !x + | x >= 25 = 25 + | test mask x = x + | otherwise = fnd test (x+step) + +noLeftIslands :: Mask -> Bool +noLeftIslands mask = noLineIslands mask 0 20 5 +noRightIslands mask = noLineIslands mask 4 24 5 + +noIslands :: Mask -> Bool +noIslands board = noisles board (count1s board) + +noisles :: Mask -> Int -> Bool +noisles _ 30 = True +noisles board ones + | (ones' - ones) `rem` n_elem /= 0 = False + | otherwise = noisles board' ones' + where board' = fill board (coordOf (first0 board)) + ones' = count1s board' + +fill :: Mask -> CellCoord -> Mask +fill m cc@(x, y) + | x < 0 || x >= n_col = m + | y < 0 || y >= 6 = m + | testBit m i = m + | otherwise = foldl (\m d -> fill m (move d cc)) (setBit m i) + [E, NE, NW, W, SW, SE] + where i = cellAt x y + +--- More Mask Generation --- +masksForColor :: Color -> [(Row, Mask)] +masksForColor c = concatMap atCell cells + where + (evens, odds) = templatesForColor c + atCell n + | even y = [(y, retag (m `shiftL` x) c) | m <- evens , isok m x y] + | odd y = [(y, retag (m `shiftL` x) c) | m <- odds , isok m x y] + where (x, y) = coordOf n + +isok :: Mask -> Row -> Col -> Bool +isok mask x y = + isValid (x+width) (y+height) && + case (y == 0, y+height==9) of + (False, False) -> noLeftIslands mask' && noRightIslands mask' + (False, True) -> noIslands (mask' `shiftL` (n_col * (y - 4))) + (True, _ ) -> noIslands mask' + where (width, height) = unpackSize (tagof mask) + mask' = untag mask `shiftL` x + +masksAtCell :: Array (Row,Col) (Array Color [Mask]) +masksAtCell = trps $ map (masksAt cells . masksForColor) colors + +masksAt :: [Int] -> [(Row,Mask)]-> [[Mask]] +masksAt [] _ = [] +masksAt (n:ns) !masks = map snd t : masksAt ns f + where + (t, f) = partition test masks + test (r, m) = n' >= 0 && n' < 25 && m `testBit` n' + where n' = n - (n_col * r) + +trps :: [[[Mask]]] -> Array (Row, Col) (Array Color [Mask]) +trps !a = array ((0,0),(9,4)) $ concatMap (uncurry (map . first . (,))) $ + zip [0..9] [copy !! y | y <- [1,0,1,0,1,2,3,4,5,6]] + where + copy = [ [(x,copy' (cellAt x y)) | x <- [0..n_col-1]] | + y <- [1,2,5,6,7,8,9]] + copy' cell = array (0,9) $ map (\clr -> (clr,a !! clr !! cell)) colors + +--- Formatting --- +format :: Bool -> String -> String +format _ [] = "" +format isodd chars | isodd = " " ++ str | otherwise = str + where + (cur, rest) = splitAt 5 chars + str = intersperse ' ' cur ++ " \n" ++ format (not isodd) rest + +toString :: Solution -> String +toString !masks = map color cells + where + masksWithRows = withRows 0 0 (reverse masks) + withRows _ _ [] = [] + withRows board r (m:rest) = (r', m) : withRows board' r' rest + where delta = first0 board `quot` n_col + board' = board `shiftR` (delta * n_col) .|. untag m + r' = r+delta + color n = maybe '.' (("0123456789" !!) . tagof . snd) + (find matches masksWithRows) + where + matches (r, m) + | n' < 0 || n' > 30 = False + | otherwise = (untag m) `testBit` n' + where n' = n - (n_col * r) + +--- Generate the solutions --- +firstZero :: CArray Int Int +firstZero = array (0,31) $ zip [0..31] + [0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5] + +solutions :: [String] +solutions = solveCell 0 colors 0 [] [] + +solveCell :: Row -> [Color] -> Mask -> Solution -> [String] -> [String] +solveCell _ [] board soln results = let s = toString soln + in s:(reverse s):results +solveCell !row !todo !board !soln results + | top/=m_top = foldr solveMask results + [(m, c) | c <- todo, m <- masks ! c, board .&. m == 0] + | otherwise = solveCell (row+1) todo (board `shiftR` n_col) soln results + where top = board .&. m_top + masks = masksAtCell ! (row, (firstZero ! top) ) + solveMask (!m,!c) results = + solveCell row (delete c todo) (untag m .|. board) (m:soln) results + +main = do + n <- return.read.head =<< getArgs + let nsolutions = take n solutions + putStrLn $ (show $ length nsolutions) ++ " solutions found\n" + putStrLn . format False . minimum $ nsolutions + putStrLn . format False . maximum $ nsolutions diff --git a/tests/meteor-contest-u.hs b/tests/meteor-contest-u.hs new file mode 100644 --- /dev/null +++ b/tests/meteor-contest-u.hs @@ -0,0 +1,268 @@ +{-# OPTIONS -O2 -fbang-patterns -optc-O3 #-} + +-- The Computer Language Benchmarks Game +-- http://shootout.alioth.debian.org/ +-- +-- Sterling Clover's translation of Tim Hochberg's Clean implementation + +module Main where +import System.Environment +import Data.Bits +import Data.List +import Data.Array.Unboxed +import Control.Arrow + +--- The Board --- +n_elem = 5 +n_col = 5 +n_row = 10 + +m_top :: Mask +m_top = 0x1F + +cells :: [Cell] +cells = [0..49] + +colors :: [Color] +colors = [0..9] + +cellAt x y = x + n_col * y +coordOf i = snd &&& fst $ i `quotRem` n_col +isValid x y = 0 <= x && x < n_col && 0 <= y && y < n_row + +--- Piece Operations --- +data Direction = E | SE | SW | W | NW | NE deriving (Enum, Eq, Ord) +type Piece = [Direction] +type CellCoord = (Int, Int) +type Mask = Int; type Color = Int; type Row = Int; +type Col = Int; type Tag = Int; type Cell = Int +type Solution = [Mask] + +pieces :: Array Int Piece +pieces = array (0,9) $ zip [0..9] $ + [[E, E, E, SE], + [SE, SW, W, SW], + [W, W, SW, SE], + [E, E, SW, SE], + [NW, W, NW, SE, SW], + [E, E, NE, W], + [NW, NE, NE, W], + [NE, SE, E, NE], + [SE, SE, E, SE], + [E, NW, NW, NW]] + +permutations :: Piece -> [Piece] +permutations p = take 12 (perms p) + where + perms p = p:(flip p) : perms (rotate p) + rotate piece = map r piece + where r E = NE + r NE = NW + r NW = W + r W = SW + r SW = SE + r SE = E + flip piece = map f piece + where f E = W + f NE = NW + f NW = NE + f W = E + f SW = SE + f SE = SW + +--- Mask Operations ---- +untag :: Mask -> Mask +untag mask = mask .&. 0x1ffffff + +retag :: Mask -> Tag -> Mask +retag mask n = untag mask .|. n `shiftL` 25 + +tagof :: Mask -> Tag +tagof mask = mask `shiftR` 25 + +tag :: Mask -> Tag -> Mask +tag mask n = mask .|. n `shiftL` 25 + +count1s :: Mask -> Int +count1s i + | i == 0 = 0 + | i .&. 1 == 1 = 1 + count1s (i `shiftR` 1) + | otherwise = count1s (i `shiftR` 1) + +first0 :: Mask -> Int +first0 i + | i .&. 1 == 0 = 0 + | otherwise = 1 + first0 (i `shiftR` 1) + +--- Making the Bitmasks --- +mod2 x = x .&. 1 +packSize a b = a*5+b +unpackSize n = quotRem n 5 + +move :: Direction -> CellCoord -> CellCoord +move E (x, y) = (x+1, y) +move W (x, y) = (x-1, y) +move NE (x, y) = (x+(mod2 y), y-1) +move NW (x, y) = (x+(mod2 y)-1, y-1) +move SE (x, y) = (x+(mod2 y), y+1) +move SW (x, y) = (x+(mod2 y)-1, y+1) + +pieceBounds :: Piece -> Bool -> (Int, Int, Int, Int) +pieceBounds piece isodd = bnds piece 0 y0 0 y0 0 y0 + where + y0 | isodd = 1 | otherwise = 0 + bnds [] _ _ xmin ymin xmax ymax = (xmin, ymin, xmax, ymax) + bnds (d:rest) x y xmin ymin xmax ymax = + bnds rest x' y' (min x' xmin) (min y' ymin) (max x' xmax) (max y' ymax) + where (x', y') = move d (x, y) + +pieceMask :: Piece -> (Mask, Mask) +pieceMask piece + | odd y1 = (tag (msk piece x2 y2 0) (packSize w2 h2), + tag (msk piece x1 (y1+1) 0 `shiftR` n_col) (packSize w1 h1)) + | otherwise = (tag (msk piece x1 y1 0) (packSize w1 h1), + tag (msk piece x2 (y2+1) 0 `shiftR` n_col) (packSize w2 h2)) + where + (xmin, ymin, xmax, ymax) = pieceBounds piece False + (x1, y1) = (-xmin, -ymin) + w1 = xmax - xmin + h1 = ymax - ymin + (xmin', ymin', xmax', ymax') = pieceBounds piece True + (x2, y2) = (-xmin', (-ymin')+1) + w2 = xmax' - xmin' + h2 = ymax' - ymin' + msk :: Piece -> Col -> Row -> Mask -> Mask + msk [] x y m = m `setBit` cellAt x y + msk (d:rest) x y m = msk rest x' y' (m `setBit` cellAt x y) + where (x', y') = move d (x, y) + +templatesForColor :: Color -> ([Mask], [Mask]) +templatesForColor c = (unzip . map pieceMask) perms + where perms | c == 5 = take 6 ps | otherwise = ps + ps = permutations $ pieces ! c + +--- Looking for Islands --- +noLineIslands :: Mask -> Cell -> Cell -> Int -> Bool +noLineIslands mask start stop step + | (fnd testBit . fnd ((not .) . testBit) . fnd testBit) start > stop = True + | otherwise = False + where + fnd test !x + | x >= 25 = 25 + | test mask x = x + | otherwise = fnd test (x+step) + +noLeftIslands :: Mask -> Bool +noLeftIslands mask = noLineIslands mask 0 20 5 +noRightIslands mask = noLineIslands mask 4 24 5 + +noIslands :: Mask -> Bool +noIslands board = noisles board (count1s board) + +noisles :: Mask -> Int -> Bool +noisles _ 30 = True +noisles board ones + | (ones' - ones) `rem` n_elem /= 0 = False + | otherwise = noisles board' ones' + where board' = fill board (coordOf (first0 board)) + ones' = count1s board' + +fill :: Mask -> CellCoord -> Mask +fill m cc@(x, y) + | x < 0 || x >= n_col = m + | y < 0 || y >= 6 = m + | testBit m i = m + | otherwise = foldl (\m d -> fill m (move d cc)) (setBit m i) + [E, NE, NW, W, SW, SE] + where i = cellAt x y + +--- More Mask Generation --- +masksForColor :: Color -> [(Row, Mask)] +masksForColor c = concatMap atCell cells + where + (evens, odds) = templatesForColor c + atCell n + | even y = [(y, retag (m `shiftL` x) c) | m <- evens , isok m x y] + | odd y = [(y, retag (m `shiftL` x) c) | m <- odds , isok m x y] + where (x, y) = coordOf n + +isok :: Mask -> Row -> Col -> Bool +isok mask x y = + isValid (x+width) (y+height) && + case (y == 0, y+height==9) of + (False, False) -> noLeftIslands mask' && noRightIslands mask' + (False, True) -> noIslands (mask' `shiftL` (n_col * (y - 4))) + (True, _ ) -> noIslands mask' + where (width, height) = unpackSize (tagof mask) + mask' = untag mask `shiftL` x + +masksAtCell :: Array (Row,Col) (Array Color [Mask]) +masksAtCell = trps $ map (masksAt cells . masksForColor) colors + +masksAt :: [Int] -> [(Row,Mask)]-> [[Mask]] +masksAt [] _ = [] +masksAt (n:ns) !masks = map snd t : masksAt ns f + where + (t, f) = partition test masks + test (r, m) = n' >= 0 && n' < 25 && m `testBit` n' + where n' = n - (n_col * r) + +trps :: [[[Mask]]] -> Array (Row, Col) (Array Color [Mask]) +trps !a = array ((0,0),(9,4)) $ concatMap (uncurry (map . first . (,))) $ + zip [0..9] [copy !! y | y <- [1,0,1,0,1,2,3,4,5,6]] + where + copy = [ [(x,copy' (cellAt x y)) | x <- [0..n_col-1]] | + y <- [1,2,5,6,7,8,9]] + copy' cell = array (0,9) $ map (\clr -> (clr,a !! clr !! cell)) colors + +--- Formatting --- +format :: Bool -> String -> String +format _ [] = "" +format isodd chars | isodd = " " ++ str | otherwise = str + where + (cur, rest) = splitAt 5 chars + str = intersperse ' ' cur ++ " \n" ++ format (not isodd) rest + +toString :: Solution -> String +toString !masks = map color cells + where + masksWithRows = withRows 0 0 (reverse masks) + withRows _ _ [] = [] + withRows board r (m:rest) = (r', m) : withRows board' r' rest + where delta = first0 board `quot` n_col + board' = board `shiftR` (delta * n_col) .|. untag m + r' = r+delta + color n = maybe '.' (("0123456789" !!) . tagof . snd) + (find matches masksWithRows) + where + matches (r, m) + | n' < 0 || n' > 30 = False + | otherwise = (untag m) `testBit` n' + where n' = n - (n_col * r) + +--- Generate the solutions --- +firstZero :: UArray Int Int +firstZero = array (0,31) $ zip [0..31] + [0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5] + +solutions :: [String] +solutions = solveCell 0 colors 0 [] [] + +solveCell :: Row -> [Color] -> Mask -> Solution -> [String] -> [String] +solveCell _ [] board soln results = let s = toString soln + in s:(reverse s):results +solveCell !row !todo !board !soln results + | top/=m_top = foldr solveMask results + [(m, c) | c <- todo, m <- masks ! c, board .&. m == 0] + | otherwise = solveCell (row+1) todo (board `shiftR` n_col) soln results + where top = board .&. m_top + masks = masksAtCell ! (row, (firstZero ! top) ) + solveMask (!m,!c) results = + solveCell row (delete c todo) (untag m .|. board) (m:soln) results + +main = do + n <- return.read.head =<< getArgs + let nsolutions = take n solutions + putStrLn $ (show $ length nsolutions) ++ " solutions found\n" + putStrLn . format False . minimum $ nsolutions + putStrLn . format False . maximum $ nsolutions diff --git a/tests/nsieve-bits-c.hs b/tests/nsieve-bits-c.hs new file mode 100644 --- /dev/null +++ b/tests/nsieve-bits-c.hs @@ -0,0 +1,42 @@ +{-# OPTIONS -O2 -optc-O -fbang-patterns #-} +-- +-- The Computer Language Shootout +-- http://shootout.alioth.debian.org/ +-- +-- Contributed by Don Stewart +-- nsieve over an ST monad Bool array +-- + +import Data.Array.IOCArray +import Data.Array.Base +import Data.Array.CArray.Base +import System.IO.Unsafe (unsafePerformIO) +import System +import Control.Monad +import Data.Bits +import Text.Printf + +main = do + n <- getArgs >>= readIO . head :: IO Int + mapM_ (sieve . (10000 *) . (2 ^)) [n, n-1, n-2] + +sieve n = do + let r = unsafePerformIO (do a <- newArray (2,n) True :: IO (IOCArray Int Bool) + go a n 2 0) + printf "Primes up to %8d %8d\n" (n::Int) (r::Int) :: IO () + +go !a !m !n !c + | n == m = return c + | otherwise = do + e <- unsafeRead a n + if e then let loop !j + | j < m = do + x <- unsafeRead a j + when x $ unsafeWrite a j False + loop (j+n) + + | otherwise = go a m (n+1) (c+1) + in loop (n `shiftL` 1) + else go a m (n+1) c + + diff --git a/tests/nsieve-bits-s.hs b/tests/nsieve-bits-s.hs new file mode 100644 --- /dev/null +++ b/tests/nsieve-bits-s.hs @@ -0,0 +1,42 @@ +{-# OPTIONS -O2 -optc-O -fbang-patterns #-} +-- +-- The Computer Language Shootout +-- http://shootout.alioth.debian.org/ +-- +-- Contributed by Don Stewart +-- nsieve over an ST monad Bool array +-- + +import Control.Monad +import Data.Array.Storable +import Data.Array.Base +import System.IO.Unsafe (unsafePerformIO) +import System +import Control.Monad +import Data.Bits +import Text.Printf + +main = do + n <- getArgs >>= readIO . head :: IO Int + mapM_ (sieve . (10000 *) . (2 ^)) [n, n-1, n-2] + +sieve n = do + let r = unsafePerformIO (do a <- newArray (2,n) True :: IO (StorableArray Int Bool) + go a n 2 0) + printf "Primes up to %8d %8d\n" (n::Int) (r::Int) :: IO () + +go !a !m !n !c + | n == m = return c + | otherwise = do + e <- unsafeRead a n + if e then let loop !j + | j < m = do + x <- unsafeRead a j + when x $ unsafeWrite a j False + loop (j+n) + + | otherwise = go a m (n+1) (c+1) + in loop (n `shiftL` 1) + else go a m (n+1) c + + diff --git a/tests/nsieve-bits-u.hs b/tests/nsieve-bits-u.hs new file mode 100644 --- /dev/null +++ b/tests/nsieve-bits-u.hs @@ -0,0 +1,41 @@ +{-# OPTIONS -O2 -optc-O -fbang-patterns #-} +-- +-- The Computer Language Shootout +-- http://shootout.alioth.debian.org/ +-- +-- Contributed by Don Stewart +-- nsieve over an ST monad Bool array +-- + +import Control.Monad.ST +import Data.Array.ST +import Data.Array.Base +import System +import Control.Monad +import Data.Bits +import Text.Printf + +main = do + n <- getArgs >>= readIO . head :: IO Int + mapM_ (sieve . (10000 *) . (2 ^)) [n, n-1, n-2] + +sieve n = do + let r = runST (do a <- newArray (2,n) True :: ST s (STUArray s Int Bool) + go a n 2 0) + printf "Primes up to %8d %8d\n" (n::Int) (r::Int) :: IO () + +go !a !m !n !c + | n == m = return c + | otherwise = do + e <- unsafeRead a n + if e then let loop !j + | j < m = do + x <- unsafeRead a j + when x $ unsafeWrite a j False + loop (j+n) + + | otherwise = go a m (n+1) (c+1) + in loop (n `shiftL` 1) + else go a m (n+1) c + + diff --git a/tests/runtests.sh b/tests/runtests.sh new file mode 100644 --- /dev/null +++ b/tests/runtests.sh @@ -0,0 +1,35 @@ +#!/bin/sh + +compile () { + ghc --make -O2 $1 -o $2 +} + +time_run () { + arg=$1 + shift + for e in $* ; do + time ./$e $arg > $e.out + done + diffn $* +} + +diffn () { + ref=$1 + ret=0 + shift + for f in $* ; do + diff $ref.out $f.out + ret=$(( $ret + $?)) + done + echo '########' Failures: $ret +} + +compile nsieve-bits-u.hs nsU +compile nsieve-bits-c.hs nsC +compile nsieve-bits-s.hs nsS + +compile meteor-contest-u.hs mcU +compile meteor-contest-c.hs mcC + +time_run 8 nsU nsC nsS +time_run 2098 mcU mcC diff --git a/tests/tests.hs b/tests/tests.hs new file mode 100644 --- /dev/null +++ b/tests/tests.hs @@ -0,0 +1,132 @@ +{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} +{-# LANGUAGE FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies, NoMonomorphismRestriction #-} +import Control.Arrow +import Test.QuickCheck +import Text.Show.Functions +import Data.Array.CArray +import Data.Array.CArray.Base (shapeToStride) +import Data.Array.Unboxed +import Data.List +import Foreign.Storable +import Text.Printf +import System.Environment (getArgs) +import System.IO +import System.Random + +instance (Ix i, Arbitrary i, Storable e, Arbitrary e) => Arbitrary (CArray i e) where + arbitrary = do + a <- arbitrary + b <- arbitrary + let l = min a b + u = max a b + es <- vector (rangeSize (l,u)) + return $ listArray (l,u) es + coarbitrary a = coarbitrary (assocs a) + +instance (Ix i, Arbitrary i, Arbitrary e, IArray UArray e) => Arbitrary (UArray i e) where + arbitrary = do + a <- arbitrary + b <- arbitrary + let l = min a b + u = max a b + es <- vector (rangeSize (l,u)) + return $ listArray (l,u) es + coarbitrary a = coarbitrary (assocs a) + +class Model a b where model :: a -> b + +instance (Ix i, IArray a e, Model i i', Model e e') => Model (a i e) ((i',i'),[e']) where + model = (model . bounds &&& map model . elems) +instance (Model i i', Model e e', Ix i', IArray a e') => Model ((i,i),[e]) (a i' e') where + model = uncurry listArray . (model *** map model) +instance (Ix i, Ix i', Model i i', Model e e', Storable e, IArray UArray e') + => Model (CArray i e) (UArray i' e') where + model = uncurry listArray . (model . bounds &&& map model . elems) +instance (Ix i, Ix i', Model i i', Model e e', Storable e', IArray UArray e) + => Model (UArray i e) (CArray i' e') where + model = uncurry listArray . (model . bounds &&& map model . elems) + +-- Types are trivially modeled by themselves +instance Model Bool Bool where model = id +instance Model Int Int where model = id +instance Model Float Float where model = id +instance Model Double Double where model = id +instance (Model a a', Model b b') => Model (a,b) (a',b') where + model (a,b) = (model a, model b) +instance (Model a a', Model b b', Model c c') => Model (a,b,c) (a',b',c') where + model (a,b,c) = (model a, model b, model c) +instance (Model a a', Model b b', Model c c', Model d d') => Model (a,b,c,d) (a',b',c',d') where + model (a,b,c,d) = (model a, model b, model c, model d) + +f =|= g = \a -> + model (f a) == g (model a) +f =||= g = \a b -> + model (f a b) == g a (model b) +infix 1 =|= +infix 1 =||= + +f =|||= g = \a b c -> + model (f a b c) == g a (model b) c +eq4 f g = \a b c d -> + model (f a b c d) == g (model a) (model b) (model c) (model d) +eq5 f g = \a b c d e -> + model (f a b c d e) == g (model a) (model b) (model c) (model d) (model e) + +(===) :: (Eq b) => (a -> b) -> (a -> b) -> a -> Bool +(f === g) x = f x == g x +infixl 1 === + +transposeArray a = ixmap ((swap *** swap) (bounds a)) swap a + where swap = (\(i,j) -> (j,i)) + +prop_flatten_flatten = flatten . flatten === flatten +prop_reshape_flatten a = reshape (0, size a - 1) a == flatten a +prop_rank = length . shape === rank +prop_shape_size = product . shape === size +prop_size = size === rangeSize . bounds +prop_shape_stride_last = last . shapeToStride . shape === const 1 +prop_transpose = transposeArray . transposeArray === id + +ca_tests :: [(String, CArray (Int,Int) Double -> Bool)] +ca_tests = [ ("flatten flatten" , prop_flatten_flatten) + , ("reshape flatten" , prop_reshape_flatten) + , ("rank" , prop_rank) + , ("shape size" , prop_shape_size) + , ("size" , prop_size) + , ("shape stride last" , prop_shape_stride_last) + , ("transpose^2" , prop_transpose) + ] + +prop_amap = (amap :: (Int -> Double) -> CArray Int Int -> CArray Int Double) + =||= (amap :: (Int -> Double) -> UArray Int Int -> UArray Int Double) + +prop_slice_all :: (Int -> Double) -> CArray (Int,Int) Int -> Property +prop_slice_all f a = size a > 0 ==> sliceWith (bounds a) (bounds a) f a == amap f a +prop_ixmapWithInd_amap :: (Int -> Double) -> CArray (Int,Int) Int -> Property +prop_ixmapWithInd_amap f a = size a > 0 ==> ixmapWithInd (bounds a) id (\_ e _ -> f e) a == amap f a + +type Acc = Int +prop_accum f a ies = all (inRange (bounds a) . fst) ies + ==> ( (accum :: (Int -> Acc -> Int) -> CArray Int Int -> [(Int, Acc)] -> CArray Int Int) + =|||= (accum :: (Int -> Acc -> Int) -> UArray Int Int -> [(Int, Acc)] -> UArray Int Int)) f a ies + +prop_composeAssoc f g h = (f . g) . h === f . (g . h) + where types = [f,g,h] :: [CArray Int Int -> CArray Int Int] + +main = do + x <- getArgs + let n = if null x then 100 else read . head $ x + conf = Config { configMaxTest = n + , configMaxFail = 1000 + , configSize = (+ 3) . (`div` 2) + , configEvery = \n args -> let s = show n in s ++ [ '\b' | _ <- s] + } + mycheck (s,a) = printf "%-25s: " s >> check conf a + mapM_ mycheck ca_tests + mapM_ mycheck [ ("amap" , prop_amap) ] + mapM_ mycheck [ ("accum" , prop_accum) ] + mapM_ mycheck [ ("composeAssoc", prop_composeAssoc) ] + mapM_ mycheck [ ("slice all" , prop_slice_all) + , ("ixmapWithInd amap" , prop_ixmapWithInd_amap) ] + +-- arb n k = generate n (mkStdGen k) arbitrary