Hi guys. I just wrote this prime number generator. It's... unusual. It seems to go increadibly fast though. Indeed, I asked it to write the first 1,000,000 prime numbers to a text file, and scanning the bit array seemed to take longer than constructing the bit array! Odd... Anyway, for your amusement: module Primes2 (compute_primes) where import Data.Word import Data.Array.IO seive :: IOUArray Word32 Bool -> Word32 -> Word32 -> IO () seive grid size p = mapM_ (\n -> writeArray grid n False) [p, 2*p .. size] next_prime :: IOUArray Word32 Bool -> IO Word32 next_prime grid = work grid 2 where work grid n = do p <- readArray grid n if p then return n else work grid (n+1) copy :: Word32 -> IOUArray Word32 Bool -> Word32 -> IO (IOUArray Word32 Bool) copy p grid size = do let size' = size * p grid' <- newArray (1,size') False mapM_ (\n -> do b <- readArray grid n if b then mapM_ (\x -> writeArray grid' (n + size*x) True) [0..p-1] else return () ) [1..size] return grid' init_grid :: IO (IOUArray Word32 Bool) init_grid = do grid <- newArray (1,6) False writeArray grid 1 True writeArray grid 5 True return grid compute_primes :: Word32 -> IO [Word32] compute_primes n = do g0 <- init_grid ps <- work n g0 6 return (2:3:ps) where work top grid size | size >= top = do p <- next_prime grid putStrLn $ "Found prime: " ++ show p if p > top then return [p] else do seive grid size p ps <- work top grid size return (p:ps) | otherwise = do p <- next_prime grid putStrLn $ "Seiving prime: " ++ show p let size' = p*size grid' <- copy p grid size seive grid' size' p ps <- work top grid' size' return (p:ps) -- Debug... show_grid :: IOUArray Word32 Bool -> IO () show_grid grid = do (b0,b1) <- getBounds grid mapM_ (\n -> do x <- readArray grid n; if x then putChar '-' else putChar '#') [b0..b1] putStrLn ":"
Andrew Coppin wrote:
copy :: Word32 -> IOUArray Word32 Bool -> Word32 -> IO (IOUArray Word32 Bool) copy p grid size = do let size' = size * p grid' <- newArray (1,size') False
mapM_ (\n -> do b <- readArray grid n if b then mapM_ (\x -> writeArray grid' (n + size*x) True) [0..p-1] else return () ) [1..size]
return grid'
Actually, thinking about this... for most kinds of arrays (whether boxed or unboxed, mutable or immutable) there's probably a more efficient way to copy the data then this. Maybe we should add something to the various array APIs to allow efficient copying of arrays / large chunks of arrays? (In the case of an unboxed array of bits, you can probably copy whole 32-bit or 64-bit words with a few machine instructions, for example.)
Andrew Coppin wrote:
Andrew Coppin wrote:
copy :: Word32 -> IOUArray Word32 Bool -> Word32 -> IO (IOUArray Word32 Bool) copy p grid size = do let size' = size * p grid' <- newArray (1,size') False
mapM_ (\n -> do b <- readArray grid n if b then mapM_ (\x -> writeArray grid' (n + size*x) True) [0..p-1] else return () ) [1..size]
return grid'
Actually, thinking about this... for most kinds of arrays (whether boxed or unboxed, mutable or immutable) there's probably a more efficient way to copy the data then this. Maybe we should add something to the various array APIs to allow efficient copying of arrays / large chunks of arrays?
(In the case of an unboxed array of bits, you can probably copy whole 32-bit or 64-bit words with a few machine instructions, for example.)
For GHC 6.6 I created
foreign import ccall unsafe "memcpy" memcpy :: MutableByteArray# RealWorld -> MutableByteArray# RealWorld -> Int# -> IO ()
{-# INLINE copySTU #-} copySTU :: (Show i,Ix i,MArray (STUArray s) e (ST s)) => STUArray s i e -> STUArray s i e -> ST s () copySTU (STUArray _ _ msource) (STUArray _ _ mdest) = -- do b1 <- getBounds s1 -- b2 <- getBounds s2 -- when (b1/=b2) (error ("\n\nWTF copySTU: "++show (b1,b2))) ST $ \s1# -> case sizeofMutableByteArray# msource of { n# -> case unsafeCoerce# memcpy mdest msource n# s1# of { (# s2#, () #) -> (# s2#, () #) }}
To allow efficient copying of STUArrays.
ChrisK wrote:
For GHC 6.6 I created
foreign import ccall unsafe "memcpy" memcpy :: MutableByteArray# RealWorld -> MutableByteArray# RealWorld -> Int# -> IO ()
{-# INLINE copySTU #-} copySTU :: (Show i,Ix i,MArray (STUArray s) e (ST s)) => STUArray s i e -> STUArray s i e -> ST s () copySTU (STUArray _ _ msource) (STUArray _ _ mdest) = -- do b1 <- getBounds s1 -- b2 <- getBounds s2 -- when (b1/=b2) (error ("\n\nWTF copySTU: "++show (b1,b2))) ST $ \s1# -> case sizeofMutableByteArray# msource of { n# -> case unsafeCoerce# memcpy mdest msource n# s1# of { (# s2#, () #) -> (# s2#, () #) }}
To allow efficient copying of STUArrays.
So... that copies the entire array into another array of the same size? (I'm having a lot of trouble understanding the code...)
Andrew Coppin wrote:
ChrisK wrote:
For GHC 6.6 I created
foreign import ccall unsafe "memcpy" memcpy :: MutableByteArray# RealWorld -> MutableByteArray# RealWorld -> Int# -> IO ()
{-# INLINE copySTU #-} copySTU :: (Show i,Ix i,MArray (STUArray s) e (ST s)) => STUArray s i e -> STUArray s i e -> ST s () copySTU (STUArray _ _ msource) (STUArray _ _ mdest) = -- do b1 <- getBounds s1 -- b2 <- getBounds s2 -- when (b1/=b2) (error ("\n\nWTF copySTU: "++show (b1,b2))) ST $ \s1# -> case sizeofMutableByteArray# msource of { n# -> case unsafeCoerce# memcpy mdest msource n# s1# of { (# s2#, () #) -> (# s2#, () #) }}
To allow efficient copying of STUArrays.
So... that copies the entire array into another array of the same size? (I'm having a lot of trouble understanding the code...)
Yes, that is what it does. The "STUArray" data type has the STUArray constructor which I import and pattern match against. The imports are:
import Data.Array.Base(unsafeRead,unsafeWrite,STUArray(..)) import GHC.Prim(MutableByteArray#,RealWorld,Int#,sizeofMutableByteArray#,unsafeCoerce#)
in 6.6.1 this is defined as
data STUArray s i a = STUArray !i !i (MutableByteArray# s) in 6.8.1 this is defined as data STUArray s i a = STUArray !i !i !Int (MutableByteArray# s)
I use sizeofMutableByteArray# to get the source size, n#. I have lost track of how unsafeCoerce# and s1# are being used...oops. It is similar to data-dependency tricks used inside Data.Array.Base, though. -- Chris
ChrisK wrote:
For GHC 6.6 I created
foreign import ccall unsafe "memcpy" memcpy :: MutableByteArray# RealWorld -> MutableByteArray# RealWorld -> Int# -> IO ()
{-# INLINE copySTU #-} copySTU :: (Show i,Ix i,MArray (STUArray s) e (ST s)) => STUArray s i e -> STUArray s i e -> ST s () copySTU (STUArray _ _ msource) (STUArray _ _ mdest) = -- do b1 <- getBounds s1 -- b2 <- getBounds s2 -- when (b1/=b2) (error ("\n\nWTF copySTU: "++show (b1,b2))) ST $ \s1# -> case sizeofMutableByteArray# msource of { n# -> case unsafeCoerce# memcpy mdest msource n# s1# of { (# s2#, () #) -> (# s2#, () #) }}
To allow efficient copying of STUArrays.
How does this guarantee that it doesn't overflow the buffer of the destination array? Reinier
Reinier Lamers wrote:
ChrisK wrote:
For GHC 6.6 I created
foreign import ccall unsafe "memcpy" memcpy :: MutableByteArray# RealWorld -> MutableByteArray# RealWorld -> Int# -> IO ()
{-# INLINE copySTU #-} copySTU :: (Show i,Ix i,MArray (STUArray s) e (ST s)) => STUArray s i e -> STUArray s i e -> ST s () copySTU (STUArray _ _ msource) (STUArray _ _ mdest) = -- do b1 <- getBounds s1 -- b2 <- getBounds s2 -- when (b1/=b2) (error ("\n\nWTF copySTU: "++show (b1,b2))) ST $ \s1# -> case sizeofMutableByteArray# msource of { n# -> case unsafeCoerce# memcpy mdest msource n# s1# of { (# s2#, () #) -> (# s2#, () #) }}
To allow efficient copying of STUArrays.
How does this guarantee that it doesn't overflow the buffer of the destination array?
As is, the above is a very unsafe operation. The check is commented out. You can uncomment the code above that says:
do b1 <- getBounds s1 b2 <- getBounds s2 when (b1/=b2) (error ("\n\nWTF copySTU: "++show (b1,b2)))
Which checks the high-level boundary matches, not just the actual length. I only have a single size of STUArray haning around, so I use the unsafe and fast version. -- Chris
Andrew Coppin wrote:
Andrew Coppin wrote:
copy :: Word32 -> IOUArray Word32 Bool -> Word32 -> IO (IOUArray Word32 Bool) copy p grid size = do let size' = size * p grid' <- newArray (1,size') False
mapM_ (\n -> do b <- readArray grid n if b then mapM_ (\x -> writeArray grid' (n + size*x) True) [0..p-1] else return () ) [1..size]
return grid'
Actually, thinking about this... for most kinds of arrays (whether boxed or unboxed, mutable or immutable) there's probably a more efficient way to copy the data then this. Maybe we should add something to the various array APIs to allow efficient copying of arrays / large chunks of arrays?
Ideally we'd have the compiler generate optimal code anyway, then you wouldn't need such things. So, let's hope those compiler/optimiser guys keep working hard! Jules
Jules Bean wrote:
Andrew Coppin wrote:
Andrew Coppin wrote:
copy :: Word32 -> IOUArray Word32 Bool -> Word32 -> IO (IOUArray Word32 Bool) copy p grid size = do let size' = size * p grid' <- newArray (1,size') False
mapM_ (\n -> do b <- readArray grid n if b then mapM_ (\x -> writeArray grid' (n + size*x) True) [0..p-1] else return () ) [1..size]
return grid'
Actually, thinking about this... for most kinds of arrays (whether boxed or unboxed, mutable or immutable) there's probably a more efficient way to copy the data then this. Maybe we should add something to the various array APIs to allow efficient copying of arrays / large chunks of arrays?
Ideally we'd have the compiler generate optimal code anyway, then you wouldn't need such things.
So, let's hope those compiler/optimiser guys keep working hard!
True - but "copyArray x y" is surely a lot clearer than some longwinded mapM thing, no? ;-)
participants (4)
-
Andrew Coppin -
ChrisK -
Jules Bean -
Reinier Lamers