Counting Coprime Pairs Using Inclusion Exclusion

The principle of inclusion exclusion is a very powerful idea that can be used to solve alot of problems that arise in combinatorics. In this blog post, I will show how it can be used to solve one problem in particular, and I will also provide a very elegant recursive implementation of it. Let's say we want to count the number of positive integers \( k \) that are less than or equal to \( m \) and relatively prime to \( n \). Using set theoretic notation, week seek to find the cardinality of the set \( \{ p \mid p \in \{ 1, 2, \dots, m \}, \gcd(n,p) == 1 \} \) where \( gcd \) is the greatest common divisor function. Let's denote the function that counts the cardinality of this set as \( f(n,m) \). Now that the problem has been defined, I will show how to solve it using brute force methods. We can map the mathematical definition of the problem to a haskell program very naturally using a list comprehension.

f :: Int -> Int -> Int
f n m = length [p | p <- [1..m], gcd n p == 1]

While this definition is certainly nice in its simplicity and readability, there are more efficient approaches. This approach is linear time in \( m \), as it requires \( m \) calls to the gcd function. I will now start talking about other more efficient methods for evaluating \( f \). Note that in the special case of \( n = m \), we can evaluate \( f \) by using Euler's totient function, which states $$ f(n,n) = \phi(n) = n \cdot (1 - \frac{1}{p_1}) \cdot \dots \cdot (1 - \frac{1}{p_k}) $$ where \( p_i \) are the prime divisors of \( n \). Intuitively, this result holds because \( 1 - \frac{1}{p_i} \) represents the fraction of integers that are not divisible by \( p_i \). Taking the product over all \( p_i \) gives the fraction of integers that are not divisble by any \( p_i \). Assuming you accept that justification, you might raise the question: "Why doesn't the formula work when \( n \neq m \), then?". Here is one way to show why the definition doesn't make sense for \( m \neq n \). Let's assume for a moment that we can extend the defintion in the following way: $$ f(n,m) = m \cdot (1 - \frac{1}{p_1}) \cdot \dots \cdot (1 - \frac{1}{p_k}) $$ where again \( p_i \) are all the prime divisors of \( n \). Writing this out in an equivalent form, we can see that $$ f(n,m) = m \cdot \frac{p_1-1}{p_1} \cdot \dots \cdot \frac{p_k-1}{p_k} $$ Unless \( n \mid m \), this runs the risk of not being an integer. However, by definition of \( f \), it is clear that \( f(n,m) \) is always an integer. This definition is suitable whenever \( n \mid m \) however. Here is a more general way to compute \( f(n,m) \). First let \( n = p_1^{a_1} \dots p_k^{a_k} \) as we did for computing the euler totient function. Then note that the numbers which are relatively prime to \( n \) and less than or equal to \( m \) must not be divisible by any of \( p_1, p_2, \dots, p_k \). Then the principle of inclusion exclusion gives us a nice way of counting how many of these things there are. First note that it is trivial to figure out how many number less than \( m \) are divisible by \( r \) -- it's just \( \lfloor \frac{m}{r} \rfloor \) since the numbers are of the form \( r, 2r, \dots, t r \) where \( t =\lfloor \frac{m}{r} \rfloor \). From here on out, I will assume that the floor is implicit, and so \( \frac{m}{r} \) is taken to be integer division/division with floor. The principle of exclusion principle tells us that we can determine the number of integers \( \leq m \) not divisble by any \( p_i \) by computing the following sum: $$ f(n,m) = m - \sum_{i=1}^k \frac{m}{p_i} + \sum_{i=1}^k \sum_{j=i+1}^k \frac{m}{p_i \cdot p_j} - \sum_{i=1}^k \sum_{j=i+1}^k \sum_{l=j+1}^k \frac{m}{p_i \cdot p_j \cdot p_l} + \dots + (-1)^k \frac{m}{\prod_{i=1}^k p_i} $$ While this sum looks pretty nasty, there is a very nice way to express it recursively. The function that I will define to compute this, which I will call \( g \), consumes a set of prime factors of \( n \) and \( m \) as defined in the problem. So that $$ f(n,m) = g(primeFactors(n), m) $$ Here is the general definition of \( g(factors, m) \) $$ \begin{align*} g(\emptyset, m) &= m \\ g(P, m) &= g(P \setminus p, m) - g(P \setminus p, \frac{m}{p}) \quad where \: p \in P \\ \end{align*} $$ First observe that \( \frac{\frac{m}{p_1}}{p_2} = \frac{m}{p_1 p_2} \), which is obvious when working in the real numbers, and less obvious but equally true when we are working with integer operations. Now let's work out a small example. Suppose we want to compute the number of integers less than or equal to 1000 that are relatively prime to 15. That is, we want to compute \( g(\{ 3, 5 \}, 1000) \), which expands to $$ \begin{align*} g(\{ 3, 5 \}, 1000) &= g(\{ 5 \}, 1000) - g(\{ 5 \}, \frac{1000}{3}) \\ g(\{ 3, 5 \}, 1000) &= [g(\emptyset, 1000) - g(\emptyset, \frac{1000}{5})] - [g(\emptyset, \frac{1000}{3}) - g(\emptyset, \frac{1000}{3 \cdot 5})] \\ g(\{ 3, 5 \}, 1000) &= 1000 - \frac{1000}{3} - \frac{1000}{5} + \frac{1000}{3 \cdot 5} \end{align*} $$ This expansion is exactly equivalent to the formula from the inclusion exclusion method, although it's much easier to implement. The answer is \( f(15, 1000) = 533 \), which matches the brute force program defined earlier. This recursive function works because once the set of prime factors have been exhausted, every combination of primes have been produced, and the associated term in the inclusion exclusion formula has the correct sign for the contribution. By that I mean all terms that have an even number of primes have a coefficient of \( 1 \), and all terms that have an odd number of terms have a \( -1 \) coefficient. Now that I've demonstrated how the method works, I will provide a Haskell implementation. If you want to run this code yourself, you will have to get the cabal package called arithmoi by running the command: cabal install arithmoi. This package provides a bunch of cool number theory functions. I will be using their function factorise :: Integer -> [(Integer, Int)] which produces the list of prime factors and associated multiplicities of a given number.

import Math.NumberTheory.Primes

factors :: Integer -> [Integer]
factors = map fst . factorise

g :: [Integer] -> Integer -> Integer
g [] m = m
g (p:ps) m = g ps m - g ps (m `div` p)

f :: Integer -> Integer -> Integer
f n m = g (factors n) m

In this blog post, I have shared a nice method for counting number of relatively prime pairs using recursion and prime factorisations. However, the method I described can be generalized to count other types things that are related to primes. I'll leave it to you to figure out how to use this method to implement a prime counting function that can count the number of primes less than \( n \) by only generating primes up to \( \sqrt{n} \).