Project Euler 問題25（1000桁のフィボナッチ数）を解いてみる（続き）

前回は、Project Euler で紹介されている問題25を解いてみました。前回は、フィボナッチ数列の定義に従い計算をしましたが、今回は、一般項を求めて計算します。

再掲載　Project Euler 問題25（1000桁のフィボナッチ数）

考察

フィボナッチ数列の漸化式は、以下となります。

$$a_{n + 2} = a_{n + 1} + a_n,\; a_1 = 1,\; a_2 = 1$$

フィボナッチ数列の一般項 $a_n$ は、以下の式となります。

$$a_n = \frac{1}{\sqrt{5}}\left\{
\left(\frac{1 + \sqrt{5}}{2}\right)^n – \,
\left(\frac{1 – \sqrt{5}}{2}\right)^n \right\}$$

一般項の2番目の項は、$(1-\sqrt{5})/2 = -0.61803 \dots$ なので、$n$ 乗していけば、小さくなっていきます。例えば $10$ 乗すると、$((1-\sqrt{5})/2)^{10} = 0.00813 \dots$ となりほとんど無視できる値となります。つまり以下となります。

$$a_n \fallingdotseq \frac{1}{\sqrt{5}}
\left(\frac{1 + \sqrt{5}}{2}\right)^n$$

近似式のように見えますが、右辺の式を整数単位で切り上げれば正確な値になります。

$a_n$ が $10^{999}$ 以上となれば、1000桁になります。

$$a_n \fallingdotseq \frac{1}{\sqrt{5}}
\left(\frac{1 + \sqrt{5}}{2}\right)^n \geqq 10^{999}$$

この式の両辺の $\log_{10}$ を取ると以下のように計算できます。

$$\log_{10} \left( \frac{1}{\sqrt{5}}
\left(\frac{1 + \sqrt{5}}{2}\right)^n \right) \geqq \log_{10} 10^{999}$$

一般的に $\log (a/b) = \log a – \log b$、$\log a^n = n \log a$ が成り立ちます。また $\log$ の定義から $log_{10} 10 = 1$ となります。よって、以下のように式を変形することができます。

$$n \log_{10} \frac{1 +\sqrt{5}}{2}\, -\, \log_{10} \sqrt{5} \geqq 999$$

$$n \geqq \frac{999 + \log_{10} \sqrt{5}}
{\log_{10} \frac{1 +\sqrt{5}}{2}}$$

$n$ を切り上げて整数にすれば、求める解答となります。複雑ですが、この式は $O(1)$ で求めることができます。関数電卓があれば、手でも計算できます。

解答例

一般項から、直接 $n$ を計算することができました。求めた式を Python で計算してみましょう。プログラムは以下となります。

"""Project Euler Problem025"""
import math

N = 999
sqrt5 = math.sqrt(5.0)
answer = math.ceil((N + math.log10(sqrt5))/math.log10((1.0 + sqrt5)/2.0))

print("Problem025: " + str(answer))

2行目で桁数を変更できます。$n = 999$ でも、$9999$ でも、$99999$ でも計算結果はすぐに出力されます。

C では、以下となります。値を大きくできるように関係する変数は long long int 型にしました。math.h を使っているため、コンパイル時（リンク時）に、「-lm」オプションが必要となります。またオンライン実行環境では、コンパイルできない環境があるかもしれません（過去に紹介した paiza.iO では動作しました）。

#include <stdio.h>
#include <math.h>

int main(void)
{
	long long int n = 999;
	long long int term;
	double sqrt5;

	sqrt5 = sqrt(5.0);
	term = (long long int)ceil(((double)n + log10(sqrt5))/log10((1.0 + sqrt5)/2.0));

	printf("Problem025: %lld\n", term);

	return 0;
}

C 版も6行目で桁数を変更できます。

結果を振り返る

前回は、問題の1000桁を増やしていきました。100万桁を超えると計算量が大きくなり計算をあきらめました。今回は、計算量が一定です。追加で結果を記します。

求めた式を近似値で記すと以下となります。

$$n \geqq \frac{(X – 1) + 0.8047 \dots}{0.4812 \dots}$$

不等式の係数は複雑ですが一次式です。そのため、X が10倍になれば、それを超えるフィボナッチ数列の項も10倍になります。

倍精度の浮動小数点は、10進数で15桁ほどの精度を持ちます。その程度までは、この式で正確な項数を求めることができます。

最後に

フィボナッチ数列の一般項を求めるのは、高校数学の履修範囲となっています（数学B「数列」）。この知識を使って、問題の計算量を $O(N)$ から $O(1)$ にすることができました。

今回で、当初の目標であった Project Euler の最初の25問を解くことができました。この25問の難易度を振り返ります。