Haskell monads: sum of primes

Question

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 Control.Monad.ST
--import Control.Monad
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

• Liviu: 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 ?

Answer 1

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 Control.Monad.ST
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

• Liviu: Thank you, I replaced the unsafe functions as mentioned.
• Liviu: 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?
• dfeuer: @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.
• Liviu: sorry for being rude, I understood what you were saying
• dfeuer: 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.

Answer 2

@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

• Liviu: This is NOT really an answer, but the latest version: PLEASE do not vote it!
• dfeuer: It's safe now! However, a brief look suggests it's probably not nearly as fast as it could be, and a lot of the code seems a bit strange. I'll update my answer to suggest further improvements.
• Liviu: @dfeuer added an optimization (for every prime we have to start from prime*prime, not from 3*prime)
• Liviu: Another problem is that ghci seems to cache the results: "*Main> main 142913828922 (5.40 secs, 1386982216 bytes) *Main> main 142913828922 (0.00 secs, 0 bytes)"
• dfeuer: I've never heard of that music site. My answer shows how, when, and why to use Control.Exception.assert. I know that QuickCheck, which started out as Haskell-only, has since been ported to a variety of other languages. I've never used HUnit, but my general understanding is that it's similar to other unit testing frameworks.