· 開発

Project Euler 6: 和の2乗と2乗の和の差を閉じた式でO(1)で求める

1から100までの和の2乗と2乗の和の差を、ループではなく n(n + 1) / 2 と n(n + 1)(2n + 1) / 6 の2つの公式から直接計算します。2つの公式をペアの組み合わせと望遠鏡和で一から導き、整数の範囲、C++17とPython 3の実装まで扱います。

Project Eulerアルゴリズム数学級数

ネタバレ注意 · この記事には解き方と答えがそのまま書かれています。先に自分で解いてみたい方は、Project Euler の問題 6からどうぞ。

Project Euler シリーズの6問目です。1から100までを足したり2乗したりするのはループ1つで済むので、この問題は計算そのものより、ループを使わずに答えを直接出す式を導く練習といえます。1問目で使った等差数列の和をもう一度取り出し、そこに2乗の和の公式を望遠鏡和(telescoping sum)で一から導いて加えます。どちらも後の問題で再び使う道具です。

ネタバレ注意: この記事には解法の全体と答えが含まれます。先に自分で解きたい場合は、元の問題から始めてください。このシリーズの記事には広告を載せません。問題文のライセンス(CC BY-NC-SA 4.0)に非営利の条件があるためです。

問題

2つの値を比べます。1つは1から n までの和を求めてから2乗した値、もう1つは1から n までをそれぞれ2乗してから足した値です。問題文の例では n = 10 のとき、2乗の和は385、和の2乗は3025で、その差は2640です。同じ差を最初の100個の自然数、つまり n = 100 について求めるのが問題です。

入力はなく、n は100に固定されています。求めるのは (和の2乗) − (2乗の和) という整数1つです。

前提知識

記号。 この記事では2つの和を次のように呼びます。

S(n) = 1 + 2 + … + n             (和)
Q(n) = 1² + 2² + … + n²          (2乗の和)
D(n) = S(n)² − Q(n)              (求める差)

n = 0 なら足すものがないので S(0) = Q(0) = 0 とします(空の和は0)。

和を2乗すると出てくるもの。 いくつかの数の和を2乗すると、各数の2乗と、異なる2数の積が2回ずつ現れます。たとえば

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

です。かっこを展開すると、a, b, c から1つ選び、さらにもう1つ選ぶすべての場合が1回ずつ現れます。同じ数を2回選べば2乗になり、異なる2数は (a, b) と (b, a) のように2通りの順序で1回ずつ現れるからです。この事実は「解法の導出」で差が何を表すかを読むのに使います。

ペアで求める等差数列の和。 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³ として1項 (k + 1)³ − k³ を展開しておくと、2乗の和を取り出せます。これが下のステップ4の核心です。

素朴な方法

1から n までループしながら和と2乗の和を同時に積み上げ、最後に和を2乗して引けば済みます。

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¹² のように大きくなれば反復回数もそれだけ増え、固定幅の整数を使う言語では、その前に2乗の和があふれるおそれがあります。逆に、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

問題文の3つの値と一致します。以下の各ステップでは、これらの数で手計算による検算をします。

ステップ2: 差が何を表すかを読む。 前提知識の展開を 1, 2, …, n に当てはめると

S(n)² = (1 + 2 + … + n)² = Q(n) + 2 × (i < j となるすべての組の積 i × j の和)

なので、D(n) は異なる2数の積をすべて足したものの2倍です。小さな例 n = 3 で確かめると S(3)² = 36、Q(3) = 14、D(3) = 22 で、2 × (1 × 2 + 1 × 3 + 2 × 3) = 2 × 11 = 22 と一致します。ここから2つのことがわかります。積はすべて正なので D(n) ≥ 0 であり、組が1つもない 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 で、問題文の例と一致します。共通因数をくくって1つの式にまとめることもできます。

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まですべてで閉じた式とループ(素朴な方法)の結果が一致しなければなりません。1つでも違えば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です。
  • 2つの割り算は上で示したとおり割り切れるので、整数除算(C++ の /、Python の //)は数学的な商と正確に一致します。
  • したがって、値が型の範囲に収まる限り、プログラムが返す s × s − q は定義どおりの D(n) です。範囲の条件は次の節で扱います。

計算量

時間と空間。 n に関係なく掛け算が数回、割り算が2回、引き算が1回なので、時間は 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) は終端の値を含まないので、終端を1つ上げています。小さな例が失敗すれば raise SystemExit(1) で終了コード1を返します。assert は python -O で実行すると取り除かれるので使っていません。

答えと検証

最初の100個の自然数について、和の2乗と2乗の和の差は 25164150 です。

手でたどれる計算は次のとおりです。

  • 和。 S(100) = 100 × 101 / 2 = 5050、S(100)² = 5050² = 25000000 + 2 × 5000 × 50 + 50² = 25502500。
  • 2乗の和。 Q(100) = 100 × 101 × 201 / 6 = 2030100 / 6 = 338350。
  • 差。 25502500 − 338350 = 上の答え。
  • 因数分解の形でもう一度。 99 × 100 × 101 × 302 / 12 = 999900 × 302 / 12 = 301969800 / 12 = 上の答え。形の違う2つの式が同じ値を出します。

検証は三重です。

  1. 小さな例と全範囲の照合。 2つのプログラムはどちらも答えを出力する前に Q(10) = 385、S(10)² = 3025、D(10) = 2640、D(1) = 0 を確認し、n = 1 から100まで閉じた式とループの結果がすべて一致するかを確認します。1つでも違えば0以外の終了コードで終わります。n = 10 は「解法の導出」の各ステップで手でもたどりました。
  2. 2つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。2つの出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則上、この検査を通ったコードだけを記事に載せます。
  3. 証明。 2つの公式と割り切れることは「解法の導出」と「正しさの証明」で示しました。

境界のケース

  • n = 1: S = 1、Q = 1 なので D = 0 です。異なる2数の組がないので、ステップ2の読み方とも一致します。小さな例に D(1) = 0 が入っています。
  • n = 0: 問題にはありませんが、S = Q = 0 なので D = 0 で、空の和を0とする慣習と合っています。コードも同じ値を返します。
  • n = 2: S² = 9、Q = 5、D = 4 = 2 × (1 × 2)。組がちょうど1つの、0でない最小のケースです。
  • 割り算の順序: n(n + 1)(2n + 1) をすべて掛けてから6で割るので、途中で端数は出ません。先に n / 6 のように割ると、整数除算で余りが捨てられて間違った値になります。
  • 範囲外の n: 「整数の範囲」のとおり、n ≥ 92682 では C++ の実装の S² が64ビットからあふれます。このプログラムは n = 100 しか扱わないので問題ありませんが、関数を別の n に再利用するときは注意が必要です。

発展

  • 3乗の和。 ステップ4と同じ方法で (k + 1)⁴ − k⁴ = 4k³ + 6k² + 4k + 1 を足し、1³ + 2³ + … + n³ を求めてみてください。結果が S(n)² に等しいこと、つまり「和の2乗 = 3乗の和」であることが確かめられます。では、この問題の D(n) は何と何の差としても読めるでしょうか。
  • 差を直接数える。 ステップ2によれば D(n) = 2 × (i < j となる組の積の和) です。この和を j ごとにまとめ、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 に従いました。
  • 解法、証明、コードは自分で書いたものです。和と2乗の和の公式は外部の資料に頼らず、上で直接導きました。Project Euler の解説(overview)、フォーラム、会員限定の資料や他の解答サイトは参照していません。Motitory は Project Euler と提携・後援の関係になく、Project Euler がこの記事を確認したこともありません。

ブログへ