First of all, it's great. However, I came across a situation where my benchmarks turned up weird results. I am new to Haskell, and this is first time I've gotten my hands dirty with mutable arrays and Monads. The code below is based on this example.

I wrote a generic monadic for function that takes numbers and a step function rather than a range (like forM_ does). I compared using my generic for function (Loop A) against embedding an equivalent recursive function (Loop B). Having Loop A is noticeably faster than having Loop B. Weirder, having both Loop A and B together is faster than having Loop B by itself (but slightly slower than Loop A by itself).

Some possible explanations I can think of for the discrepancies. Note that these are just guesses:

• Something I haven't learned yet about how Haskell extracts results from monadic functions.
• Loop B faults the array in a less cache efficient manner than Loop A. Why?
• I made a dumb mistake; Loop A and Loop B are actually different.
  • Note that in all 3 cases of having either or both Loop A and Loop B, the program produces the same output.

Here is the code. I tested it with ghc -O2 for.hs using GHC version 6.10.4 .

import Control.Monad
import Control.Monad.ST
import Data.Array.IArray
import Data.Array.MArray
import Data.Array.ST
import Data.Array.Unboxed

for :: (Num a, Ord a, Monad m) => a -> a -> (a -> a) -> (a -> m b) -> m ()
for start end step f = loop start where
    loop i
        | i <= end   = do
            f i
            loop (step i)
        | otherwise  = return ()

primesToNA :: Int -> UArray Int Bool
primesToNA n = runSTUArray $ do
    a <- newArray (2,n) True :: ST s (STUArray s Int Bool)
    let sr = floor . (sqrt::Double->Double) . fromIntegral $ n+1

    -- Loop A
    for 4 n (+ 2) $ \j -> writeArray a j False

    -- Loop B
    let f i
        | i <= n     = do
            writeArray a i False
            f (i+2)
        | otherwise  = return ()
        in f 4

    forM_ [3,5..sr] $ \i -> do
        si <- readArray a i
        when si $
            forM_ [i*i,i*i+i+i..n] $ \j -> writeArray a j False
    return a

primesTo :: Int -> [Int]
primesTo n = [i | (i,p) <- assocs . primesToNA $ n, p]

main = print $ primesTo 30000000