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

Question

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.

Answer 1

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 + 0``1] == 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