## 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

### 7 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