Three multiplications, no loop: Project Euler #1 in closed form

A walk-through of Project Euler problem #1 that skips the obvious one-line solution, derives the O(1) closed form from Gauss and inclusion–exclusion, and times all three implementations against each other so the cost of being clever is on the page rather than in folklore.

Why Python

Python is the wrong language if you want raw throughput and the right one when the interesting object is a number that might overflow a 64-bit register. The Project Euler "one-minute rule" is generous enough that a naive O(N) Python loop solves problem #1 in microseconds, so speed is not the choice driver here. What pushes me to Python is int: arbitrary precision by default, no u128, no BigInteger, no i65535 ceremony. The closed-form answer for N = 10**18 has 36 decimal digits, and I want to write n * (n + 1) // 2 without thinking about whether the intermediate n * (n + 1) overflows. timeit and perf_counter are in the standard library; the REPL is one keystroke away. For a problem whose interest is mathematical, not performance-engineering, the language should disappear, and Python disappears better than C, Rust, Zig, or Go for this shape of work.

The only language I seriously considered as an alternative was Lean or Coq, where I could state and prove the inclusion–exclusion identity. That would have been a different post.

The problem, restated

Project Euler problem #1 asks for a single integer: the sum of every positive integer strictly less than 1000 that is divisible by 3, by 5, or by both. The warm-up the problem ships with is the case N = 10, where the qualifying integers are 3, 5, 6, 9 and the sum is 23. The published task is the same question with N = 1000. (Per Project Euler's request, I'm not going to put the numeric answer for N = 1000 in this post; the code below reproduces it on any machine in well under a millisecond.)

A few things deserve to be pinned down before writing any code:

• Strictly below, not "up to and including." For N = 10 the multiple 10 is not counted. This bites you the first time you write range(1, N + 1) instead of range(1, N).
• Multiples of 3 or 5, inclusive-or. 15 qualifies because it is a multiple of 3 (and also a multiple of 5). The trap is in the sum: if you add sum-of-multiples-of-3 and sum-of-multiples-of-5 you have counted 15, 30, 45, … twice.
• Positive integers. 0 is technically a multiple of every integer; whether you include it is irrelevant because adding zero changes nothing.

That last subtlety is small but worth saying out loud, because once you start writing closed forms you stop iterating, and any wrong-by-one error becomes a wrong-by-O(N) error in the result.

What "obvious" looks like

The two obvious-and-slow solutions are both worth writing, because they're what we will measure against:

def naive_loop(N):
    total = 0
    for k in range(N):
        if k % 3 == 0 or k % 5 == 0:
            total += k
    return total

This is O(N) and does the modulo work twice (sometimes once, with short-circuiting) per integer. It is the version I'd write first if the question were ever in front of me on a whiteboard. For N = 1000 it is plenty.

A small step up:

def step_sum(N):
    s3  = sum(range(3,  N, 3))
    s5  = sum(range(5,  N, 5))
    s15 = sum(range(15, N, 15))
    return s3 + s5 - s15

Three passes, but each one steps over only the integers that matter. The s15 subtraction is already inclusion–exclusion in disguise: it removes the multiples of 15 that we counted in both s3 and s5. Same O(N), smaller constant — about an order of magnitude faster in practice.

The interesting solution is the one where we never iterate at all.

Where Gauss comes in

Every multiple of d strictly less than N can be written as d, 2d, 3d, …, kd where k = ⌊(N − 1) / d⌋. The sum of those multiples is

d + 2d + 3d + … + kd  =  d · (1 + 2 + 3 + … + k)

and the inner sum is the k-th triangular number. Gauss's schoolboy trick — pair the first and last term, the second and second-to-last, each pair sums to k + 1, there are k/2 pairs (or extend to the standard induction proof for arbitrary k) — gives

1 + 2 + 3 + … + k  =  k(k + 1) / 2

So the sum of positive multiples of d strictly less than N is

S(d, N) = d · k · (k + 1) / 2,    where k = ⌊(N − 1) / d⌋

This is O(1). One floor-division, one multiplication of moderately-sized integers, one shift. No loop, no list, no allocation.

For PE #1 we need the sum over multiples-of-3 or multiples-of-5. If the two sets were disjoint we could just add S(3, N) + S(5, N). They are not disjoint — they meet at the multiples of lcm(3, 5) = 15 — so we have to subtract the intersection once to avoid double-counting:

answer = S(3, N) + S(5, N) − S(15, N)

This is the two-set case of the inclusion–exclusion principle.

Inclusion–exclusion, generalised

Inclusion–exclusion is one of the few pieces of combinatorics that earns its keep in everyday programming, and it's worth stating in full because PE #1 is the smallest non-trivial instance.

If A_1, A_2, …, A_n are finite sets, then

|A_1 ∪ A_2 ∪ … ∪ A_n|
    = Σ |A_i|
    − Σ |A_i ∩ A_j|        (i < j)
    + Σ |A_i ∩ A_j ∩ A_k|  (i < j < k)
    − …
    + (−1)^(n+1) |A_1 ∩ A_2 ∩ … ∩ A_n|

The pattern: alternating signs, with a + on odd-cardinality intersections and a − on even-cardinality ones. The same identity holds when you replace cardinality with any function that is additive over disjoint unions — and summing the elements is exactly such a function. So if A_d = {k ∈ [1, N) : d | k} and we want Σ {k : k ∈ A_3 ∪ A_5}, we get

sum(A_3 ∪ A_5) = sum(A_3) + sum(A_5) − sum(A_3 ∩ A_5)

and the intersection of "multiples of 3" with "multiples of 5" is exactly the multiples of lcm(3, 5) = 15. That's the key bridge: intersection of divisibility constraints is divisibility by the LCM. Two integers d and e both divide k iff lcm(d, e) divides k.

For more divisors the formula has more terms but the structure is identical:

sum_union({d_1, …, d_n}, N)
    = Σ_{S ⊆ D, S ≠ ∅}  (−1)^(|S|+1) · S(lcm(S), N)

That is 2^n − 1 terms — exponential in the number of divisors but completely independent of N. For PE #1, n = 2, so we have 2² − 1 = 3 terms, exactly the three Gauss sums we wrote down above.

The implementation

def triangular(n):
    """1 + 2 + … + n  via Gauss's formula."""
    return n * (n + 1) // 2


def sum_multiples_below(d, N):
    """Sum of positive multiples of d strictly less than N."""
    k = (N - 1) // d
    return d * triangular(k)


def closed_form(N):
    return (
        sum_multiples_below(3, N)
        + sum_multiples_below(5, N)
        - sum_multiples_below(15, N)
    )

Three things worth saying about this code:

1. Floor division (N - 1) // d. I want k = ⌊(N − 1) / d⌋, the number of positive multiples of d that are strictly less than N. The N − 1 accounts for the strict inequality: if N were itself a multiple of d, we should not count it, which is exactly what (N − 1) // d gives us. Using N // d would be correct for "≤ N" and wrong for "< N", and the test case at N = 6 makes the bug visible: multiples of 3 below 6 are {3} (one element); (6 − 1) // 3 = 1 but 6 // 3 = 2.
2. Integer division everywhere. n * (n + 1) is always even — one of n and n + 1 is even — so n * (n + 1) // 2 is exact. No floats involved at any point. Python's int is arbitrary precision, so the intermediate value can be as large as we like.
3. The lcm(3, 5) = 15 is hard-coded. That's fine for PE #1 because the divisors are baked into the question. The general version below pulls the LCM out into math.lcm.

What it costs to run

I wrote a small driver that times the three implementations on the same N = 1000, then runs only the closed form on far larger N to show the time stays flat. This is the entire solver.py:

from time import perf_counter


def naive_loop(N):
    total = 0
    for k in range(N):
        if k % 3 == 0 or k % 5 == 0:
            total += k
    return total


def step_sum(N):
    s3  = sum(range(3,  N, 3))
    s5  = sum(range(5,  N, 5))
    s15 = sum(range(15, N, 15))
    return s3 + s5 - s15


def triangular(n):
    return n * (n + 1) // 2


def sum_multiples_below(d, N):
    return d * triangular((N - 1) // d)


def closed_form(N):
    return (
        sum_multiples_below(3, N)
        + sum_multiples_below(5, N)
        - sum_multiples_below(15, N)
    )


def time_it(fn, N, repeats=1):
    t0 = perf_counter()
    for _ in range(repeats):
        ans = fn(N)
    return ans, (perf_counter() - t0) / repeats


if __name__ == "__main__":
    N = 1000
    for name, fn, reps in [
        ("naive_loop",  naive_loop,  10_000),
        ("step_sum",    step_sum,    10_000),
        ("closed_form", closed_form, 1_000_000),
    ]:
        ans, sec = time_it(fn, N, reps)
        print(f"{name:12s}  mean={sec*1e6:9.3f} us  ({reps} reps)")

    print("All three agree:",
          naive_loop(N) == step_sum(N) == closed_form(N))

    for N_big in [10**6, 10**9, 10**12, 10**15, 10**18]:
        ans, sec = time_it(closed_form, N_big, 100_000)
        print(f"  N=10^{len(str(N_big))-1:2d}   "
              f"answer has {len(str(ans)):2d} digits   "
              f"mean={sec*1e6:6.3f} us")

Running it inside the sandbox (Python 3.13.12, x86_64 Linux, no warm-up tricks beyond perf_counter averaging):

N = 1000
------------------------------------------------------------
naive_loop    mean=  116.599 us  (10000 reps)
step_sum      mean=    9.756 us  (10000 reps)
closed_form   mean=    0.916 us  (1000000 reps)

All three return the same value: True

Closed-form on huge N (only this one survives):
  N=10^ 6   answer has 12 digits   mean= 1.038 us
  N=10^ 9   answer has 18 digits   mean= 1.036 us
  N=10^12   answer has 24 digits   mean= 1.364 us
  N=10^15   answer has 30 digits   mean= 1.403 us
  N=10^18   answer has 36 digits   mean= 1.403 us

Three observations:

• For the published N = 1000, every implementation is well inside the one-minute rule. The naive loop is about 100 microseconds; the closed form is under one microsecond. Even the slow version solves the problem 600,000× faster than the rule allows.
• The closed form is roughly 100× faster than the naive loop and 10× faster than step_sum. That ratio is bigger than it sounds: at N = 1000 the work the loop does is small, so most of those 116 microseconds is interpreter overhead per iteration, which the closed form skips entirely.
• As N grows by twelve orders of magnitude — from 10^6 to 10^18 — the closed-form runtime grows by less than 50%. The growth that does exist is the cost of multiplying ever-bigger Python ints; the algorithm itself is O(1), so the only thing that scales is bigint arithmetic, and for 36-digit numbers it scales gently.

The naive loop and step_sum do not appear in the second table because at N = 10^18 they would not finish in this universe.

Sanity-checking the closed form

Two cheap checks to make sure I'm not lying to myself.

Check 1: the warm-up. The problem hands you a known case: multiples of 3 or 5 below 10 are 3, 5, 6, 9, summing to 23. Plugging N = 10 into the closed form:

• S(3, 10): k = ⌊9/3⌋ = 3, sum = 3 · 3 · 4 / 2 = 18. (That's 3 + 6 + 9.)
• S(5, 10): k = ⌊9/5⌋ = 1, sum = 5 · 1 · 2 / 2 = 5. (That's 5.)
• S(15, 10): k = ⌊9/15⌋ = 0, sum = 15 · 0 · 1 / 2 = 0. (No multiples of 15 below 10.)

Total: 18 + 5 − 0 = 23. ✓

Check 2: density argument. Asymptotically, the fraction of integers below N that are multiples of 3 or 5 is

1/3 + 1/5 − 1/15 = 5/15 + 3/15 − 1/15 = 7/15

The mean value of an integer in [1, N) is about N/2, and there are about N of them, so the sum should be roughly (7/15) · (N/2) · N = 7 N² / 30 ≈ 0.2333 · N². Plotting actual against approximation:

N=     100  approx=          2333.3  ratio=0.99343
N=    1000  approx=        233333.3  ratio=0.99929
N=   10000  approx=      23333333.3  ratio=0.99993
N=  100000  approx=    2333333333.3  ratio=0.99999
N= 1000000  approx=  233333333333.3  ratio=0.99999

The ratio of actual to asymptotic approaches 1 from below as N grows, exactly what you'd expect when the floor function is rounding off a small fraction of each progression. (I'm omitting the actual column from this table on purpose — the N = 1000 row is the published answer to PE #1.)

Generalising: any set of divisors

The closed form has nothing to do with the specific divisors 3 and 5. Pull them out as a parameter and the same pattern handles arbitrary divisor sets:

from itertools import combinations
from math import lcm
from functools import reduce


def triangular(n):
    return n * (n + 1) // 2


def sum_multiples_below(d, N):
    return d * triangular((N - 1) // d)


def sum_union(D, N):
    """Sum of integers k in [1, N) divisible by at least one d in D."""
    D = list(D)
    total = 0
    for size in range(1, len(D) + 1):
        sign = -1 if size % 2 == 0 else 1
        for subset in combinations(D, size):
            l = reduce(lcm, subset)
            total += sign * sum_multiples_below(l, N)
    return total

The structure mirrors the inclusion–exclusion identity exactly: iterate over non-empty subsets of D; for each subset, take the LCM, compute the closed-form sum of multiples of that LCM, and add it with sign (−1)^(|S|+1).

Two cross-checks confirm it agrees with brute force where brute force is feasible:

D={3,5}, N=1000              match: True
D=[2, 3, 5, 7, 11], N=100000 match: True  digits: 10
D=[2, 3, 5, 7, 11], N=10^18  digits of answer: 36

The brute force for D = {2, 3, 5, 7, 11}, N = 100000 is 100000 modulo checks per element, fine on a laptop. At N = 10^18 brute force is gone but sum_union still returns instantly because the work is O(2^|D|) regardless of N.

A subtlety worth flagging: this generalisation has 2^n − 1 terms, which is fine for the small divisor sets that appear in early Project Euler problems but blows up if you hand it 30 divisors. For n = 30 you have over a billion subsets — at which point a Mobius-function approach over the lattice of divisors, or simply giving up and using a sieve, is the better engineering call. PE #1 doesn't push you anywhere near that threshold; PE #1 is the case n = 2.

The two endpoints I keep getting wrong

Every time I revisit this kind of problem after a few years away, I rederive the same off-by-one error and have to swat it down. So, for the future me reading this: the question asks for multiples strictly below N. The number of multiples of d in the open interval [1, N) is

k = ⌊(N − 1) / d⌋

The number of multiples of d in the closed interval [1, N] is

k = ⌊N / d⌋

These differ exactly when d divides N. For PE #1 with d ∈ {3, 5, 15} and N = 1000, the divisor 5 divides 1000 and so does the divisor … no, neither 3 nor 15 divide 1000. But 5 does. So if you accidentally compute the inclusive sum, you over-count by 1000 because you've added 5 · 200 = 1000 (the multiple of 5 at the boundary) and then subtracted a non-existent multiple of 15. The result is wrong by 1000, and the error is the kind that survives a casual cross-check because both wrong and right answers are plausible six-digit numbers. The cure is to write (N - 1) // d and to test against the warm-up case N = 10.

Why the closed form is "right"

A subjective claim, but I'll stand on it. For PE #1 specifically, the closed form is overkill — the naive loop is fast enough, easy to read, and harder to get wrong. If a junior engineer wrote naive_loop in a PR I would not ask them to change it.

But "right" in the Project Euler sense is not just "produces the answer." It's "produces the answer in a way that scales." Project Euler is a series, and problem #1 is a 5-percent warm-up; the later problems are not generous. The habit you want to build on the warm-ups is the one that survives:

• Naive loop: dies at N ≈ 10^9 because Python iterates one billion times in tens of seconds.
• step_sum (range-with-step): dies at N ≈ 10^10 for the same reason, with a smaller constant.
• Closed form: alive at N = 10^18, alive at N = 10^100 (Python int is happy to keep multiplying), and the change of mental model that buys you that survival is the same change of mental model that lets you tackle PE #401 (sum of squares of divisors), PE #436 (Lah numbers), and a hundred others where direct iteration is hopeless from the first character of the problem.

Spending an extra ten minutes on a five-line problem to derive the closed form is a tax you pay once and harvest interest on for years. The five-line solution is correct; the three-multiplication solution is the one that teaches you anything.

Edge cases I did not implement

A few I considered and decided were out of scope for the post but worth listing so the reader can fill them in:

• N ≤ 0. My closed_form returns 0 for N = 1 (correct: empty range) and a negative number for N = 0 because (0 - 1) // d = -1. The cleanest fix is to clamp: k = max(0, (N - 1) // d). I left it out because the problem promises N = 1000.
• Negative divisors. sum_multiples_below(-3, N) is meaningless as written. The natural extension is "multiples of |d|", which the brute-force loop handles transparently and the closed form would need a d = abs(d) line for.
• Overflow. Python doesn't have any. C, Rust, and Zig solutions of this problem need to think about it: at N = 10^18, the intermediate k * (k + 1) is on the order of 10^36, which needs 128-bit integers in C and a BigInt in Rust. This is the place where the language choice would actually matter.
• Modular variants. "Sum of multiples of 3 or 5 below N, modulo 10^9 + 7" appears in competitive programming. The closed form survives intact — every +, -, *, // becomes its modular counterpart, and the // 2 becomes "multiply by the modular inverse of 2."

References and reproduction

• Project Euler problem #1: https://projecteuler.net/problem=1. Account-walled to see the answer-checker; problem text is public.
• Inclusion–exclusion principle: standard combinatorics, but Concrete Mathematics (Graham, Knuth, Patashnik) chapter 5 has the cleanest treatment I know.
• Gauss's triangular-number formula: famously rediscovered by him at age 9 or so. The induction proof is a few lines and worth working out by hand once.
• All code in this post is solver.py and general.py in the artefacts tarball; both are also reproduced in full in the sections above. They run on stock Python 3.10+ (I tested 3.13). No third-party dependencies.

To reproduce: copy solver.py from this post into a file, run python3 solver.py, and the last line of the table is the answer to PE #1. The answer is left for the reader to compute, on Project Euler's request.

[exit 0]
<div class="problem_content" role="problem">
<p>If we list all the natural numbers below $10$ that are multiples of $3$ or $5$, we get $3, 5, 6$ and $9$. The sum of these multiples is $23$.</p>
<p>Find the sum of all the multiples of $3$ or $5$ below $1000$.</p>
</div><br>
<br><script src="js/info_copy.1776805016.js"></script>
<div id="snackbar">Copied to Clipboard</div>
</div> <!--end_content-->

</div> <!--end_container-->

<div id="footer" class="noprint">
Project Euler: <a href="copyright">Copyright Information</a> | <a href="privacy">Privacy Policy</a>
</div>

<div id="modal_window">
   <div id="modal_content" class="message_body">
   <p>The page has been left unattended for too long and that link/button is no longer active. Please refresh the page.</p>
   </div>
</div>

<script src="js/mathjax_config.1766083028.js"></script>
<script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@4/tex-mml-svg.js?config=newcm"></script><script type="application/json" id="page_expire">86400</script>
<script src="js/general.1681854148.js"></script>

</body>
</html>

[exit 0]
N = 1000
------------------------------------------------------------
naive_loop    ans=<redacted>  mean=  116.599 us  (10000 reps)
step_sum      ans=<redacted>  mean=    9.756 us  (10000 reps)
closed_form   ans=<redacted>  mean=    0.916 us  (1000000 reps)

All three return the same value: True

Closed-form on huge N (only this one survives):
  N=10^ 6   answer has 12 digits   mean= 1.038 us
  N=10^ 9   answer has 18 digits   mean= 1.036 us
  N=10^12   answer has 24 digits   mean= 1.364 us
  N=10^15   answer has 30 digits   mean= 1.403 us
  N=10^18   answer has 36 digits   mean= 1.403 us

[exit 0]
N=10 -> 23 23 23

[exit 0]
D={3,5}, N=1000  match: True
D=[2, 3, 5, 7, 11], N=100000  match: True  digits: 10
D=[2, 3, 5, 7, 11], N=10^18  digits of answer: 36

[exit 0]
Python 3.13.12
Linux 24e369ce4c89 6.8.0-90-generic #91-Ubuntu SMP PREEMPT_DYNAMIC Tue Nov 18 14:14:30 UTC 2025 x86_64 GNU/Linux

[exit 0]
N=     100  actual=            2318  approx=            2333.3  ratio=0.993429
N=    1000  actual=          233168  approx=          233333.3  ratio=0.999291
N=   10000  actual=        23331668  approx=        23333333.3  ratio=0.999929
N=  100000  actual=      2333316668  approx=      2333333333.3  ratio=0.999993
N= 1000000  actual=    233333166668  approx=    233333333333.3  ratio=0.999999