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.

About these ads

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] + [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 )
    
  7. Haskellian said

    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]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 611 other followers

%d bloggers like this: