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$ の定義から $lo{g}_{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問の難易度を振り返ります。