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.
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][…] 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 […]
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).
I have worked on this problm http://www.spoj.pl/problems/PRIME1/ about this,but this semms too hard for me!!
[…] 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 […]
[…] found a solution here which says the function should be $(-1/2 × (L + 1 + Pk))mod Pk$ , where L:start of range, Pk: […]
[…] 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 […]
[…] 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 […]
[…] 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 […]