Sum, square, subtract, slowly — then notice it's a polynomial

Project Euler problem #6 asks for an eight-digit number that two short loops would happily compute in microseconds. That's the boring answer. The interesting answer is the one-line polynomial it collapses to, why the polynomial is integer-valued, and what stops working when you ask the loop for n equal to a quintillion.

The problem, restated

Let

Project Euler problem #6 (projecteuler.net/problem=6) asks for

evaluated at $n = 100$. The challenge gives you $D(10) = 3025 - 385 = 2640$ as a warmup so you can verify a small case; the published task is just the same shape at a bigger $n$.

I'm not going to print $D(100)$ in this post — Project Euler asks solvers not to publish answers, and the point here is the derivation, not the digits. You'll run the code at the bottom and read the answer off your own terminal.

Why Python

This problem fits cleanly in three or four lines of any language, so "what language is fastest at integer arithmetic" is a non-question. The real choice was driven by what I wanted to be able to show:

1. Native arbitrary-precision integers. The naive sum at $n=100$ fits in a 32-bit int; at $n=10^{6}$ it doesn't; at $n=10^{18}$ it doesn't even fit in a 128-bit int. I want one closed-form function I can hand any $n$ without thinking about overflow, because the entire point of the closed form is that going from $n=100$ to $n=10^{18}$ should cost zero extra implementation effort. C's __int128 runs out around $n \approx 5\times 10^{9}$ for the squared sum-of-integers term $S(n)^{2}$ (and around $n \approx 10^{13}$ for $Q(n)$ on its own); Python never does.
2. fractions.Fraction for an independent algebraic check. I want to compute the closed form a second time using exact rationals, not relying on the integer-divisibility shortcut, and assert the result's denominator collapses to 1. That's a one-line existence proof that my factored formula has no hidden floor-division bug.
3. time.perf_counter_ns resolution. The closed form clocks in at ~500 ns. Anything coarser than nanoseconds and the comparison plot is flat for all interesting $n$.

If this were a heavy-arithmetic Project Euler problem (a 200-digit modular exponentiation, or a sieve to $10^{10}$), I'd reach for C with GMP or Rust with num-bigint. For a four-line closed-form polynomial problem whose interesting story is algebraic, Python is the right tool.

The obvious approach, and why it's obviously fine here

The most direct translation of the definition is:

def naive_diff(n: int) -> int:
    """Direct definition. No identities, just two loops collapsed into one."""
    s = 0          # running sum:   1 + 2 + ... + k
    sq = 0         # running sum of squares: 1 + 4 + ... + k**2
    for k in range(1, n + 1):
        s += k
        sq += k * k
    return s * s - sq

This is $O(n)$ in arithmetic operations, and at $n=100$ it returns in about ten microseconds. That's roughly $6\times 10^{6}$× faster than the Project Euler "one-minute rule" allows. We are done by the rules of the site.

So why keep going? Because the closed form is a much cleaner mathematical object than the loop, and Project Euler's whole pedagogy is to nudge you off the loop and onto the identity. The problem is rated 5% on Project Euler's difficulty scale specifically because it's the introduction to that pattern.

Closed forms: Gauss and Faulhaber

Two identities do all the work.

Gauss's identity (allegedly from his schooldays):

The standard proof is the "fold" trick: write $1 + 2 + \dots + n$ on one line and $n + (n-1) + \dots + 1$ underneath. Add column-wise; every column sums to $n+1$, and there are $n$ of them, so twice the sum is $n(n+1)$.

The sum-of-squares formula (known to Archimedes via On Spirals (Prop. 10) and On Conoids and Spheroids, roughly 1,850 years before Faulhaber tabulated such sums up to $p=17$, and well before Jakob Bernoulli (Ars Conjectandi, 1713) gave the general formula for $\sum_{k} k^{p}$ via the Bernoulli numbers):

There are at least four standard proofs (telescoping $k^{3} - (k-1)^{3}$; induction, which is rote once you guess the formula; finite differences, since $k^{2}$ is a degree-2 polynomial and the sum is therefore degree 3 with coefficients pinned by three samples; generating functions, via $\sum k^2 x^k = x(1+x)/(1-x)^{3}$ summed by an extra factor of $1/(1-x)$). The cleanest in my opinion is telescoping: $k^{3} - (k-1)^{3} = 3k^{2} - 3k + 1$, sum both sides from $k=1$ to $n$, the left telescopes to $n^{3}$, the right gives $3Q(n) - 3S(n) + n$, then solve for $Q(n)$ using Gauss.

Both formulas give integer outputs for integer $n$, which is the kind of fact I keep wanting to check until I've checked it: $n(n+1)$ is always even (consecutive integers), so the $/2$ is exact. $n(n+1)(2n+1)$ is always divisible by 6, because $n(n+1)$ covers the factor of 2 and one of $\{n, n+1, 2n+1\}$ covers the factor of 3 (case split on $n \bmod 3$).

Simplifying the difference

The lazy thing is to stop here. Both pieces are closed-form, so

is already $O(1)$. Done.

But this expression hides structure. Factor $n(n+1)$ out of both terms:

Put the bracket over a common denominator of 12:

The numerator $3n^{2} - n - 2$ factors. Discriminant: $1 + 24 = 25$, so roots at $n = (1 \pm 5)/6 = 1$ or $-2/3$. Therefore $3n^{2} - n - 2 = (n - 1)(3n + 2)$. Substituting:

That's a much prettier object. It's manifestly zero at $n=0$ and $n=1$ (so $D(1) = 0$, which agrees with the definition: $1^{2} - 1^{2} = 0$), and it grows like $\tfrac{1}{4}n^{4}$ for large $n$ (the leading term is $3n^{4}/12 = n^{4}/4$).

It also has a sanity check built in: the three consecutive integers $n-1, n, n+1$ in the numerator immediately tell you the product is divisible by both 2 and 3, so by 6. The remaining factor of 2 in the denominator wants one more even factor — which arrives via case analysis. If $n$ is odd, of two consecutive evens, exactly one is divisible by 4 and the other only by 2, so together they contribute 2·4 = 8 — already past what we need. If $n$ is even, $n$ itself is even and $3n+2$ is also even, so we get at least 4. Either way the product is divisible by 12. I confirmed this empirically up to $n=1000$:

$ python3 extra_checks.py
recurrence D(n)-D(n-1)=n^2(n-1) holds for n=2..500: OK
n*(n-1)*(n+1)*(3n+2) ≡ 0 (mod 12) for n=1..1000: OK
explicit vs factored at n=10**18: agree (72 digits)

That last line — agreeing at $n=10^{18}$ — is the one I care about most. If my factored form silently lost a factor somewhere, a 72-digit number would have caught it.

A nicer angle: Nicomachus

There's a side identity worth pointing out, because once you see it, problem 6 stops being a problem at all and becomes a question about the gap between two faces of the same object.

Nicomachus's theorem: the sum of the first $n$ cubes equals the square of the sum of the first $n$ integers.

The visual proof is the classic one: arrange unit squares into an L-shape for each $k$, and the L's tile the square of side $S(n)$. The algebraic check at small $n$:

n=  1: (sum k)^2 = sum k^3 =          1, D(n) = sum k^3 - sum k^2 = 0
n=  5: (sum k)^2 = sum k^3 =        225, D(n) = sum k^3 - sum k^2 = 170
n= 10: (sum k)^2 = sum k^3 =       3025, D(n) = sum k^3 - sum k^2 = 2640
n= 50: (sum k)^2 = sum k^3 =    1625625, D(n) = sum k^3 - sum k^2 = 1582700

The right-most column at $n=10$ is the problem-statement warmup of $2640$.

Once you accept Nicomachus, problem 6 reduces to

That last equality is dropping the trivially zero $k=1$ term and rewriting $k^{3} - k^{2} = k^{2}(k-1)$. Now $D(n)$ is itself a Faulhaber-style sum, just one with a $(k-1)$ shift. And it gives us our recurrence for free:

I checked this for $n = 2, \dots, 500$ (above), and the chain of identities works through it cleanly. Try it at $n=10$: $D(10) - D(9)$ should equal $100 \cdot 9 = 900$. My naive loop gives $D(10) = 2640$ and $D(9) = 1740$, and $2640 - 1740 = 900$. ✓

What the closed form actually looks like as a curve

Plotting $D(x) = x(x-1)(x+1)(3x+2)/12$ as a function of a continuous real $x$ shows the polynomial's shape directly. The integer $D(n)$ values are the integer-point samples of this curve.

The two roots at $x=0$ and $x=1$ are the only places the polynomial is zero on the non-negative reals, and that matches the fact that for $n \in \{0, 1\}$ the sum and sum-of-squares are equal (trivially when $n=0$, and at $n=1$ both are 1, so $S(1)^2 - Q(1) = 1 - 1 = 0$).

The growth rate of $n^{4}/4$ is what makes the loop version slower-than-necessary in an interesting way: the answer grows as $n^{4}$, but the loop's cost only grows as $n$. The bottleneck for huge $n$ is big-integer arithmetic in the running sums, whose word-length grows with $\log n$, not the iteration count alone.

The whole solver, top to bottom

"""Project Euler #6 — Sum Square Difference.

Compute  D(n) = (1 + 2 + ... + n)^2  -  (1^2 + 2^2 + ... + n^2)

Two implementations:
  1. naive_diff(n)       — O(n) loop, no math identities used.
  2. closed_form_diff(n) — O(1) using Faulhaber / Gauss closed forms.

We use Python because:
  * native arbitrary-precision int means the same formula works for n=100
    and n=10**18 without overflow noise distracting from the math;
  * time.perf_counter_ns gives nanosecond resolution, perfect for showing
    O(n) vs O(1) on a log axis;
  * fractions.Fraction lets us double-check the closed form algebraically
    against the loop for small n without any float fuzz.
"""

from __future__ import annotations
import time
from fractions import Fraction

def naive_diff(n: int) -> int:
    """Direct definition. No identities, just two loops collapsed into one."""
    s = 0          # running sum:   1 + 2 + ... + k
    sq = 0         # running sum of squares: 1 + 4 + ... + k**2
    for k in range(1, n + 1):
        s += k
        sq += k * k
    return s * s - sq

def closed_form_diff(n: int) -> int:
    """O(1) using:
         sum_{k=1..n} k       = n(n+1)/2          (Gauss)
         sum_{k=1..n} k^2     = n(n+1)(2n+1)/6    (Faulhaber, p=2)

       Difference simplifies to  n(n-1)(n+1)(3n+2)/12  — see post for the
       algebra.  We compute the *factored* form because every factor is
       integer-divisible at the right place, and Python's // never has
       to fall back to a Fraction.
    """
    # n(n-1)(n+1)(3n+2) is always divisible by 12 for n >= 1 — proof in post.
    return n * (n - 1) * (n + 1) * (3 * n + 2) // 12

def closed_form_diff_explicit(n: int) -> int:
    """Same as above but written in the un-factored form so a reader can
    map each piece back to its symbolic origin.  Useful as a sanity peer."""
    sum_n = n * (n + 1) // 2
    sum_sq_n = n * (n + 1) * (2 * n + 1) // 6
    return sum_n * sum_n - sum_sq_n

def verify_against_naive(upto: int = 200) -> None:
    """Cross-check both closed forms against the naive loop for n=1..upto.
    A single mismatch raises AssertionError; silent return means OK."""
    for n in range(1, upto + 1):
        a = naive_diff(n)
        b = closed_form_diff(n)
        c = closed_form_diff_explicit(n)
        assert a == b == c, (n, a, b, c)

def time_call(fn, n, repeats=5):
    """Return median wall-clock time in nanoseconds across `repeats` runs."""
    samples = []
    for _ in range(repeats):
        t0 = time.perf_counter_ns()
        fn(n)
        samples.append(time.perf_counter_ns() - t0)
    samples.sort()
    return samples[len(samples) // 2]

def fmt_ns(ns: int) -> str:
    if ns < 1_000:
        return f"{ns} ns"
    if ns < 1_000_000:
        return f"{ns/1_000:.2f} us"
    if ns < 1_000_000_000:
        return f"{ns/1_000_000:.2f} ms"
    return f"{ns/1_000_000_000:.3f} s"

def main():
    # 1.  Algebraic sanity: closed form matches naive for the first 200 n.
    verify_against_naive(200)
    print("verify_against_naive(200): OK")

    # 2.  The challenge target itself.  We DO NOT print the answer.
    n_target = 100
    target_naive = naive_diff(n_target)
    target_closed = closed_form_diff(n_target)
    assert target_naive == target_closed
    digits = len(str(target_naive))
    print(f"n=100:  naive == closed_form,  result has {digits} digits")

    # 3.  Runtime comparison across a sweep of n values.
    print()
    print(f"{'n':>16} | {'naive O(n)':>14} | {'closed O(1)':>12}")
    print("-" * 50)
    for n in [10, 100, 1_000, 10_000, 100_000, 1_000_000]:
        t_n = time_call(naive_diff, n)
        t_c = time_call(closed_form_diff, n)
        print(f"{n:>16} | {fmt_ns(t_n):>14} | {fmt_ns(t_c):>12}")

    # 4.  Show what closed-form unlocks: huge n.
    print()
    for n in [10**9, 10**12, 10**18]:
        t_c = time_call(closed_form_diff, n, repeats=11)
        digits = len(str(closed_form_diff(n)))
        print(f"closed_form_diff(10**{len(str(n))-1}):  {fmt_ns(t_c)}  "
              f"({digits}-digit answer)")

    # 5.  Independent check via Fraction (no integer-divisibility shortcuts).
    n = 100
    sum_n = Fraction(n * (n + 1), 2)
    sum_sq_n = Fraction(n * (n + 1) * (2 * n + 1), 6)
    diff = sum_n * sum_n - sum_sq_n
    assert diff.denominator == 1
    assert int(diff) == target_closed
    print()
    print("Fraction cross-check at n=100: OK (denominator collapses to 1)")

if __name__ == "__main__":
    main()

A few choices worth highlighting:

• Why // 12 instead of / 12? Because the proof above shows the numerator is divisible by 12 for every $n \ge 1$, floor division is exact and never silently loses precision. If I'd used /, Python would coerce to float and the result at $n=100$ (an eight-digit integer) would still be exact, but at $n=10^{18}$ (a 72-digit integer) it would catastrophically round to a float that has at most 16 significant digits. The wrong tool would be invisibly wrong.
• time.perf_counter_ns rather than time.time. The closed form runs in roughly 500 nanoseconds. time.time() has microsecond resolution at best, would round most calls to zero, and would also include monotonicity-skew on systems where the wall clock drifts. perf_counter_ns is monotonic and nanosecond-resolved.
• Median-of-five timing. A single timing call gets contaminated by GC pauses, JIT-like effects in the Python interpreter, or the kernel scheduling another process onto the core. The median across five repeats removes the worst outliers without inflating the budget.

Runtime, measured

This is what the sandbox prints when I run python3 sumsqdiff.py:

verify_against_naive(200): OK
n=100:  naive == closed_form,  result has 8 digits

               n |     naive O(n) |  closed O(1)
--------------------------------------------------
              10 |        1.02 us |       271 ns
             100 |        9.25 us |       378 ns
            1000 |      106.50 us |       400 ns
           10000 |        1.04 ms |       454 ns
          100000 |       12.04 ms |       591 ns
         1000000 |      134.67 ms |       536 ns

closed_form_diff(10**9):  487 ns  (36-digit answer)
closed_form_diff(10**12):  613 ns  (48-digit answer)
closed_form_diff(10**18):  658 ns  (72-digit answer)

Fraction cross-check at n=100: OK (denominator collapses to 1)

A few observations about that table:

1. The naive column scales linearly. Going from $n=10$ to $n=10^{6}$ is a 100,000× increase in $n$ and the runtime grows by about 130,000× — close enough, given that the per-iteration cost goes up slightly as the running sums turn into multi-word big integers around $n \sim 10^{4}$. The published target is $n=100$, which finishes in nine microseconds; we are 6,600,000× under the Project Euler one-minute rule.

2. The closed-form column is essentially flat at around 400–600 nanoseconds across six orders of magnitude of $n$. The slow drift upward (271 ns at $n=10$, 658 ns at $n=10^{18}$) is Python's big-integer multiplication going from one-word to four-word operands. It's not algorithmic; it's the cost of the answer being a 72-digit number.

3. At $n = 10^{18}$, the closed form gives a 72-digit answer in 658 ns. The equivalent naive loop would take roughly $10^{18}/10^{6} \cdot 135\text{ ms} \approx 1.35\times 10^{11}$ seconds, or about 4,300 years — and that's a lower bound, since the per-iteration big-integer cost grows further as the running sums climb past 16-word operands. The closed form isn't just faster; it's the only way you can ask the question.

Anchoring runtime in real Python timing is important here because Project Euler's one-minute rule is the only externally-imposed correctness criterion. The naive loop at $n=100$ takes 9 microseconds, six million times faster than required. The closed form takes 378 nanoseconds — for a problem this small, you cannot tell them apart by eye on any human screen. The whole point of the closed-form work is what it unlocks, not what it speeds up at the published $n$.

Three independent sanity checks

When the answer is a single number, "I ran my code and got X" is a thin evidence base. I want three of these, and each one has to test a different thing.

Check 1 — naive vs closed agree on small $n$. verify_against_naive(200) runs the loop and the two closed-form variants for every $n$ from 1 to 200 and asserts the three answers agree. If my algebra is wrong, this catches it.

Check 2 — the recurrence $D(n) - D(n-1) = n^{2}(n-1)$. A completely different derivation gives a relationship between consecutive $D$ values:

Use $S(n) - S(n-1) = n$ and $S(n) + S(n-1) = n(n+1)/2 + (n-1)n/2 = n^{2}$:

check_recurrence (defined in extra_checks.py, bundled separately — it's not in sumsqdiff.py's main() flow) walks $n$ from 2 to 500, computes $D(n) - D(n-1)$ via the naive loop, and asserts it equals $n^{2}(n-1)$:

def check_recurrence(upto=500):
    for n in range(2, upto + 1):
        delta = naive_diff(n) - naive_diff(n - 1)
        assert delta == n * n * (n - 1), (n, delta, n*n*(n-1))

If my naive loop has an off-by-one bug, this catches it.

Check 3 — Fraction denominator collapses. The block below evaluates the unfactored expression $S(n)^{2} - Q(n)$ with exact rationals — Fraction(n(n+1), 2) and Fraction(n(n+1)(2n+1), 6) — and asserts the denominator collapses to 1. That confirms the explicit closed form is integer-valued at $n=100$ without leaning on Python's //. The factored-/-12 derivation is then cross-checked by the next line, int(diff) == target_closed: if I'd miscounted and the 12 should have been a 24, the Fraction denominator would still be 1, but this equality against closed_form_diff(100) would fail.

n = 100
sum_n = Fraction(n * (n + 1), 2)
sum_sq_n = Fraction(n * (n + 1) * (2 * n + 1), 6)
diff = sum_n * sum_n - sum_sq_n
assert diff.denominator == 1

The denominator does collapse to 1 at $n=100$, confirming the explicit form is integer-valued there, and int(diff) matches closed_form_diff(100), which is what actually catches a wrong constant in the factored form.

If I had not derived the factored form, I would have skipped checks 2 and 3 and only had check 1. Check 1 alone is weak — the naive loop and the closed form share a common ancestor (my keyboard), so if I'd made a subtle off-by-one in both, the cross-check would silently pass. Checks 2 and 3 are valuable precisely because they test the closed form via mechanisms my naive loop wasn't involved in writing.

What I didn't bother with

A few alternatives I considered and rejected:

• Pre-computing the answer in NumPy. np.arange(1, 101).sum() ** 2 - (np.arange(1, 101)**2).sum() is one line and runs in microseconds. It's not faster than the pure-Python loop at $n=100$, and at large $n$ it silently overflows int64 and gives a wrong answer with no warning. Closed form wins on every axis.
• Memoizing across calls. functools.lru_cache would let naive_diff(100) return instantly on the second call. It does nothing for a single-shot script and complicates timing.
• C extension. A ctypes call into a hand-rolled C function would be faster on a per-call basis if I needed to call it billions of times. The closed form is already 500 ns; I do not need it to be 50 ns. The right answer at this scale is to use the better algorithm, not the faster language.
• Computing $D(n)$ symbolically with sympy. SymPy can simplify $S(n)^2 - Q(n)$ to the factored form mechanically. I considered including this as a fourth check, but sympy is a heavyweight dependency, and the algebra by hand fits in eight lines. It also has more pedagogical value when written out.

The closed form, point by point

For readers who want to see $D(n)$ before they compute it, the first few values lay out like this:

You can verify the recurrence in the $D(n)$ column: $4 - 0 = 4 = 2^{2}\cdot 1$, $22 - 4 = 18 = 3^{2}\cdot 2$, $70 - 22 = 48 = 4^{2}\cdot 3$, $170 - 70 = 100 = 5^{2}\cdot 4$. The differences ladder up as $n^{2}(n-1)$, exactly as derived.

The discrete points sit exactly on the continuous polynomial from the earlier plot — they have to, because the closed form is just the polynomial evaluated at integer $n$. This is the whole story of the simplification: a problem stated as a discrete combinatorial sum turned out to be the restriction of a degree-4 polynomial to the integers.

Running it yourself

The script is self-contained and depends only on the standard library. Save it as sumsqdiff.py and run:

$ python3 sumsqdiff.py

Inside the sandbox I'm running on, the wall-clock time for the whole script — including the verification loop over $n \in [1, 200]$, the timing sweep, the huge-$n$ tail, and the Fraction check — was under two seconds end-to-end. The actual answer for $n=100$ is in the last digit count line ("result has 8 digits"). To see the number itself, change one line near the bottom of main() from a digit-count print to a value print. Project Euler asks not to publish the digits, so I won't.

References

• Project Euler problem 6, Sum Square Difference: https://projecteuler.net/problem=6. Rated 5% on the site's difficulty scale; first published 2001-12-14.
• Faulhaber's formula for $\sum k^p$: J. Faulhaber, Academia Algebræ (1631); modern exposition in Knuth, Johann Faulhaber and Sums of Powers (1993), Math. Comp. 61, pp. 277-294. arXiv: https://arxiv.org/abs/math/9207222.
• Nicomachus's theorem ($\sum k^3 = (\sum k)^2$): standard combinatorial proof via L-shaped gnomons (squared triangular numbers); see the visual derivation at https://en.wikipedia.org/wiki/Squared_triangular_number.
• Python time.perf_counter_ns: https://docs.python.org/3/library/time.html#time.perf_counter_ns.
• Python fractions.Fraction: https://docs.python.org/3/library/fractions.html.

Artefacts

The solver (sumsqdiff.py) and the extra-checks script (extra_checks.py) are in the download bundle attached to this post. Only sumsqdiff.py is pasted in full above; extra_checks.py appears only as a small snippet plus its console output, and the full file lives in the bundle.