The following prime number sieve is based on the representation of an odd integer as the difference of two squares . The algorithm uses a fact that all odd composite numbers can be expressed as difference $x^2-y^2$ for some integers $x,y$ such that $0 \le y \le x-3$ . Note that prime numbers cannot be expressed in this way .
Pseudocode :
Input : an integer g > 1
if mod(g,2) == 1 then g := g+1
A := array of length ceiling((g-1)/2)
for n from 0 to floor((g-1)/2) :
A[n] := 2n+1
A[0] := 2
for x from 3 to floor((g+9)/6) :
for y from x-3 to 0 by step -2 :
if x^2-y^2 > g then break
A[(x^2-y^2-1)/2] := 0
Output : all nonzero elements of array A
PARI/GP implementation :
Sieve(g)=
{
A=vector(floor((g-1)/2));
for(n=1,floor((g-1)/2),
A[n]=2*n+1);
for(x=3,floor((g+9)/6),
forstep(y=x-3,0,[-2],
if(x^2-y^2>g,break);
A[(x^2-y^2-1)/2]=0));
A=concat(2,A);
for(j=0,floor((g-1)/2),
if(!(A[j+1]==0),print(A[j+1])))
}
You can run this code here .
Question :
When g is big enough you may notice that some elements of array are equalized to zero more than once . For example $45$ , because $45=7^2-2^2$ and $45=9^2-6^2$ . Is there something I could change in the algorithm to achieve that each element of array which corresponds to composite number is equalized to zero exactly once ?