Number theory

An integer is either a perfect square or its square root is irrational. Said a different way, when you compute the square root of an integer, there are either no figures to the right of the decimal or there are an infinite number of figures to right of the decimal and they don’t repeat. There’s no middle ground. You can’t hope, for example, that the decimal expansion might stop or repeat after a hundred or so terms.

This theorem came up recently when I was talking to one of my kids about her math class, so I decided to look up the proof. It’s easier than I expected, not much harder than the familiar proof that the square root of 2 is irrational. Here goes.

Suppose a/b is a fraction in lowest terms, i.e. a and b are relatively prime, and a/b a solution to xⁿ = c where n > 0 is an integer and c is an integer. Then

(a/b)ⁿ = aⁿ / bⁿ = c

and so

aⁿ / b = c b^n-1.

Now the right side of the equation above is an integer, so the left side must be an integer as well. But b is relatively prime to a, and so b is relatively prime to aⁿ. The only way aⁿ / b could be an integer is for b to equal 1 or -1. And so a/b must be an integer.

This proof shows that what we said about square roots extends to cube roots and in fact to all integer roots. For example, the fifth root of an integer is either an integer or an irrational number.

Update: Note that there’s another proof in the comments, one that I believe is easier to follow.

Quadratic congruences are much more complex than linear congruences. If you’re interested in the details, see these notes on quadratic congruences.

When I started looking into this, I thought my number theory class years ago covered the details but that I’d forgotten them. Looking back, I see that my instructor and my textbook neglected some cases. Then I looked at a few other references and saw that they also left out some cases. I put together some notes by pulling together information from several incomplete but complementary references.

The Chinese Remainder Theorem (CRT) is a tool for solving problems involving modular arithmetic. The theorem is called the “Chinese” remainder theorem because the Chinese mathematician Sun Tsu stated a special case of the theorem sometime between 280 and 473 A.D. In its simplest form, the theorem says that if you want to solve a congruence (mod α β), you can solve it (mod α) and (mod β) separately, then combine the solutions to these sub-problems into a solution to the original problem, provided α and β are relatively prime.

More generally, the CRT says that the system of equations x ≡ a_i (mod m_i) has a common solution if the numbers m₁, m₂, …, m_r are relatively prime in pairs (i.e. m_i and m_j have no factors in common if i ≠ j). The solution is unique mod m where m is the product of the m_i’s. In a typical application, you may start with a problem stated (mod m) and then factor m into prime powers. Each of these prime powers is then an m_i in the CRT. For example, to solve a congruence (mod 1400), you might set m₁ = 8, m₂ = 25, m₃ = 7.

You can find an existence proof of the Chinese Remainder Theorem in any number theory book. Here I present an explicit solution given in Donald Knuth’s book Seminumerical Algorithms. For other algorithms and much more information on the CRT, see the book Chinese Remainder Theorem: Applications in Computing, Coding, and Cryptography. (In case it seems “coding” and “cryptography” are redundant, the book uses “coding” to refer to error-correcting codes, such as the way music is encoded on a CD, not secret codes, which would fall under cryptography.)

First we find numbers y_i such that y_i ≡ 1 (mod m_i) and y_i ≡ 0 (mod m_j) for i ≠ j. Then the sum a₁ y₁ + a₂ y₂ + … + a_r y_r is the solution we’re looking for.

But how to find the y_i? We can let y_i = (m/m_i)^φ(m_i).

But what is φ and why does this work? The function φ is called the Euler phi function. For a positive integer n, φ(n) is the number of positive integers less than n and relatively prime to n. (I’ll explain how to compute it in a minute.) Euler proved that if two numbers a and n are relatively prime, then a^φ(n) ≡ 1 (mod n). Because m/m_i is relatively prime to m_i, y_i ≡ 1 (mod m_i). For j ≠ i, m/m_i is a multiple of m_j and so y_i ≡ 0 (mod m_j).

So how do you compute φ(n)? This is easy once you know two facts about this function. First, if α and β are relatively prime, then φ(α β) = φ(α) φ(β). Second, for a prime p, φ(p) = p-1 and φ(p^k) = p^k – p^k-1 for k > 1. Putting these facts together, you can easily calculate φ(n) once you have its factorization into prime powers.

The solution given here is not be the fastest way to compute a solution to the problem posed by the CRT, but it is interesting because it’s explicit. On the other hand, it’s not terribly inefficient. The term (m/m_i)^φ(m_i) can be calculated using no more than 2 log₂φ(m_i) multiplications using fast exponentiation.

Related posts:

How to solve linear conguences
How to solve quadratic congruences

How do you efficiently compute aⁿ for integer n? You could multiply a*a*a*…*a, n-1 multiplications, but there are much more efficient ways. For example, a⁸ could be computed by ((a²)²)² in three multiplications. This example suggests it may be worthwhile to think about the binary representation of the exponent.

Here’s an algorithm. Write the exponent n in binary. Read the binary representation from left to right, starting with the second bit from the left. Start with the number a, and every time you read a 0 bit, square what you’ve got. Every time you read a 1 bit, square what you’ve got and multiply by a. It follows that aⁿ can be computed using no more than 2 log₂(n) multiplications.

As an example, suppose we want to calculate 7¹³. In binary, 13 is written as 1101. We skip over the leading 1 then read 1, 0, and 1. So we calculate ((((7²)7)²)²)7. In C notation, we would compute as follows.

x = 7;
x *= x;
x *= 7;
x *= x;
x *= x;
x *= 7;

In number theory calculations, such as arise in cryptography, it’s often necessary to compute aⁿ (mod m) for very large integers a, n, and m. These numbers may be hundreds or thousands of digits long and so of course the calculations cannot be carried out using ordinary machine integers. You could compute aⁿ first then reduce the result (mod m), but it would be more efficient to reduce (mod m) after every multiplication. This puts an upper bound on the size of numbers (and hence the amount of memory) needed for the calculation.

The same algorithm works when a is not an integer but a matrix. For example, consider a large matrix M that gives the connections between notes in a graph. M contains a 1 in position (i, j) if the i^th and j^th nodes are connected by an edge and contains a 0 otherwise. The n^th power of M tells which nodes are connected by paths of length less than or equal to n. If M is large, we don’t want to do unnecessary multiplications in computing Mⁿ.

For some exponents, the method given above can be improved. The smallest exponent for which you can improve on the binary algorithm is 15. Using the binary algorithm, computing a¹⁵ requires six multiplications. But you could compute the same quantity in five multiplications as follows.

b = a*a*a;
c = b*b;
c *= c;
c *= b;

For more examples, see Seminumerical Algorithms section 4.6.3.

A linear congruence is the problem of finding an integer x satisfying

ax ≡ b (mod m)

for specified integers a, b, and m. This problem could be restated as finding x such that

1. the remainder when ax is divided by m is b,
2. ax % b == m in C code,
3. ax – b is divisible by m, or
4. in base m, ax ends in b.

Two solutions are considered the same if they differ by a multiple of m. (It’s easy to see that x is a solution if and only if x + km is a solution for all integers k.)

For example, we may want to solve 7x ≡ 3 (mod 10). It turns out x = 9 will do, and in fact that is the only solution. However, linear congruences don’t always have a unique solution. For example, 8x ≡ 3 (mod 10) has no solution; 8x is always an even integer and so it can never end in 3 in base 10. For another example, 8x ≡ 2 (mod 10) has two solutions, x = 4 and x = 9.

In this post, we answer two questions.

1. How many solutions does ax ≡ b (mod m) have?
2. How can we compute them?

The answer to the first question depends on the greatest common divisor of a and m. Let g = gcd(a, m). If b is not divisible by g, there are no solutions. If b is divisible by g, there are g solutions.

So if g does divide b and there are solutions, how do we find them? The brute force solution would be to try each of the numbers 0, 1, 2, … m-1 and keep track of the ones that work. That works in theory, but it is impractical for large m. Cryptography applications, for example, require solving congruences where m is extremely large and brute force solutions are impossible.

First, suppose a and m are relatively prime. That is, assume g = gcd(a, m) = 1. We first put the congruence ax ≡ b (mod m) in a standard form. We assume a > 0. If not, replace ax ≡ b (mod m) with -ax ≡ -b (mod m). Also, we assume a < m. If not, subtract multiples of m from a until a < m.

Now solve my ≡ -b (mod a). This is progress because this new problem is solving a congruence with a smaller modulus since a < m. If y solves this new congruence, then x = (my + b)/a solves the original congruence. We can repeat this process recursively until we get to a congruence that is trivial to solve. The algorithm can be formalized into a procedure suitable for programming. The result is closely related to the Euclidean algorithm.

Now what if the numbers a and m are not relatively prime? Then first solve the congruence (a/g)y ≡ (b/g) (mod (m/g)) using the algorithm above. Then the solutions to ax ≡ b (mod m) are x = y + tm/g where t = 0, 1, 2, …, g-1.

Here are a couple examples.

First, let’s solve 7x ≡ 13 (mod 100). Since 7 and 100 are relatively prime, there is a unique solution. The algorithm says we should solve 100y ≡ -13(mod 7). Since 100 ≡ 2 (mod 7) and -13 ≡ 1 (mod 7), this problem reduces to solving 2y ≡ 1 (mod 7), which is small enough to solve by simply sticking in numbers. We find y = 4. Then x = (100*4 + 13)/7 = 59. You can verify that 7*59 = 413 so 7*59 ≡ 13 (mod 100). In general, we may have to apply the algorithm multiple times until we get down to a problem small enough to solve easily.

Now let’s find all solutions to 50x ≡ 65 (mod 105). Since gcd(50, 105) = 5 and 65 is divisible by 5, there are 5 solutions. So we first solve 10x ≡ 13 (mod 21). The algorithm above says we can solve this by first solving 21y ≡ -13 (mod 10), which reduces immediately to y ≡ 7 (mod 10), and so we take y = 7. This says we can take x = (105*7 + 65)/50 = 16. The complete set of solutions to our original congruence can be found by adding multiples of 105/5 = 21. So the solutions are 16, 37, 58, 79, and 100.

I intend to write posts in the future about how to solve simultaneous systems of linear congruences and how to solve quadratic congruences.

Update: Here are the posts I intended to write: systems of congruences, quadratic congruences.

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(n) 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( p ):
			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 hexidecimal 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.

* W. D. Stangl, “Counting Squares in Z_n”, Mathematics Magazine, pp. 285-289, Vol. 69 No. 4 October 1996.