Aside from increasing the ranges of y which is between -1000 to 1000, what other improvements should I make to these codes so that when I run it on any pari gp compiler online, it returns very valid answers?
When I ran these codes from within the ranges of -1000 to 1000 and from to -10000 to 10000, it returned some answers but nothing beyond these ranges. Though I am aware that between ranges of y: -10^12 to 10^12 for instance, there are some desirable solutions, like when y = 22 741 480 906, I do not know what else to do to arrive at such results systematically aside from increasing ranges of y as mentioned above.
{
for(y=-1000, 1000,
H= hyperellratpoints((16224+'w^3-624*'w*y-3*'w^2*y^2+1248*y^3+3*'w*y^4-3*y^6)/6, 10^7);
for(i=1, #H,
z= H[i][2];
if(z>0, if(z==floor(z),
w= H[i][1];
k= (w^2-624*y+3*y^4)/144;
if(k==floor(k),
x= (3*y^2-w)/12;
if(x==floor(x),
if(x!=0, if(12*x^2-6*y^2*x+y*(y^3-52)-12*k==0 && -z^2+(y^3-52)^2-288*k*x==0,
print("("x", "y", "z", "k")")
))
)
)
))
)
)
}
From a cursory look at your code, it seems you are looking at integer solutions of an equation 6*z^2 = P(y,w) for some explicit polynomial P with integer coefficients, and you further require two auxiliary quantities to be integral. You're doing this by trying all integral values of y in a certain range, computing for each such y all rational solutions of bounded height for the equation z^2=P/6, then filtering out using various tests (in particular z > 0).
To answer your main question: there's no ready-made function for finding integer points on general surfaces in PARI/GP.
What you can and should do is
y,w and x = (3y^2-w)/12 are integral, then 3 divides w = 3W and you now want x = (y^2 - W) / 4 to be integral which implies that 'y', 'w' and 'W' share the same parity (and in fact 4 | W when even, congruent to 1 mod 4 otherwise).You may want to split the problem in two cases at this point: y odd or even, etc. Once done as much as possible, use the new variables. All integrality tests can then be removed: if you have integral solutions in the new variables, they will automatically satisfy the divisibility conditions (and divisions by constants simplify).
I'm not sure why you're checking 12*x^2-6*y^2*x+y*(y^3-52)-12*k==0 and -z^2+(y^3-52)^2-288*k*x==0 as both follow formally from the definitions of k and x, hence are always true.
the special case x = 0 (or y = 3*w^2) provides trivial solutions since your equation becomes z = |y^3 - 52| (I understand you're excluding those).
It all boils down to finding integral solutions of an equation a^2 = Q(b,c) where everybody in sight is integral. You will either have to use more about the geometry of the surface (for instance understand why x = 0 leads to all these trivial solutions), or simply loop over b,c integral in some range and check whether the result is a square. This test is exceedingly fast and you can certainly check about 10^16 pairs (b,c) in reasonable time on a single computer. (Make that 10^20 if you have access to a computing cluster.)
A final remark: to test whether x is an integer in PARI/GP (which is no longer needed in the version I recommend) use either type(x) == "t_INT" or denominator(x) == 1 if you already know that x is a rational number.