· 開発

Project Euler 2: 偶数のフィボナッチ数だけを足す

フィボナッチ数列では3項ごとに偶数が現れることを証明し、偶数項どうしの漸化式 E(k) = 4E(k−1) + E(k−2) を導いて、400万以下の偶数項だけを足します。C++17とPython 3の実装、正しさの証明、整数の範囲まで扱います。

Project Eulerアルゴリズム数学フィボナッチ数

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

Project Euler シリーズの2問目です。フィボナッチ数列を順に作りながら偶数だけを足しても解けますが、先に偶数の項がどこに現れるかを調べておくと、偶数項から次の偶数項へ直接進む漸化式が得られます。偶奇(2で割った余り)を追う方法も、漸化式を変形する方法も、後の問題でよく使います。

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

問題

1と2から始まり、以降の項は直前の2項の和で作るフィボナッチ数列を考えます。最初のいくつかの項は 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, … です。このうち値が 400万(4,000,000)を超えない項だけを対象にして、その中の偶数の項をすべて足した値を求めます。

境界は「超えない」、つまり4,000,000 以下です。4,000,000自体はフィボナッチ数ではないので、この問題では「未満」と読んでも答えは同じですが、関数として一般化するときは「以下」を正確に実装しておきます。

前提知識

フィボナッチ数列と番号の付け方。 この記事では問題の初期値どおり F₁ = 1、F₂ = 2 とし、n ≥ 3 では Fₙ = Fₙ₋₁ + Fₙ₋₂ と定義します。すると F₃ = 3、F₄ = 5、F₅ = 8、F₆ = 13、F₇ = 21、F₈ = 34 です。よく見る 1, 1, 2, 3, … とは番号が一つずれるので、ほかの資料の公式を使うときは番号を合わせる必要があります。

数列は狭義単調増加です。 F₂ > F₁ で、n ≥ 3 ではすべての項が正なので Fₙ = Fₙ₋₁ + Fₙ₋₂ > Fₙ₋₁ です。したがって同じ値は2回現れず、ある項が上限を超えればそれ以降の項もすべて超えます。「上限を超えたら止める」ループが正しい根拠はこれです。

偶奇の足し算。 偶数 + 偶数 = 偶数、偶数 + 奇数 = 奇数、奇数 + 奇数 = 偶数です。2で割った余りどうしを足して、もう一度2で割った余りを見ればよいからです(0 + 0 = 0、0 + 1 = 1、1 + 1 = 2 → 0)。フィボナッチ数列の各項は直前の2項の和なので、直前の2項の偶奇がわかれば次の項の偶奇が決まります。

指数的な増加。 項が指数的に大きくなると、上限 N 以下の項の数は N ではなく log N に比例します。後の「計算量」で、偶数項について具体的な不等式として示します。

素朴な方法

まず思いつく方法は二つあります。

  1. 1から4,000,000までのすべての整数について、フィボナッチ数かどうか、偶数かどうかを調べる。上限 N に対して O(N) 回以上の確認が必要で、フィボナッチ数かどうかの判定にも別にコストがかかります。
  2. フィボナッチ数列を実際に作りながら、上限以下の項のうち偶数のものだけを足す。
a, b = 1, 2
total = 0
while a <= 4000000:
    if a % 2 == 0:
        total += a
    a, b = b, a + b

二つ目の方法でもすでに十分速いです。フィボナッチ数は指数的に増えるので、4,000,000以下の項は32個しかなく、ループも32回で終わります。時間計算量は O(log N)、空間は O(1) です。ただし32個のうち実際に足す偶数項は11個で、残りの21回は偶奇を調べて捨てるだけです。以下では偶数項がどこに現れるかを正確に証明し、偶数項だけを直接作る漸化式を導きます。漸近的な計算量は同じですが、ループの回数はおよそ3分の1になり、偶奇の判定もなくなります。何より「なぜそうなるか」を式で説明できます。

解法の導出

ステップ1: 小さな例でパターンを見る。 最初の10項と偶奇を並べて書きます。

n   :  1  2  3  4  5  6  7  8  9 10
Fₙ  :  1  2  3  5  8 13 21 34 55 89
偶奇:  奇 偶 奇 奇 偶 奇 奇 偶 奇 奇

偶数項は F₂ = 2、F₅ = 8、F₈ = 34 で、番号は2、5、8と3つおきです。89以下の偶数項の和は 2 + 8 + 34 = 44 で、この値は後でコードの小さな例の確認に使います。

ステップ2: 偶数項がちょうど3項ごとに現れることを証明する。 連続する2項の偶奇の組 (Fₙ, Fₙ₊₁) を見ます。次の組は (Fₙ₊₁, Fₙ + Fₙ₊₁) なので、今の組だけで決まります。

(奇, 偶) → (偶, 奇 + 偶 = 奇) → (奇, 偶 + 奇 = 奇) → (奇, 奇 + 奇 = 偶) → (奇, 偶) → …

最初の組 (F₁, F₂) = (1, 2) は (奇, 偶) で、3回進むと再び (奇, 偶) に戻ります。次の組は今の組だけで決まるので、この周期3の循環は永遠に繰り返されます。したがって Fₙ が偶数になるのは、n を3で割った余りが2のとき、つまり n = 3k − 1(k = 1, 2, 3, …)のときだけです。以下、偶数項を Eₖ = F₃ₖ₋₁ と呼びます。E₁ = F₂ = 2、E₂ = F₅ = 8、E₃ = F₈ = 34 です。

ステップ3: 3つ離れた項の関係を作る。 n ≥ 7 なら、定義 Fₘ = Fₘ₋₁ + Fₘ₋₂ を m = n, n − 1, n − 2, n − 4 に使えます(どれも3以上)。1段階ずつ展開します。

Fₙ = Fₙ₋₁ + Fₙ₋₂
   = (Fₙ₋₂ + Fₙ₋₃) + Fₙ₋₂              ← Fₙ₋₁ を展開
   = 2Fₙ₋₂ + Fₙ₋₃
   = 2(Fₙ₋₃ + Fₙ₋₄) + Fₙ₋₃             ← Fₙ₋₂ を展開
   = 3Fₙ₋₃ + 2Fₙ₋₄
   = 3Fₙ₋₃ + Fₙ₋₄ + (Fₙ₋₅ + Fₙ₋₆)       ← Fₙ₋₄ の一つを展開
   = 3Fₙ₋₃ + (Fₙ₋₄ + Fₙ₋₅) + Fₙ₋₆
   = 3Fₙ₋₃ + Fₙ₋₃ + Fₙ₋₆               ← Fₙ₋₄ + Fₙ₋₅ = Fₙ₋₃
   = 4Fₙ₋₃ + Fₙ₋₆

最後の Fₙ₋₄ + Fₙ₋₅ = Fₙ₋₃ も定義を m = n − 3 に使ったもので、n ≥ 7 なら n − 3 ≥ 4 なので使えます。

ステップ4: 偶数項どうしの漸化式にする。 n = 3k − 1 を代入すると n − 3 = 3(k − 1) − 1、n − 6 = 3(k − 2) − 1 なので

Eₖ = 4Eₖ₋₁ + Eₖ₋₂    (k ≥ 3、このとき n = 3k − 1 ≥ 8 ≥ 7)

です。手で検算すると E₃ = 4 · 8 + 2 = 34、E₄ = 4 · 34 + 8 = 144 で、実際に F₁₁ = 144 です(F₉ = 55、F₁₀ = 89、F₁₁ = 144)。

ステップ5: 初期値を整える。 漸化式は k ≥ 3 で証明したので、E₁ = 2、E₂ = 8 を初期値にすれば足ります。コードでは便宜上 E₀ = 0 を一つ加えます。4 · E₁ + E₀ = 4 · 2 + 0 = 8 = E₂ なので、(E₀, E₁) = (0, 2) から始めても同じ式で E₂ 以降が正しく続きます。E₀ は足さない補助の値にすぎません。

最終的なアルゴリズム

上限を limit とすると:

  1. previous = 0(E₀)、current = 2(E₁)、sum = 0 から始めます。
  2. current ≤ limit の間、sum に current を足し、(previous, current) を (current, 4 · current + previous) に置き換えます。
  3. current > limit になったら sum を返します。
  4. 答えを出力する前に、小さな上限で結果を確かめます。1 → 0(偶数項なし)、2 → 2、8 → 2 + 8 = 10(上限の値そのものを含むかの確認)、89 → 44(ステップ1の例)。一つでも違えば0以外の終了コードで終わります。

C++ の実装には、非常に大きな上限で次の項を計算するときにあふれないようにする判定がもう一つあります。後の「整数の範囲」で説明します。

正しさの証明

ループが j 回目の判定(j = 1, 2, …)をするとき、次が成り立ちます(不変条件)。

  • previous = Eⱼ₋₁、current = Eⱼ
  • sum = E₁ + E₂ + … + Eⱼ₋₁

最初(j = 1)は previous = 0 = E₀、current = 2 = E₁、sum = 0(空の和)なので成り立ちます。current ≤ limit なら sum に Eⱼ を足して E₁ + … + Eⱼ となり、新しい current は 4Eⱼ + Eⱼ₋₁ = Eⱼ₊₁(j = 1 ならステップ5、j ≥ 2 ならステップ4の漸化式)、新しい previous は Eⱼ なので、j + 1 でも成り立ちます。

E はフィボナッチ数列の一部なので狭義単調増加で、いつか必ず limit を超えます。初めて Eⱼ > limit になったときにループが終わり、そのとき sum = E₁ + … + Eⱼ₋₁ です。単調増加なので E₁, …, Eⱼ₋₁ はすべて limit 以下で、Eⱼ 以降はすべて limit を超えます。ステップ2で偶数のフィボナッチ項はちょうど E₁, E₂, … であることを示し、数列が単調増加なので各値は1回ずつしか現れません。したがって sum は limit 以下の偶数のフィボナッチ項の和とちょうど一致します。

計算量

ループの回数。 k ≥ 3 では Eₖ = 4Eₖ₋₁ + Eₖ₋₂ > 4Eₖ₋₁、また E₂ = 8 = 4E₁ なので、k ≥ 2 で Eₖ ≥ 4Eₖ₋₁ です。よって Eₖ ≥ 2 · 4^(k−1) となり、limit 以下の偶数項の個数 K は 2 · 4^(K−1) ≤ limit から K ≤ 1 + log₄(limit / 2) です。limit = 4,000,000 なら log₄(2,000,000) ≈ 10.47 なので K ≤ 11 で、実際に11個です。ループは O(log limit) 回で、1回は固定幅の整数では定数時間なので時間計算量は O(log limit)、保持する値は3つだけなので空間計算量は O(1) です。

整数の範囲。 この問題の値はどれも小さく、最後に足す偶数項は400万以下、ループを終えるときに作る次の偶数項も1500万程度なので、32ビット整数にも収まります。それでも関数をほかの上限で使い回せるように範囲を確認しておきます。

  • 和の大きさ: 上の Eⱼ₋₁ ≤ Eⱼ / 4 を繰り返し使うと、最後に足した項を E_K として、和は E_K (1 + 1/4 + 1/16 + …) < (4/3) E_K ≤ (4/3) limit です。和は常に上限の4/3倍より小さくなります。
  • C++ は unsigned long long を使います。C++ 標準はこの型が少なくとも64ビットであることを保証し(C 標準の ULLONG_MAX の最小値に従う)、符号なし整数の演算は 2^N を法として定義されているので未定義動作になりません。ただしあふれると黙って誤った値になるので、あふれないようにする必要があります。
  • 次の項 4 · current + previous は、上限が 2^64 に近いとあふれることがあります。あふれて小さな値に回り込むと、ループが終わらなくなったり誤った項を足したりします。そこで C++ の実装は次の項を計算する前に current > (limit − previous) / 4 かどうかを見ます。整数では 4 · current > limit − previous は current > ⌊(limit − previous) / 4⌋ と同じなので、この条件が真なら次の項は limit を超え、ループはどのみち終わります。計算せずに止めても結果は同じです。ループの中では previous ≤ current ≤ limit なので、limit − previous も0を下回って回り込むことはありません。
  • この判定により、次の項は limit 以下のときにしか計算されないのであふれず、和は (4/3) limit より小さいので limit ≤ 2^63 なら 2^64 未満です。つまり C++ の関数は limit ≤ 2^63 で正確です。
  • Python の int は精度に上限がなくあふれないので、この判定は置いていません。アルゴリズムは同じです。

境界のケース

  • limit < 2: 最小の偶数項2より小さいので偶数項はなく、結果は0です。ループの条件が最初から偽なので特別扱いは要りません。Python では0や負の値を渡しても0になります。C++ は符号なしの型なので負の値は受け取りません。
  • limit が偶数項と等しいとき: 「超えない」は「以下」なので、その項も足します。limit = 2 なら2、limit = 8 なら10でなければならず、コードの小さな例の確認にこの二つが入っています。
  • 初期値が 1, 1 のとき: 1, 1, 2, 3, … で始める別の定義では1が2回現れますが、1は奇数なので偶数項の和には影響しません。1と2から始めるという条件は、偶数項の和に関しては結果を変えません。
  • 非常に大きな limit(C++): 上の整数の範囲のとおり、limit ≤ 2^63 で正確です。それより大きいと和があふれることがあります。

C++17 による実装

#include <iostream>

// Sum of the even Fibonacci terms (1, 2, 3, 5, 8, ...) that do not exceed limit.
// The even terms are 2, 8, 34, ..., and each one is 4 * previous + the one before.
unsigned long long sum_even_fibonacci_up_to(unsigned long long limit) {
    unsigned long long previous = 0;
    unsigned long long current = 2;
    unsigned long long sum = 0;
    while (current <= limit) {
        sum += current;
        // previous <= current <= limit, so limit - previous does not wrap.
        // If 4 * current + previous would exceed limit, stop before computing it.
        if (current > (limit - previous) / 4) {
            break;
        }
        const unsigned long long next = 4 * current + previous;
        previous = current;
        current = next;
    }
    return sum;
}

int main() {
    if (sum_even_fibonacci_up_to(1) != 0 || sum_even_fibonacci_up_to(2) != 2 ||
        sum_even_fibonacci_up_to(8) != 10 || sum_even_fibonacci_up_to(89) != 44) {
        return 1;
    }
    std::cout << sum_even_fibonacci_up_to(4000000) << '\n';
    return 0;
}

小さな例の確認には assert ではなく if を使いました。assert は NDEBUG を定義したビルドで消えるからです。確認が一つでも外れれば return 1 で0以外の終了コードを返します。

Python 3 による実装

def sum_even_fibonacci_up_to(limit: int) -> int:
    """Sum of the even Fibonacci terms (1, 2, 3, 5, 8, ...) that do not exceed limit.

    The even terms are 2, 8, 34, ..., and each one is 4 * previous + the one before.
    """
    previous, current = 0, 2
    total = 0
    while current <= limit:
        total += current
        previous, current = current, 4 * current + previous
    return total


def main() -> None:
    if (
        sum_even_fibonacci_up_to(1) != 0
        or sum_even_fibonacci_up_to(2) != 2
        or sum_even_fibonacci_up_to(8) != 10
        or sum_even_fibonacci_up_to(89) != 44
    ):
        raise SystemExit(1)
    print(sum_even_fibonacci_up_to(4000000))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ小さな例です。previous, current = current, 4 * current + previous は右辺をすべて評価してからまとめて代入するので、一時変数は要りません。小さな例が合わなければ raise SystemExit(1) で終了コード1を返します。assert は Python を -O 付きで実行すると取り除かれるので使っていません。

答えと検証

漸化式で400万以下の偶数項を順に作ると次のようになります(かっこ内はそこまでの和)。

E₁  = 2                              (2)
E₂  = 8                              (10)
E₃  = 4 · 8      + 2      = 34       (44)
E₄  = 4 · 34     + 8      = 144      (188)
E₅  = 4 · 144    + 34     = 610      (798)
E₆  = 4 · 610    + 144    = 2584     (3382)
E₇  = 4 · 2584   + 610    = 10946    (14328)
E₈  = 4 · 10946  + 2584   = 46368    (60696)
E₉  = 4 · 46368  + 10946  = 196418   (257114)
E₁₀ = 4 · 196418 + 46368  = 832040   (1089154)
E₁₁ = 4 · 832040 + 196418 = 3524578
E₁₂ = 4 · 3524578 + 832040 = 14930352 > 4,000,000 → 停止

したがって答えは 1089154 + 3524578 = 4613732 です。

検証は三つの方法で行います。

  1. 小さな例。 二つのプログラムはどちらも、答えを出力する前に上限1、2、8、89の結果が0、2、10、44かどうかを確かめ、違えば0以外の終了コードで終わります。上限89の値44は「解法の導出」のステップ1で、最初の10項から手でも確認しました。
  2. 二つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。両方の出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則として、この検査に通ったコードだけを記事に載せます。
  3. 手計算と別の方法。 上の表の各項が本当にフィボナッチ数か(たとえば 610 = F₁₄、2584 = F₁₇)は、2項ずつ足していけば確かめられます。「素朴な方法」のループで32個の項をすべて作り、偶数だけを足して同じ値になるかを比べてみるのもよい練習です。

発展

  • 奇数項の和。 同じ上限で奇数項の和はどう求めるでしょうか。全項の和から偶数項の和を引けば求まります。全項の和 F₁ + F₂ + … + Fₙ には Fₙ₊₂ を使った閉じた式があるので、この記事の番号の付け方(F₁ = 1、F₂ = 2)に合わせて自分で導いてみてください。
  • ほかの余り。 3の倍数のフィボナッチ項は何項ごとに現れるでしょうか。ステップ2と同じように (Fₙ mod 3, Fₙ₊₁ mod 3) の組を追えば周期がわかります。それらの項の間にも、ステップ3と同じ方法で漸化式が作れるか考えてみてください。
  • 3項ではなく m 項。 ステップ3の展開を一般化すると、Fₙ を Fₙ₋ₘ と Fₙ₋₂ₘ で表せます。m = 2 のときの係数は何か、m = 3 のときになぜ4が出てくるのかを確かめてみてください。
  • オーバーフローの実験。 C++ のコードからあふれの判定を外し、limit に ULLONG_MAX に近い値を渡すと何が起きるでしょうか。limit が 2^63 より大きいとき、和が実際にあふれるかどうかも Python の結果と比べてみてください。
  • このシリーズの後の問題。 フィボナッチ数と漸化式は、1〜100番の中の別の問題にもまた出てきます。項が固定幅の整数に収まらないほど大きくなることもあるので、今回の整数の範囲の分析を覚えておくと役に立ちます。

出典とライセンス

  • Project Euler, Problem 2 — 元の問題。この記事の「問題」の節はその内容を独自に要約したもので、原文と同じ CC BY-NC-SA 4.0 で共有します。問題の帰属表示は Project Euler です。
  • Project Euler, Copyright — 問題文のライセンスとリンクの方法。非営利の条件に従い、この記事には広告を載せません。
  • 解法、証明、コードは独自に書いたものです。Project Euler の解説(overview)、フォーラム、会員限定の資料やほかの解答サイトは参照していません。Motitory は Project Euler と提携・スポンサー関係になく、Project Euler がこの記事を確認したこともありません。
  • C++17 最終作業ドラフト N4659 — §6.9.1 基本型(符号なし整数の 2^N を法とする演算)、§21.3.6 <climits>(範囲は C 標準に従う)、§22.3 assert(NDEBUG)
  • Python 標準ライブラリ: Numeric Types — int の無制限の精度
  • Python 言語リファレンス: The assert statement — 最適化モードでは assert が取り除かれる

ブログへ