In a previous post, fast way to tell whether a number is a square, the question came up of how many integers are squares modulo an integer m. That is, for how many numbers y in the list 0, 1, 2, …, m-1 does the equation x2 ≡ y (mod m) have a solution. The application in the earlier post was only concerned with m equal to a power of 2, but it turns out there is a completely general solution.
Michael Lugo pointed to a paper (*) that gives a solution for all values of m. I will summarize the paper here and give a Python script implementing the solution.
The paper defines the function s(n) as the number of distinct squares modulo n. The first observation is that s(n) is a multiplicative function. That means that if m and n are relatively prime, s(mn) = s(m) s(n). Since every integer can be factored into prime powers, knowing s() is multiplicative reduces the problem to computing s(pn) for primes p. The paper then shows how to compute s() for powers of odd primes and powers of 2. The full solution is summarized in the code below.
Assume we have factored our number m. For example, suppose we want to know how many possible last three digits there are for squares base 10. That is equivalent to computing s(1000), and so we factor 1000 into 23 53. We represent that factorization as a list of pairs: [ (2, 3), (5, 3) ]. The function squares below tells us that s(1000) is 159.
def isOdd( n ):
return n%2
def s( factor ):
(p, n) = factor
if isOdd( p ):
if n == 1:
r = (p+1)/2
elif n == 2:
r = (p*p - p + 2)/2
elif isOdd( n ):
r = (p**(n+1) + 2*p + 1)/(2*p + 2)
else:
r = (p**(n+1) + p + 2)/(2*p + 2)
else:
if isOdd( n ):
r = (2**(n-1) + 5)/3
else:
r = (2**(n-1) + 4)/3
return r
def squares( factorization_list ):
product = 1
for factor in factorization_list:
product *= s(factor);
return product
# test cases to exercise each branch
assert squares( [(5, 1)] ) == 3
assert squares( [(5, 2)] ) == 11
assert squares( [(5, 3)] ) == 53
assert squares( [(5, 4)] ) == 261
assert squares( [(2, 3)] ) == 3
assert squares( [(2, 4)] ) == 4
assert squares( [(2,3), (5,3)] ) == 159
The importance of this result to the original problem is that looking at the last several bits of a number can screen out no more than 5/6 of the non-squares. Looking at the last four bits, i.e. the last hexadecimal digit, screens 3/4 of the non-squares, so there’s not room to do much better.
Update: So about 1 out of every 6 integers is a square mod 2n for large n. What about mod pn for odd prime p? The corresponding ratio is 1/2.
Update (21 April 2010): Fixed a couple typos. Added syntax coloring to the source code. See an implementation of the algorithm in J by Tracy Harms.
* W. D. Stangl, “Counting Squares in Z_n”, Mathematics Magazine, pp. 285-289, Vol. 69 No. 4 October 1996.
Well. In computers, we use fixed-length integers, especially if we want speed. You have either 32 bits or 64 bits. So your statement is false (from a computer programming point of view). We can always tell whether a number is a square from the 32 or 64 bits… ;-) This is not just a common routine, this is a fundamental part of our CPU architecture (which is register-based).
I find it fascinating that one integer out of 6 is a square.
As I understand we have pairs (p, n), where p — factor of our module, that is prime, n — it’s power. Out of all prime numbers only 2 is even. Why not just check p if it equals 2 instead? Or it doesn’t matter and was used for generality?
Of course, you can replace isOdd(p) by p != 2. The semantics of the resulting program don’t change (assuming correct input), i.e. it does not really matter.
There is a small difference, though: depending on how p is stored, isOdd(p) can be faster to determine than p != 2.
isOdd(p) can be faster to determine than p != 2.
Ah, that’s what I was interested about. Thank you!
Is {lfloor log nrfloor}^2 a square mod n?
Arun, yes, it is always a square. Every integer square is also a square modulo n for any integer n.
To add one more function:
def squares_integer(n):
from sympy import factorint
return squares(factorint(n).items())
This is an important problem. There are highly composite moduli with much lower ratios of quadratic residues than your remarks would suggest. E.G..: 24 has 6 (= 1/4); 5040 has 192 = (1/26.25). Using Python 3.x as a test platform, I found empirically, working with 20 digit numbers (base 10), that sqrt(x) is about 40 times slower than x % m, so testing modulo 5040, then calculating sqrt(x) when necessary improves performance by abut one order of magnitude. There are some diminishing returns related to the magnitude of the test modulus, but I can envision using a modulus of order say 1E10 in an ordered array of residues and searching the list using a binary algorithm. What am I missing? Thanks
Solution of quadratic equations over rings is challenging.
Let’s take the ring of 2×2 matrices over R or C. It can be easily verified that the equation X^2 = [[0,1],[0,0]] has not solutions.