3
$\begingroup$

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.

$\endgroup$

2 Answers 2

2
$\begingroup$
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.

$\endgroup$
1
  • 2
    $\begingroup$ Do instead of Table would probably be faster. $\endgroup$
    – Roman
    Commented Apr 11, 2023 at 18:31
2
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.