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