# Prime number sieve using difference of two squares

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 ?

• I presume you mean apart from the obvious of pre-testing the array element to see if it has been zeroed already? The number of different ways of calculating a composite number from the difference of two squares depends directly on the number of prime factors the number has. To discover the number of prime factors you either have to factorise the composite number or to try to factorise the difference of two squares found. – James Arathoon Aug 24 '17 at 9:00
• If you put two primes into their difference of squares form and multiply them to find the difference of two squares for the resultant composite number you will see that along the way information that would quickly allow you to undo what has just been done is lost. – James Arathoon Aug 24 '17 at 9:13

Not an answer just some elementary results to help clarify my comments above.

Multiplying two numbers using difference of two squares notation.

Let $N=n_1n_2$, where $n_1=p_1q_1$ and $n_2=p_2q_2$

Also let $p_1=(m_1+n_1)$, $q_1=(m_1-n_1)$, $p_2=(m_2+n_2)$ and $q_2=(m_2-n_2)$

$N=n_1 \times n_2$ can be written in three different ways using different pairings of $p$ and $q$, thus

$N=p_1p_2 \times q_1q_2=p_1q_1 \times p_2q_2=p_1q_2 \times p_2q_1$

a) First pairing $N=p_1p_2 \times q_1q_2$

\begin{align} N&=p_1p_2 \times q_1q_2\\ &= \left( \frac{p_1p_2+q_1q_2}{2}\right)^2-\left(\frac{p_1p_2-q_1q_2}{2} \right)^2\\ &=\left( m_1m_2+n_1n_2 \right)^2-\left( m_1n_2+n_1m_2 \right)^2 \end{align}

b) Second pairing $N=p_1q_1 \times p_2q_2$

\begin{align} N&=p_1q_1 \times p_2q_2\\ &= \left( \frac{p_1q_1+p_2q_2}{2}\right)^2-\left(\frac{p_1q_1-p_2q_2}{2} \right)^2\\ &=\left( \frac{\left(m_1^2-n_1^2 \right)+\left(m_2^2-n_2^2\right)}{2} \right)^2-\left( \frac{\left(m_1^2-n_1^2 \right)-\left(m_2^2-n_2^2 \right)}{2} \right)^2 \end{align}

c) Third pairing $N=p_1q_2 \times p_2q_1$

\begin{align} N&=p_1q_2 \times p_2q_1\\ &= \left( \frac{p_1q_2+p_2q_1}{2}\right)^2-\left(\frac{p_1q_2-p_2q_1}{2} \right)^2\\ &=\left( m_1m_2-n_1n_2 \right)^2-\left( n_1m_2-m_1n_2 \right)^2 \end{align}

From these calculations it can be seen that $$(m_1m_2 \times n_1n_2)=(n_1m_2 \times m_1n_2)$$