Number Theory
An brief introduction to modular arithmetic, prime numbers, and factoring in computer science

Modular Arithmetic
Let’s say that you tried to compute the 100,000th, or 10⁵, Fibonacci number. Soon you realize that the number is completely off, however.
This is because the number in question has 63 digits, a lot more than even the widest bit integer in C++ (128) can support. Therefore, many competitive programming problems will instead ask you for the result mod some number, usually either 10⁹ + 7 or 998353244. A number, such as x, mod another number, such as y, simply means the remainder of x divided by y. For example, 5 mod 2 would be 1. Most programming languages use the % operator to denote mod.
We use modular arithmetic to deal with such numbers and find the result in mod efficiently. Instead of doing calculations with extremely large integers directly, we do calculations with their remainder, or mod some number, according to the following rules:
(a + b) mod c = ((a mod c) + (b mod c)) mod c
(a - b) mod c = ((a mod c) - (b mod c)) mod c
(a * b) mod c = ((a mod c) * (b mod c)) mod c
(Note that the pattern does not apply to division)
With negative numbers, the modulo operation may give a negative output, resulting in inconsistent behavior in the code. To combat this, we can do ((a mod c) + c) mod c to ensure the result is positive.
Modular Inverse
A modular inverse of a number can be thought of as its reciprocal. More formally, define a modular inverse of x as x⁻¹ such that x * x⁻¹ 1 mod c for some divisor c. If c is prime (which it usually is), then it can be proven that x⁻¹ always exists for any x between 1 and c - 1, inclusive. Finding the modular inverse of a number can be done by finding xc - 2 and using binary exponentiation in the process. This can only be done if c is prime, but the vast majority of the time it will be. This process is also encapsulated in the templates following with the / operator.
As the operator suggests, modular inverses are useful for division calculations. For example, if we wanted to compute 3 divided by 5 modulo 11, we’d calculate 3 * (modular inverse of 5), which would be 3 * 9 = 27 ≡ 5 mod 11.
Modnum Templates
Modnum templates, such as Ecnerwala’s, can make your code much more readable and less prone to bugs compared to coding add, multiply functions etc. yourself. If that’s not an option, however, for contests such as USACO which don’t allow templates, we recommend trying to memorize or understand part of Benq’s short template.
Prime Numbers
Prime numbers are the building blocks of all integers. The fundamental theorem of arithmetic tells us that every number has a unique prime factorization.
Check if one number is prime
To check that n is prime, we iterate from 2 to √(n), checking if they divide n. If none of them divide n, n is prime, otherwise, it’s composite. This takes O(√(n)) time.
Code:
def is_prime(x):
    if x <= 1:
        return False
    i = 2
    while i * i <= x:
        if x % i == 0:
        return False 
        i+=1
    return True
                              
Code:
static boolean isPrime(int x) {
  if (x <= 1) {
    return false;
  }
  for (int i = 2; i * i <= x; i++) {
    if (x % i == 0) {
      return false;
    }
  }
  return true;
}
Code:
bool is_prime(int x) {
  if (x <= 1) {
    return false;
  }
  for (int i = 2; i * i <= x; i++) {
    if (x % i == 0) {
      return false;
    }
  }
  return true;
}
Finding primes from 1 through n
Say we are trying to find all prime numbers from 1 to n. We keep an boolean array prime[n]. Let prime[i]=1 indicate that i is a prime and prime[i]=0 indicate that i is composite. We iterate from 2 to n. If our current number i is marked as prime, we go through the array and mark all the other multiples of i as composite. This way, after iterating through all the numbers, all the primes will have value 1 and all composites will have value 0 (we ignore 1, since it is neither prime nor composite). The algorithm takes O(n log n) time, and is known as the Sieve of Eratosthenes.
Code:
prime = [True for i in range(N + 1)]
prime[0] = False
prime[1] = False
for i in range(2, N + 1):
  if not prime[i]:
    continue 
  for j in range(i * 2, N + 1, i): 
    prime[j] = False
Code:
boolean[] prime = new boolean[N + 1];
prime[0] = false;
prime[1] = false;
for (int i = 2; i <= N; i++) {
  if (!prime[i]) {
    continue;
  }
  for (int j = i * 2; j <= N; j++) {
    prime[j] = false;
  }
}
Code:
vector<bool> prime(N + 1);
prime[0] = false;
prime[1] = false;
for (int i = 2; i <= N; i++) {
  if (!prime[i]) {
    continue;
  }
  for (int j = i * 2; j <= N; j++) {
    prime[j] = false;
  }
}
Factors

Quick Prime Factorization
Suppose we want to compute the prime factorization for every number between 2…N. Usually this would take around O(√(N)) per query, but we can speed this up immensely by precomputing the highest factor of every number to a runtime of O(log N) per query with O(N log N) preprocessing.
Code:
hfac = [1 for i in range(N + 1)]
for i in range(2, N + 1):
    if hfac[i] == 1:
        for j in range(i, N + 1, i):
        hfac[j] = i

# 60 -> [2, 2, 3, 5]
def factors(x):
    ret = []
    while x != 1:
        ret.append(hfac[x])
        x = int(x / hfac[x])
    ret.reverse()
    return ret      
Code:
int[] hfac = new int[N + 1];
Arrays.fill(hfac, 1);
for (int i = 2; i <= N; i++) {
    if (hfac[i] == 1) {
        for (int j = i; j <= N; j += i) {
            hfac[j] = i;
        }
    }
}

// 60 -> [2, 2, 3, 5]
static List factors(int x) {  
    List ret = new ArrayList<>();
    while (x != 1) {
        ret.add(hfac[x]);
        x /= hfac[x];
    } 
    Collections.reverse(ret);
    return ret;
}
                                
Code:
vector hfac(N + 1, 1);
for (int i = 2; i <= N; i++) {
    if (hfac[i] == 1) {
        for (int j = i; j <= N; j += i) {
            hfac[j] = i;
        }
    }
}

// 60 -> [2, 2, 3, 5]
vector factors(int x) {  
    vector ret; 
    while (x != 1) {
        ret.push_back(hfac[x]);
        x /= hfac[x];
    } 
    reverse(ret.begin(), ret.end());
    return ret;
}                                
You can play around with all the code we've used in this article on Replit:
Copyright ©2023 Howard County Hour of Code