[seqfan] Testing potential squares

[Keith F. Lynch]

I recently wanted to test a large number of large numbers to see
if they were squares.  I wanted to do so as quickly as possible.
The odds of a positive integer around N being square are about
1/(2*sqrt(N)).  So if N is around a trillion (10^12), the odds
are about 1 in 2 million.

The obvious way to test is to take the square root.  But that's slow.

To quickly rule out most non-squares, I created, just once, at the
beginning of the program, a table of all quadratic residues mod 2^16
(65536).  Then I tested every potential square by ANDing it with
2^16-1 (65535), and looking up the result in that table, both of which
are very fast.  That rules out nearly 5/6 of all non-squares.  The
proportions of quadratic residues modulo higher and higher powers of 2
quickly approach 1/6.  (2^16 already gets 1/5.999268.  See A023105.)

To rule out most of the remaining non-squares, I searched for numbers
that had a lower proportion of quadratic residues than any smaller
number (e.g. 1441440, the first number such that fewer than one
percent of numbers are quadratic residues modulo it), intending to
create a table of residues mod such a number and test every potential
square by calculating it mod that number and seeing if it was in the
table.  After calculating the first few numbers that have a record
low proportion of quadratic residues, I looked them up in OEIS, but
couldn't find them.  Later I found they're at A085635; I had missed
it because it lists the first term as 2 when I think it should be 1.

Should I correct that, or is there a reason it's 2?  Thanks.  I also
have more terms to add to it and to its companion sequence A084848.
They currently show 44 terms, the last of which is around a hundred
thousand and has a ratio of about 50.  I've computed 193 terms so
far, of which the last is about two thirds of a trillion and has a
ratio of about 1588.

Later I decided that, having already weeded most non-squares out mod
2^16, I think I should instead look for odd numbers that have a lower
proportion of quadratic residues than any smaller odd number.  That
list is not in OEIS.  I'll add it as soon as someone tells me whether
I should include 1 as its first term, and if not, why not.  Thanks.

Those two tests (mod 65536 and mod 740025 (the first odd number to
have fewer than 3% of numbers modulo it to be residues)), together
leave fewer than 1/200 of non-squares.  So instead of one in two
million, squares are one in ten thousand.  So I only need try taking
the square root of ten thousand numbers for every number that turns
out to actually be a square.  Still not very good.

If I have one gigabyte to devote to the task, I can instead test my
prospective square mod 7661478825, a larger record-setting odd number,
which has about 1/233 residues.  That requires one bit per residue
to get it to fit in one gigabyte.  If N is my prospective square mod
7661478825, I shift N right 3 bits to get the byte to look up in the
table.  There will be about a 1/29 chance of a hit, i.e. at least one
residue in the block of 8 containing the correct bit.  Only if there
is a hit do I check the bit, which I do by ANDing N with 7, and using
that as the index of which bit in that byte to look at.  So I've ruled
out all but about 1/1400 non-squares.

I can do better.  Suppose I look at the record low proportions of
residues, not for all odd numbers, but separately for numbers all of
whose prime factors are congruent to 1 mod 4 and for numbers all of
whose prime factors are congruent to 3 mod 4.  Since numbers in one of
those sets will all be relatively prime to all numbers in the other of
the two sets, and since for any two relatively prime numbers A and B,
the number of residues mod AB will always be the product of the number
of residues mod A with the number of residues mod B, the result will
be the same by testing mod A followed by testing mod B as by testing
mod A*B.  (A table of size A*B would require a prohibitive amount
of memory.)

If, once again, I have one gigabyte for the purpose, I can test mod
3877273323 and (if it passes that test) mod 3833254945.  The former
is the largest record-holding number less than 4*10^9 all of whose
prime factors are congruent to 3 mod 4, and the latter is the largest
record-holding number less than 4*10^9 all of whose prime factors are
congruent to 1 mod 4.  The factorizations are 3^2 7^2 11 19 23 31 and
5 13 17 29 37 53 61 respectively.  And the proportions of quadratic
residues are about 1/127 and about 1/85 respectively, or about
1/10800.  Combined with the first test, mod 65536, I've ruled out all
but about 1/65000 squares.  And testing two mod tables rather than one
doesn't slow things down much, since only if it passes the first one
do I try the second one.  So I only need try taking the square root
of about 30 numbers for every number that turns out to actually be a
square (assuming the prospective squares are around a trillion).

I could do even better by having four tables, one for numbers all of
whose prime factors are mod 1 mod 8, another for numbers all of whose
prime factors are mod 3 mod 8, and likewise for 5 and 7.  Or eight
tables, one for each possible remainder mod 16.  (Again, these are for
all of the prime factors of the large number with the record low (for
its size) proportion of quadratic residues, not for that number itself.)

I looked into instead checking mod some power of each successive
prime.  After checking mod 2^16 (with an AND) I could check mod 3^16,
then 5^16, then 7^16, etc.  But that would require more taking of mods,
since in the limit of high powers, powers of 2 have 1/6 residues, and
powers of odd primes have P/(2(P+1)) residues, i.e. 3/8, 5/12, 7/16,
etc.  That's 1/2 in the limit of large power of large primes.

That's not a ridiculously bad approach, especially if you're short
on memory.  The average number of mods, assuming that the potential
square passed the 65536 test, would be the sum of the partial products
of each odd prime over two more than the prime, i.e.
1 + 3/8 + (3/8 * 5/12) + (3/8 * 5/12 * 7/16) + (3/8 * 5/12 * 7/16 * 11/24) + ...
which comes to about 1.658653437258581054.  For comparison, my two odd
mods approach would, with the numbers given, require 1.0079 mods on
average.

(Of course in every algorithm you'd stop doing mods and take the square
root once the probability of the number being square got high enough.
So only the first few digits of my 1.658... number above are meaningful
unless you're testing truly stupendous potential squares.  But anything
worth doing is worth overdoing.  Especially since that presumably-
transcendental prime-related number doesn't seem to have ever appeared
anywhere before.)

I'm curious whether anyone can think of a faster way to test large
numbers for squareness.

I'm also curious what the ultimate behavior of the A085635/A084848
proportions is.  I know they grow without limit.  But at what rate
do the record-setting ratios increase?  Thanks.