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 *x*^{2} ≡ *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*(*p ^{n}*) 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 2^{3} 5^{3}. 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 2* ^{n}* for large

*n*. What about mod

*p*for odd prime

^{n}*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.

Hi John,

I think there’s a typo in the code: the second “isOdd( p )” in the function “s” should have been an “isOdd( n )” if I’m right.

-Felix

Felix: Thanks. I fixed the bug you found.

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