Project Euler 5: Building the Smallest Number Divisible by 1 to 20 from Least Common Multiples

This is the fifth problem in the Project Euler series. “Divisible by every number” is exactly the definition of a common multiple, and the smallest one is the least common multiple. So this time, instead of hunting for the answer, we build it directly. Along the way we prove from scratch why Euclid’s algorithm and the two-number least common multiple formula are correct. Both come back often in later problems.

Spoiler warning: this post contains a full solution and the answer. If you want to try it yourself first, start with the original problem. Posts in this series never carry ads, because the problem text is licensed under CC BY-NC-SA 4.0, which is non-commercial.

Problem

The problem’s example is that 2520 is the smallest number that every integer from 1 to 10 divides with no remainder. The task is to find the smallest positive integer that is divisible by every integer from 1 to 20.

There is no input; the upper limit is fixed at 20. The output is a single number: the smallest positive integer meeting that condition.

Background

Divisibility and the least common multiple. When an integer m is a multiple of a positive integer k, meaning m = k × q for some integer q, we say “k divides m” and write k | m. A common multiple of several numbers is a multiple of all of them, and their least common multiple (lcm) is the smallest positive common multiple. Multiplying a few positive integers together always gives one common multiple, so the least common multiple always exists.

Euclid’s algorithm. The greatest common divisor (gcd) of two numbers comes from the following property. If dividing a by b (> 0) gives quotient q and remainder r (0 ≤ r < b), then

gcd(a, b) = gcd(b, r),    r = a − q × b
gcd(a, 0) = a

The reason is that the two pairs have the same common divisors. If d divides a and b, it also divides a − q × b = r; conversely, if d divides b and r, it also divides q × b + r = a. With the same common divisors, the greatest of them is the same too. The remainder strictly decreases every step, so it eventually reaches 0, and the number left at that point is the greatest common divisor. For example, gcd(840, 9) is 3, because 840 = 93 × 9 + 3 and 9 = 3 × 3 + 0.

The gcd is an integer combination of the two numbers (Bézout’s identity). Each remainder r = a − q × b in the algorithm is an integer combination of the two numbers before it, and those are integer combinations of the original a and b, so by following the steps every remainder has the form s × a + t × b with integers s and t. The last nonzero remainder is the gcd, so there are integers s and t with gcd(a, b) = s × a + t × b. In the example above, 3 = 1 × 840 − 93 × 9. The proof below uses this fact once.

The least common multiple through prime factorization (for hand checks). Comparing the exponent of each prime, the least common multiple takes the largest exponent for every prime. For example 2520 = 2³ × 3² × 5 × 7, where 8 = 2³ and 9 = 3² are the largest powers of 2 and 3 not above 10. The final algorithm does not factor anything; we use factorization only to check results by hand.

Brute force

Try m = 1, 2, 3, … in turn and take the first m that every k from 2 to 20 divides.

m = 1
while not all(m % k == 0 for k in range(2, 21)):
    m += 1

It is guaranteed to finish: the product of 1 through 20 meets the condition, so the loop stops at or before it. The problem is that the number of candidates it tests is the answer itself. Each candidate costs up to 19 remainder operations, so the time is O(answer × n), where n is the upper limit. And each time the limit passes a prime p, the answer is multiplied by p (no smaller number has p as a factor), so the answer and the number of candidates grow quickly with n. For n = 10 there are only 2520 candidates, but for n = 20 there are far more; exactly how many we will see in “Answer and checks”.

Stepping through candidates 20 at a time (the answer is a multiple of 20) or 2520 at a time (the answer is a common multiple of 1 to 10, and so, by Lemma A below, a multiple of 2520) cuts the count by that factor. But it still looks at a number of candidates proportional to the answer. We want a method whose running time does not depend on how large the answer is.

Derivation

Step 1: reread the problem. “A positive integer divisible by every integer from 1 to n” is a positive common multiple of 1, 2, …, n, and the smallest one is, by definition, lcm(1, 2, …, n). Write it L(n). There is nothing to search for: we just compute L(20).

Step 2 (Lemma A): every common multiple is a multiple of the least common multiple. Let l be the least common multiple of some numbers and m any positive common multiple. Dividing m by l gives m = q × l + r with 0 ≤ r < l. Each of the numbers divides both m and l, so it divides r = m − q × l as well; r is a common multiple. If r > 0 it would be a positive common multiple smaller than l, contradicting the minimality of l, so r = 0, that is, l | m. Conversely, any multiple of l is clearly a common multiple. So the common multiples are exactly the multiples of l.

Step 3: fold the numbers in one at a time. L(k) = lcm(L(k − 1), k). For a positive integer m,

• m is a common multiple of 1, …, k
• ⇔ m is a common multiple of 1, …, k − 1, and k | m
• ⇔ (Lemma A) L(k − 1) | m and k | m
• ⇔ m is a common multiple of L(k − 1) and k

so the two sets of common multiples are the same, and therefore so are their smallest elements. All that remains is the least common multiple of two numbers.

Step 4: for two numbers, lcm(a, b) = a / gcd(a, b) × b. For positive integers a and b, let g = gcd(a, b) and write a = g × a′, b = g × b′.

• gcd(a′, b′) = 1. If some d > 1 divided both a′ and b′, then g × d would be a common divisor of a and b, contradicting the maximality of g.
• l = g × a′ × b′ is a common multiple, because l = a × b′ = a′ × b. Also l = (a / g) × b, and a / g is an integer since g | a.
• l is the smallest. Write any positive common multiple as m = a × x. Then m = g × a′ × x, and since b = g × b′ divides m, b′ | a′ × x. Take integers s and t with s × a′ + t × b′ = 1 from Bézout’s identity. Then x = s × a′ × x + t × b′ × x, and b′ divides both terms on the right, so b′ | x. Hence x ≥ b′ and m = a × x ≥ a × b′ = l.

Checking a small example by hand. Follow the problem’s example, n = 10. Start from L = 1 and fold in k = 2, …, 10 in turn.

k    g = gcd(L, k)    L / g × k
2    gcd(1, 2)   = 1     1 / 1 × 2  = 2
3    gcd(2, 3)   = 1     2 / 1 × 3  = 6
4    gcd(6, 4)   = 2     6 / 2 × 4  = 12
5    gcd(12, 5)  = 1    12 / 1 × 5  = 60
6    gcd(60, 6)  = 6    60 / 6 × 6  = 60
7    gcd(60, 7)  = 1    60 / 1 × 7  = 420
8    gcd(420, 8) = 4   420 / 4 × 8  = 840
9    gcd(840, 9) = 3   840 / 3 × 9  = 2520
10   gcd(2520, 10) = 10   2520 / 10 × 10 = 2520

The result, 2520, matches the problem’s example. In terms of prime factors, 4 raises the exponent of 2 from 1 to 2, 8 raises it from 2 to 3, and 9 raises the exponent of 3 from 1 to 2. 6 and 10 already divide L, so g = k and L does not change. This agrees with the final 2520 = 2³ × 3² × 5 × 7.

Final algorithm

1. Start with L = 1.
2. For k = 2, 3, …, n in order, compute g = gcd(L, k) with Euclid’s algorithm and replace L by L / g × k.
3. Return L. If n < 2 the loop never runs and the result is 1.
4. Before printing the answer, check small inputs: gcd(12, 18) = 6, gcd(840, 9) = 3, lcm(4, 6) = 12, L(1) = 1 and L(10) = 2520, and the brute-force search for n = 10 must also give 2520. If any of these differ, exit with a nonzero status.

Correctness proof

Right after the k-th iteration, the following holds (the invariant).

L = lcm(1, 2, …, k).

Before the loop (k = 1), L = 1 = lcm(1), so it holds. Suppose it holds up to k − 1. In the k-th iteration, g = gcd(L, k) is computed exactly by Euclid’s algorithm from the background section, the new value L / g × k is lcm(L(k − 1), k) by step 4, and that is L(k) by step 3. So when the loop finishes at k = n, L = L(n), which by step 1 is the smallest positive integer divisible by every integer from 1 to n. The loop runs exactly n − 1 times, and each call of Euclid’s algorithm ends because the remainder decreases every step.

Complexity

Time. gcd is called n − 1 times. gcd(L, k) becomes (k, L mod k) after the first division and only deals with numbers smaller than k from then on. In Euclid’s algorithm, two steps (a, b) → (b, r) → (r, ·) starting from a ≥ b leave the first number below half of a: if b ≤ a / 2 then r < b ≤ a / 2, and if b > a / 2 the quotient is 1, so r = a − b < a / 2. So one call of gcd(L, k) takes O(log k) divisions, and the whole computation takes O(n log n) arithmetic operations (while the values fit in a machine integer). Only a few values are stored, so the extra space is O(1). Unlike the O(answer × n) of the brute force, this does not depend on how large the answer is.

Integer range.

• Every value during the loop is at most the final result. L only grows or stays the same (L(k − 1) | L(k)), the intermediate L / g is at most L, and L / g × k is that step’s result L(k). So a type that holds the final result never overflows along the way.
• Divide before multiplying. Writing L × k / g makes L × k up to k times larger than the result. That does not matter for n = 20, but take L(43) = 2⁵ × 3³ × 5² × 7 × 11 × 13 × 17 × 19 × 23 × 29 × 31 × 37 × 41 × 43 = 9419588158802421600 (about 9.42 × 10¹⁸) and fold in k = 44: L × 44 exceeds 2⁶⁴ − 1 (about 1.84 × 10¹⁹), the largest 64-bit unsigned value, while L / 44 × 44 = L does not (4 and 11 divide L, so g = 44).
• The C++ code uses unsigned long long, which the standard guarantees to be at least 64 bits wide. By the calculation above, 44, 45 and 46 are not powers of a new prime, so L(46) = L(43) fits in 64 bits, while L(47) = 47 × L(46) does not. That is why the code comment says limit <= 46. On overflow, unsigned integers wrap around modulo 2ᴺ (N being the type’s width in bits) instead of causing undefined behavior, which means the program does not stop but silently produces a wrong value, so staying within range matters.
• Python’s int has unlimited precision, so it cannot overflow. But / produces a floating-point result, so we use integer division //. Since g divides L, // gives the exact quotient.

C++17

#include <iostream>

// Greatest common divisor by Euclid's algorithm: gcd(a, b) = gcd(b, a mod b).
unsigned long long gcd(unsigned long long a, unsigned long long b) {
    while (b != 0) {
        const unsigned long long r = a % b;
        a = b;
        b = r;
    }
    return a;
}

// Least common multiple of two positive numbers. Dividing before multiplying
// keeps every intermediate value no larger than the result.
unsigned long long lcm(unsigned long long a, unsigned long long b) {
    return a / gcd(a, b) * b;
}

// lcm(1, 2, ..., limit), folding in one number at a time. Returns 1 when limit < 2.
// Meant for limit <= 46, where the result fits in 64 bits.
unsigned long long lcm_up_to(unsigned long long limit) {
    unsigned long long result = 1;
    for (unsigned long long k = 2; k <= limit; ++k) {
        result = lcm(result, k);
    }
    return result;
}

// Naive search: the first positive m divisible by every k in 1..limit.
// Used only to cross-check the small example.
unsigned long long smallest_by_search(unsigned long long limit) {
    for (unsigned long long m = 1;; ++m) {
        bool divisible = true;
        for (unsigned long long k = 2; k <= limit; ++k) {
            if (m % k != 0) {
                divisible = false;
                break;
            }
        }
        if (divisible) {
            return m;
        }
    }
}

int main() {
    if (gcd(12, 18) != 6 || gcd(840, 9) != 3 || lcm(4, 6) != 12 ||
        lcm_up_to(1) != 1 || lcm_up_to(10) != 2520 ||
        smallest_by_search(10) != 2520) {
        return 1;
    }
    std::cout << lcm_up_to(20) << '\n';
    return 0;
}

smallest_by_search is the brute force as written above; it is used only to check that it agrees with the improved method on the official example, n = 10. The small examples are checked with if, not assert. assert does nothing in builds where NDEBUG is defined, whereas the comparisons in the if and the return 1 are part of the program’s observable behavior, so -O2 optimization cannot remove them. If any check fails, the program exits with a nonzero status.

Python 3

def gcd(a: int, b: int) -> int:
    """Greatest common divisor by Euclid's algorithm: gcd(a, b) = gcd(b, a mod b)."""
    while b != 0:
        a, b = b, a % b
    return a


def lcm(a: int, b: int) -> int:
    """Least common multiple of two positive numbers.

    Dividing before multiplying keeps every intermediate value no larger than the result.
    """
    return a // gcd(a, b) * b


def lcm_up_to(limit: int) -> int:
    """lcm(1, 2, ..., limit), folding in one number at a time. Returns 1 when limit < 2."""
    result = 1
    for k in range(2, limit + 1):
        result = lcm(result, k)
    return result


def smallest_by_search(limit: int) -> int:
    """Naive search: the first positive m divisible by every k in 1..limit.

    Used only to cross-check the small example.
    """
    m = 1
    while True:
        if all(m % k == 0 for k in range(2, limit + 1)):
            return m
        m += 1


def main() -> None:
    if (
        gcd(12, 18) != 6
        or gcd(840, 9) != 3
        or lcm(4, 6) != 12
        or lcm_up_to(1) != 1
        or lcm_up_to(10) != 2520
        or smallest_by_search(10) != 2520
    ):
        raise SystemExit(1)
    print(lcm_up_to(20))


if __name__ == "__main__":
    main()

Same algorithm, same small examples as the C++ version. a, b = b, a % b evaluates the whole right-hand side before assigning, so it does the job of the temporary r in C++. Because range(2, limit + 1) excludes its end value, the end is one higher. If a small example fails, raise SystemExit(1) gives exit status 1. assert is not used because running with python -O removes it.

Answer and checks

The smallest positive integer divisible by every integer from 1 to 20 is 232792560.

Continuing from the n = 10 result 2520 and folding in k = 11, …, 20 gives:

k    g = gcd(L, k)    L
11   1      27720
12   12     27720
13   1      360360
14   14     360360
15   15     360360
16   8      720720
17   1      12252240
18   18     12252240
19   1      answer above
20   20     answer above

The parts that can be checked by hand:

• Prime factorization. The primes up to 20 are 2, 3, 5, 7, 11, 13, 17 and 19, and the largest powers not above 20 are 16 = 2⁴, 9 = 3² and the remaining primes themselves. 16 × 9 = 144, × 5 = 720, × 7 = 5040, × 11 = 55440, × 13 = 720720, × 17 = 12252240, × 19 = the answer above.
• Divisible by all of them. In the factorization of any k ≤ 20, each prime power pᵉ is at most k, hence at most 20, so its exponent does not exceed the one above. Therefore k divides the answer.
• The smallest. Lowering any one of the exponents for 16, 9, 5, 7, 11, 13, 17, 19 makes that number fail to divide it. For example, reducing 2⁴ to 2³ means 16 no longer divides it.
• The link to 2520. The answer above = 2520 × 92378, and 92378 = 2 × 11 × 13 × 17 × 19. Going from n = 10 to n = 20 adds only one more 2 (for 16) and four new primes.

This also shows that the brute force would have to test that many candidates to reach this value. I did not measure its running time for this post.

The checks come in three layers.

1. Small examples. Before printing the answer, both programs check gcd(12, 18) = 6, gcd(840, 9) = 3, lcm(4, 6) = 12, L(1) = 1 and L(10) = 2520, and also run the brute-force search for n = 10 and check that it equals 2520. If any differ, they exit with a nonzero status. The n = 10 case was also traced by hand in the derivation.
2. Agreement between the two implementations. In this series’ repository a check script compiles the C++ code with g++ -std=c++17 -O2 -Wall -Wextra -Werror and runs it (warnings count as failures), and runs the Python code with python3 -I -B. Both outputs must match the stored answer.txt byte for byte, and a nonzero exit status or a timeout (1 second for C++, 10 seconds for Python) is a failure. By the series rules only code that passes this check is published.
3. The proof. The correctness proof showed that folding in least common multiples one at a time gives exactly L(n).

Boundary cases

• n = 1: the loop never runs and the function returns 1. 1 is the smallest positive integer divisible by 1, so that is right. The small examples include L(1) = 1. n = 0 is not part of the problem, but for the same reason it also returns 1, matching the convention that the least common multiple of an empty set is 1.
• k that already divides L (6, 10, 12, 14, 15, 18, 20): g = gcd(L, k) = k, so L / k × k = L stays the same.
• A new prime k (7, 11, 13, 17, 19): every prime factor of L is smaller than k, so g = 1 and L is multiplied by k.
• A prime power k (4, 8, 9, 16): only one exponent goes up by one. For 16, for example, g = 8, so L only doubles.
• Multiplying before dividing: as in “Integer range”, L × k can overflow once n is large. Dividing first keeps every intermediate value no larger than the result.

Going further

• Compute it straight from prime powers. As used in the answer check, L(n) is the product, over the primes p ≤ n, of the largest power pᵉ not above n. Change the program to find the primes with the sieve of Eratosthenes and build this product. Compared with calling gcd n − 1 times, which one does fewer operations? Working with primes comes up again in problems 7 and 10 of this series.
• Does the order matter? The argument in step 3 does not depend on going “from 1 to k”. Show in the same way that folding the numbers in any order gives the same least common multiple.
• Finding the Bézout coefficients. Extending Euclid’s algorithm also yields the s and t in gcd(a, b) = s × a + t × b. Try it by hand for 840 and 9.
• The product of gcd and lcm. Step 4 gives gcd(a, b) × lcm(a, b) = a × b right away. Does the same relation hold for three or more numbers? Look for a small counterexample.

Sources and license

• Project Euler, Problem 5 — the original problem. The “Problem” section of this post is my own summary of it and is shared under the same CC BY-NC-SA 4.0 license as the original. Attribution for the problem goes to Project Euler.
• Project Euler, Copyright — the license of the problem text and how to link to problems. Because of the non-commercial condition, this post carries no ads (ads: false).
• The minimum width of C++ integer types and unsigned arithmetic follow the C++ working draft [basic.fundamental], and assert and NDEBUG follow [assertions.assert]. For Python, simultaneous assignment follows Assignment statements, assert and -O follow The assert statement, and the precision of int and the results of / and // follow Numeric Types.
• The solution, proofs and code are my own; the properties of Euclid’s algorithm and the least common multiple are proved directly above. I did not consult Project Euler’s overviews, forum, member-only material or other solution sites. Motitory is not affiliated with or sponsored by Project Euler, and Project Euler has not reviewed this post.