## 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 formi+j+ 2ijwhere integersiandjrange from 1 ≤i≤jandi+j+ 2ij≤n. For each remaining integerk, 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.

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]