Sieve of Sundaram (for finding prime numbers)

Question

The Challenge

Implement the Sundaram sieve for finding prime numbers below n. Take an input integer, n, and output the prime numbers below n. You can assume that n will always be less than or equal to one million.

---

Sieve

1. Start with a list of the integers from 1 to n.

2. Remove all numbers that are in the form i + j + 2ij where:

   • i and j are less than n. j is always greater than or equal to i, which is greater than or equal to 1.
     
     

   • i + j + 2ij is less than or equal to n
     
     

3. Multiply the remaining numbers by 2, and add 1.

This will yield all the prime numbers (except 2, which should be included in your output) less than 2n + 2.

---

Here is an animation of the sieve being used to find primes below 202.

---

Output

Your output should be every prime integer ≤ n (in ascending order) followed by a newline:

2
3
5

Where n is 5.

---

Examples

> 10
2
3
5
7

> 30
2
3
5
7
11
13
17
19
23
29

Inputs are denoted by >.

Answer 1

Pyth, 23 bytes

2j@JSQmhyd-Jm+sdy*Fd^J2

Demonstration

Really just implements the algorithm as given.

Answer 2

Haskell, 93 90 bytes

import Data.List
g n=unlines[show$2*x+1|r<-[[1..n]],x<-2:(r\\[i+j+2*i*j|j<-r,i<-r]),2*x<n]

How it works: [i+j+2*i*j|j<-r,i<-r] are all i+j+2ij which are removed (\\) from [1..n]. Scale to 2x+1 and turn them into a string (show). Join with NL (unlines).

Answer 3

Scala, 115 124 122 115 114 bytes

n=>{println(2);for{m<-1 to n;if !(for{j<-1 to n;i<-1 to j}yield i+j+2*i*j).contains(m);o=2*m+1;if o<=n}println(o)}

An anonymous function; takes n as an argument and prints the result to stdout.

Answer 4

JavaScript (ES7), 107 105 bytes

Array comprehensions are awesome! But I wonder why JS has no range syntax (e.g. [1..n])...

n=>{for(a=[i=1];i<n;a[i++]=i);for(i=0;i++<n;)for(j=0;j<n;a[i+j+++2*i*j]=0);return[for(i of a)if(i)i*2+1]}

This was tested successfully in Firefox 40. Breakdown:

n=>{
  for(a=[i=1];i<n;a[i++]=i); // fill a list with 1..n
  for(i=0;i++<n;)            // for each integer i in 0..n
    for(j=0;j<n;)            //   for each integer j in 0..n
      a[i+j+++2*i*j-1]=0;    //     set the corresponding item of the list to 0
  return[for(i of a)         // filter the list by:
          if(i)              //   item != 0 AND item != undefined
           i*2+1]            // and return each result * 2 + 1
}

Alternative, ES6-friendly solution (111 bytes):

n=>{for(a=[i=1];i<n;a[i++]=i);for(i=0;i++<n;)for(j=0;j<n;a[i+j+++2*i*j]=0);return a.filter(x=>x).map(x=>x*2+1)}

Suggestions welcome!

Answer 5

MATLAB, 98

n=1:input('');m=n;for p=m for i=1:p j=i:p;for k=i+j+2*i*j n(n==k)=[];end;end;end;disp(2*n'+1);

And in a readable form

n=1:input(''); %Ask for the input number (e.g. 100) and form a range
m=n; %Back up the range as we will be editing 'n', but need 'm' as a loop list
for p=m %For each number between 1 and n inclusive
    for i=1:p %'i' is all numbers greater than or equal to 1 up to p
        j=i:p; %'j' is all numbers greater than or equal to i up to p
        for k=i+j+2*i*j %Calculate the numbers to remove, and loop through them
            n(n==k)=[]; %Remove that value from the 'n' array
        end
    end
end
disp([2;2*n'+1]); %An display the list including the number 2 seperated by a new line.

Answer 6

Java8 : 168 165 bytes

N->{int[]A=new int[N*N];int i=1,j;N=N/2;for(;i<N;i++)for(j=i;j<N;)A[i+j+2*i*j++]=1;System.out.println(N>1?2:\"\");for(i=1;i<N;i++)if(A[i]<1)System.out.println(2*i+1);}

For more bigger number data type with wide range can be used. We do not need to iterate for whole N indexes N/2 is sufficient.

To understand properly following is the equivalent method.

static void findPrimeSundar(int N){
    int[] A = new int[N*N];
    int i=1,j;
    N=N/2;
    for(;i<N;i++)
      for(j=i;j<N;)
        A[i+j+2*i*j++]=1;
    System.out.println(N>1?2:"");
    for(i=1;i<N;i++)
        if(A[i]<1)System.out.println(2*i+ 1);
}

Answer 7

CJam, 35 bytes

2li:V,:)__2m*{_:+\:*2*+}%m2f*:)&+N*

Try it online

This seems somewhat lengthy relative to isaacg's Pyth solution, but it's... what I have.

Explanation:

2       Push a 2, will be part of final output.
li      Get input and convert to integer n.
:V      Save in variable V for later use.
,       Generate list [0 ... n-1].
:)      Increment list elements to get list [1 ... n].
__      Create two copies, one for sieve, and for clamping results.
2m*     Cartesian power, generating all i,k pairs.
{       Loop over all i,j pairs.
  _     Copy pair.
  :+    Calculate sum i + j.
  \     Swap copy of pair to top.
  :*    Calculate product i * j.
  2*    Multiply by 2, to get 2 * i * j.
  +     Add both values, to get i + j + 2 * i * j.
}%      End loop over all i,j pairs.
m       Sieve operation, remove the calculated values from the list of all values.
2f*     Multiply the remaining values by 2...
:)      ... and add 1 to the. We now have the list of all primes up to 2 * n + 2.
&       Intersect with [1 ... n] list, because output is only values <= n.
+       Concatenate with the 2 we pushed at the start.
N*      Join with newlines.

Answer 8

Perl 6, 96 bytes

If I strictly follow the description the shortest I managed to get is 96 bytes.

->\n {$_=@=1..n;for 1..n {for $^i..n {.[$i+$^j+2*$i*$j-1]=0}};2,|.[0..n].map(* *2+1).grep(3..n)}

->\n {
  $_=@=1..n; # initialize array
  for 1..n { # $i
    for $^i..n { # $j
      .[$i+$^j+2*$i*$j-1]=0 # remove value
    }
  };
  2,|.[0..n].map(* *2+1).grep(3..n)
}

---

If I could do the 2n + 1 on initialization of the array, pre-inserting 2, and limiting that to only the values less than or equal to n; it can be reduced to 84 bytes.

->\n {$_=@=2,{++$*2+1}...^*>n;for 1..n {for $^i..n {.[$i+$^j+2*$i*$j]=$}};.grep(?*)}

If I also ignore that j is supposed to be at least i, I can reduce it to 82 bytes.

->\n {$_=@=2,{++$*2+1}...^*>n;for 1..n X 1..n ->(\i,\j){.[i+j+2*i*j]=$};.grep(?*)}

---

Example usage:

my $code = ->\n {...} # insert one of the lambdas from above

say $code(30).join(',');
# 2,3,5,7,11,13,17,19,23,29

my &code = $code;
say code 11;
# (2 3 5 7 11)

Answer 9

PHP, 126 Bytes

$r=range(1,$n=$argn/2-1);for(;++$i**2<=$n;)for($j=$i;$n>=$d=$j+$i+2*$i*$j++;)unset($r[$d-1]);foreach($r as$v)echo 1+2*$v."\n";

Online Version

Answer 10

Julia 0.6, 65 bytes

n->[2;(p=setdiff(1:n,[i+j+2i*j for i=1:n for j=i:n])*2+1)[p.<=n]]

Try it online!

Not a great challenge in terms of golfing, but I just had to do it for the name. :)

Answer 11

Husk, 24 bytes

¶↑≤¹:2m+1mDṠ-om§+oDΠΣπ2ḣ

Try it online!

Ok, I'm not sure why increment → isn't working properly here. Once that works, this will be 23 (tied with Pyth).