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))
78498
> (length (primes #e2e6))
148933
> (- 148933 78498)
70435
> (length (primes-range #e1e6 #e2e6 #e1e5))
70435

You can run the example at http://programmingpraxis.codepad.org/ZQA6ldK8.

Pages: 1 2

9 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)
      where
        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
      where
        middle (_, v, _) = v
        stop (sb, _, _) = sb == r
        start = (l, [], map (firstQ l) primes)
        incr (startblock, primessofar, q) = (nextblock, nextprimes, nextq)
          where
    	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 http://www.spoj.pl/problems/PRIME1/ 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: […]

  7. […] can see a simple implementation of a segmented sieve here. Note that a segmented sieve will be very much faster than O’Neill’s priority-queue […]

  8. […] ver una implementación simple de un tamiz segmentado aquí. Tenga en cuenta que un tamiz segmentado será mucho más rápido que el tamiz de cola de prioridad […]

  9. […] can see a simple implementation of a segmented sieve here. Note that a segmented sieve will be very much faster than O’Neill’s priority-queue […]

Leave a comment