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² ≡ 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ⁿ) 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³ 5³. 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ⁿ for large n. What about mod pⁿ 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.