## Sieve Of Sundaram

### October 7, 2011

In 1934, an Indian math student, S. P. Sundaram, invented a new method of sieving primes:

From a list of integers from 1 to n, remove all integers of the form i + j + 2ij where integers i and j range from 1 ≤ ij and i + j + 2ijn. For each remaining integer k, the integer 2k+1 is prime, and the list gives all the odd primes (thus excluding the prime 2).

Your task is to implement the sieve of Sundaram and calculate the list of primes to a million. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

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*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 )
```