PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 624 / Two heads are better than one

コイン投げで表が連続して2回出るまでに M 回かかるとする。M が n で割り切れる確率を P(n) であらわす。
P(2)=3/5, P(3)=9/31 である。

有理数 a/b と素数p に対して a≡bq (mod p) をみたす最小の自然数 q を Q(a/b, p) であらわす。

たとえば Q(P(2), 109)=Q(3/5, 109)=66 である。
5x66 =330 ≡ 3 (mod 109)
同様に Q(P(3), 109)=46 である。

Q(P(10^18), 10^9+9) を求めよ。

本家:Problem 624 - Project Euler

下準備

M-1 回目までは「裏は連続していいが,表は連続しない」というルールで表と裏が並びます。これは入試数学でもよくある設定で,フィボナッチ数 f[M-1] 通りあります(証明略)。

「M が n で割り切れる」=「M は n の倍数」から M=kn とおけて,P(n) は次のように表せます。

 P(n)=\displaystyle\sum_{k=1}^{\infty} \frac{f_{kn-1}}{2^{kn}}

フィボナッチ数列の一般項を使ってこの和を求め,P(2)~P(10) を求めました。

この形の P(n) で計算できる Q はせいぜい n=100 程度までなので,もっとうまい方法で P(n) を求めないといけません。

行列を使う

上で求めた P(n) は分母,分子とも指数関数です。底は

  • 分子は -1, 1-√5, 1+√5
  • 分母は -1, 4, 1-√5, 1+√5

これらを解にもつ特性方程式を作って,対応する漸化式を求めました。
P(n)=a(n)/b(n) とおきます。

 a_1=1,\, a_2=3,\, a_3=9
 a_{n+3}=a_{n+2}+6 a_{n+1}+4 a_{n}

 b_1=1,\, b_2=5,\, b_3=31,\, b_4=145
 b_{n+4}=5 b_{n+3}+2 b_{n+2}-20 b_{n+1}-16 b_{n}

これを行列に直して P(n) の分母,分子を求めます。

a≡bq (mod p) は

 q =( (a \bmod p) \cdot (b^{-1} \bmod p) \bmod p)

でOK。あっさり終わります*1

おまけ

MatrixPowerMod は昔は普通に使えたような気がしますが,現行バージョンのマニュアルには載っていません。stackexchange に助けられました。

mathematica.stackexchange.com

*1:最近の問題なので答えは載せないことにします