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 + 2ijn implies j < (ni)/(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.

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 613 other followers

%d bloggers like this: