Introduction

Here are my notes for techniques for solving Project Euler problems. It is a continuous work in progress. I tend to update this when there is a technique that applies to more than one problem. I keep a list of candidate topics that I might write about in the future here:

• Seiving
• Modular arithmetic
• Dynamic Programming
• Lucy’s Algorithm
• Diophantine equations
• closed form solutions
• summing multiplicative functions
• groups, adjuncts, etc

GCD and extended GCD algorithms

Greatest Common Divisor (gcd)

Say we have two integers \(a\) and \(b\), the gcd is the biggest number that cleanly divides both \(a\) and \(b\) without a remainder. That is, the gcd is the biggest number \(m\) such that \(a = mx\) and \(b = my\) for some \(x\) and \(y\). There is a simple algorithm to find this number:

Subtract the smaller number from the bigger number until both numbers are equal. When they are equal, we have the gcd

For example, if we want to find the gcd of 64 and 26:

So the gcd of 64 and 26 is 2. This is called the Euclidean algorithm. More information, and more efficient algorithms can be found on Wikipedia

When \(gcd(a, b) = 1\), then \(a\) and \(b\) are said to be coprime. that is, they have no factors in common.

Optimising the GCD

If one number is much larger the other, then we spend a lot of time repeatedly subtracting it. Repeated subtraction is just division. What we are left with is the remainder. So instead of repeatedly subtracting one number from the other, we can take the remainder of one number from the other and save a lot of time:

pub fn gcd(mut a: u32, mut b: u32) -> u32 {
    while a != 0 {
        (a, b) = (b % a, a);
    }

    b
}

fn main() {
    println!("gcd({}, {}) = {}", 64, 26, gcd(64, 26));
}

Lowest Common Multiple (lcm)

The lowest common multiple is the smallest number that is a multiple of both \(a\) and \(b\). The lcm is equal to:

\[\frac{ab}{gcd(a, b)}\]

Say we want to find the smallest solution to the equation \(ax + by = 0\). Then we have:

\[ \begin{align*} a\frac{b}{gcd(a, b)} - b\frac{a}{gcd(a, b)} &= 0 \\ (x, y) = \left(\frac{b}{gcd(a, b)}, -\frac{a}{gcd(a, b)}\right) \end{align*} \]

Bézout’s Identity

Let \(a\) and \(b\) be any two integers, then an important middle step for many problems is to find an \(x\) and \(y\) such that \(ax + by = gcd(a, b)\) This is called Bézout’s identity. We start with two equations:

\[ \begin{align*} 1a + 0b &= a \\ 0a + 1b &= b \end{align*} \]

At each step, just like the gcd, we subtract the larger equation from the smaller one until both are the same. Since we start with \(a\) and \(b\) just like the Euclidean algorithm, we will end up with the gcd at the end. Using 64 and 26 the same example as above:

So now we have two solutions to our equation: \(x = 11, y = -27\) and \(x = -2, y = 5\). This process is called the extended Euclidean algorithm. Using the gcd optimisation above, we get the following algorithm:

// Solve ax + by = gcd(a, b). Returns (x, y, gcd(a, b))
fn extended_gcd(mut a: i32, mut b: i32) -> (i32, i32, i32) {
    let (mut x_1, mut y_1, mut x_2, mut y_2) = (1, 0, 0, 1);

    while a != 0 {
        let quotient = b / a;
        // subtract the first equation from the second and swap
        (x_1, x_2) = (x_2 - x_1 * quotient, x_1);
        (y_1, y_2) = (y_2 - y_1 * quotient, y_1);
        (a, b) = (b % a, a);
    }

    (x_2, y_2, b)
}

fn main() {
    let (x, y, gcd) = extended_gcd(64, 26);
    println!("Solution: {} * 64 + {} * 26 = {}", x, y, gcd);
}

We note that in the first column of the table above, the value of \(x\) always increases, whilst the value of \(y\) always decreases. The opposite is true for the second column. These solutions are somewhat special, the left column contains the smallest positive value of \(a\) that satisfies the identity, whilst the right hand column contains the smallest positive value of \(b\) that satisfies the identity.

There are infinitely many other solutions that we can make by adding or subtracting zero to the above equation:

\[ \begin{align*} ax + by &= 0 \\ 13a - 32b &= 0 \tag*{from “lcm” section above} \\ 11a - 27b &= 2 \tag*{from table above} \\ (11a - 27b) + (13a - 32b) &= 2 \\ 24a - 59a &= 2 \end{align*} \]

Chinese remainder theorem (CRT)

The Chinese remainder theorem states that if we have the solution to an equation modulo \(m_1\) and modulo \(m_2\), then we can find the solution modulo \(m_1m_2\). In other words, if we have two equations of the form:

\[ \begin{align*} x &= a_1 \mod m_1 \\ x &= a_2 \mod m_2 \\ \end{align*} \]

then we can find \(x \mod m_1m_2\).

When to use this

Often we can find the answer to an equation modulo some small prime, but it is difficult to find the full solution. In this case, we can calculate the answer modulo a series of small primes, and then combine the answers using CRT.

Sometimes a question will directly require the use of CRT to find a solution.

Calculating a CRT with coprime \(m_1\) and \(m_2\):

To start, we want to find a pair of numbers \(n_1\), \(n_2\) such that:

\[ \begin{align*} n_1 &= 1 \mod m_1 \\ &= 0 \mod m_2 \\ n_2 &= 1 \mod m_2 \\ &= 0 \mod m_1 \end{align*} \]

Then \(a_1n_1 + a_2n_2\) will be a solution. To find \(n_1\) we can rearrange:

\[ \begin{align*} n_1 &= 0 \mod m_2 \\ &= m_2k_2 \\ &= 1 \mod m_1 \\ &= 1 - m_1k_1 \\ 1 &= m_2k_2 + m_1k_1 \\ \end{align*} \]

We can use Bézout’s identity to find a solution to \(1 = m_2k_2 + m_1k_1\) as long as \(m_1\) and \(m_2\) are coprime. In the case that they are not coprime, then we can make them coprime by dividing one of them by the gcd. In this case, solutions only exist if \(a_1 = a_2 \mod gcd(m_1,m_2)\)

// solve ax + by = gcd(a, b). return (x, y, gcd(a, b))
pub fn extended_gcd(a: i64, b: i64) -> (i64, i64, i64) {
    let mut r = (a, b);
    let mut s = (1, 0);

    while r.1 != 0 {
        let quotient = r.0 / r.1;
        r = (r.1, r.0 - quotient * r.1);
        s = (s.1, s.0 - quotient * s.1);
    }

    (s.0, (r.0 - s.0 * a) / b, r.0)
}

// Returns x such that x % m_1 = a_1 and x % m_2 = a_2. Returns None if there are no solutions
fn crt(a_1: i64, m_1: i64, mut a_2: i64, mut m_2: i64) -> Option<i64> {
    let (mut k_1, mut k_2, gcd) = extended_gcd(m_1, m_2);

    if gcd > 1 {
        // No solution
        if a_1 % gcd != a_2 % gcd {
            return None;
        }

        m_2 = m_2 / gcd;
        a_2 = a_2 % m_2;
        (k_1, k_2, _) = extended_gcd(m_1, m_2);
    }

    Some((m_2 * k_2 * a_1 + m_1 * k_1 * a_2).rem_euclid(m_1 * m_2))
}

fn main() {
    assert_eq!(crt(1, 2, 3, 5), Some(3));
    assert_eq!(crt(1, 10, 0, 2), None);
    assert_eq!(crt(6, 10, 1, 5), Some(6));

    let result = crt(20, 23, 7, 11).unwrap();
    assert_eq!(result % 11, 7);
    assert_eq!(result % 23, 20);
    assert!(result < 11 * 23);
}

Stern-Brocot Tree

The Stern-Brocot tree allows us to iterate over all fractions in a given range. It has a few main properties we care about:

• All fractions in the tree are fully reduced
• It will eventually iterate through all fractions
• Fractions with smaller denominators are higher in the tree
• Fractions in the tree are ordered from smallest to largest

We start at \(\left(\frac{0}{1}\frac{1}{0}\right)\). For any given node \(\left(\frac{a}{c}\frac{b}{d}\right)\), the left child is given by \(\left(\frac{a}{c}\frac{a + b}{c + d}\right)\) and the right child is given by \(\left(\frac{a + b}{c + d}\frac{b}{d}\right)\)

This can be implemented with relatively few lines of code:

fn stern_brocot(left: (u64, u64), right: (u64, u64)) {
    let mediant = (left.0 + right.0, left.1 + right.1);

    // terminate if the denominator is greater than 3
    if mediant.1 > 3 {
        return;
    }

    // terminate if the fraction increases past 2
    if mediant.0 > mediant.1 * 2 {
        return;
    }

    stern_brocot(left, mediant);
    println!("Visited fraction: {}/{}", mediant.0, mediant.1);
    stern_brocot(mediant, right);
}

fn main() {
    stern_brocot((0, 1), (1, 0));
}

Best approximation of anything

The Stern-Brocot tree can be used to find the best approximation of anything with a fraction. Say we want to find the best approximation for \(\sqrt{13}\) where the denominator is less than \(1000\). As we go deeper into the tree, the denominators increase, and the upper and lower bounds become tighter. Once the denominator exceeds 1000, then we stop, since going down the tree will only make the denominator larger than it already is. The tree is ordered left to right, so we only need to check the branch that contains the number we want to approximate.

fn distance((num, den): (u32, u32)) -> f64 {
    (num as f64 / den as f64 - 13.0_f64.sqrt()).abs()
}

fn best_approximation() -> (u32, u32) {
    let mut left = (0, 1);
    let mut right = (1, 0);

    loop {
        let mediant = (left.0 + right.0, left.1 + right.1);

        // Stop when the denominator is > 1000
        if mediant.1 > 1000 {
            if distance(left) < distance(right) {
                return left;
            } else {
                return right;
            }
        }

        // If sqrt(13) is less than the mediant, go left, otherwise go right
        if (mediant.0 * mediant.0) / (mediant.1 * mediant.1) >= 13 {
            right = mediant;
        } else {
            left = mediant;
        }
    }
}

fn main() {
    let fraction = best_approximation();
    println!("Best approximation for sqrt(13) is {}/{}", fraction.0, fraction.1);
}

Pell’s Equation

We can solve basically every quadratic diophantine equation of two variables by reducing it to the Pell equation

Solve Pell’s equation

Pell’s equation has the form \(x^2 - ny^2 = 1\) where \(n\) is a positive integer that is not a square¹. It has a trivial solution \(x = \pm1, y = 0\), but our task is to find non trivial ones. Rearranging gives the following:

\[ \begin{align*} \sqrt{n + \frac{1}{y^2}} &= \frac{x}{y} \end{align*}
\]

Looking at this equation, we can see that as \(y\) increases, \(1/y^2\) decreases, and \(x/y\) becomes a better and better approximation for \(\sqrt{n}\).

Probably the easiest way to solve this equation is to use the Stern Brocot tree², which iterates through all reduced fractions in a range. Any solution \(x/y\) must be between \(\sqrt{n}\) and \(\sqrt{n + \frac{1}{y^2}}\) inclusive. We simply ignore branches outside this bound. As we progress down the tree \(y\) increases so the upper bound becomes closer and closer to \(\sqrt{n}\). We will find every solution because the Stern-Brocot tree iterates through every reduced fraction.

use std::cmp::Ordering;

fn find_pell_solutions(n: i64, left: (i64, i64), right: (i64, i64)) {
    let mediant = (left.0 + right.0, left.1 + right.1);

    // terminate if the denominator is greater than 100
    if mediant.1 > 1000 {
        return;
    }

    if mediant.0.pow(2) == 1 + mediant.1.pow(2) * n {
        println!("Found solution! {}^2 - {}*{}^2 = 1", mediant.0, n, mediant.1);
    }

    // if x/y < sqrt(n), then the left branch of the tree is completely out of
    // bounds
    if mediant.0.pow(2) > mediant.1.pow(2) * n {
        find_pell_solutions(n, left, mediant);
    }

    // if x/y > sqrt(n + 1/y^2), then the right branch of the tree is
    // completely out of bounds
    if mediant.0.pow(2) < mediant.1.pow(2) * n + 1 {
        find_pell_solutions(n, mediant, right);
    }
}

fn main() {
    // Find solutions for n = 2
    find_pell_solutions(2, (0, 1), (1, 0));
}

A non recursive solution

Recursive functions can be refactored to use a loop if you manage the stack yourself:

fn find_pell_solutions(n: i64) {
    let mut stack: Vec<((i64, i64), (i64, i64))> = vec![((0, 1), (1, 0))];

    while let Some((left, right)) = stack.pop() {
        let mediant = (left.0 + right.0, left.1 + right.1);

        if mediant.1 > 1000 {
            continue;
        }

        if mediant.0.pow(2) == 1 + mediant.1.pow(2) * n {
            println!("Found solution! {}^2 - {}*{}^2 = 1", mediant.0, n, mediant.1);
        }

        if mediant.0.pow(2) > mediant.1.pow(2) * n {
            stack.push((left, mediant));
        }

        if mediant.0.pow(2) < mediant.1.pow(2) * n + 1 {
            stack.push((mediant, right));
        }
    }
}

fn main() {
    find_pell_solutions(2);
}

Negative pell equation

The negative pell equation has the form \(x^2 - ny^2 = -1\) we rearrange to get

\[ \begin{align*} \sqrt{n - \frac{1}{y^2}} &= \frac{x}{y} \end{align*}
\]

\(x/y\) is between \(\sqrt{n - \frac{1}{y^2}}\) and \(\sqrt{n}\), and we can use the Stern-Brocot tree again:

fn find_pell_solutions(n: i64, left: (i64, i64), right: (i64, i64)) {
    let mediant = (left.0 + right.0, left.1 + right.1);

    if mediant.1 > 1000 {
        return;
    }

    if mediant.0.pow(2) + 1 == mediant.1.pow(2) * n {
        println!("Found solution! {}^2 - {}*{}^2 = 1", mediant.0, n, mediant.1);
    }

    // if x/y < sqrt(n - 1/y^2), then the left branch of the tree is completely
    // out of bounds
    if mediant.0.pow(2) + 1 > mediant.1.pow(2) * n {
        find_pell_solutions(n, left, mediant);
    }

    // if x/y > sqrt(n), then the right branch of the tree is completely out of
    // bounds
    if mediant.0.pow(2) < mediant.1.pow(2) * n {
        find_pell_solutions(n, mediant, right);
    }
}

fn main() {
    find_pell_solutions(2, (0, 1), (1, 0));
}

Generalised pell equation

The generalised pell equation has the form \(x^2 - ny^2 = a\) we rearrange to get

\[ \begin{align*} \sqrt{n - \frac{a}{y^2}} &= \frac{x}{y} \end{align*}
\]

\(x/y\) is between \(\sqrt{n - \frac{a}{y^2}}\) and \(\sqrt{n}\), and we can solve it like the others.

An aside: solving \(x^2 - n^2y^2 = a\)

What happens when \(n\) is a square in the generalised Pell equation? Then we have \(x^2 - (ny)^2 = a\), which is a difference of squares. Factorising gives:

\[ a = (x - ny)(x + ny) \]

We can solve this equation by factorising \(a = pq\). Suppose \(p = x + ny\) and \(q = x - ny\), then we get the following:

\[ x = \frac{p + q}{2}, y = \frac{p - q}{2n} \]

A simple algorithm is to factorise \(a\) to find candidates for \(x\) and \(y\). Then filter out all solutions where \(x\) and \(y\) are fractions.

Solve \(ax^2 + bx + c = dy^2\)

We apply the quadratic equation directly which gives:

\[ x = \frac{-b \pm \sqrt{b^2 - 4a(c - dy^2)}}{2a} \]

We note that \(x\) is only rational if the part under the square root is a square. Expanding gives:

\[ \begin{align*} z^2 &= b^2 - 4a(c - dy^2) \\ z^2 -4ady^2 &= b^2 - 4ac \\ z^2 -Ny^2 &= A \end{align*} \]

Now we can solve it using any technique for the generalised Pell equation. If \(N\) is a square, we can solve using factorisation instead Once we have solutions, we can substitute back to get \(x\), filtering out any solutions where \(x\) is a fraction.

Generalised quadratic diophantine equations

Say we have a diophantine equation for which we want to find integer solutions for \(x, y\):

\[ax^2 + bxy + cy^2 + dx + ey + f = 0\]

We apply the quadratic equation directly to solve for \(y\):

\[ y = \frac{-(bx + e) \pm \sqrt{(bx + e)^2 - 4c(ax^2 + dx + f)}}{2c} \]

We note that \(y\) is only rational if the part under the square root is a square. Expanding gives:

\[ \begin{align*} z^2 &= (bx + e)^2 - 4c(ax^2 + dx + f) \\ &= (b^2 + 4ac)x^2 + (2be -4cd)x + (e^2 - 4cf) \\ &= Ax^2 + Bx + C \end{align*} \]

And now we can solve it using the technique above. Once we have candidate solutions, filter out any solutions where \(y\) is a fraction.

1. If \(n\) is a square, then we have \(x^2 - m^2y^2 = 1\). This is a difference of squares. The only solutions are the trivial solutions \(x = \pm1, y = 0\). ↩

2. Most literature online says that we can find a solution using continued fractions, however finding good guidance or working code online is surprisingly hard. Solving the generalized Pell equation by John P. Robertson is a good resource that is terse, but comprehensive. ↩

Lucy’s algorithm

Lucy’s algorithm can help us find sums of multiplicative functions. I will discuss several variants in this article

Not Quite Lucy Algorithm

Say we have the following problem:

• We have a function \(f(n)\)
• \(f(n)\) is completely multiplicative, i.e, \(f(1) = 1$ and $f(a)f(b) = f(ab)\)
• We have an easy way to calculate \(F(l) = \sum_{n=1}^{l}{f(n)}\)
• We want to find the sum of \(f(p)\) for every prime \(p\) less than a limit \(l\)

Then we can use Lucy’s algorithm to calculate that sum.

We denote \(lpf(n)\) as the lowest prime factor of \(n\). We define the function \(S(l, p)\) as the sum of all \(f(n)\) less than or equal to \(l\) where the lowest prime factor of \(n\) is greater than \(p\).¹

We note that \(S(l, \lfloor \sqrt{l} \rfloor)\) is the sum of \(f(n)\) for all primes between \(\sqrt{l}\) and \(l\). \(1\) is also included in this sum, so we will need to subtract that out near the end.

The set of \(n\) with \(lpf(n)\) greater than \(p\) is set of \(n\) with \(lpf(n)\) greater than \(p - 1\) minus the set of \(n\) with \(lpf(n)\) exactly equal to \(p\). We can write this as follows:

\[S(l, p) = S(l, p - 1) - f(p) \cdot S(\lfloor l/p \rfloor, p - 1)\]

This formula relies on \(f(n)\) being completely multiplicative. We also note that in order for the lowest prime factor to be exactly \(p\), \(p\) must be a prime. If \(p\) is composite, \(S(\lfloor l/p \rfloor, p - 1)\) is zero. This gives us a recursive way to find \(S(l, p)\). The base case is \(S(l, 1) = F(l)\).

This is just enough information to write code for the algorithm. In this example, \(f(n) = 1\), which will allow us to count the number of primes less than or equal to \(l\):

use std::collections::{HashMap, HashSet};

fn capital_s(
    cache: &mut HashMap<(u64, u64), u64>,
    primes: &HashSet<u64>,
    l: u64,
    p: u64)
-> u64 {
    if let Some(result) = cache.get(&(l, p)) {
        return *result;
    }

    if p <= 1 {
        // return sum of f(n) for all n
        l
    } else if p >= l {
        // smallest prime factor is greater than limit, return 0
        1
    } else if !primes.contains(&p) {
        // p is composite, keep going
        capital_s(cache, primes, l, p - 1)
    } else {
        // normal recursive case
        let result = capital_s(cache, primes, l, p - 1) - capital_s(cache, primes, l / p, p - 1);
        cache.insert((l, p), result);
        result
    }
}

fn main() {
    let l = 1000;
    // Find these primes using a sieve
    let mut small_primes = HashSet::from([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]);
    let mut cache = HashMap::new();
    let result = capital_s(&mut cache, &small_primes, l, 31) + small_primes.into_iter().count() as u64 - 1;

    println!("There are {} primes less than or equal to {}", result, l);
}

Multiplicative Lucy

\[\begin{align*} S(l, p) &= \sum_{\substack{n \leq l \\ lpf(n) > p}}{f(n)} \\ &= \sum_{\substack{n \leq l\\lpf(n) > p - 1}}{f(n)} - \sum_{\substack{n \leq l\\lpf(n) = p}}{f(n)}\\ &= S(l, p - 1) - \sum_{\substack{pn \leq l \\ lpf(n) > p - 1}}{f(pn)} \\ &= S(l, p - 1) - \left(\sum_{\substack{pn \leq l \\ lpf(n) > p - 1 \\ p \nmid n}}{f(pn)} + \sum_{\substack{pn \leq l \\ lpf(n) > p - 1 \\ p \mid n}}{f(pn)} \right)\\ &= S(l, p - 1) - \left(f(p)\sum_{\substack{n \leq \lfloor l/p \rfloor \\ lpf(n) > p}}{f(n)} + \sum_{\substack{p^2n \leq l \\ lpf(n) > p - 1}}{f(p^2n)} \right)\\ &= S(l, p - 1) - \left(f(p)S(\lfloor l / p \rfloor, p) + f(p^2)S(\lfloor l/p^2 \rfloor, p) + … \right)\\ &= S(l, p - 1) - \sum_{n \geq 1}{f(p^n)S(\lfloor l / p^n \rfloor, p)}\\ \end{align*} \]

We denote \(lpf(n)\) as the lowest prime factor of \(n\). We define the function \(S(l, p)\) as the sum of all \(f(n)\) less than or equal to \(l\) where the lowest prime factor of \(n\) is greater than \(p\).¹

\[S(l, p) = \sum_{\substack{n \leq l \\ lpf(n) > p}}{f(n)}\]

We note that for any composite less than \(l\) the lowest prime factor is always smaller than \(\sqrt{l}\). Therefore, any number between \(\sqrt{l}\) and \(l\) that has a lowest prime factor greater than \(\sqrt{l}\) must be prime. Another way of writing this is as follows:

\[S(l, \lfloor \sqrt{l}\rfloor) = \sum_{n\in \mathbb{P},\ \lfloor \sqrt{l}\rfloor<n\leq l} n\]

So we can sum all primes less than \(l\) by calculating \(S(l, \lfloor \sqrt{l}\rfloor)\), and then adding the primes below \(\sqrt{l}\). We can find these small primes with a simpler method such as the sieve of Eratosthenes.

We want to calculate \(S(l, p)\) Recursively. The base case is \(S(l, 1)\) which we can calculate easily:

\[ \begin{align*} S(l, 1) &= \sum_{n \leq l}{n} \\ &= \frac{n(n + 1)}{2} \end{align*} \]

\(S(l, p)\) has the following recurrence relation² which is all the information we need to get started:

\[ \begin{align*} S(l, p) &= S(l, p - 1) - p \cdot S(\lfloor l/p \rfloor, p - 1) \end{align*} \]

We also note that when \(p\) is composite then \(S(l, p) = S(l, p - 1)\), because it is impossible to have a lowest prime factor equal to a composite.

We can write an algorithm as follows:

• Calculate

1. This is not the usual definition of \(S\) that you find for lucy’s algorithm, but it simplifies the recursive formulas a lot, which is why I use it. ↩ ↩2

2. The recurrence relation is proved as follows: \[\begin{align*} S(l, p) &= \sum_{\substack{n \leq l\\lpf(n) > p - 1}}{n} - \sum_{\substack{n \leq l\\lpf(n) = p}}{n}\\ &= S(l, p - 1) - \sum_{\substack{pn \leq l \\ lpf(n) > p - 1}}{pn} \\ &= S(l, p - 1) - p\sum_{\substack{n \leq \lfloor l/p \rfloor \\ lpf(n) > p - 1}}{n} \\ &= S(l, p - 1) - p \cdot S(\lfloor l/p \rfloor, p - 1) \end{align*} \] ↩

Principle of Inclusion Exclusion

Mobius inversion

Say we have the following problem:

• We have two related functions \(f\) and \(g\) such that \[ \begin{align*} g(1) &= f(1) + f(2) + f(3) … \\ g(2) &= f(2) + f(4) + f(6) … \\ g(3) &= f(3) + f(6) + f(9) … \\ \vdots \\ g(n) &= \sum_{m = 1}^{\infty}f(mn) \end{align*}
  \]
• We want to calculate \(f(1)\)
• We have a good way to find \(g(n)\) but no easy way to find \(f(n)\)

Then we can use a technique called Mobius inversion.

We start with a simpler problem. Let’s say we want to calculate the sum of all \(f(n)\) when \(n\) is divisible by \(2\) or \(3\). We can calculate this using \(g(2) + g(3) - g(6)\). We subtract \(g(6)\) because we counted the multiples of \(6\) twice.

Say we add another prime and we want to find the sum of all \(f(n)\) when \(n\) is divisible by \(2, 3\) or \(5\). Then we can calculate this with \(g(2) + g(3) + g(5) - g(6) - g(10) - g(15) + g(30)\). We subtract \(g(6), g(10)\) and \(g(15)\) because we counted them twice. We add \(g(30)\) because we add it three times for \(g(2), g(3)\) and \(g(5)\) and then subtract it three times for \(g(6), g(10)\) and \(g(15)\). This leaves us counting it zero times, so we have to add it back.

Adding more primes we can see a pattern emerge:

• Add all \(g(n)\) when \(n\) is divisible by any prime
• Subtract all \(g(n)\) when \(n\) is divisible by at least two primes
• Add all \(g(n)\) when \(n\) is divisible by at least three primes
• Subtract all \(g(n)\) when \(n\) is divisible by at least four primes
• And so on

To simplify the process we can introduce the Mobius function \(\mu(n)\) which allows us to know when to add or subtract a value of \(g(n)\). It has the following definition:

\[ \mu(n) = \begin{cases} -1 & \text{if n is the product of an odd number of distinct primes} \\ 1 & \text{if n is the product of an even number of distinct primes} \\ 0 & \text{otherwise} \\ \end{cases} \]

Knowing this, we can calculate \(f(1)\) using \(g(1)\) and then subtracting all the multiples of \(3, 5, 7, 11 …\) until we are left with just \(f(1)\). We assume that \(f(n) = 0\) when \(n\) is larger than some big number. The final expression for \(f(1)\) is as follows:

\[\boxed{f(1) = \mu(1)g(1) + \mu(2)g(2) + \mu(3)g(3) … = \sum_{m = 1}^{\infty}\mu(m)g(m)}\]

Variations

We can usually define some \(g’(x)\) and \(f’(x)\) to solve a range of related problems. Here are some common examples, but other substitutions also work.

Calculate \(f(a)\)

Say you want to calculate \(f(a)\) instead of \(f(1)\). Then define \(f’(m) = f(ma)\) and \(g’(m) = g(ma)\) then we have \(g’(n) = \sum_{m = 1}^{\infty}f’(mn)\) which is our original problem. Solving and substituting back gives the general solution:

\[f(n) = \sum_{m = 1}^{\infty}\mu(m)g(mn)\]

Calculate \(f(N)\) where \(g(n) = \sum_{m = 1}^{\infty}f\left(\left\lfloor\frac{N}{mn}\right\rfloor\right)\)

We define \(f’(x) = f\left(\left\lfloor\frac{N}{x}\right\rfloor\right)\). Then we have \(g(n) = \sum_{m = 1}^{\infty}f’(mn)\) which is our original problem. Solving and substituting back gives the solution:

\[f(N) = \sum_{m = 1}^{\infty}\mu(m)g(m)\]

This also works similarly for other scenarios, such as \(g(n) = \sum_{m = 1}^{\infty}f\left(\frac{N}{mn}\right)\)

Calculate \(f(N)\) where \(g(n) = \sum_{m = 1}^{\infty}f\left(\left\lfloor\frac{n}{m}\right\rfloor\right)\)?

We define \(f’(x) = f\left(\left\lfloor\frac{N}{x}\right\rfloor\right)\), and \(g’(x) = g\left(\left\lfloor\frac{N}{x}\right\rfloor\right)\) Then we have the following:

\[ \begin{align*} g’(n) &= \sum_{m = 1}^{\infty}f\left(\left\lfloor\frac{\left\lfloor\frac{N}{n}\right\rfloor}{m}\right\rfloor\right) \\ &= \sum_{m = 1}^{\infty}f\left(\left\lfloor\frac{N}{mn}\right\rfloor\right) \\ &= \sum_{m = 1}^{\infty}f’(mn) \end{align*}
\]

Which is our original problem. Solving and substituting gives:

\[ \begin{align*} f’(1) &= \sum_{m = 1}^{\infty}\mu(m)g’(m) \\ f(N) &= \sum_{m = 1}^{\infty}\mu(m)g\left(\left\lfloor\frac{N}{m}\right\rfloor\right) \end{align*}
\]

Which is a more “traditional” looking Mobius inversion.

What if \(g(n)\) doesn’t sum to infinity?

Say you have \(g(n) = \sum_{m = 1, mn \leq k}f(mn)\) instead of \(g(n) = \sum_{m = 1}^{\infty}f(mn)\). Then define \(f’(x)\) as follows:

\[ f’(n) = \begin{cases} f(x) & x \leq k \\ 0 & \text{otherwise} \\ \end{cases} \]

Then you have \(g(n) = \sum_{m = 1}^{\infty}f’(mn)\) which is our original problem. We also have \(g(n) = 0\) when \(n > k\) by definition. Therefore the solution is given by:

\[f(1) = \sum_{m = 1}^{k}\mu(m)g(m)\]

Diophantine Equations

\(a = xy\)

Say we want to find \(x\), \(y\), and we are given \(a\). This is just factorization. If we have to do this repeatedly for many numbers, the most efficient way is to use a sieve, which will factorise many numbers at once.

\(a = x^2 - y^2\)

This is a difference of squares. In this case, we can factorise out \(a = (x + y)(x - y)\). Let \(m = x + y\) and \(n = x - y\), then we have \(a = mn\) which we can solve by factorisation. Once we have \(m\) and \(n\), we can find \(x\) and \(y\) again.

Modular Arithmetic

This is the big one. Modular arithmetic occurs in almost every difficult competitive programming question, since numbers will almost certainly overflow a fixed width integer. Modular arithmetic forms a finite group. This is good for us since a finite group can be represented with a finite amount of memory.

We note the following:

\[ \begin{align*} (a \mod m) + (b \mod m) &\equiv a + b \mod m \\ (a \mod m) - (b \mod m) &\equiv a - b \mod m \\ (a \mod m) \cdot (b \mod m) &\equiv a \cdot b \mod m \\ (a \mod m)^b &\equiv a^b \mod m \end{align*} \]

fn main() {
    let (a, b, m) = (5, 9, 11);

    assert_eq!(((a % m) + (b % m)) % m, (a + b) % m);
    assert_eq!(((a % m) - (b % m)) % m, (a - b) % m);
    assert_eq!(((a % m) * (b % m)) % m, (a * b) % m);
}

This allows us to reduce the values of \(a\) and \(b\) before we do any calculations, and still end up with the same result, allowing us to work with much smaller numbers than we would normally for most calculations.

Exponentiation is a little different. In the case we have coprime \(a\) and \(p\), then \(a^{\phi(p)} \equiv 1\mod p\) which is given by Euler’s theorem. Therefore:

\[ \begin{align*} a^{b \mod \phi(p)} \equiv a^b \mod p \end{align*} \]

Division

The first step to calculate \(\frac{b}{a} = ba^{-1}\) is to first calculate \(a^{-1}\). We know that \(aa^{-1} \equiv 1 \mod p\). which can be rewritten as:

\[ \begin{align*} aa^{-1} &= 1 \mod p \\ aa^{-1} &= py + 1 \\ aa^{-1} - py &= 1 \\ \end{align*} \]

We know how to solve the similar equation \(ax + by = gcd(a, b)\) using the Bézout’s identity. If \(a\) and \(p\) are coprime, then \(gcd(a, p) = 1\) and we can find the solution to \(aa^{-1} - py = 1 = gcd(a, p)\). If they are not coprime, then there is no inverse, and we cannot perform modular division.

So the process to find \(\frac{b}{a} \mod p\) is as follows:

• Solve for \(x\) in \(ax - py = 1\) using Bézout’s identity. A solution only exists if \(gcd(a, p) = 1\)
• \(x = a^{-1}\) is the inverse of \(a\)
• The solution is \(ba^-1\)

Example

Let’s find \(4/3 \mod 7\). Firstly, we solve \(3x + 7y = gcd(3, 7)\)

Which gives us \(3x + 7y = 3\cdot5 - 7\cdot2 = 1 = gcd(3, 7)\). The gcd is \(1\) so we have a solution \((x, y) = (5, 2)\). \(x\) is the inverse of \(a\). The final step is to multiply by \(b\)

\[ \begin{align*} 4/3 &= 4\cdot3^{-1} \mod 7 \\ &= 4\cdot5 \mod 7 \\ &= 6 \mod 7\\ \end{align*} \]

To verify this is correct we can multiply by \(3\)

\[ \begin{align*} 4/3 &= 6 \mod 7 \\ 4 &= 6\cdot3 \mod 7 \\ &= 18 \mod 7 \\ &= 4 \mod 7 \\ \end{align*} \]

Square roots

Dynamic Programming

Dynamic programming is a method of solving problems that can be expressed recursively. Note that, we can always rewrite a non recursive function recursively. Basically, we save the result of intermediate computations. This is called memoization. Consider the Fibbonacci sequence:

\[ \begin{align*} F_0 &= 0 \\ F_1 &= 1 \\ F_n &= F_{n - 1} + F_{n - 2} \end{align*} \]

And a naive implementation:

fn fibonacci(n: u64) -> u64 {
    match n {
        0 => 0,
        1 => 1,
        n => fibonacci(n - 1) + fibonacci(n - 2)
    }
}

assert_eq!(fibonacci(5), 5);
assert_eq!(fibonacci(10), 55);

This takes \(O(F_n)\) time to calculate, since every addition we perform adds \(1\) to the result. This is because in the process of calculating our final results, we perform many duplicate calls.

Instead of calculating these intermediate values many times over, we can store them in a cache and retrieve them when necessary. This allows us to calculate much larger versions of \(F_n\)

use std::collections::HashMap;

fn fibonacci(cache: &mut HashMap<u64, u64>, n: u64) -> u64 {
    // We have already calculated this before, skip the calculation and return
    // the value
    if let Some(result) = cache.get(&n) {
        return *result;
    }

    // Perform the actual calculation
    let result = match n {
        0 => 0,
        1 => 1,
        n => fibonacci(cache, n - 1) + fibonacci(cache, n - 2)
    };

    // Store the result of the calculation for later
    cache.insert(n, result);
    return result;
}

let mut cache = HashMap::new();
assert_eq!(fibonacci(&mut cache, 60), 1548008755920);
assert_eq!(fibonacci(&mut cache, 78), 8944394323791464);