[seqfan] Re: Testing potential squares

[Brendan McKay]

All modern general-purpose CPU chips have a single instruction
for the square root of a floating-point number.  If you convert N
to floating point, take the square-root, and round to an integer,
you will have the only integer whose square might be N provided
N has less than twice the number of bits in a floating-point mantissa.
This will happen within several machine instructions provided you
use a language (C, C++, probably Java) that optimizes floating-point
arithmetic.  My Linux and Mac computers both have "long double"
type with at least 80 bits of mantissa, which will get you to well
over 10^40 (use the sqrtl() library function).

Brendan.

On 21/8/18 10:09 pm, Keith F. Lynch wrote:
> 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.
>
> --
> Seqfan Mailing list - http://list.seqfan.eu/