PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 258 / ラグ付フィボナッチ数列

数列を以下のように定義する。

  • 0 ≤ k ≤ 1999 に対して gk = 1
  • k ≥ 2000 に対して gk = g(k-2000) + g(k-1999)

k = 10^18 に対して gk mod 20092010 を求めよ。

Problem 258 - Project Euler


10^18 という大きな添字,20092010という大きな法をどうするかがポイントです。添字は特性方程式で処理し,法は中国式剰余定理で処理しました。

■10^18 の処理

特性方程式 x^(2000)-x-1=0 を使って x^(10^18) を次数下げします。普通に

f[x_]:=PolynomialRemainder[x^(10^18), x^2000 - x - 1];
f[1]

とすると,かなり時間がかかった後で「再帰の回数が限界を超えた」エラーが出てしまいます。ここでしばらくハマりましたが,回避法がわかりました。

  • 係数がかなり大きくなるので,Modulus を指定してあらかじめ割っておく
  • 代入は「/. x->1」を使った ReplaceAll で行う
■20092010 の処理

法の 20092010 を素因数分解します。

 20092010=2\cdot 5\cdot 859\cdot 2339

これらを PolynomialRemainder の Modulus に指定して,最後に ChineseRemainder でまとめます。似たような計算を4回行うことになりますが,意外に時間はかかりませんでした。

PolynomialRemainder は素数しか法に指定できないので,20092010 が各素数を1個ずつしか含まなくて助かりました。