3

Well, I have written the following code (using the fast square root test found in this answer):

Clear["Global`*"];
sQ[n_] := FractionalPart@Sqrt[n + 0``1] == 0;
r = 31265;
a = 2;
Monitor[Parallelize[
  While[True, 
   If[sQ[9 (-4 + r)^2 + 12 a (1 + a) (5 + a (-2 + r) - r) (-2 + r)], 
    Print[a]]; a++]; a], a]

Is there a quick and smart way to adjust this code, so that it chooses a value of r (which is given by a list {...,...,...}) and runs through the values of a which is given by a lower and upper bound?

So, for example I have for r: $\left\{16420, 19605, 31265, 31368, 83135\right\}$ and it needs to check $2\le a\le 10^9$. And when it finds a solution it need to print r and a.

Thanks a lot.

Jan Eerland
  • 2,001
  • 10
  • 17

2 Answers2

2
expr = 9 (-4 + r)^2 + 12 a (1 + a) (5 + a (-2 + r) - r) (-2 + r);

rlist = {16420, 19605, 31265, 31368, 83135};

Last@Reap@
  Table[If[sQ[expr], Sow[{a, r}], Nothing], {a, 2, 10^5}, {r, rlist}]

For this brief test run up to 10^5:

{{{259, 31265}, {1191, 19605}, {1310, 83135}, {6936, 16420}, {14858, 
   31368}}}

I am guessing that the task will become harder with the magnitude of a.

Syed
  • 52,495
  • 4
  • 30
  • 85
2

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 + ac3
    mpz_mul(n2,a,n2); mpz_add(n2,c1,n2); // n2 = c1 + a(c2 + ac3)
    mpz_mul(n2,a,n2); mpz_add(n2,c0,n2); // n2 = c0 + a(c1 + a(c2 + ac3))
    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
Roman
  • 47,322
  • 2
  • 55
  • 121