Prime number sieve using difference of two squares Here is a working pari/gp implementation of prime number sieve. Is there something that I could change in code to achieve a better running time?

Sieve(g)=
{
A=vector(floor((g-1)/2));
for(n=1,floor((g-1)/2),
   A[n]=2*n+1);
for(x=3,floor(sqrt(g)),
   if(Mod(x,2)==1,i=0,i=1);
   forstep(y=i,x-3,[2],
      A[(x^2-y^2-1)/2]=0));
for(x=ceil(sqrt(g)),ceil((g+9)/6),
   i=ceil(sqrt(x^2-g));
   if(Mod(x+i,2)==0,i=i+1);
   forstep(y=i,x-3,[2],
      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])))
}
Good comments make code more maintainable. I had to reverse engineer what this code was doing, when a few comments could have made it clear.



For my taste the code could use a lot more whitespace, but I'm going to try to stick to your style in the proposed changes.



Separation of responsibilities: it would be cleaner to return a vector rather than to "return" with print.



There's a lot in common between the two main loops:


for(x=3,floor(sqrt(g)),
   if(Mod(x,2)==1,i=0,i=1);
   forstep(y=i,x-3,[2],
      A[(x^2-y^2-1)/2]=0));
for(x=ceil(sqrt(g)),ceil((g+9)/6),
   i=ceil(sqrt(x^2-g));
   if(Mod(x+i,2)==0,i=i+1);
   forstep(y=i,x-3,[2],
      A[(x^2-y^2-1)/2]=0));



It might clarify things to refactor them into one:

for(x=3,floor((g+9)/6),
  forstep(y=x-3,0,-2,
    if(x^2-y^2&gt;g,break);
    A[(x^2-y^2-1)/2]=0
  ));


Note that this fixes an out-by-one in the upper bound of x.

Now I'm ready to look at performance. The low-hanging fruit here is x^2-y^2: it's usually faster to calculate a constant product or square of a loop index by addition. I'm also going to use bitwise operations for multiplication and division by 2 because that means not having to worry about how clever the compiler is.

for(x=3,floor((g+9)/6),
  d=(x&lt;&lt;1)-1;
  forstep(y=x-3,0,-2,
    d+=(y+1)&lt;&lt;2;
    if(d&gt;g,break);
    A[d&gt;&gt;1]=0
  ));


You should easily be able to prove the invariant that d == x^2 - y^2, and since it's odd (d-1)/2 == d&gt;&gt;1. In fact, we can instead work with (x^2-y^2)&gt;&gt;1:

cap=(g-1)&gt;&gt;1;
for(x=3,floor((g+9)/6),
  d=x-1;
  forstep(yp=(x-2)&lt;&lt;1,2,-4,
    d+=yp;
    if(d&gt;cap,break);
    A[d]=0
  ));


With Try It Online! and pushing towards the size limit for allocation of A I see a speedup from 4.4s (user) to 3.5s. (Output removed for timing purposes. Timings available by expanding the "Debug" section).



Other optimisations:


If the vector is just initialised as vector((g-1)&gt;&gt;1) then it will initially contain 0 at all indices. Then change A[d]=0 to A[d]=1 and change the logic in the output loop.
Ditch the concat and output 2 as a special case.


Beyond that, you need to look at different algorithms. The binary quadratic form x^2-y^2 doesn't seem to be a common one in prime sieving. You could see how it compares with Sundaram's sieve (BQF x+y+2xy) or failing that, if you're working with large enough ranges, Atkin-Bernstein (e.g. a combination of 3x^2+y^2, 4x^2+y^2, 5x^2+3y^2, 15x^2+y^2, and 3x^2-y^2 for different totients of 60).