Segmented Sieve Of Eratosthenes

February 5, 2010

We use primes from the previous Sieve of Eratosthenes exercise and isqrt from the Standard Prelude:

(define (primes-range l r b)
  (let* ((ps (cdr (primes (+ (isqrt r) 1))))
         (qs (map (lambda (p) (modulo (* -1/2 (+ l 1 p)) p)) ps))
         (zs '()) (z (lambda (p) (set! zs (cons p zs)))))
    (do ((t l (+ t b b))
         (qs qs (map (lambda (p q) (modulo (- q b) p)) ps qs)))
        ((= t r) (reverse zs))
      (let ((bs (make-vector b #t)))
        (do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs))
          (do ((j (car qs) (+ j (car ps)))) ((<= b j))
            (vector-set! bs j #f)))
        (do ((j 0 (+ j 1))) ((= j b))
          (if (vector-ref bs j) (z (+ t j j 1))))))))

Ps and qs are as in the problem description, zs accumulates the resulting primes, and t is the beginning of each sieving block, initially l. The outer loop, on t, establishes the sieving blocks and updates qs from one block to the next. The inner loop on qs and its nested loop on j perform the sieving. The final loop on j collects the primes from the current block. Here is a sample use:

> (length (primes #e1e6))
> (length (primes #e2e6))
> (- 148933 78498)
> (length (primes-range #e1e6 #e2e6 #e1e5))

You can run the example at

Pages: 1 2

6 Responses to “Segmented Sieve Of Eratosthenes”

  1. Mithrandir said

    My first try

    firstQ :: Int -> Int -> Int
    firstQ l pk = (-(1+l+pk) `div` 2) `mod` pk
    nextQ :: Int -> Int -> Int -> Int
    nextQ b pk qk = (qk - b) `mod` pk
    sieveBlock :: [Int] -> [Int] -> [Int] -> [Int]
    sieveBlock [] [] blk = blk
    sieveBlock (p:ps) (o:ofs) blk = filter f (sieveBlock ps ofs blk)
        f x = (x < o) || ((x - o) `mod` p /= 0)
    getValsfromBlock :: Int -> [Int] -> [Int]
    getValsfromBlock l bl = map (\x->l+1+2*x) bl
    sieve :: Int -> Int -> Int -> [Int] -> [Int]
    sieve l r b primes = middle $ until stop incr start
        middle (_, v, _) = v
        stop (sb, _, _) = sb == r
        start = (l, [], map (firstQ l) primes)
        incr (startblock, primessofar, q) = (nextblock, nextprimes, nextq)
    	array = [0 .. b-1]
    	nextblock = startblock + 2 * b
    	nextprimes = primessofar ++ (getValsfromBlock startblock bl)
    	nextq = zipWith (nextQ b) primes q
    	bl = sieveBlock primes q array
    main = print $ sieve 100 200 10 [3, 5, 7, 11, 13]
  2. […] Of Eratosthenes February 5, 2010 Mithrandir Leave a comment Go to comments Today’s Programming Praxis problem describes an algorithm which may be used to find the prime numbers in a large interval, so large […]

  3. BMeph said

    Won’t doing “(primes (+ (isqrt r) 1))” sometimes add an extra prime that you won’t need?

    For instance, try it with l = 200, and r = 325 (basically, any time r is at least the value of
    an even square, but less than an odd square).

  4. garretthh07 said

    I have worked on this problm about this,but this semms too hard for me!!

  5. […] I would probably choose the first alternative given your problem size. Implementations of both the segmented Sieve of Eratosthenes and Lehmer’s prime-counting function are available at my […]

  6. […] found a solution here which says the function should be $(-1/2 × (L + 1 + Pk))mod Pk$ , where L:start of range, Pk: […]

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: