They say you cannot get out of the monad, but I need to use the result (at least for assertions / unit tests). But print how does it do it?

Otherwise, I worked very hard for the next code with my new limited Haskell abilities (especially when talking about monads), so it needs a review.

import qualified Data.Vector as Vec
import qualified Data.Vector.Mutable as MutableVec

assert :: Monad a => Bool -> a ()
assert False = error "assertion failed!"
assert _     = return ()

get_primes_limit :: Int -> Int
get_primes_limit limit = 1 + floor (sqrt $fromIntegral limit) update_sieve4_step_M 2 limit arr = do mapM_ (\i -> MutableVec.write arr i 0) [4, 6 .. limit-1] update_sieve4_step_M prime limit arr = do mapM_ (\i -> MutableVec.write arr i 0) [3*prime, 5*prime .. limit-1] update_sieve4_step limit arrIn = do arr <- Vec.unsafeThaw arrIn MutableVec.write arr 1 0 update_sieve4_step_M 2 limit arr mapM_ (\p -> update_sieve4_step_M p limit arr) [3, 5 .. (get_primes_limit limit)] Vec.unsafeFreeze arr update_sieve4 limit = update_sieve4_step limit (Vec.generate limit (id)) --TODO: cannot get this out of the monad (not m Int); print does it somehow primes_sum4_1 limit = do sieve <- update_sieve4 limit return (Vec.sum sieve) --primes_sum4 :: Control.Monad.Primitive.PrimMonad m => Int -> m Int primes_sum4 limit = do result <- primes_sum4_1 limit --print result return result primes_sum4_validation = do result <- primes_sum4_1 9 --assert (17 == result) print "Nothing to validate!" validate = do --assert (17 == primes_sum4 9) --assert (17 == primes_sum4 10) --assert (17 == primes_sum4 11) --assert (28 == primes_sum4 12) --assert (100 == primes_sum4 29) --assert (129 == primes_sum4 30) --assert (1060 == primes_sum4 100) --assert (76127 == primes_sum4 1000) --assert (5736396 == primes_sum4 10000) --assert (32405717 == primes_sum4 25000) --assert (454396537 == primes_sum4 100000) print "Nothing to validate!" main = do primes_sum4 2000000 • Well, I discovered runST: primes_sum4 limit = runST$ do primes_sum4_1 limit Immediately after posting, when closing the browser tabs :)))) What do I do now ? Jan 4 '15 at 17:05

In Haskell, unsafe can mean various things, but it almost always means that you have to have a very solid understanding of what the unsafe thing is about before you use it if you want to have any hope of staying out of trouble.

unsafeFreeze and unsafeThaw in particular are low-level functions exposed to allow programmers to dig into the internal guts of vectors to build new abstractions offering a safe interface. For example, you might use a pure function to build a new vector, then use unsafeThaw to produce a mutable vector. You'd then have a new function hiding away the unsafety and exposing a plain old "make a new vector like this within ST/IO/PrimMonad". Or you might create a mutable vector, do a bunch of mutation, and then unsafeFreeze it to make a pure one. You'd then have a new function hiding away the unsafety and exposing a plain old "make an immutable vector like this".

What you've written is something you should never see: an action that takes what, to the outside world, is a pure value, and actually changes it. Such a "function" breaks all the usual rules of Haskell; it's as strange as something that changes every 3 in the world into a 4.

What you should do instead, if you want to keep mutating the same vector, is put that whole thing within the monad. Make the vector, do all your mutations (holding on to whatever you need to return) and then produce a result at the end.

Side note: the general rule about unsafe also applies to any function or constructor imported from a module named GHC.*, *.Internal or *.Private, as well as certain functions whose names begin with unchecked or end in Descending or Ascending. The primitive member of the PrimMonad class is also deeply unsafe.

### Update

You should read The Genuine Sieve of Eratosthenes by Melissa O'Neill, which shows how to lazily produce a list of primes in an efficient manner. It's a very accessible and interesting article.

I've pasted my own version of a fairly simple vector-based sieve below, with extensive comments. To run the testing code you will need the arithmoi and hspec packages. You can install those (and vector) in a Cabal sandbox if you wish.

module Primes where

import qualified Data.Vector.Unboxed as Vec
import qualified Data.Vector.Unboxed.Mutable as MutableVec
import Control.Exception (assert)

-- For testing
import Test.Hspec
import Test.QuickCheck
import qualified Math.NumberTheory.Primes.Sieve as MNPS

-- Approximate square root of an Int. Note that this may not work properly
-- when the Int is too large to be represented precisely as a Double.
intSqrt :: Int -> Int
intSqrt n = ceiling (sqrt $fromIntegral n) -- We make a vector indicating which *odd* numbers are prime. An odd number n -- is prime if (primalityTable n) ! getIndex n is True. primalityTable :: Int -> Vec.Vector Bool primalityTable 0 = Vec.empty primalityTable 1 = Vec.singleton False primalityTable 2 = Vec.fromList [False, False] primalityTable upto = runST (do let limit = let root = intSqrt upto in root - fromEnum (even root) -- Make a mutable vector and fill it with True arr <- MutableVec.replicate (getIndex (upto - fromEnum (even upto)) + 1) True -- 1 is not prime assert (MutableVec.length arr > 1)$ MutableVec.write arr (getIndex 1) False

-- For each element of the vector up to the limit, check if the element has
-- been crossed off. If not, use it to cross off other elements.
flip mapM_ [0 .. getIndex limit] $\i -> do -- Does the index i represent a prime? isPrime <- assert (i < MutableVec.length arr)$ MutableVec.read arr i
if isPrime
then crossOff i arr
else return ()
Vec.freeze arr)

-- Since we're skipping all the even numbers, we use these functions to convert
-- numbers and the vector indices they correspond to. That way, we don't have
-- to try to keep track of the conversions in our poor heads, except where we
-- want to.  Since we're using Control.Exception.assert, we'll get useful
-- information about where in the source the assertion failed, and the
-- assertion will go away when we compile with optimization enabled (unless we
-- use -fno-ignore-asserts).
getIndex :: Int -> Int
getIndex n = assert (odd n) $n quot 2 getValue :: Int -> Int getValue n = n + n + 1 -- Cross off odd multiples of the prime represented by the index passed in. How -- do we calculate this efficiently? We are working with the prime p = getValue -- i = 2 * i + 1. We first cross off the value p*p. Instead of monkeying around -- with getValue and getIndex in the inner loop, we will use thte fact that -- moving up in value by 2*p corresponds to moving up in index by p. crossOff :: Int -> MutableVec.MVector s Bool -> ST s () crossOff i arr = mapM_ (\q -> assert (q < MutableVec.length arr) MutableVec.write arr q False) [startingIndex, startingIndex+p .. MutableVec.length arr - 1] where p = getValue i startingIndex = getIndex (p*p) -- Note that we can reasonably expect the intermediate vector here to be fused away by -- the fancy compiler rewrite rules in the vector library. primesSum :: Int -> Int primesSum uptoNotEqual | uptoNotEqual <= 2 = 0 | otherwise = 2 + Vec.sum (Vec.imap (\i isPrime -> if isPrime then getValue i else 0) (primalityTable (uptoNotEqual-1))) primesSumMNPS :: Int -> Int primesSumMNPS upto = sum . takeWhile (< upto) . map fromInteger$ MNPS.primes

staticTests = [(0,0),(1,0),(2,0),(3,2),(4,5),(5,5),(6,10),(9,17),(10,17),(11,17),
(12,28),(29,100),(30,129),
(100,1060),(1000,76127),(10000,5736396),
(25000,32405717), (100000,454396537),
(2000000,142913828922)]

testOnList summer lst =
flip mapM_ lst $\(n,s) -> it ("works for "++show n)$ summer n shouldBe s

main :: IO ()
main = do
hspec $do -- Since we're using primesMNPS to check primesSum, -- we want to be sure it passes the static tests. describe "Primes.primesMNPS"$
testOnList primesSumMNPS staticTests
describe "Primes.primesSum" $do testOnList primesSum staticTests it "works just like primesSumMNPS"$
property $\n -> let upto = abs (n rem 30000000) in primesSum upto == primesSumMNPS upto • Thank you, I replaced the unsafe functions as mentioned. Jan 17 '15 at 17:00 • Is this English "an action that takes what, to the outside world, is a pure value, and actually changes it." ? I don't get the 3 and 4 tale - is this some kind of rhetoric? Jan 17 '15 at 17:02 • @Liviu, it is English, and the tale is some kind of rhetoric, yes. In Java, say, a vector is a pointer to some boxes in memory, and an immutable vector is a pointer to some boxes in memory that you're not allowed to change. In Haskell, mutable and immutable vectors are radically different sorts of things from each other, although under the hood they are implemented similarly. An immutable vector in Haskell is not (conceptually) a pointer to boxes at all; it is a compound value much like a tuple. Jan 17 '15 at 19:18 • sorry for being rude, I understood what you were saying Jan 17 '15 at 19:30 • By the way, the reason the code above has so many assertions is that the original version was buggy. I fixed the bugs, but kept getting segmentation faults. Eventually I narrowed the problem down to a couple functions in arithmoi that I was using; I switched to some others and filed a bug report. Jan 19 '15 at 7:45 @dfeuer's review: replaced unsafe functions with their "safe" counterparts import Control.Monad.ST import qualified Data.Vector as Vec import qualified Data.Vector.Mutable as MutableVec assert :: Monad a => Bool -> a () assert False = error "assertion failed!" assert _ = return () get_primes_limit :: Int -> Int get_primes_limit limit = 1 + floor (sqrt$ fromIntegral limit)

update_sieve4_step_M 2 limit arr = do
mapM_ (\i -> MutableVec.write arr i 0) [4, 6 .. limit-1]
update_sieve4_step_M prime limit arr = do
mapM_ (\i -> MutableVec.write arr i 0) [prime*prime, (prime+2)*prime .. limit-1]

update_sieve4_step limit arrIn = do
arr <- Vec.thaw arrIn
MutableVec.write arr 1 0
update_sieve4_step_M 2 limit arr
mapM_ (\p -> update_sieve4_step_M p limit arr) [3, 5 .. (get_primes_limit limit)]
Vec.freeze arr

update_sieve4 limit = update_sieve4_step limit (Vec.generate limit (id))

primes_sum4_1 limit = do
sieve <- update_sieve4 limit
return (Vec.sum sieve)

primes_sum4 limit = runST $do primes_sum4_1 limit validate = do assert (17 == primes_sum4 9) assert (17 == primes_sum4 10) assert (17 == primes_sum4 11) assert (28 == primes_sum4 12) assert (100 == primes_sum4 29) assert (129 == primes_sum4 30) assert (1060 == primes_sum4 100) assert (76127 == primes_sum4 1000) assert (5736396 == primes_sum4 10000) assert (32405717 == primes_sum4 25000) assert (454396537 == primes_sum4 100000) assert (False || (142913828922 == primes_sum4 2000000)) print "Validation done!" main = do print$ primes_sum4 2000000