# Writing a program that finds for what $(x,y)$ a function gives a perfect square number

The overal question I am trying to answer is:

For what $$(x,y)$$, which are positive integers, is the following number a perfect square number?

$$9 \left(x^3 (y-2)^2+3 x^2 (y-2)-2 x (y-45) (y-2)+7 (y-1)^2\right)\tag1$$

Now, I am using the following code:

ParallelTable[
If[IntegerQ[
Sqrt[9 (x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) +
7 (y - 1)^2)]], {x, y}, Nothing], {x, 1, 100}, {y, 1,
100}] /. {} -> Nothing


But that is a bit slow for bigger values of $$x$$ and $$y$$. Is there a faster way to test when the number is a perfect square.

• Note that your overall factor of 9 is redundant - it will not affect whether the result is square or not. Commented Sep 12, 2020 at 13:25
• Try FindInstance[7353 + 9 y (-1841 + 718 y) == z^2, {z, y}, PositiveIntegers]. You can get extremely large numbers this way. The equation is just what you get when x -> 9. Fixing a variable like this appears to make it a bit easier for FindInstance. I got {x->9, y -> 959638328, z -> 77142029727}. Here is a solution with x->9 and enormous y and z pastebin.com/jWvWBUpx which you can get by asking for more results from the above. Commented Sep 12, 2020 at 14:47

You can do

FindInstance[n^2 == 9 (x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2),
{x, y, n}, PositiveIntegers]

(*    {{x -> 3, y -> 5, n -> 102}}    *)


to find an exemplary instance.

If you want all solutions up to $$x,y\le s$$, you can do

With[{s = 20},
Solve[{n^2 == 9 (x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2),
1 <= x <= s && 1 <= y <= s}, {x, y, n}, PositiveIntegers]]

(*    {{x -> 1, y -> 8, n -> 87},
{x -> 3, y -> 5, n -> 102},
{x -> 3, y -> 8, n -> 159},
{x -> 9, y -> 8, n -> 537}}    *)


Faster: using the squareness test of this answer and a Sow/Reap combination, and eliminating the prefactor of 9 (see @mikado's comment):

sQ[n_] := FractionalPart@Sqrt[n + 01] == 0
With[{s = 1000},
Reap[Do[If[
sQ[x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2],
Sow[{x, y}]], {x, s}, {y, s}]][[2, 1]]]

(*    {{1, 8}, {1, 128}, {3, 5}, {3, 8}, {9, 8}, {11, 1}, {47, 8}}    *)


Further, using the parallelization trick of this Q&A,

SetSharedFunction[ParallelSow];
ParallelSow[expr_] := Sow[expr]

With[{s = 10^4},
Reap[ParallelDo[
If[sQ[x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2],
ParallelSow[{x, y}]], {x, s}, {y, s}]][[2, 1]]]

(*    {{1, 8}, {1, 128}, {1, 1288}, {3, 5}, {3, 8}, {9, 8}, {11, 1}, {47, 8}}    *)


Other than that, this kind of calculation is done much more efficiently in a low-level language like C. Here's my attempt to use pure C with 128-bit integers (because for $$x=y=10^6$$ we overflow 64-bit integers), going up to $$s=10^6$$ in about two hours:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <inttypes.h>

typedef __int128 myint;

static myint
compute_isqrt(const myint x)
{
myint r = sqrt(x);
while (r*r <= x) {
if (r*r == x)
return r;
r++;
}
return -1;
}

static myint
isqrt(const myint x)
{
if (x < 0)
return -1;
switch(x & 0xf) {
case 0:
case 1:
case 4:
case 9:
return compute_isqrt(x);
default:
return -1;
}
}

#define M 1000000

int main() {
for (myint x=1; x<=M; x++)
for (myint y=1; y<=M; y++) {
myint z = x*x*x*(y-2)*(y-2)+3*x*x*(y-2)-2*x*(y-45)*(y-2)+7*(y-1)*(y-1);
myint n = isqrt(z);
if (n >= 0) {
printf("%" PRId64 " %" PRId64 " %" PRId64 "\n",
(int64_t)x, (int64_t)y, (int64_t)n);
}
}
return EXIT_SUCCESS;
}


Save as perfectsquare.c, compile with

gcc -Wall -O3 perfectsquare.c -o perfectsquare


and run with

time ./perfectsquare


Here are all the solutions $$\{x,y,n/3\}$$ up to $$s=10^6$$:

1 8 29
1 128 329
1 1288 3171
1 13168 32271
1 126848 310729
3 5 34
3 8 53
3 42680 225859
3 61733 326678
3 476261 2520154
3 688856 3645101
9 8 179
11 1 0
47 8 1949
15577 8 11664979

• Using your method with $s=10^6$ takes already way to long. Commented Sep 12, 2020 at 13:43
• @Jan try again, please. Commented Sep 12, 2020 at 13:56