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.
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
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
*.Private, as well as certain functions whose names begin with
unchecked or end in
primitive member of the
PrimMonad class is also deeply unsafe.
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
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 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
then crossOff i arr
else return ()
-- 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
| uptoNotEqual <= 2 = 0
| otherwise = 2 + Vec.sum (Vec.imap (\i isPrime -> if isPrime then getValue i else 0)
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),
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