## Sieve Of Sundaram

### October 7, 2011

We use n for the maximum prime and m for the half-size list of integers:

```(define (sundaram n)   (let* ((m (quotient n 2)) (pv (make-vector (+ m 1) #t)))     (do ((i 1 (+ i 1))) ((< (quotient m 4) i))       (do ((j i (+ j 1))) ((< (quotient (- m i) (+ i i 1)) j))         (vector-set! pv (+ i j (* 2 i j)) #f)))     (let loop ((i 1) (ps (list 2)))       (cond ((= i m) (reverse ps))             ((vector-ref pv i) (loop (+ i 1) (cons (+ i i 1) ps)))             (else (loop (+ i 1) ps))))))```

We employ two optimizations. Since the condition i + j + 2ijn implies j < (ni)/(2i+1), the j loop stops early. Similarly, i can never be more than one-quarter of n, allowing the i loop to stop early.

You can run the program at http://programmingpraxis.codepad.org/O6J0rkDn. Shown there are several other sieves, for comparison. The sieves of Euler and Sundaram are slowest, the simple sieve of Eratosthenes is fastest, and the sieve of Atkin is between. We also show a segmented version of the sieve of Eratosthenes, which is a little bit slower than the simple sieve of Eratosthenes (due to the time spent resetting the sieve from one segment to the next) but has the advantages of using much less space and providing a lower bound other than zero.

Pages: 1 2

### 7 Responses to “Sieve Of Sundaram”

1. khigia said

In Go:

package main

import (
“fmt”
“math”
)

func main() {
n := 2000

// initilize the sieve as array of boolean set to true
sieve := make([]bool, n)
for i, _ := range sieve {
sieve[i] = true
}

// since j >= i then i+j+2ij >= i+i+2ii
// since i+j+2ij < n then i+i+2ii < n
// ii+2i+1 < n+1
// (i+1)(i+1) i i < sqrt(n+1)
boundi := int(math.Sqrt(float64(n+1)))
for i := 1 ; i <= boundi ; i++ {
// i+j+2ij < n j(1+2i) j < (n-i)/(1+2i)
boundj := (n-i)/(1+2*i)
for j := i ; j <= boundj; j++ {
remidx := i + j + 2 * i * j
if remidx < n {
sieve[remidx] = false
}
}
}
for i, v := range sieve {
if v {
fmt.Println(2 * i + 1)
}
}
}

2. Mike said

Python 3

For any value of i, minimum k to be sieved is 2*i*i + 2*i, and the
step between k’s to be sieved is 2*i + 1. Stop when minimum k too big.

```from itertools import count

def sundaram(n):
nk = (n-1)//2
ks = list(range(nk+1))

for i in count(1):
step  = 2*i+1
start = i*(step + 1)
if start > nk:
break

ks[start::step] = (0 for _ in range(start, nk+1, step))

return [2] + [2*k+1 for k in filter(None, ks)]
```
3. murkey said
```#!/usr/bin/ruby

# Sieve of Sundaram

def primes(m)
n = (m - 1) / 2

ints = Array(1..n)

for j in 1..(n/3).to_i
i = 1
while i + j + 2*i*j <= n
notprime = i + j + 2*i*j
ints[notprime - 1] = nil
i += 1
end
end

ints.compact!
ints.each do |k|
prime = 2*k + 1
print "#{prime} "
end
end

primes(1000000)

```
4. CyberSpace17 said

My C++ Solution
http://ideone.com/m6MAW

5. deang said

Ideas as to how prime still returns 15,21,25,27,33,35,… ?
The removal loop appears to work and the first two loops do
appear to produce l=1, j=2 simultaneously, which would exclude the 7 in N and therefore the 15 in prime..

n<-100
N<-1:n

for (l in 1:floor(1/2 *((2*n+1)^.5-1))) {
for (j in l:((n-l)/(1+2*l))) {
for (b in 1:length(N)) {
if (N[b]==(l+j+2*l*j)) {
N<-N[-(b)]
}
}
}
}
prime<-2*N+1
prime

6. markehammons said

Scala version:

```  val n = (1000000 - 2)/2
val nonPrimes = for (i <- 1 to n; j <- i to (n - i) / (2 * i + 1)) yield i+j+(2*i*j)
val primes = 2 +:((1 to n) diff (nonPrimes) map (2*_+1))
println( primes )
```
```cartProd :: [a] -> [(a, a)]