PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 297 / ゼッケンドルフ表記

すべての正整数はフィボナッチ数列の連続しない項の和で一意に表せる。たとえば 100=3+8+89 である。 このような表し方はゼッケンドルフ表記と呼ばれる。

整数 n>0 に対して z(n) を n のゼッケンドルフ表記中の項数とする。

z(5)=1, z(14)=2, z(100)=3

0< n< 10^6 に対して Σz(n)=7894453 である。

0< n< 10^17 に対して Σz(n) を求めよ。

Problem 297 - Project Euler

実験

「Zeckendorf mathematica」でググると,mathematica による実装が何件もヒットします。stackexchange のコードを利用しました。

mathematica.stackexchange.com

Σz(n) も求められますが,10^6までで40秒以上かかります。

z(n) のスピードアップ

n のゼッケンドルフ表記は必ず n 以下の最大のフィボナッチ数(fib(i)とする)を含みます。すべての自然数はゼッケンドルフ表記できるので,z(n) を求めるのに n 以下の最大のフィボナッチ数をどんどん引いていっても大丈夫。証明は省略しますが,「連続しない項の和」の条件も問題ありません。

 z(n)=z(n-fib(i))+1

この方法で 10^6 までの和を求めると約4秒。10倍のスピードになりましたが,10^17を相手にするのはまだ無理です。

Σz(n) のスピードアップ

n=a から n=b までの z(n) の和を S(a, b)で表します。

フィボナッチ数までの和

 S(1,\, fib(i))
 =S(1,\, fib(i-1))+S(fib(i-1)+1,\, fib(i)-1)+z\left(fib(i)\right)

第2項の添字をずらします。
fib(i-1)+1 ≦ n ≦ fib(i)-1=fib(i-1)+fib(i-2)-1 のとき

 z(n)=z(n-fib(i-1))+1

このような n の個数は fib(i-2)-1 個。添字を fib(i-1) ずらすと n の範囲は 1 ~ fib(i-2)-1 に変わります。

 S(fib(i-1)+1,\, fib(i)-1)
 =S(1,\, fib(i-2)-1)+fib(i-2)-1
 =S(1,\, fib(i-2))+fib(i-2)-2

はじめの式にもどして  z\left(fib(i)\right)=1 を使うと

 S(1,\, fib(i))
 =S(1,\, fib(i-1))+S(1,\, fib(i-2))+fib(i-2)-1

非フィボナッチ数までの和

非フィボナッチ数 N 未満の最大のフィボナッチ数を fib(i) とします。

 S(1,\, N)=S(1,\, fib(i))+S\left(fib(i)+1,\, N\right)

第2項の添字をずらします。

 S(1,\, N)=S(1,\, fib(i))+S\left(1,\, N-fib(i)\right)+N-fib(i)

本番

上で導いた2式を両方使うのがスジですが,実は片方だけですみます。
たとえば 10^17-1 のゼッケンドルフ表記は次のようになっていて,N-fib(i) がフィボナッチ数になるのは最後の2のときだけです。これは他の数でも同様です。

「非フィボナッチ数までの和」の右辺第1項の添字を1つずらすと,添字がずっと非フィボナッチ数のまま計算できます。

 S(1,\, fib(i))=S(1,\, fib(i)-1)+z\left(fib(i)\right)
 =S(1,\, fib(i)-1)+1

 \therefore S(1,\, N)
 =S(1,\, fib(i)-1)+S\left(1,\, N-fib(i)\right)+N-fib(i)+1

この式を使うと 0.003 秒で答えが出ます。