I think Mathematica is the wrong tool for this job. Here's a code written in C that is MUCH faster (about one minute per $r$-value up to $a=10^9$). I'm using 128-bit integers, which limit the ranges of $a$ and $r$ somewhat; see GMP solution further down.
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.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;
}
}
static int check3(const myint a, const myint coeff[4]) {
if (isqrt(coeff[0] + a * (coeff[1] + a * (coeff[2] + a * coeff[3]))) == -1)
return EXIT_FAILURE;
return EXIT_SUCCESS;
}
int main(int argc, char *argv[]) {
long long int rr;
if ((argc != 2) || (sscanf(argv[1], "%lld", &rr) != 1)) {
fprintf(stderr, "usage: %s <r>\n", argv[0]);
return EXIT_FAILURE;
}
myint r = rr;
myint coeff[4] = {9 * (r - 4) * (r - 4), -12 * (r - 2) * (r - 5),
36 * (r - 2), 12 * (r - 2) * (r - 2)};
for (myint a = 2; a <= 1000000000; a++)
if (check3(a, coeff) == EXIT_SUCCESS)
printf("%" PRId64 " %" PRId64 "\n", (uint64_t)r, (uint64_t)a);
return EXIT_SUCCESS;
}
Here I saved the above code as perfectsquare.c
, compiled it with GCC, and ran it for several values of $r$:
--------% gcc -Ofast perfectsquare.c -o perfectsquare
--------% time ./perfectsquare 16420
16420 6936
./perfectsquare 16420 46.55s user 0.13s system 99% cpu 46.888 total
--------% time ./perfectsquare 19605
19605 1191
./perfectsquare 19605 49.36s user 0.17s system 99% cpu 49.568 total
--------% time ./perfectsquare 31265
31265 259
31265 325791721
./perfectsquare 31265 62.85s user 0.16s system 99% cpu 1:03.03 total
--------% time ./perfectsquare 31368
31368 14858
./perfectsquare 31368 62.76s user 0.22s system 99% cpu 1:03.01 total
--------% time ./perfectsquare 83135
83135 1310
./perfectsquare 83135 111.60s user 0.28s system 99% cpu 1:51.91 total
update: GMP version without speed penalty
For using larger numbers (going beyond 128-bit integers) we can use the GNU Multiple Precision Arithmetic Library. It turns out that it is not slower than GCC's built-in 128-bit integer type!
#include <stdio.h>
#include <gmp.h>
#include <stdlib.h>
int main(int argc, char *argv[]) {
if (argc != 3) {
fprintf(stderr, "usage: %s <r> <amax>\n", argv[0]);
return EXIT_FAILURE;
}
mpz_t r, amax, r2, r4, r5, c0, c1, c2, c3, a, n2;
mpz_init_set_str(r, argv[1], 10);
mpz_init_set_str(amax, argv[2], 10);
mpz_init(r2); mpz_sub_ui(r2,r,2); // r2 = r-2
mpz_init(r4); mpz_sub_ui(r4,r,4); // r4 = r-4
mpz_init(r5); mpz_sub_ui(r5,r,5); // r5 = r-5
mpz_init(c0); mpz_mul(c0,r4,r4); mpz_mul_si(c0,c0,9); // c0 = 9*(r-4)^2
mpz_init(c1); mpz_mul(c1,r2,r5); mpz_mul_si(c1,c1,-12); // c1 = -12*(r-2)*(r-5)
mpz_init(c2); mpz_mul_si(c2,r2,36); // c2 = 36*(r-2)
mpz_init(c3); mpz_mul(c3,r2,r2); mpz_mul_si(c3,c3,12); // c3 = 12*(r-2)^2
mpz_init(n2);
for (mpz_init_set_si(a,2); mpz_cmp(a,amax) <= 0; mpz_add_ui(a,a,1)) {
mpz_mul(n2,a,c3); mpz_add(n2,c2,n2); // n2 = c2 + a*c3
mpz_mul(n2,a,n2); mpz_add(n2,c1,n2); // n2 = c1 + a*(c2 + a*c3)
mpz_mul(n2,a,n2); mpz_add(n2,c0,n2); // n2 = c0 + a*(c1 + a*(c2 + a*c3))
if (mpz_perfect_square_p(n2)) {
mpz_out_str(stdout, 10, a);
printf("\n");
}
}
return EXIT_SUCCESS;
}
Save as perfectsquare_gmp.c
and compile with
gcc -I/opt/local/include -Ofast perfectsquare_gmp.c -L/opt/local/lib -lgmp -o perfectsquare_gmp
(on macOS with MacPorts; linking will be different on different OS types), and now very large numbers are possible (limited by memory & patience):
--------% time ./perfectsquare_gmp 31265 1000000000
259
325791721
./perfectsquare_gmp 31265 1000000000 61.81s user 0.10s system 99% cpu 1:01.95 total