[Haskell-cafe] A tale of Project Euler

[Sebastian Sylvan]

On Nov 29, 2007 6:43 PM, Andrew Coppin <andrewcoppin at btinternet.com> wrote:
> Daniel Fischer wrote:
> > One thing: since You check the array bounds, the system needn't check them
> > again, use unsafeWrite and unsafeRead. And use Int for the index, that would
> > be MUCH faster.
> >
>
> I can't find the functions you're talking about. I have however changed
> the index type. (Make little or no noticable speed difference. But then,
> it's already pretty damn fast in the first place...)
>
> > Another thing: you can stop sieving when p*p > size, another speedup
> >
>
> Saves a few hundred milliseconds.
>
> > Fifth thing: better use an STUArray, don't drag IO in if it's not necessary.
> >
>
> I don't understand the ST monad.
>


There's not a whole lot to understand if you just want to use it
(though it's all very cool from a theoretical standpoint too). Here
are my minor changes to your program.

import Data.Array.ST
import Control.Monad.ST


calc_primes :: [Word64]
calc_primes = runST ( do
   grid <- newArray (2,size) True
   seive 2 grid )
 where
   seive :: Word64 -> STUArray s Word64 Bool -> ST s [Word64]
   seive p g = do
     mapM_ (\n -> writeArray g n False) [p, 2*p .. size]
     mp' <- next (p+1) g
     case mp' of
       Nothing -> return [p]
       Just p' -> do
         ps <- seive p' g
         return (p:ps)

   next :: Word64 -> STUArray s Word64 Bool -> ST s (Maybe Word64)
   next p g = do
     if p == size
       then return Nothing
       else do
         t <- readArray g p
         if t
           then return (Just p)
           else next (p+1) g


The benefit should be obvious: No pesky IO type, so you can use it in
your pure code. You just need to give a type signature somewhere to
show the type system that you're using STUArray, but the rest just
uses the same type class as you already used for the IOUArrays.

And here are the modifications to use the unsafe reads and writes (42%
speedup for me):

import Data.Array.Base

size = 1000001 :: Word64

calc_primes :: [Word64]
calc_primes = runST ( do
   grid <- newArray (2,size) True
   seive 2 grid )
 where
   seive :: Word64 -> STUArray s Word64 Bool -> ST s [Word64]
   seive p g = do
     mapM_ (\n -> unsafeWrite g (fromIntegral n) False) [p, 2*p .. size]
     mp' <- next (p+1) g
     case mp' of
       Nothing -> return [p]
       Just p' -> do
         ps <- seive p' g
         return (p:ps)

   next :: Word64 -> STUArray s Word64 Bool -> ST s (Maybe Word64)
   next p g = do
     if p == size
       then return Nothing
       else do
         t <- unsafeRead g (fromIntegral p)
         if t
           then return (Just p)
           else next (p+1) g

-- 
Sebastian Sylvan
+44(0)7857-300802
UIN: 44640862