Re: question on trying to use quadratic residues to eliminate needless checks

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

* Bill Allombert [2025-03-24 22:10]:
> > I am trying to determine if using quadratic residues will shorten down the
> > time to check two numbers, a,b such that a^2 + b^2 = square (say s^2)
> > 
> > In regard to body cuboids, I looked at the side = 288807105787200
> > 
> > Using the divisors of side^2, we find that 1,492,627 integers or fractions
> > exist, such that x^2 + 288807105787200^2 = s^2 where x is in Z or Q
> > (actually n/2)
> > 
> > There are 1,262,992 integers with this property and 229,635 rationals (all
> > have denominator 2) with this property.
> > 
> > I am looking for two numbers r,s from this pile of 1,492,627 values such
> > that r^2 + s^2 = t^2 (t either integer or rational)
> 
> Assuming you do the computation in C, not in GP:
> 
> Pick k a small integer as large as possible so that the computation does not
> use too much memory (try k=7).
> 
> Take N = 4*vecprod(primes(k)).
> Compute a bit array B so that B[i]==1 iff i is a square modulo N.
> 
> Compute the array of r^2 mod N using C arithmetic for all r.
> 
> Compute i = r^2 + s^2 mod N using C arithmetic and check whether B[i]==1.
> If yes do the full check.
> 
> This strategy reduces the number of full tests by a factor 1/2^(k+1).

Presumably a bit less that that: a «full test» starts by successive
modular checks (more or less the ones you're emulating through B) each
would eliminate (heuristically) about half of the non squares. Hence,
most of the time, you're not paying the «full price» because the «full test»
is not run to completion.

Of course, doing it in GP is going to be a disaster.

Here's a different idea: your number 288807105787200 is highly composite
(it was chosen that way as a product of tiny primes <= 31 and that's why its
square has so many divisors)

- replace all your rational inputs r by 4*r^2; we obtain a set S of integers
  and want to test the pairs s < s' in S such that s + s' is a perfect square.
  (Of course 2*s = 8*r^2 is not a perfect square). Note that the complete
  factorization of all those s are already known from divisors(,1).

- whenever the valuations of s and s' at one p <= 31 differ, then the smallest
  one must be even else their sum can't be a square. In fact you can fix s,
  list the primes at which it has an odd valuation (in fact 1 for any
  prime p >= 7 by your choice of input number) and eliminate all s'
  which have a higher valuation at any such prime (meaning 2 for any
  prime p >= 7)

You can combine this with Bill's proposition (but not in GP).

But in fact, I don't believe it's worth the effort. All the integers s + s'
you want to test have less than 128 bits.

? b = [2^127, 2^128-1]; N = 10^7;
? for (i = 1, N, random(b)) \\ doing nothing, for reference
time = 1,087 ms.
? for(i = 1, N, issquare(random(b)))
time = 1,598 ms.

So about 0.5 seconds to process 10^7 random integers a bit larger than
your actual target size.

I would expect your 1113966934251 checks to require

? 1113966934251 / 10^7 * 0.5 / 3600
%1 = 15.471762975708333333333333333333333333

About 15 hours, on a single 2GHz processor, on my low-end aging laptop.

If you came up with a different estimate (days...), I assume something
inefficient is being done in your code. Maybe handling rational numbers
(squaring, addition, issquare) instead of my integers 4*r^2.

Cheers,

    K.B.
-- 
Pr. Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/