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 ≤ i ≤ j and i + j + 2ij ≤ n. 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.

Posted by programmingpraxis

Filed in Exercises

7 Comments »

7 Responses to “Sieve Of Sundaram”

khigia said
October 7, 2011 at 4:29 PM
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)
}
}
}

Mike said

October 7, 2011 at 7:44 PM

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)]

murkey said

October 8, 2011 at 6:14 AM

#!/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)

CyberSpace17 said
October 21, 2011 at 8:53 PM
My C++ Solution
http://ideone.com/m6MAW
deang said
October 23, 2011 at 5:10 PM
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

markehammons said

November 30, 2012 at 4:42 AM

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 )

Haskellian said

May 5, 2014 at 10:36 PM

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]

Programming Praxis