数列を以下のように定義する。
- 0 ≤ k ≤ 1999 に対して gk = 1
- k ≥ 2000 に対して gk = g(k-2000) + g(k-1999)
k = 10^18 に対して gk mod 20092010 を求めよ。
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 を素因数分解します。
これらを PolynomialRemainder の Modulus に指定して，最後に ChineseRemainder でまとめます。似たような計算を4回行うことになりますが，意外に時間はかかりませんでした。
PolynomialRemainder は素数しか法に指定できないので，20092010 が各素数を1個ずつしか含まなくて助かりました。