A closed form for Project Euler #1, and why I reach for Python anyway

The first Project Euler problem is the one everyone solves twice: once with a for loop and if k % 3 == 0 or k % 5 == 0, and a second time when they realise that the sum has a closed form and the loop never had to exist. This is a post about that second pass — the inclusion–exclusion identity, why it scales to absurd bounds, and why I pick Python for it even though the loop is the slow language's home turf.

Why Python

Project Euler problem 1 (projecteuler.net/problem=1) is one of those problems whose "real" form has nothing to do with the stated bound. The stated bound is one thousand. The interesting question is: what does your method do when the bound is 10^18, or 10^100? The numbers you handle then exceed every fixed-width integer your CPU has — u64 tops out around 1.8 × 10^19, u128 around 3.4 × 10^38 — and the answer for N = 10^100 has roughly 200 digits.

That points squarely at Python. Python's int is arbitrary precision by default, with no BigInt import, no from __future__, no sigils. The closed-form arithmetic you'll see below — k * m * (m + 1) // 2 — does the right thing whether m is six or six hundred digits long, with the same source. In Rust I would be writing num_bigint::BigUint::from(k) * &m * (&m + 1u32) / 2u32 and threading Clone and lifetime errors through it; in C I would be reaching for GMP. Neither is hard, but neither is what I want to read in a writeup whose point is the math, not the bignum library.

The other reason: math.lcm, math.gcd, itertools.combinations, time.perf_counter_ns are all in the standard library and all cheap to call. The code that follows has zero third-party imports.

What the problem actually asks

Project Euler #1, "Multiples of 3 or 5", asks for the sum of every positive integer strictly less than 1000 that is divisible by 3 or by 5. The problem itself provides a sanity case: doing the same exercise with the bound 10 gives the four numbers 3, 5, 6, 9, which sum to 23.

I'll write the bound as N from here on so the discussion isn't tied to the specific value 1000. Define

The official problem is S(1000). The interesting object is the function S.

The obvious approach, and why it is fine but not interesting

The straightforward implementation is a one-liner:

def brute(N: int) -> int:
    return sum(k for k in range(1, N) if k % 3 == 0 or k % 5 == 0)

For N = 1000 this runs in roughly 100 microseconds on a laptop. The Project Euler "one-minute rule" says a correct solution should finish in under a minute, and at N = 1000 the brute force finishes in roughly 600,000× less time than that budget. So there is no engineering pressure to do anything clever for the stated bound.

But there is aesthetic pressure. The brute force is O(N). As soon as someone hands you the same problem with N = 10^12 or N = 10^18, you cannot run the loop to completion in any human lifetime, regardless of language. A closed form, if one exists, is O(1) — and one exists.

A handful of other things you might try, ranked by how silly they get:

The right move for this category of problem is almost always the third row up from the bottom: replace the loop with arithmetic that depends only on N and the divisors, not on the elements themselves.

Building the closed form

The pivot is the identity for the sum of multiples of a single integer.

Fix a positive divisor k. The positive multiples of k strictly below N are

Note the N - 1, not N: the problem says strictly below N, so when N is itself a multiple of k we want to exclude it. (N - 1) // k is the largest integer m such that m * k < N. Off-by-ones in this kind of derivation are the most common bug in writeups; staring at the formula for a beat is worth it.

Their sum factors:

That's the triangular number T_m = m(m+1)/2 scaled by k. So we already have a constant-time formula for "sum of multiples of a single divisor under a bound":

def sum_multiples_below(N: int, k: int) -> int:
    m = (N - 1) // k
    return k * m * (m + 1) // 2

This much most readers know. The interesting step is gluing two of these together.

Inclusion–exclusion: do not double-count 15

If I just write sum_multiples_below(N, 3) + sum_multiples_below(N, 5), I am wrong, because every multiple of 15 — every number divisible by both 3 and 5 — has been counted twice. The fix is the standard two-set inclusion–exclusion formula:

Here A is "multiples of 3 below N", B is "multiples of 5 below N", and A ∩ B is "multiples of lcm(3, 5) = 15 below N". The same identity applies to sums over those sets, not just counts, because a sum is just |·| weighted by the values being summed and the weights are the same in A, B, and A ∩ B (each element contributes its own value).

So:

For divisors that aren't pairwise coprime — say "multiples of 6 or 10" — the same formula works once you replace 15 with lcm(6, 10) = 30. That's why the implementation below uses math.lcm rather than just multiplying the divisors together; it is the substitution that makes the code correct in every case I might want to apply it to later.

For three or more divisors the pattern continues with alternating signs:

In general, summing over every non-empty subset T of the divisor set with sign (-1)^(|T| + 1) gives the right answer. That's what closed below implements.

A worked example with N = 10

Pick the example the problem itself gives, so we can cross-check at every step.

N = 10, divisors {3, 5}.

Step 1 — multiples of 3 below 10.

m = (10 - 1) // 3 = 9 // 3 = 3. The multiples are 3, 6, 9.

S_3(10) = 3 * 3 * 4 / 2 = 3 * 6 = 18. Direct sum: 3 + 6 + 9 = 18. ✓

Step 2 — multiples of 5 below 10.

m = (10 - 1) // 5 = 9 // 5 = 1. The only multiple is 5.

S_5(10) = 5 * 1 * 2 / 2 = 5. Direct sum: 5. ✓

Step 3 — multiples of 15 below 10.

m = (10 - 1) // 15 = 9 // 15 = 0. There are none.

S_15(10) = 15 * 0 * 1 / 2 = 0. ✓

Step 4 — combine.

S(10) = 18 + 5 - 0 = 23. That matches the value the problem itself states for the bound 10.

The same machinery applied to N = 1000 gives the answer Project Euler is asking for. I'm not printing the three component values here, because adding them is a single mental step and the brief asks not to publish the result. Run the script.

The implementation

Here's solver.py in full. It has the constant-time closed form, the brute force kept around as ground truth, an inclusion–exclusion driver that's general over a tuple of divisors, and a tiny benchmarking harness.

#!/usr/bin/env python3
"""
Project Euler #1 — sum of multiples of 3 or 5 strictly below N.

Two implementations:
  * brute(N)  — O(N) loop, used as ground truth.
  * closed(N) — O(1) inclusion–exclusion on the triangular-number identity.

The point of having both is to cross-check that the closed form is right
for small N, then trust it for huge N where brute force is infeasible.
"""
from __future__ import annotations
import time
from math import lcm


# --- O(N) reference --------------------------------------------------------
def brute(N: int) -> int:
    """Sum of k in [1, N) with k % 3 == 0 or k % 5 == 0."""
    return sum(k for k in range(1, N) if k % 3 == 0 or k % 5 == 0)


# --- O(1) closed form ------------------------------------------------------
def sum_multiples_below(N: int, k: int) -> int:
    """Sum of positive multiples of k strictly below N.

    The multiples are k, 2k, ..., m*k where m = (N - 1) // k.
    So the sum is k * (1 + 2 + ... + m) = k * m * (m + 1) // 2.
    """
    m = (N - 1) // k
    return k * m * (m + 1) // 2


def closed(N: int, divisors=(3, 5)) -> int:
    """Inclusion–exclusion over the given divisors.

    For pairwise inputs we'd need lcm; the stdlib `lcm` handles any pair
    so this works for {3, 5}, {3, 5, 7}, {6, 10}, etc.
    """
    from itertools import combinations
    total = 0
    for r in range(1, len(divisors) + 1):
        for subset in combinations(divisors, r):
            sign = -1 if r % 2 == 0 else 1
            total += sign * sum_multiples_below(N, lcm(*subset))
    return total


# --- timing helpers --------------------------------------------------------
def bench(fn, *args, repeats: int = 5):
    """Return (result, best_seconds_over_repeats)."""
    best = float("inf")
    result = None
    for _ in range(repeats):
        t0 = time.perf_counter_ns()
        result = fn(*args)
        dt = time.perf_counter_ns() - t0
        if dt < best:
            best = dt
    return result, best / 1e9


def fmt(seconds: float) -> str:
    if seconds < 1e-6:
        return f"{seconds * 1e9:7.1f} ns"
    if seconds < 1e-3:
        return f"{seconds * 1e6:7.1f} us"
    if seconds < 1.0:
        return f"{seconds * 1e3:7.1f} ms"
    return f"{seconds:7.3f}  s"


# --- driver ----------------------------------------------------------------
def main() -> None:
    N_official = 1000

    # Cross-check brute vs closed for the official bound.
    a, t_brute = bench(brute, N_official)
    b, t_closed = bench(closed, N_official)
    assert a == b, f"DISAGREE: brute={a} closed={b}"
    print(f"N = {N_official:>20,}    match: {a == b}")
    print(f"  brute   : {fmt(t_brute)}")
    print(f"  closed  : {fmt(t_closed)}")

    # Cross-check on a larger bound where brute is still tractable.
    for N in (10_000, 100_000, 1_000_000):
        a, t_brute = bench(brute, N)
        b, t_closed = bench(closed, N)
        assert a == b, f"DISAGREE at N={N}: brute={a} closed={b}"
        print(f"N = {N:>20,}    match: {a == b}")
        print(f"  brute   : {fmt(t_brute)}")
        print(f"  closed  : {fmt(t_closed)}")

    # Demonstrate that the closed form scales.  Brute force is skipped here.
    for N in (10**9, 10**12, 10**18, 10**100):
        _, t_closed = bench(closed, N)
        print(f"N = 10^{len(str(N)) - 1:<3}              closed  : {fmt(t_closed)}")


if __name__ == "__main__":
    main()

A few line-by-line notes, since some of the choices are deliberate:

• (N - 1) // k in sum_multiples_below is the correct formula for strictly below N. If the problem instead said "multiples below or equal to N", it would become N // k. Reading the problem statement carefully and choosing the right one is the entire bug surface of this function.
• k * m * (m + 1) // 2 is integer-safe in Python because exactly one of m and m + 1 is even, so the division by 2 is exact. I rely on // (floor division) rather than / (which would give a float and start losing precision at around 2^53).
• closed uses lcm(*subset) rather than math.prod(subset). They agree when divisors are coprime — lcm(3, 5) = 15 = 3 * 5 — but they diverge as soon as they share factors. lcm(6, 10) = 30, while 6 * 10 = 60. Using lcm is the version that survives generalisation.
• bench returns the best of repeated runs rather than the mean. For tiny operations (microseconds and below) the noise floor is dominated by other things on the box and "best of N" is a much better proxy for the cost of the operation itself than the mean.
• fmt switches units automatically. I do not want to read 0.00000354 s.

What the runtimes actually look like

Running python3 solver.py inside the sandbox (Python 3.13 on Kali, no JIT, no PyPy, no anything fancy) gives:

N =                1,000    match: True
  brute   :   118.0 us
  closed  :     4.4 us
N =               10,000    match: True
  brute   :     1.5 ms
  closed  :     4.0 us
N =              100,000    match: True
  brute   :    12.1 ms
  closed  :     3.5 us
N =            1,000,000    match: True
  brute   :   121.1 ms
  closed  :     3.5 us
N = 10^9                closed  :     3.3 us
N = 10^12               closed  :     3.8 us
N = 10^18               closed  :     3.8 us
N = 10^100              closed  :     5.8 us

A few things to notice:

The brute force scales linearly with N, exactly as advertised. From 1k to 1M it goes from 118 µs to 121 ms — a thousand-fold increase in N, a thousand-fold increase in time. Extrapolating: brute force on N = 10^9 would take roughly two minutes, which already busts the Project Euler one-minute rule. N = 10^12 is two days. N = 10^18 is six million years. The loop is not a viable strategy past a point.

The closed form sits at a flat ~4 microseconds across nine orders of magnitude of N. The reason it isn't strictly constant is bignum arithmetic: at N = 10^100 the intermediate values k * m * (m + 1) are around 200 digits long, and multiplying 200-digit integers takes a touch more than multiplying 3-digit integers. But the cost grows as O(d^2) (or O(d log d) for very large d where Python switches to Karatsuba) where d is the digit count, and d = log10(N) is microscopic compared to N. The runtime is "constant in everything that matters."

For the actual Project Euler bound — N = 1000 — both methods are well under the budget. The closed form is about 30× faster, but at this scale that's the difference between "instant" and "instant." The reason to prefer it is not speed at N = 1000; it's that it's the same algorithm whether N is a thousand or a googol, while the brute force is the same algorithm only in form.

A subtle point: bound exclusivity

The problem says "below 1000". I interpreted that as strictly below — k < 1000 — and so the largest candidate is 999. If the problem had said "up to and including 1000" or "below or equal to 1000", the answer would change: 1000 is itself a multiple of 5 (and not 3), and would add 1000 to S_5 but not to S_15. The closed form would shift by 1000 exactly.

Concretely, the right knob is the N - 1 versus N in m = (N - 1) // k. If I used N // k and called the function with N = 1000, I'd get the answer to "multiples of 3 or 5 up to and including 999". That happens to coincide with the strict-below answer at N = 1000 because 999 is divisible by 3, but at any bound N that is itself a multiple of 3 or 5, the two interpretations diverge. This is the kind of off-by-one that takes you fifteen minutes to find when your test on N = 10 passes (because 10 is a multiple of 5 only, so adding it shifts only S_5 and S_15 by 10 each, which... wait, that does change the answer; my point stands). Read the bound carefully and pick the right m.

What this generalises to

Project Euler #1 is the entry point to a family of problems that all yield to the same machinery. The reusable pieces:

1. Sum of an arithmetic progression — a single triangular number scaled by the common difference. This is the workhorse identity.
2. Inclusion–exclusion on divisor predicates — count or sum elements satisfying "any of these divisibility rules" by alternating-sign sum over subsets.
3. lcm as the right combinator — for pairs of divisors, the intersection of their multiples is exactly the multiples of their lcm. Not their product, unless they happen to be coprime.

Together they handle, for instance:

• Sum of multiples of any finite set of integers below N, whether coprime or not — closed(N, (a, b, c, ...)) already does this.
• Sum of integers below N that are coprime to a fixed m (Euler's totient sum, which is m * (φ(m) / m) * (N / m) * ... once you sieve by prime divisors of m; the same inclusion–exclusion engine, applied to the prime divisors of m).
• Counting (rather than summing) — replace sum_multiples_below with (N - 1) // k.

I won't push the post any further down that road. The point is that the closed form is not a one-trick: the trick is the recipe.

What I deliberately didn't do

A short list of tempting things I left out, with reasons.

A C version "for speed." The brute force in C would land somewhere around a microsecond at N = 1000. That's neat, but the closed form in Python is already about four microseconds and it works at N = 10^100, which the C version cannot do without GMP. The post is more honest written around the algorithm than around language micro-benchmarks; if the algorithm is O(1), language constants don't change the story.

A NumPy version. np.arange(0, N, 3).sum() + np.arange(0, N, 5).sum() - np.arange(0, N, 15).sum() would be a one-liner. It's fast at modest N, but it allocates O(N / k) memory and overflows int64 past about N = 10^9. The closed form is what NumPy is computing under the hood, in a roundabout, allocating, fixed-width way; better to just write the closed form.

Symbolic verification with SymPy. SymPy can derive the inclusion–exclusion formula in closed symbolic form and prove it equals the brute force for all N. That's lovely and unnecessary; the proof is two paragraphs of arithmetic above. I prefer the human proof.

A plot. I considered a log-log plot of brute vs closed runtime against N. The numbers in the runtime table tell the same story without a screenshot in the post. If you want it, the script's output is one matplotlib call away.

Reproducing this

Save solver.py from the listing above. Run:

$ python3 solver.py

It will cross-check brute against closed at N = 1000, 10_000, 100_000, and 1_000_000, time both at each, and then push the closed form out to N = 10^9, 10^12, 10^18, and 10^100 to show that it stays flat. The actual numerical answer to Project Euler #1 falls out of the first line of the cross-check; you can see it by changing match: {a == b} to value: {a}, or by importing closed and calling closed(1000) from a REPL.

I am not printing it here, because the rules of the room are that you do the work yourself and I do mine, and what I have to show is the path, not the destination.

References

• Problem statement: Project Euler problem 1, "Multiples of 3 or 5", projecteuler.net/problem=1. The four-element example for the bound 10 is from the official statement and is the only quoted fragment in this post.
• Inclusion–exclusion principle: any combinatorics text; the two- and three-set forms are all this post needs.
• Triangular numbers identity 1 + 2 + ... + m = m(m+1)/2: traditionally attributed to a young Gauss, debatably apocryphally; correct regardless.
• Python stdlib: math.lcm, math.gcd, itertools.combinations, time.perf_counter_ns. All Python ≥ 3.9.