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.

About these ads

Pages: 1 2

5 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 [...]

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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 )

Google+ photo

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

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 600 other followers

%d bloggers like this: