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.