Project Euler 304 / Primonacci - 砂場で遊ぼう

問題概略
どう解くか
- 素数は NextPrime で処理
- フィボナッチ数は MatrixPowerMod で処理
実装
- 並列処理
- 前の項を使い回す

問題概略

$n$ より大きい最小の素数を $\verb|next_prime|(n)$ であらわす。
これを用いて $a_n$ を次のように定義する。

$a_1=\verb|next_prime|(10^{14})$

$a_n=\verb|next_prime| (a_{n-1}) \quad (n>1)$

また，次の式をみたすフィボナッチ数 $f(n)$ を用いて $b_n=f(a_n)$ と定義する。

$f(0)=0,\, f(1)=1$

$f(n)=f(n-1)+f(n-2) \quad (n>1)$

$\displaystyle\sum_{n=1}^{10^5} b_n=\sum_{n=1}^{10^5} f(a_n)$ を $\bmod\, 1234567891011$ で求めよ。
https://projecteuler.net/problem=304

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

どう解くか

この問題のポイントは「大きな素数を求めること」「大きなフィボナッチ数を求めること」です。これらをうまく処理できる言語にとっては，この問題は「やるだけ」です。

mathematica はこの条件をみたしていて，問題文をほぼそのままコードに落とすだけで解けました。

素数は NextPrime で処理

大きな素数には組み込み関数の NextPrime が使えます。例として $a_{100000}$ を求めてみました。

    In[]:= a[n_] := NextPrime[10^14, n];
        a[10^5]
    
    Out[]= 100000003235443

フィボナッチ数は MatrixPowerMod で処理

フィボナッチ数 $f_n$ は行列を使って求めます。
まず漸化式 $f_{n+2}=f_{n+1}+f_n$ を行列のかけ算に直します。

$\begin{pmatrix} f_{n+2}\\f_{n+1} \end{pmatrix} =\begin{pmatrix} f_{n+1}+f_n\\f_{n+1} \end{pmatrix} =\begin{pmatrix} 1&1\\1&0 \end{pmatrix} \begin{pmatrix} f_{n+1}\\f_{n} \end{pmatrix}$

$f_n$ は次のベクトルの第 1 成分です。

$\begin{pmatrix} f_{n}\\f_{n-1} \end{pmatrix} =\begin{pmatrix} 1&1\\1&0 \end{pmatrix}^{n-1} \begin{pmatrix} f_{1}\\f_{0} \end{pmatrix} =\begin{pmatrix} 1&1\\1&0 \end{pmatrix}^{n-1} \begin{pmatrix} 1\\0 \end{pmatrix}$

この計算を組み込み関数の MatrixPowerMod で処理すると，同時に $\bmod$ の計算もできて楽です。この関数は現行のマニュアルには載っていませんが，何の問題もなく使えます。StackExchange を参考にしました。

mathematica.stackexchange.com

$b_n$ の $\bmod$ をとったものを $c_n$ とおいて，次のように書きました。

    modnum = 1234567891011;
    mat = {{1, 1}, {1, 0}}; 
    c[n_] := Mod[First@#, modnum] &@
        (Algebra`MatrixPowerMod[mat, a@n - 1, modnum] . {1, 0});

実装

並列処理

上にあげたコードを使って計算したらちょっと遅かったので並列処理しました。計算時間は約 200 秒です。

    nmax = 10^5;
    a[n_] := NextPrime[10^14, n];
    modnum = 1234567891011;
    mat = {{1, 1}, {1, 0}}; 
    c[n_] := Mod[First@#, modnum] &@
        (Algebra`MatrixPowerMod[mat, a@n - 1, modnum] . {1, 0});
    ans = Mod[Total[ParallelMap[c, Range@nmax]], modnum]

前の項を使い回す

もう少し速くしたかったので， $b_n$ を毎回一から計算するかわりに $b_{n-1}$ を再利用する方法を試してみました。

\begin{align*}
\begin{pmatrix}
b_n\\ \square
\end{pmatrix}
=\begin{pmatrix}
1&1\\1&0
\end{pmatrix}^{a_n-1}
\begin{pmatrix}
1\\0
\end{pmatrix}\quad \therefore
\begin{pmatrix}
b_n\\ \square
\end{pmatrix}
=\begin{pmatrix}
1&1\\1&0
\end{pmatrix}^{a_n-a_{n-1}}
\begin{pmatrix}
b_{n-1}\\ \square
\end{pmatrix}
\end{align*}

計算回数はかなり減りますが，この方法だと並列化できないのでかなり遅くなります。

$n=10^3$ までの和をとる時間を比較しました。元のコードが約 4 秒で，新しく書いたコードは約 12 秒です。

結論として，mathematica では「できるだけ組み込み関数を使って，並列化する」のが最適解だと思います。

    nmax = 10^3;
    a[n_] := a[n] = NextPrime[10^14, n];
    modnum = 1234567891011;
    mat = {{1, 1}, {1, 0}};
    v[1] = Algebra`MatrixPowerMod[mat, a@1 - 1, modnum] . {1, 0};
    v[n_] := v[n] = 
      Mod[#, modnum] &@
        (Algebra`MatrixPowerMod[mat, a[n] - a[n - 1], modnum] . v[n - 1]);
    ans = First@Mod[Total[v /@ Range@nmax], modnum]