Return to Answer

###Update

Update

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
intSqrt :: Int -> Int
intSqrt n = ceiling (sqrt $ fromIntegral n)
-- I believe ceil makes more sense here than floor+1.
-- We make a vector indicating which *odd* numbers are prime. An odd number n
-- is prime if (primalityTable n)!n is True.
primalityTable :: Int -> Vec.Vector Bool
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
 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 <- 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 -> 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
-- Note that MNPS.primeSieve gives too many results for very small
-- arguments, so we need the takeWhile
primesSumMNPS upto = sum . takeWhile (< upto) . map fromInteger $ MNPS.primeList (MNPS.primeSieve (toInteger upto))
staticTests = [(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 for "++show n) $ summer n `shouldBe` s
main :: IO ()
main = hspec $ do
 -- Since we're using primesMNPS to check primesSum, we want to be sure
 -- it passes the static tests.
 describe "Primes.primesMNPS" $ do
 testOnList primesSumMNPS staticTests
 describe "Primes.primesSum" $ do
 testOnList primesSum staticTests
 it "works just like primesSumMNPS" $
 property $ \n -> let upto = abs (n `rem` 4000000)
 in primesSum upto == primesSumMNPS upto

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

###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
intSqrt :: Int -> Int
intSqrt n = ceiling (sqrt $ fromIntegral n)
-- I believe ceil makes more sense here than floor+1.
-- We make a vector indicating which *odd* numbers are prime. An odd number n
-- is prime if (primalityTable n)!n is True.
primalityTable :: Int -> Vec.Vector Bool
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
 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 <- 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 -> 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
-- Note that MNPS.primeSieve gives too many results for very small
-- arguments, so we need the takeWhile
primesSumMNPS upto = sum . takeWhile (< upto) . map fromInteger $
 MNPS.primeList (MNPS.primeSieve (toInteger upto))
staticTests = [(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 for "++show n) $ summer n `shouldBe` s
main :: IO ()
main = hspec $ do
 -- Since we're using primesMNPS to check primesSum, we want to be sure
 -- it passes the static tests.
 describe "Primes.primesMNPS" $ do
 testOnList primesSumMNPS staticTests
 describe "Primes.primesSum" $ do
 testOnList primesSum staticTests
 it "works just like primesSumMNPS" $
 property $ \n -> let upto = abs (n `rem` 4000000)
 in primesSum upto == primesSumMNPS upto