偶数のフィボナッチ数

C Project Euler

けさは Project Euler の Problem 2 (Oct 19, 2001) をやってみました. ただし, 厳密には, けさになってはじめてやってみたというわけではありません. 数日前から考えてはいたのですが, 期待していたほど美しい解き方を思いつくことができず, 結局いまに至ってしまったというのが実情です.

出題内容は次のとおりです.

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

つまり, 4000000 未満のフィボナッチ数のうち, 偶数であるものをすべて足し合わせると, 合計はいくつになりますか, という問題です.

最も初歩的な解き方

おそらく最も初歩的な解き方は, フィボナッチ数列の漸化式を使ってフィボナッチ数を 1 項ずつ求め, それが偶数であれば合計に足していく, というものになろうかと思います.

なお, 問題文ではフィボナッチ数を 1, 2, 3, 5, 8, 13, ... と列挙していますが, 以下では 0, 1, 1, 2, 3, 5, 8, 13, ... と列挙することにします.

int get_sum0(int sup) {
  int sum = 0;
  int f0 = 0;
  int f1 = 1;
  int f2 = f1 + f0;
  while (f2 < sup) {
    if (f2 % 2 == 0)
      sum += f2;
    f0 = f1;
    f1 = f2;
    f2 = f1 + f0;
  }
  return sum;
}

引数名の sup というのは, 一応 「上限」 を意味するらしい英単語 supremum を略記したつもりです. 引数に 4000000 という値を渡すと, この関数は 4613732 という答えを返します. これが正解のようです.

少しループの回数を減らした解き方

前述の解き方は, フィボナッチ数を 1 項ずつ調べていますが, 少し考えれば, そもそもフィボナッチ数は 3 項ごとに偶数になることが明らかです. フィボナッチ数列の並びは, ... 奇数, 奇数, 偶数, 奇数, 奇数, 偶数, ... となるからです. したがって, できればループの回数も 3 分の 1 に抑えたい気がします.

そこで, 第 0 項 (= 0) と 第 1 項 (= 1) を出発点として, 3 項ごとにフィボナッチ数を求める方法はないものかと考えてみました.

まず, k を 0 以上の整数とすると, fk+3 は fk + 2fk+1 とあらわすことができます.

同様に, fk+6 は fk+3 + 2fk+4 とあらわすことができます. すでに fk+3 は 1 つ前のステップで得られていますので, あとは何らかの方法で fk+4 を得ておく必要があります.

ところで, fk+4 は 2fk + 3fk+1 とあらわすことができます.

したがって, これらを繰り返すことで, 3 項ごとに現れるフィボナッチ数を次々に得ることができそうです. たとえば, 次のようになろうかと思います.

int get_sum1(int sup) {
  int sum = 0;
  int a = 0;
  int b = 1;
  int f = a + 2 * b;
  while (f < sup) {
    sum += f;
    b = 2 * a + 3 * b;
    a = f;
    f = a + 2 * b;
  }
  return sum;
}

実際に, 引数に 4000000 という値を渡すと, この関数は 4613732 という答えを返します.

ちょっと乱暴かもしれない解き方

隣り合う 2 つのフィボナッチ数の比は黄金比に収束します. これを利用すれば, 黄金比の 3 乗 (= 2 + √5) をかけていくことで, 3 項ごとに現れるフィボナッチ数を次々に得ることができるかもしれません.

#include <math.h>

int get_sum2(int sup) {
  double r = 2 + sqrtl(5);
  int sum = 0;
  int f = 2;
  while (f < sup) {
    sum += f;
    f = round(f * r);
  }
  return sum;
}

ただし, これは乱暴な解き方かもしれません. というのも, 黄金比の 3 乗をかけたときの小数点以下の取り扱いについては単純に四捨五入してしまってよいものなのだろうか, ならびに, 浮動小数点数で計算を続けていたらやがては誤差が発生して計算が不正確になってしまうことがありうるのではないだろうか, といった点について, ほとんど深く考えていないからです.

実際に, わたしの手元の環境で, 引数に 4000000 という値を渡すと, この関数は 4613732 という答えを返しましたが, たぶんたまたま誤差が発生せずに済んだか, 発生した誤差がたまたま合計としては相殺されて隠れてしまったか, なんじゃないかという気がしています. (具体的に検証はしていませんが...)

(なお, 蛇足になりますが, 上記のコードで使った round という関数は, C89/C90 の標準ライブラリには含まれていません. C99 で追加された関数のようです.)

処理速度の比較

第 1 の解き方と第 2 の解き方を比較すると, 後者は約半分の処理時間で済むようです. ループの回数が 3 分の 1 になったとはいえ, さすがに処理時間も 3 分の 1 になるわけではありませんでしたが.

第 3 の解き方の場合は, 第 1 の解き方の 2.5 倍くらいの処理時間がかかりました.

その他

インターネット上で Project Eulerに挑戦: 問題2 (Oct 16, 2011) という記事を見かけました. わたしの理解を超える内容を含むもので, ただただ腰を抜かすばかりですが, ひとつ関心を覚える点がありました. すなわち, 偶数のフィボナッチ数を順々に足し上げていくことをしなくても, 合計を得ることができる方法が示唆されていた点です. たとえば, f3 (= 2), f6 (= 8), f9 (= 34) の合計 (= 44) であれば, それは f11 (= 89) を使って, (f11 - 1) / 2 とあらわすことができるようです. うまく一般化すれば, もう少し処理速度の速い解き方ができるのかもしれません. またの機会に考えてみたいと思います.