Smallest Multiple — the LCM of 1..20 from three angles

Project Euler problem #5 is a one-liner if you reach for math.lcm and a multi-hour wall clock if you reach for brute force. The interesting part isn't the answer; it's that the three obvious ways of solving it differ by nearly six orders of magnitude (~600,000×) in runtime and one of them lays the structure of the integers bare in a way the others hide.

What the problem asks

Project Euler #5 — "Smallest Multiple" — asks for the smallest positive integer that every integer from 1 to 20 divides without a remainder. In standard notation, that's lcm(1, 2, 3, ..., 20). I'll restate it once mathematically and then leave the official phrasing on Colin Hughes' page where it belongs:

Find the smallest positive integer $M$ such that $k \mid M$ for every $k \in \{1, 2, \ldots, 20\}$.

The problem is published at difficulty level 0 on PE, which means it's an early-archive warm-up. That doesn't mean it's uninteresting — it means the mathematical content is just below the surface, and the value of writing it up is in seeing the three methods side by side.

I deliberately won't print the numerical answer in this post (PE asks solvers not to spoil it, and that's a reasonable ask). The artefacts at the bottom will compute it in a few microseconds on whatever machine you run them on.

Why Python

Two reasons. First, the Project Euler "one-minute rule" — a solution should run in under a minute — is so loose for this problem that any language wins; performance is not a constraint, so we should optimise for clarity. Second, Python ships math.gcd and math.lcm in the standard library (since 3.5 and 3.9 respectively), and CPython's int is arbitrary-precision out of the box. For an LCM problem where the answer can grow exponentially (as e^N) with N (more on that below), arbitrary precision matters: in C you'd be choosing between uint64_t and dragging in GMP after about N = 40, and that's a distraction from the arithmetic.

A C or Rust implementation would be three times the lines for no insight gained. A SymPy implementation would be one line (sympy.lcm_list(range(1, 21))) with too much black box for a teaching post. Plain Python is the right granularity.

The naïve approach, and why it has to scale badly

Here is the first thing anyone writes when they see this problem:

"Start at 1. Keep going up. Return the first number every k in 1..20 divides."

That's correct. It is also intolerably slow, and the reason it's slow tells you what the right approach has to look like.

Think about it: the answer must be a positive integer that 20 divides, so it must be a multiple of 20 — start the search at 20, step by 20, and immediately you've cut the search space by a factor of 20. You can extend this to "the candidate must be a multiple of lcm(1..N)" — but the moment you've said that out loud, you've answered the problem; brute force is no longer brute. So as a genuine "look at every integer in turn" algorithm, you'll be testing ~lcm(1..N) / N candidates, and each test is O(N) modulus operations. The total work is O(lcm(1..N)) modular operations, and by the Prime Number Theorem, log lcm(1..N) ~ N (Chebyshev's $\psi$ function), so brute force is exponential in N in the worst case.

For N = 20, that's "a couple of seconds" — annoying but tolerable. For N = 40, you'd be waiting over a year. The brute-force code I wrote terminates in 2.87 s for N = 20 on the sandbox machine (timing below).

The brute-force routine, from smallest_multiple.py:

def smallest_multiple_brute(n: int) -> int:
    candidate = n
    while True:
        for k in range(2, n + 1):
            if candidate % k:
                break
        else:
            return candidate
        candidate += n

Two micro-optimisations are baked in already — start at n and step by n — and even with those it's the slowest of the three by a factor of ~600,000. That's the cost of refusing to use structure.

The LCM identity — the textbook approach

The textbook fact is

with the proviso that both sides are non-negative integers (the formula on the right is automatically an integer because gcd(a, b) divides a). The proof is one sentence: write a and b in their prime factorisations $a = \prod p_i^{\alpha_i}$, $b = \prod p_i^{\beta_i}$; then lcm = \prod p_i^{\max(\alpha_i, \beta_i)} and gcd = \prod p_i^{\min(\alpha_i, \beta_i)}; and max + min = α + β.

Because LCM is associative — $\mathrm{lcm}(a, \mathrm{lcm}(b, c)) = \mathrm{lcm}(\mathrm{lcm}(a, b), c)$ — you can compute lcm(1..N) by folding the binary operation from the left:

Each step costs one gcd (Euclidean algorithm, O(log min(a,b)) division steps) and one multiplication on big integers. With CPython's arbitrary-precision arithmetic, that's a few microseconds total for N = 20:

def smallest_multiple_lcm(n: int) -> int:
    return reduce(math.lcm, range(1, n + 1), 1)

This is the answer most people would commit to git. It's correct, it's fast, and it doesn't tell you anything about why the answer has the value it has. For that, you need the third approach.

Prime-power product — the structural answer

Step back from the gcd–lcm identity and look at the prime factorisations directly. Every integer k in 1..N has a unique factorisation $k = \prod p^{e_p(k)}$ over primes p. The LCM of a set of integers takes the max exponent across the set, prime by prime:

What's $\max_{k \le N} e_p(k)$? It's the largest exponent $e$ such that $p^e \le N$, because the integer $p^e$ itself appears in $\{1, \ldots, N\}$. (Any $k \le N$ with $e_p(k) = e$ is at least $p^e$, so $p^e \le N$.) That is to say,

This is a complete, closed-form structural description of the answer in terms of the primes below N. The algorithm writes itself: sieve the primes up to N, and for each prime keep multiplying it into itself until one more multiplication would push past N.

def primes_upto(n: int) -> list[int]:
    """Plain sieve of Eratosthenes."""
    if n < 2:
        return []
    sieve = bytearray(b"\x01") * (n + 1)
    sieve[0] = sieve[1] = 0
    for i in range(2, int(n ** 0.5) + 1):
        if sieve[i]:
            sieve[i * i :: i] = b"\x00" * len(sieve[i * i :: i])
    return [i for i, v in enumerate(sieve) if v]

def smallest_multiple_primes(n: int) -> int:
    result = 1
    for p in primes_upto(n):
        pk = p
        while pk * p <= n:
            pk *= p
        result *= pk
    return result

(From smallest_multiple.py.) This is the same asymptotic complexity class as the LCM fold for small N, but it has two virtues the LCM fold lacks. First, the intermediate value never decreases — there's no gcd "correction" being applied per step — so you can see the answer accumulating prime by prime. Second, you get a free decomposition: the loop emits exactly the table that constitutes the structural answer.

For N = 20, that table is:

(Produced by python3 smallest_multiple.py 20, output included in the timing section below.)

The answer is the product of the right-hand column. Eight numbers, six of them prime, two of them prime powers (16 and 9). I'm not going to multiply them out in the post body, but you can do it on the back of a napkin in about ten seconds — and that's the entire point of this approach. The LCM of the first 20 positive integers is a tiny, comprehensible mathematical object.

Why exponent 4 on 2 and not 5? Because $2^4 = 16 \le 20$ but $2^5 = 32 > 20$. Why exponent 2 on 3 and not 3? Because $3^2 = 9 \le 20$ but $3^3 = 27 > 20$. Why exponent 1 on 5 through 19? Because $5^2 = 25 > 20$, so no $k \le 20$ has more than one factor of 5, and similarly for the larger primes. The structure makes the answer feel inevitable once you see it.

The whole solver, in one file

Here is smallest_multiple.py in full — the three approaches plus the timing harness:

#!/usr/bin/env python3
"""
Project Euler #5 — Smallest Multiple.

Three solutions to lcm(1..N), wired up for benchmarking and cross-checking.

Run:
    python3 smallest_multiple.py            # N = 20
    python3 smallest_multiple.py 40         # any N you like
"""

from __future__ import annotations
import math
import sys
import time
from functools import reduce

# ---------------------------------------------------------------------------
# Approach A. Brute-force search.
#
# Start at N (the candidate must be at least N) and keep adding N until the
# candidate is divisible by every k in 2..N. This is the "obvious" first idea
# and it is correct, but it gets very slow as N grows because the answer
# grows roughly as exp(N) (by the prime number theorem applied to log lcm).
# ---------------------------------------------------------------------------
def smallest_multiple_brute(n: int) -> int:
    candidate = n
    while True:
        for k in range(2, n + 1):
            if candidate % k:
                break
        else:
            return candidate
        candidate += n

# ---------------------------------------------------------------------------
# Approach B. Pairwise LCM reduction.
#
# lcm(a, b) = a * b // gcd(a, b). Fold this across 1..N.
# Each gcd is O(log min(a, b)) via the Euclidean algorithm, so the whole
# reduction is essentially O(N log M) where M is the running LCM.
# Python's math.lcm already does this for you in one call, but the fold makes
# the algorithm visible.
# ---------------------------------------------------------------------------
def smallest_multiple_lcm(n: int) -> int:
    return reduce(math.lcm, range(1, n + 1), 1)

# ---------------------------------------------------------------------------
# Approach C. Prime-power product.
#
# For each prime p <= N, the largest power of p that appears in any k in 1..N
# is p^floor(log_p N). Multiply those across all primes <= N and you have the
# LCM directly, no gcd loop, no growing intermediate.
#
#   lcm(1..N) = product over primes p <= N of  p^floor(log_p N)
#
# This is the structural answer: it tells you WHY the LCM looks the way it
# does. For N=20 the primes are 2,3,5,7,11,13,17,19 with exponents
# 4,2,1,1,1,1,1,1.
# ---------------------------------------------------------------------------
def primes_upto(n: int) -> list[int]:
    """Plain sieve of Eratosthenes."""
    if n < 2:
        return []
    sieve = bytearray(b"\x01") * (n + 1)
    sieve[0] = sieve[1] = 0
    for i in range(2, int(n ** 0.5) + 1):
        if sieve[i]:
            sieve[i * i :: i] = b"\x00" * len(sieve[i * i :: i])
    return [i for i, v in enumerate(sieve) if v]

def smallest_multiple_primes(n: int) -> int:
    result = 1
    for p in primes_upto(n):
        # Largest k such that p**k <= n.
        pk = p
        while pk * p <= n:
            pk *= p
        result *= pk
    return result

def time_it(fn, *args, repeats: int = 5) -> tuple[int, float]:
    best = math.inf
    value = None
    for _ in range(repeats):
        t0 = time.perf_counter_ns()
        value = fn(*args)
        dt = time.perf_counter_ns() - t0
        if dt < best:
            best = dt
    return value, best / 1e6  # ms

def main() -> None:
    n = int(sys.argv[1]) if len(sys.argv) > 1 else 20

    print(f"Project Euler #5 — smallest multiple of 1..{n}")
    print(f"{'method':<22}{'time (ms, best of 5)':>26}  agreement")
    print("-" * 64)

    a, ta = time_it(smallest_multiple_brute, n) if n <= 22 else (None, None)
    b, tb = time_it(smallest_multiple_lcm, n)
    c, tc = time_it(smallest_multiple_primes, n)

    # We do not print the actual answer, per Project Euler's request.
    # We only confirm that the three methods agree.
    ok_bc = (b == c)
    ok_ab = (a is None) or (a == b)
    print(f"{'A brute force':<22}{(ta if ta is not None else float('nan')):>20.4f} ms      "
          f"{'(skipped, too slow)' if ta is None else 'matches B/C' if ok_ab else 'DISAGREES'}")
    print(f"{'B pairwise lcm':<22}{tb:>20.4f} ms      "
          f"{'matches C' if ok_bc else 'DISAGREES'}")
    print(f"{'C prime powers':<22}{tc:>20.4f} ms      reference")

    print()
    print("Structural decomposition (primes p <= N, exponent floor(log_p N)):")
    for p in primes_upto(n):
        k = 0
        pk = 1
        while pk * p <= n:
            pk *= p
            k += 1
        print(f"  p = {p:>3}    exponent = {k}    p^k = {pk}")
    digits = len(str(c))
    print(f"\nThe answer has {digits} decimal digits. Run this script to see it.")

if __name__ == "__main__":
    main()

The brute-force method is a guard: if you ask for N larger than 22, it's silently skipped, because for N = 23 it's already over a minute. The two fast methods are run in either case.

Verification and timing

Running on the sandbox machine (Kali container, Python 3.13, single thread):

$ python3 smallest_multiple.py 20
Project Euler #5 — smallest multiple of 1..20
method                      time (ms, best of 5)  agreement
----------------------------------------------------------------
A brute force                    2885.0630 ms      matches B/C
B pairwise lcm                      0.0047 ms      matches C
C prime powers                      0.0054 ms      reference

Structural decomposition (primes p <= N, exponent floor(log_p N)):
  p =   2    exponent = 4    p^k = 16
  p =   3    exponent = 2    p^k = 9
  p =   5    exponent = 1    p^k = 5
  p =   7    exponent = 1    p^k = 7
  p =  11    exponent = 1    p^k = 11
  p =  13    exponent = 1    p^k = 13
  p =  17    exponent = 1    p^k = 17
  p =  19    exponent = 1    p^k = 19

(Direct shell transcript.) All three methods produce the same integer. The "best of 5" timing isolates the algorithmic cost from cold-start noise; the actual call costs are so small that interpreter overhead dominates.

The headline ratio is 2.87 seconds versus 5 microseconds: brute force is ~600,000× slower than either of the structural approaches. The PE one-minute rule is comfortably satisfied by all three for N = 20, but only the latter two are usable past N ≈ 25.

A scan across small N to make the divergence visible (from benchmark.py):

  N      brute (ms)    lcm-fold (ms)    prime-pow (ms)
 10           0.068          0.00281           0.00606
 15           5.289          0.00351           0.00531
 20        2874.383          0.00482           0.00562
 25       (skipped)          0.00593           0.00848

The brute-force column grows by roughly a factor of e ≈ 2.7 per +1 in N on average (it's not exactly geometric because the multiplier is the prime-power update; the steep jump from N=15 to N=20 reflects the LCM growing by a factor of 2×17×19 = 646 across N=16..20 (2^4 replaces 2^3, then primes 17 and 19 enter)). The two structural methods are basically flat — they're bottlenecked by the sieve and the big-int multiply, not the iteration count.

A sanity check: Chebyshev's psi

The structural formula

is the Chebyshev function $\psi$ in disguise. The Prime Number Theorem is equivalent to $\psi(N) \sim N$, so $\log \mathrm{lcm}(1, \ldots, N)$ should grow linearly in N with slope 1 (natural log). Computing this from the prime-power product as a cross-check (from growth_table.py):

   N   digits    log10 lcm    log lcm / N
   5        2       1.7782         0.8189
  10        4       3.4014         0.7832
  20        9       8.3670         0.9633
  50       22      21.4912         0.9897
 100       41      40.8434         0.9405
 200       90      89.5280         1.0307
 500      218     217.8647         1.0033
1000      433     432.8530         0.9967
2000      867     866.1793         0.9972

The right-hand column — $\ln \mathrm{lcm}(1..N) / N$ — wobbles around 1, exactly as PNT predicts, with the wobble coming from the irregular spacing of prime powers (every time $N$ crosses a fresh prime power $p^k$, $\psi$ jumps by $\log p$, not $\log p^k$ — so the jumps at 4, 8, 9, 16, 25 all contribute the same amount as the parent prime). By N = 500 we're within half a percent of 1.

This is a much stronger result than the problem asks for. PE #5 wants one specific integer for N = 20; the structural method also tells you that the answer for any other N would be the exponential of something very close to N, which is why this problem class doesn't admit a polynomial-time brute force.

A note on the search-by-stepping refinement

There's a folk improvement to brute force that often shows up in PE writeups: instead of stepping by N, step by the LCM of the integers you've already proven divide your candidate. Concretely — find the smallest multiple of 2, then the smallest multiple of lcm(2, 3) divisible by 4, and so on. This effectively recovers Approach B at the price of being uglier; it's not really a separate algorithm, just the LCM fold written iteratively with extra steps. I mention it because if you wrote the brute-force code and noticed the redundant work, you've already discovered the right answer.

What I didn't do

A few things I deliberately did not pursue, in case you're wondering whether they'd help:

• A "smarter" sieve. For N = 20 the sieve runs over 21 bytes; profiling it is pointless. For very large N you'd switch to a segmented sieve or a wheel, but PE #5 isn't asking for that.
• GMP. CPython's big-int implementation is plenty for N up to a few thousand. The 41-digit answer for N = 100 comes back in 12 microseconds; you don't need a faster bignum library.
• Memoisation across N. Each call recomputes the sieve from scratch. If you were doing this in a loop over many N, you'd cache; for one call, the constant factors don't matter.

Run it yourself

Everything in this post is in two files in the artefacts tarball — smallest_multiple.py (the solver) and growth_table.py (the Chebyshev check). With Python ≥ 3.9 anywhere:

python3 smallest_multiple.py        # the answer to PE #5
python3 smallest_multiple.py 100    # any N you like
python3 growth_table.py             # psi(N)/N for a sweep of N

The number you get back is your answer to submit to Project Euler. The structural decomposition above tells you why it has the form it does.

References

• Project Euler problem #5, "Smallest Multiple": https://projecteuler.net/problem=5. Published 30 November 2001 by Colin Hughes; difficulty level 0.
• Hardy & Wright, An Introduction to the Theory of Numbers, §22.1 — Chebyshev's $\psi$ function and the equivalence with the Prime Number Theorem.
• CPython source, Modules/mathmodule.c, for the math.gcd / math.lcm implementations used by Approach B.

Artefacts

The tarball for this post contains:

• smallest_multiple.py — the three solvers and the timing harness.
• growth_table.py — the Chebyshev $\psi$ sanity check.
• benchmark.py — the cross-N timing scan.