· 개발

Project Euler 6: 합의 제곱과 제곱의 합, 그 차이를 닫힌 식으로 O(1)에 구하기

1부터 100까지의 합을 제곱한 값과 제곱들의 합의 차이를, 반복문 대신 n(n + 1) / 2와 n(n + 1)(2n + 1) / 6 두 공식으로 바로 계산합니다. 두 공식을 짝짓기와 망원합으로 직접 유도하고, 정수 범위와 C++17·Python 3 구현까지 다룹니다.

Project Euler알고리즘수학급수

스포일러 주의 · 이 글에는 풀이 과정과 정답이 그대로 나옵니다. 직접 먼저 풀어 보고 싶다면 Project Euler 6번 문제부터 보세요.

Project Euler 시리즈의 여섯 번째 문제입니다. 1부터 100까지 더하고 곱하는 일은 반복문 하나면 끝나므로, 이 문제는 계산 자체보다 반복 없이 바로 답을 내는 식을 끌어내는 연습에 가깝습니다. 1번 문제에서 쓴 등차수열의 합을 다시 꺼내고, 여기에 제곱의 합 공식을 망원합(telescoping sum)으로 직접 유도해 더합니다. 둘 다 뒤의 문제에서 다시 쓰는 도구입니다.

스포일러 주의: 이 글에는 풀이 전체와 정답이 들어 있습니다. 직접 먼저 풀어 보고 싶다면 원래 문제부터 보세요. 이 시리즈의 글에는 광고를 싣지 않습니다. 문제 본문의 라이선스(CC BY-NC-SA 4.0)가 비영리 조건을 두기 때문입니다.

문제

두 값을 비교합니다. 하나는 1부터 n까지의 합을 구한 뒤 제곱한 값이고, 다른 하나는 1부터 n까지 각각을 제곱한 뒤 더한 값입니다. 문제에 나온 예로 n = 10이면 제곱의 합은 385, 합의 제곱은 3025이고 그 차이는 2640입니다. 같은 차이를 처음 100개의 자연수, 즉 n = 100에 대해 구하는 것이 문제입니다.

따로 받는 입력은 없고 n은 100으로 정해져 있습니다. 구할 것은 (합의 제곱) − (제곱의 합)인 정수 하나입니다.

필요한 배경 지식

합의 기호. 이 글에서는 두 합을 다음처럼 부릅니다.

S(n) = 1 + 2 + … + n             (합)
Q(n) = 1² + 2² + … + n²          (제곱의 합)
D(n) = S(n)² − Q(n)              (구하는 차이)

n = 0이면 더할 것이 없으므로 S(0) = Q(0) = 0으로 둡니다(빈 합은 0).

합을 제곱하면 생기는 것. 여러 수의 합을 제곱하면 각 수의 제곱과, 서로 다른 두 수의 곱이 두 번씩 나옵니다. 예를 들어

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

입니다. 괄호를 풀면 a, b, c 중 하나를 고르고 또 하나를 고르는 모든 경우가 한 번씩 나오는데, 같은 수를 두 번 고르면 제곱이고, 서로 다른 두 수는 (a, b)와 (b, a)처럼 두 순서로 한 번씩 나오기 때문입니다. 이 사실은 “풀이 유도”에서 차이가 무엇인지 읽는 데 씁니다.

짝짓기로 구하는 등차수열의 합. 1 + 2 + … + n을 앞에서부터 쓴 줄과 뒤에서부터 쓴 줄을 위아래로 놓으면 각 열의 합이 모두 n + 1이고, 열이 n개 있습니다. 그래서 2 × S(n) = n(n + 1), 곧

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

입니다. n과 n + 1 중 하나는 짝수이므로 n(n + 1)은 언제나 짝수이고, 2로 나눈 결과는 정수입니다. 이 시리즈 1번 문제에서 쓴 것과 같은 공식입니다.

망원합. 이웃한 항의 차를 더하면 중간 항이 모두 지워지고 처음과 끝만 남습니다.

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

f(k) = k³으로 두고 한 항 (k + 1)³ − k³을 전개해 두면 제곱의 합을 끌어낼 수 있습니다. 이것이 아래 4단계의 핵심입니다.

가장 단순한 방법

1부터 n까지 돌면서 합과 제곱의 합을 함께 쌓고, 끝에서 합을 제곱해 빼면 됩니다.

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

솔직히 말하면 n = 100에서는 이것으로 충분합니다. 반복은 100번, 시간은 O(n), 추가 공간은 O(1)입니다. 한계는 n이 커질 때 드러납니다. n이 10⁹, 10¹²처럼 커지면 반복 횟수도 그만큼 늘어나고, 고정 폭 정수를 쓰는 언어에서는 제곱의 합이 그 전에 넘칠 수 있습니다. 반대로 S(n)과 Q(n)을 n에 대한 식으로 쓸 수 있으면 n이 얼마든 산술 연산 몇 번으로 끝납니다. 이 글은 그 식을 끌어내는 데 집중하고, 가장 단순한 방법은 식이 맞는지 대조하는 데 씁니다.

풀이 유도

1단계: 문제의 예 n = 10을 직접 계산합니다.

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

문제에 주어진 세 값과 같습니다. 아래 단계마다 이 숫자들로 손 검산을 합니다.

2단계: 차이가 무엇인지 읽습니다. 배경 지식의 전개를 1, 2, …, n에 쓰면

S(n)² = (1 + 2 + … + n)² = Q(n) + 2 × (i < j인 모든 쌍의 곱 i × j의 합)

이므로 D(n)은 서로 다른 두 수의 곱을 모두 더한 것의 두 배입니다. 작은 예 n = 3으로 확인하면 S(3)² = 36, Q(3) = 14, D(3) = 22이고, 2 × (1 × 2 + 1 × 3 + 2 × 3) = 2 × 11 = 22로 맞습니다. 여기서 두 가지를 얻습니다. 곱은 모두 양수이므로 D(n) ≥ 0이고, 쌍이 하나도 없는 n = 1(그리고 n = 0)에서는 D = 0입니다. 앞의 사실은 C++에서 부호 없는 뺄셈이 안전하다는 근거가 됩니다.

3단계: S(n) = n(n + 1) / 2. 배경 지식의 짝짓기 그대로입니다. n = 10이면 10 × 11 / 2 = 55이고, 1단계와 같습니다.

4단계: Q(n) = n(n + 1)(2n + 1) / 6. (k + 1)³을 전개하면

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

입니다. 이것을 k = 1부터 n까지 더합니다. 왼쪽은 망원합이라 (n + 1)³ − 1이 남고, 오른쪽은 3Q(n) + 3S(n) + n입니다.

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

n = 10에서 손으로 확인해 보면 왼쪽은 11³ − 1 = 1330, 오른쪽은 3Q + 3 × 55 + 10 = 3Q + 175이므로 3Q = 1155, Q = 385입니다. 1단계와 같습니다.

이제 일반식으로 정리합니다. 왼쪽에서 n을 옮기면 (n + 1)³ − 1 − n = (n + 1)³ − (n + 1) = (n + 1)((n + 1)² − 1) = (n + 1) × n × (n + 2)입니다. 그래서

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

입니다. n = 10이면 10 × 11 × 21 / 6 = 2310 / 6 = 385입니다.

나눗셈이 정확하다는 것. 위 식은 정수 Q(n)과 같으므로 n(n + 1)(2n + 1)은 반드시 6의 배수이지만, 코드에서 정수 나눗셈을 믿으려면 직접 확인해 두는 편이 좋습니다. 2로 나누어떨어지는 것은 n(n + 1)이 짝수이기 때문입니다. 3에 대해서는 n을 3으로 나눈 나머지로 나눕니다. 나머지가 0이면 n이, 2이면 n + 1이, 1이면 2n + 1(= 2 × (3m + 1) + 1 = 3(2m + 1))이 3의 배수입니다. 2와 3은 서로소이므로 곱은 6의 배수입니다.

5단계: 합칩니다.

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

n = 10이면 55² − 385 = 3025 − 385 = 2640으로 문제의 예와 같습니다. 원한다면 공통 인수를 묶어 한 식으로 쓸 수도 있습니다.

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

n = 10이면 9 × 10 × 11 × 32 / 12 = 31680 / 12 = 2640입니다. 이 인수분해 꼴은 n = 1에서 (n − 1) 때문에 D = 0이 바로 보인다는 장점이 있지만, 코드는 각 단계를 따라가기 쉬운 S(n)² − Q(n) 꼴을 씁니다.

최종 알고리즘

  1. s = n(n + 1) / 2를 구합니다.
  2. q = n(n + 1)(2n + 1) / 6을 구합니다.
  3. s × s − q를 돌려줍니다.
  4. 정답을 출력하기 전에 작은 입력을 확인합니다. Q(10) = 385, S(10)² = 3025, D(10) = 2640, D(1) = 0이어야 하고, n = 1부터 100까지 모두 닫힌 식과 반복문(가장 단순한 방법)의 결과가 같아야 합니다. 하나라도 다르면 0이 아닌 종료 코드로 끝냅니다.

정확성 증명

  • 3단계의 S(n) = n(n + 1) / 2는 모든 n ≥ 0에서 성립합니다. 짝짓기 논증이 n에 의존하지 않고, n = 0이면 양쪽이 0입니다.
  • 4단계의 항등식 (k + 1)³ − k³ = 3k² + 3k + 1은 전개만으로 성립하고, 이를 k = 1, …, n에 대해 더한 식과 그 뒤의 정리는 모두 등식 변형이므로 Q(n) = n(n + 1)(2n + 1) / 6도 모든 n ≥ 1에서 성립합니다. n = 0이면 양쪽이 0입니다.
  • 두 나눗셈은 위에서 보인 대로 나누어떨어지므로, 정수 나눗셈(C++의 /, Python의 //)이 수학적 몫과 정확히 같습니다.
  • 따라서 값이 타입 범위 안에 있는 한 프로그램이 돌려주는 s × s − q는 정의 그대로의 D(n)입니다. 범위 조건은 다음 절에서 다룹니다.

복잡도

시간과 공간. n에 상관없이 곱셈 몇 번, 나눗셈 두 번, 뺄셈 한 번이므로 시간은 O(1), 저장하는 값도 몇 개뿐이라 추가 공간은 O(1) 입니다(값이 기계어 정수에 들어가는 동안). 가장 단순한 방법은 O(n) 시간이었습니다.

정수 범위.

  • n = 100에서 나오는 값은 n(n + 1) = 10100, S = 5050, 나눗셈 전의 n(n + 1)(2n + 1) = 2030100, Q = 338350이고, 가장 큰 값인 S² 역시 3000만이 안 됩니다. 32비트 정수로도 충분하지만, C++ 구현은 일반화했을 때 여유를 두려고 unsigned long long(표준이 최소 64비트를 보장)을 씁니다.
  • 일반화할 때 먼저 넘치는 곳은 S²입니다. S² ≤ 2⁶⁴ − 1이려면 S ≤ 2³² − 1 = 4294967295여야 합니다. n = 92681이면 n(n + 1) = 8589860442, S = 4294930221로 이 범위 안이고, n = 92682이면 S = 4295022903으로 2³²를 넘으므로 S²이 64비트를 넘칩니다. 그래서 코드 주석에 n <= 92681이라고 적었습니다. 나눗셈 전의 n(n + 1)(2n + 1)은 대략 2n³이라 n이 200만 정도가 될 때까지 64비트에 들어가므로, 이 범위에서는 S² 쪽이 더 먼저 한계에 닿습니다.
  • 뺄셈 s × s − q는 2단계에서 D(n) ≥ 0을 보였으므로 부호 없는 정수에서도 0 아래로 돌아가지 않습니다.
  • 부호 없는 정수는 넘치면 정의되지 않은 동작 대신 2ᴺ(N은 타입의 비트 수)으로 나눈 나머지로 돌아갑니다. 프로그램이 멈추지 않고 조용히 틀린 값을 낸다는 뜻이므로, 위의 범위를 넘는 n에는 128비트 정수나 임의 정밀도 정수가 필요합니다.
  • Python의 int는 정밀도에 제한이 없어 넘칠 일이 없습니다. 다만 /는 부동소수점 결과를 내므로 정수 나눗셈 //를 씁니다. 나누어떨어짐을 보였으므로 //는 정확한 몫입니다.

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는 가장 단순한 방법을 그대로 옮긴 함수이고, n = 1부터 100까지 닫힌 식과 같은 값을 내는지 대조하는 데만 씁니다. 공식의 유도에 실수가 있다면 이 대조에서 걸립니다. 작은 예의 확인은 assert가 아니라 if로 했습니다. assert는 NDEBUG가 정의된 빌드에서 아무 일도 하지 않지만, if의 비교와 return 1은 프로그램의 관찰 가능한 동작이라 -O2 최적화로도 없어지지 않습니다.

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()

C++와 같은 알고리즘, 같은 확인입니다. 합을 담는 변수를 sum이 아니라 total로 부른 것은 Python 내장 함수 sum을 가리지 않기 위해서이고, C++에서도 같은 이름을 맞춰 썼습니다. range(1, n + 1)은 끝값을 포함하지 않으므로 하나 더 높여 두었습니다. 작은 예가 틀리면 raise SystemExit(1)로 종료 코드 1을 냅니다. assert는 python -O로 실행하면 제거되므로 쓰지 않았습니다.

정답과 검증

처음 100개의 자연수에 대해 합의 제곱과 제곱의 합의 차이는 25164150입니다.

손으로 따라갈 수 있는 계산은 다음과 같습니다.

  • 합. S(100) = 100 × 101 / 2 = 5050, S(100)² = 5050² = 25000000 + 2 × 5000 × 50 + 50² = 25502500.
  • 제곱의 합. Q(100) = 100 × 101 × 201 / 6 = 2030100 / 6 = 338350.
  • 차이. 25502500 − 338350 = 위의 정답.
  • 인수분해 꼴로 다시. 99 × 100 × 101 × 302 / 12 = 999900 × 302 / 12 = 301969800 / 12 = 위의 정답. 서로 다른 모양의 식이 같은 값을 냅니다.

검증은 세 겹입니다.

  1. 작은 예와 전 구간 대조. 두 프로그램 모두 답을 출력하기 전에 Q(10) = 385, S(10)² = 3025, D(10) = 2640, D(1) = 0을 확인하고, n = 1부터 100까지 닫힌 식과 반복문의 결과가 모두 같은지 확인합니다. 하나라도 다르면 0이 아닌 종료 코드로 끝납니다. n = 10은 “풀이 유도”의 각 단계에서 손으로도 따라가 보았습니다.
  2. 두 구현의 일치. 이 시리즈의 저장소에서는 검사 스크립트가 C++ 코드를 g++ -std=c++17 -O2 -Wall -Wextra -Werror로 컴파일해 실행하고(경고도 실패로 처리), Python 코드를 python3 -I -B로 실행합니다. 두 출력이 저장해 둔 answer.txt와 바이트 단위로 같아야 하고, 0이 아닌 종료 코드나 시간 초과(C++ 1초, Python 10초)는 실패입니다. 시리즈 규칙상 이 검사를 통과한 코드만 글에 싣습니다.
  3. 증명. 두 공식과 나누어떨어짐은 “풀이 유도”와 “정확성 증명”에서 보였습니다.

경계 사례

  • n = 1: S = 1, Q = 1이라 D = 0입니다. 서로 다른 두 수의 쌍이 없으니 2단계의 해석과도 맞습니다. 작은 예에 D(1) = 0이 들어 있습니다.
  • n = 0: 문제에는 없지만 S = Q = 0이라 D = 0이고, 빈 합을 0으로 두는 관례와 맞습니다. 코드도 같은 값을 냅니다.
  • n = 2: S² = 9, Q = 5, D = 4 = 2 × (1 × 2). 쌍이 하나뿐인 가장 작은 0이 아닌 경우입니다.
  • 나눗셈 순서: n(n + 1)(2n + 1)을 다 곱한 뒤 6으로 나누므로 중간에 소수점이 생기지 않습니다. 먼저 n / 6처럼 나누면 정수 나눗셈에서 나머지가 버려져 틀린 값이 나옵니다.
  • 범위 밖의 n: 위 “정수 범위”처럼 n ≥ 92682이면 C++ 구현의 S²이 64비트를 넘칩니다. 이 프로그램은 n = 100만 다루므로 문제가 없지만, 함수를 다른 n에 재사용할 때는 주의가 필요합니다.

더 생각해 볼 것

  • 세제곱의 합. 4단계와 같은 방법으로 (k + 1)⁴ − k⁴ = 4k³ + 6k² + 4k + 1을 더해 1³ + 2³ + … + n³을 구해 보세요. 결과가 S(n)²과 같다는 것, 즉 “합의 제곱 = 세제곱의 합”이라는 것을 확인할 수 있습니다. 그렇다면 이 문제의 D(n)은 무엇과 무엇의 차이로도 읽힐까요?
  • 차이를 직접 세기. 2단계에 따르면 D(n) = 2 × (i < j인 쌍의 곱의 합)입니다. 이 합을 j에 대해 묶어 i < j인 i의 합 S(j − 1)을 쓰면 다른 경로로 같은 공식이 나옵니다. 직접 계산해 보세요.
  • 더 높은 거듭제곱. 망원합을 (k + 1)ᵖ⁺¹ − kᵖ⁺¹로 바꾸면 1ᵖ + … + nᵖ도 n에 대한 다항식으로 얻을 수 있습니다. 이 다항식의 차수는 왜 p + 1일까요?
  • 1번 문제와의 연결. 이 시리즈 1번 문제에서는 등차수열의 합을 포함-배제와 함께 썼습니다. 같은 S(n) 공식이 여기서 다시 쓰였다는 점을 비교해 보세요.

출처와 라이선스

  • Project Euler, Problem 6 — 원래 문제. 이 글의 “문제” 절은 그 내용을 직접 요약한 것이고, 원문과 같은 CC BY-NC-SA 4.0으로 공유합니다. 문제의 저작자 표시는 Project Euler입니다.
  • Project Euler, Copyright — 문제 본문의 라이선스와 링크 방식. 비영리 조건에 따라 이 글에는 광고를 싣지 않습니다(ads: false).
  • C++ 정수 타입의 최소 폭과 부호 없는 정수의 산술은 C++ 표준 작업 초안 [basic.fundamental], assert와 NDEBUG는 [assertions.assert]를 따랐습니다. Python의 assert와 -O는 The assert statement, int의 정밀도와 /·//의 결과는 Numeric Types를 따랐습니다.
  • 풀이, 증명, 코드는 직접 쓴 것입니다. 합과 제곱의 합 공식은 외부 자료에 기대지 않고 위에서 직접 유도했습니다. Project Euler의 해설(overview), 포럼, 회원 전용 자료나 다른 풀이 사이트는 참고하지 않았습니다. Motitory는 Project Euler와 제휴나 후원 관계가 없고, Project Euler가 이 글을 검토하지도 않았습니다.

블로그로