## 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* + 2*ij* ≤ *n* implies *j* < (*n*−*i*)/(2*i*+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.

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)

}

}

}

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.

My C++ Solution

http://ideone.com/m6MAW

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

Scala version:

In Haskell:

cartProd :: [a] -> [(a, a)]

cartProd xs = [(x,y) | x <- xs, y <- xs]

sieveSundaram :: Integer -> [Integer]

sieveSundaram n = map ((+1) . (*2)) $ filter valid [1..n]

where valid x = not . any (\(i,j) -> i + j + 2*i*j == x) . filter (uncurry (<=)) $ cartProd [1..x]