This week we will learn about number theory that can be useful in contest problems. Contest problems often revolve around prime numbers and counting, in which number theory is a useful tool.
An integer \(x\) is a divisor of integer \(y\) if there exists integer \(k\) such that \( kx = y \). The greatest common divisor (GCD) of two integers is the largest integer that is a divisor of both of them. Definitions for zero can vary. We can compute GCD's using Euclid's algorithm (or the Euclidean algorithm), which is based on the following observation: Given two integers \(a \geq b\), it is true that $$ \gcd(a,b) = \gcd(a-b,b). $$ In other words, we can replace the larger number with its difference with the smaller number. By extension, by repeating this process, we can replace one of the numbers with its remainder when divided by the other. This therefore gives the following simple, elegent recursive algorithm for computing GCDs.
// Computes the greatest common divisor of a and b
function gcd(a, b) {
if (b == 0) return a
else return gcd(b, a % b)
}
GCDs are associative, so we can compute the GCD of more than two variables using the fact that $$ \gcd(a,b,c) = \gcd(a,\gcd(b,c)) = \gcd(\gcd(a,b),c) $$
The least common multiple of two integers \(x,y\) is the smallest integer \(m\) such that \( m = k_1 x = k_2 y \) for some integers \(k_1,k_2\). In other words, it is the smallest number that is both a multiple of \(x\) and \(y\). LCMs are easy to compute by using the observation $$ \text{lcm}(a,b) = \frac{|a\cdot b|}{\gcd(a,b)}. $$ This works for all inputs except when \(a = b = 0\), which would cause division by zero. In this case, we can define \( \text{lcm}(0,0) = 0 \) as a special case.
Similar to GCD, LCMs are also associative, so we can write $$ \text{lcm}(a,b,c) = \text{lcm}(a,\text{lcm}(b,c)) = \text{lcm}(\text{lcm}(a,b),c) $$
Testing whether a number is prime is a useful tool. Recall that a number is prime if it has two distinct divisors, 1 and itself. This means that 1 is not a prime number! Let's look at two common algorithms for primality checking.
The most naive algorithm for checking whethe an integer \(n\) is prime would be to check all possible divisors from \(1\) to \(n\). This takes \( O(n) \) time, though, which is very slow. A simple way to speed this up and obtain a good algorithm is to use the fact that if \( d \) is a divisor of \( n \), then \( \frac{n}{d} \) is also a divisor. This means that if \( n \) has any non-trivial divisor (a divisor other than 1 or itself) then it has one that is at most \( \sqrt{n} \). Therefore our algorithm only has to check up to \( \sqrt{n} \) when searching for divisors.
// Checks whether a number is prime in O(sqrt(n)) time
function is_prime(n) {
if (n == 1) return false;
for (i = 2; i * i <= n; i++) {
if (n % i == 0) return false;
}
return true;
}
Note the handy use of writing \( i \times i \leq n \) rather than \( i \leq \sqrt{n} \) since this avoids potentially hazardous and more expensive floating-point computations. Also, remember that 1 is not a prime number! I'll say this multiple times because this is a very common thing to forget and causes lots of bugs!
If the \( O(\sqrt{n}) \) algorithm is too slow, a faster algorithm, called the Miller-Rabin algorithm can be used. Deriving the Miller-Rabin algorithm is somewhat complicated, so we won't cover the full details, but let's get some intuition. Fermat's little theorem (not to be confused with Fermat's last theorem) says that for any integer \(a\) and prime number \(p\), $$ a^p \equiv a \mod p. $$ Therefore, if \( a^n \not\equiv a \mod p \), we immediately know that \( n \) is not prime. Of course, it is possible that \( a^n \equiv a \mod p \) just by chance, even if \( p \) is not prime. To check whether a number is "probably prime", we could therefore check this equation a bunch of times for many random \( a \) values. If it ever fails, we know for sure that \( n \) is not prime, but we might get false positives. This algorithm is called the Fermat primaility test.
The Miller-Rabin algorithm is essentially a more clever version of this that uses a few extra number theory tricks. You can find a more detailed description and some example implementations here.
Sometimes, it is useful to have a list of prime numbers handy. A trivial way to generate all of the primes up to \(n \) would of course simply be to loop over all of the numbers up to \(n\) and run a primality test on them. This is not super fast, though, and there are much better algorithms for this purpose. The most common algorithm used for this purpose is the Sieve of Eratosthenes. The algorithm is remarkably simply. Start with a list containing all of the integers from \( 0 \) to \( n \). Initially, 0 and 1 are marked as not prime, while all other numbers are marked as prime. We now loop over each number from 2 to \(n\) in order. If a number is prime, we know that all multiples of it are not prime, so when we find a prime number, we can loop over all of its multiples and mark them as not prime. When we are done, all of the composite (not prime) numbers will be marked as such, and all of the primes will be untouched.
A few optimizations can be made to speed this up. Instead of looping all the way to \(n\), notice that just like in the trial division algorithm, we can stop at \(\sqrt{n} \) since we will have hit at least one of the factors of each number by then. For each factor \( i \), we can also start striking off primes at \( i \times i \), instead of \( 2i \), since all smaller multiples will have been hit already.
Here is an example implementation that uses these ideas.
// Generates a list is_prime such that is_prime[n] = true if and only if n is prime
function sieve(n) {
is_prime = [true] * (n+1)
is_prime[0] = is_prime[1] = false
for (i = 2; i*i <= n; i++) {
if (is_prime[i]) {
for (j = i*i; j <= n; j += i)
is_prime[j] = false
}
}
return is_prime
}
The complexity of this algorithm is \( O(n \log\log(n)) \).
Facorizing a number means to compute a list containing all of its prime factors. Let's look at two algorithms for doing this, both of which are based on other algorithms we've talked about already.
The simplest method for factoring is based on the same concept of trial division as used in our primality testing algorithm. To gather the factors of a number, it suffices to simply loop over all possible factors and check whether they divide the number. As was the case before, we only have to loop up to \( \sqrt{n} \). The reason that this works is because a number can only have one factor that is larger than \( \sqrt{n} \), so we can simply factor out every other factor, and it will be what is left over at the end. If the number has no factor larger than \( \sqrt{n} \), then we will be left with 1.
// Returns the list of prime factors of n. This includes duplicate
// entries if a factor has multiplicity greater than one
function factor(n) {
factors = []
for (i = 2; i*i <= n; i++) {
while (n % i == 0) {
factors.push(i)
n /= i
}
}
if (n > 1) factors.push(n) // checks whether there is a factor > sqrt(n) left over
return factors
}
Trial division is fine for factoring a small quantity of not too large numbers, but what if we want to factor a lot of huge numbers! We will need something faster. One way to achieve this by using a modified version of the sieve algorithm.
When we perform the sieve, in addition to just figuring out which numbers are prime, notice that we can actually gather more information than that. For every composite number, the sieve also identifies a prime factor of that number (the number that was responsible for marking it as not prime). Therefore by running the sieve, we can precompute one factor of every number up to \( n \). With this information, it is easy to factor any number by simply looking up the stored factor of that number, and then dividing that factor out and repeating until we hit one.
// Precomputes a factor for every number up to n
// Assumes that "fact" is a global variable
function sieve(n) {
fact = [0,1,2,...,n]
for (i = 2; i*i <= n; i++) {
if (fact[i] == i) {
for (j = i*i; j <= n; j += i) {
fact[j] = i
}
}
}
}
// Returns the list of prime factors of n. This includes duplicate
// entries if a factor has multiplicity greater than one
function factor(n) {
factors = []
while (n > 1) {
f = fact[n]
factors.push(f)
n /= f
}
return factors
}
We just call sieve once at the beginning of our program, and use factor as many times as we'd like. Since a number \(n \) has at most \( O(\log(n)) \) factors, this factoring algorithm takes just \(O( \log(n) \) time.
There are also other more efficient variants of the sieve algorithm, some which run in \( O(n) \) time, such as this one.