Project Euler 6: The Square of the Sum Minus the Sum of the Squares in O(1)

This is the sixth problem in the Project Euler series. Adding and squaring the numbers from 1 to 100 takes a single loop, so the real exercise here is not the arithmetic but deriving a formula that gives the answer without looping at all. We bring back the arithmetic series sum from problem 1 and add a formula for the sum of squares, derived from scratch with a telescoping sum. Both are tools that come back 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

We compare two quantities. One is the sum of 1 to n, squared; the other is the sum of the squares of 1 to n. The problem’s example is n = 10: the sum of the squares is 385, the square of the sum is 3025, and the difference is 2640. The task is to find the same difference for the first one hundred natural numbers, that is, n = 100.

There is no input; n is fixed at 100. The output is a single integer: (square of the sum) − (sum of the squares).

Background

Notation. This post names the two sums as follows.

S(n) = 1 + 2 + … + n             (the sum)
Q(n) = 1² + 2² + … + n²          (the sum of squares)
D(n) = S(n)² − Q(n)              (the difference we want)

For n = 0 there is nothing to add, so S(0) = Q(0) = 0 (an empty sum is 0).

What squaring a sum produces. Squaring a sum of several numbers gives the square of each number plus every product of two different numbers, twice. For example

(a + b + c)² = a² + b² + c² + 2(ab + ac + bc)

Expanding the brackets produces every way of picking one of a, b, c and then another one. Picking the same number twice gives a square, and two different numbers appear once in each order, as (a, b) and (b, a). We use this in the derivation to see what the difference actually is.

The arithmetic series by pairing. Write 1 + 2 + … + n forwards on one line and backwards on the line below. Every column adds up to n + 1, and there are n columns. So 2 × S(n) = n(n + 1), that is,

S(n) = n(n + 1) / 2

One of n and n + 1 is even, so n(n + 1) is always even and dividing by 2 gives an integer. This is the same formula used in problem 1 of this series.

Telescoping sums. Adding up the differences of neighbouring terms cancels everything in the middle and leaves only the ends.

(f(2) − f(1)) + (f(3) − f(2)) + … + (f(n + 1) − f(n)) = f(n + 1) − f(1)

Taking f(k) = k³ and expanding one term (k + 1)³ − k³ lets us pull out the sum of squares. That is the heart of step 4 below.

Brute force

Loop from 1 to n, accumulating both the sum and the sum of squares, then square the sum and subtract.

total = 0
squares = 0
for k in range(1, 101):
    total += k
    squares += k * k
difference = total * total - squares

To be honest, for n = 100 this is plenty: 100 iterations, O(n) time, O(1) extra space. The limits show up as n grows. For n like 10⁹ or 10¹² the loop grows with it, and in a language with fixed-width integers the sum of squares may overflow before the loop even ends. If instead we can write S(n) and Q(n) as formulas in n, any n takes just a few arithmetic operations. This post focuses on deriving those formulas and uses the brute force only to check them.

Derivation

Step 1: work out the problem’s example, n = 10, directly.

S(10) = 1 + 2 + … + 10 = 55,               S(10)² = 3025
Q(10) = 1 + 4 + 9 + 16 + 25 + 36 + 49 + 64 + 81 + 100 = 385
D(10) = 3025 − 385 = 2640

These match the three values given in the problem. We check every step below against these numbers by hand.

Step 2: see what the difference is. Applying the expansion from the background section to 1, 2, …, n gives

S(n)² = (1 + 2 + … + n)² = Q(n) + 2 × (sum of i × j over all pairs i < j)

so D(n) is twice the sum of the products of all pairs of different numbers. Checking a small case, n = 3: S(3)² = 36, Q(3) = 14, D(3) = 22, and 2 × (1 × 2 + 1 × 3 + 2 × 3) = 2 × 11 = 22. Two facts follow. Every product is positive, so D(n) ≥ 0, and for n = 1 (and n = 0) there are no pairs, so D = 0. The first fact is why the unsigned subtraction in C++ is safe.

Step 3: S(n) = n(n + 1) / 2. This is the pairing argument from the background section. For n = 10, 10 × 11 / 2 = 55, matching step 1.

Step 4: Q(n) = n(n + 1)(2n + 1) / 6. Expanding (k + 1)³ gives

(k + 1)³ − k³ = 3k² + 3k + 1

Add this up for k = 1 to n. The left side telescopes to (n + 1)³ − 1, and the right side is 3Q(n) + 3S(n) + n.

(n + 1)³ − 1 = 3Q(n) + 3S(n) + n

Checking by hand at n = 10: the left side is 11³ − 1 = 1330 and the right side is 3Q + 3 × 55 + 10 = 3Q + 175, so 3Q = 1155 and Q = 385, as in step 1.

Now for the general formula. Moving n to the left gives (n + 1)³ − 1 − n = (n + 1)³ − (n + 1) = (n + 1)((n + 1)² − 1) = (n + 1) × n × (n + 2). So

3Q(n) = n(n + 1)(n + 2) − 3 × n(n + 1) / 2
      = n(n + 1) × (2(n + 2) − 3) / 2
      = n(n + 1)(2n + 1) / 2
 Q(n) = n(n + 1)(2n + 1) / 6

For n = 10, 10 × 11 × 21 / 6 = 2310 / 6 = 385.

The division is exact. Since the formula equals the integer Q(n), n(n + 1)(2n + 1) must be a multiple of 6, but before trusting integer division in code it is worth checking directly. Divisibility by 2 holds because n(n + 1) is even. For 3, split on the remainder of n modulo 3: if it is 0, n is a multiple of 3; if 2, n + 1 is; if 1, then 2n + 1 = 2 × (3m + 1) + 1 = 3(2m + 1) is. Since 2 and 3 are coprime, the product is a multiple of 6.

Step 5: put them together.

D(n) = S(n)² − Q(n) = (n(n + 1) / 2)² − n(n + 1)(2n + 1) / 6

For n = 10, 55² − 385 = 3025 − 385 = 2640, matching the problem’s example. If you like, you can pull out the common factor and write it as one expression.

D(n) = n(n + 1) × (3n(n + 1) − 2(2n + 1)) / 12
     = n(n + 1)(3n² − n − 2) / 12
     = (n − 1) n (n + 1)(3n + 2) / 12

For n = 10, 9 × 10 × 11 × 32 / 12 = 31680 / 12 = 2640. The factored form has the advantage that the factor (n − 1) makes D = 0 at n = 1 obvious, but the code uses the S(n)² − Q(n) form because each step is easier to follow.

Final algorithm

1. Compute s = n(n + 1) / 2.
2. Compute q = n(n + 1)(2n + 1) / 6.
3. Return s × s − q.
4. Before printing the answer, check small inputs: Q(10) = 385, S(10)² = 3025, D(10) = 2640 and D(1) = 0, and for every n from 1 to 100 the closed form must equal the loop (the brute force). If any of these differ, exit with a nonzero status.

Correctness proof

• S(n) = n(n + 1) / 2 from step 3 holds for every n ≥ 0: the pairing argument does not depend on n, and for n = 0 both sides are 0.
• The identity (k + 1)³ − k³ = 3k² + 3k + 1 in step 4 holds by expansion, and summing it over k = 1, …, n and rearranging are all equalities, so Q(n) = n(n + 1)(2n + 1) / 6 holds for every n ≥ 1. For n = 0 both sides are 0.
• Both divisions are exact, as shown above, so integer division (/ in C++, // in Python) gives exactly the mathematical quotient.
• Therefore, as long as the values stay within the type’s range, the s × s − q the program returns is D(n) by definition. The range condition is covered in the next section.

Complexity

Time and space. Whatever n is, the work is a few multiplications, two divisions and one subtraction, so the time is O(1), and only a few values are stored, so the extra space is O(1) (while the values fit in a machine integer). The brute force was O(n) time.

Integer range.

• For n = 100 the values are n(n + 1) = 10100, S = 5050, n(n + 1)(2n + 1) = 2030100 before division, and Q = 338350; even the largest, S², is under 30 million. A 32-bit integer would be enough, but the C++ code uses unsigned long long (which the standard guarantees to be at least 64 bits wide) to leave room when generalizing.
• When generalizing, S² is the first to overflow. For S² ≤ 2⁶⁴ − 1 we need S ≤ 2³² − 1 = 4294967295. At n = 92681, n(n + 1) = 8589860442 and S = 4294930221, which is within range; at n = 92682, S = 4295022903 exceeds 2³², so S² overflows 64 bits. That is why the code comment says n <= 92681. The pre-division product n(n + 1)(2n + 1) is roughly 2n³ and fits in 64 bits until n is around two million, so within this range S² hits the limit first.
• The subtraction s × s − q cannot go below 0 even with unsigned integers, because step 2 showed D(n) ≥ 0.
• On overflow, unsigned integers wrap around modulo 2ᴺ (N being the type’s width in bits) instead of causing undefined behavior. The program does not stop but silently produces a wrong value, so n beyond the range above needs 128-bit or arbitrary-precision integers.
• Python’s int has unlimited precision, so it cannot overflow. But / produces a floating-point result, so we use integer division //. Since the divisions were shown to be exact, // gives the exact quotient.

C++17

#include <iostream>

// 1 + 2 + ... + n = n(n + 1) / 2. One of n and n + 1 is even, so the division is exact.
unsigned long long sum_up_to(unsigned long long n) {
    return n * (n + 1) / 2;
}

// 1^2 + 2^2 + ... + n^2 = n(n + 1)(2n + 1) / 6. The product is always divisible
// by 6, so the division is exact.
unsigned long long sum_of_squares(unsigned long long n) {
    return n * (n + 1) * (2 * n + 1) / 6;
}

// (1 + 2 + ... + n)^2 - (1^2 + 2^2 + ... + n^2). The difference is never negative,
// so the unsigned subtraction does not wrap. Meant for n <= 92681, where
// (n(n + 1) / 2)^2 still fits in 64 bits.
unsigned long long sum_square_difference(unsigned long long n) {
    const unsigned long long s = sum_up_to(n);
    return s * s - sum_of_squares(n);
}

// Direct loop over 1..n. Used only to cross-check the closed forms.
unsigned long long difference_by_loop(unsigned long long n) {
    unsigned long long total = 0;
    unsigned long long squares = 0;
    for (unsigned long long k = 1; k <= n; ++k) {
        total += k;
        squares += k * k;
    }
    return total * total - squares;
}

int main() {
    if (sum_of_squares(10) != 385 || sum_up_to(10) * sum_up_to(10) != 3025 ||
        sum_square_difference(10) != 2640 || sum_square_difference(1) != 0) {
        return 1;
    }
    for (unsigned long long n = 1; n <= 100; ++n) {
        if (sum_square_difference(n) != difference_by_loop(n)) {
            return 1;
        }
    }
    std::cout << sum_square_difference(100) << '\n';
    return 0;
}

difference_by_loop is the brute force as written above; it is used only to check that the closed form gives the same value for every n from 1 to 100. A mistake in the derivation of either formula would be caught here. 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.

Python 3

def sum_up_to(n: int) -> int:
    """1 + 2 + ... + n = n(n + 1) / 2. One of n and n + 1 is even, so the division is exact."""
    return n * (n + 1) // 2


def sum_of_squares(n: int) -> int:
    """1^2 + 2^2 + ... + n^2 = n(n + 1)(2n + 1) / 6.

    The product is always divisible by 6, so the division is exact.
    """
    return n * (n + 1) * (2 * n + 1) // 6


def sum_square_difference(n: int) -> int:
    """(1 + 2 + ... + n)^2 - (1^2 + 2^2 + ... + n^2), which is never negative."""
    s = sum_up_to(n)
    return s * s - sum_of_squares(n)


def difference_by_loop(n: int) -> int:
    """Direct loop over 1..n. Used only to cross-check the closed forms."""
    total = 0
    squares = 0
    for k in range(1, n + 1):
        total += k
        squares += k * k
    return total * total - squares


def main() -> None:
    if (
        sum_of_squares(10) != 385
        or sum_up_to(10) * sum_up_to(10) != 3025
        or sum_square_difference(10) != 2640
        or sum_square_difference(1) != 0
    ):
        raise SystemExit(1)
    for n in range(1, 101):
        if sum_square_difference(n) != difference_by_loop(n):
            raise SystemExit(1)
    print(sum_square_difference(100))


if __name__ == "__main__":
    main()

Same algorithm, same checks as the C++ version. The running sum is called total rather than sum so it does not shadow Python’s built-in sum, and the C++ code uses the same name to match. Because range(1, n + 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

For the first one hundred natural numbers, the difference between the square of the sum and the sum of the squares is 25164150.

The calculation can be followed by hand:

• The sum. S(100) = 100 × 101 / 2 = 5050, and S(100)² = 5050² = 25000000 + 2 × 5000 × 50 + 50² = 25502500.
• The sum of squares. Q(100) = 100 × 101 × 201 / 6 = 2030100 / 6 = 338350.
• The difference. 25502500 − 338350 = the answer above.
• Again with the factored form. 99 × 100 × 101 × 302 / 12 = 999900 × 302 / 12 = 301969800 / 12 = the answer above. Two differently shaped formulas give the same value.

The checks come in three layers.

1. Small examples and the full range. Before printing the answer, both programs check Q(10) = 385, S(10)² = 3025, D(10) = 2640 and D(1) = 0, and check that the closed form and the loop agree for every n from 1 to 100. If any differ, they exit with a nonzero status. The n = 10 case was also traced by hand at every step of 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. Both formulas and the exactness of the divisions were shown in the derivation and the correctness proof.

Boundary cases

• n = 1: S = 1 and Q = 1, so D = 0. There are no pairs of different numbers, which matches the reading in step 2. The small examples include D(1) = 0.
• n = 0: not part of the problem, but S = Q = 0 and D = 0, matching the convention that an empty sum is 0. The code gives the same value.
• n = 2: S² = 9, Q = 5, D = 4 = 2 × (1 × 2). The smallest nonzero case, with exactly one pair.
• Order of division: n(n + 1)(2n + 1) is fully multiplied before dividing by 6, so no fraction ever appears. Dividing first, say n / 6, would drop the remainder in integer division and give a wrong value.
• n out of range: as in “Integer range”, for n ≥ 92682 the S² in the C++ code overflows 64 bits. This program only handles n = 100, so it is fine, but take care when reusing the function for other n.

Going further

• The sum of cubes. Using the method of step 4, add up (k + 1)⁴ − k⁴ = 4k³ + 6k² + 4k + 1 to find 1³ + 2³ + … + n³. You should find that it equals S(n)², that is, “the square of the sum = the sum of the cubes”. What other difference, then, could D(n) in this problem be read as?
• Counting the difference directly. By step 2, D(n) = 2 × (sum of the products of pairs i < j). Group this sum by j and use S(j − 1) for the sum of the i below j, and the same formula comes out by another route. Try working it out.
• Higher powers. Changing the telescoping sum to (k + 1)ᵖ⁺¹ − kᵖ⁺¹ also gives 1ᵖ + … + nᵖ as a polynomial in n. Why is that polynomial of degree p + 1?
• The link to problem 1. Problem 1 of this series used the arithmetic series together with inclusion-exclusion. Compare how the same S(n) formula is reused here.

Sources and license

• Project Euler, Problem 6 — 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, 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 formulas for the sum and the sum of squares are derived directly above without relying on outside material. 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.