Project Euler 329 / 素数カエル

問題概略

素数カエルは 1 から 500 までの番号がついた 500 個の正方形の上を跳びまわる。
左右のいずれかの正方形に等確率で跳ぶことができ，正方形外に跳び出すことはない。
1 か 500 に着地した場合，次は自動的に移動可能な正方形に跳ぶ。

• 素数が書かれた正方形にいる場合，次の正方形に跳ぶ直前にカエルは の確率で「P」 (PRIME) と鳴き， の確率で「N」 (NOT PRIME) と鳴く。
• 素数でない数が書かれた正方形にいる場合，次の正方形に跳ぶ直前にカエルは の確率で「P」と鳴き， の確率で「N」と鳴く。

カエルの開始位置はすべての正方形に対し等確率でランダムである。最初の 15 回の鳴き声が PPPPNNPPPNPPNPN になる確率を求めよ。

https://projecteuler.net/problem=329

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

drive.google.com

条件つき確率

まずは誤答例を紹介します。
以下，「PPPPNNPPPNPPNPN」と鳴くことを「正しく鳴く」ということにして，事象 が起きる確率を で表します。

\begin{align*}
P(\text{$n$回目に正しく鳴く}) &=\sum_{k=1}^{500} P(\text{$n$回目に$k$番で正しく鳴く})\\
&=\sum_{k=1}^{500} P(\text{$n$回目に$k$番にいる})P(\text{$k$番で正しく鳴く})
\end{align*}

「求める確率はこれの積」とやると不正解。実際に計算するとどうなるかお見せしましょう。
分母，分子の桁数がえらいことになります。言語は mathematica です。

    In[]:= AbsoluteTiming[
        p[1, k_] := 1/500;
        p[n_, 1] := p[n, 1] = 1/2 p[n - 1, 2];
        p[n_, 2] := p[n, 2] = p[n - 1, 1] + 1/2 p[n - 1, 3];
        p[n_, 499] := p[n, 499] = p[n - 1, 500] + 1/2 p[n - 1, 498];
        p[n_, 500] := p[n, 500] = 1/2 p[n - 1, 499];
        p[n_, k_] := p[n, k] = 1/2 (p[n - 1, k - 1] + p[n - 1, k + 1]);
        croaks = Characters@"PPPPNNPPPNPPNPN";
        q[n_, k_] := q[n, k] =
         If[croaks[[n]] == "P", If[PrimeQ@k, 2/3, 1/3], 
          If[PrimeQ@k, 1/3, 2/3]];
        r[n_] := Sum[p[n, k]*q[n, k], {k, 1, 500}];
        ans = Apply[Times, r /@ Range@15]]
       
    Out[]= {0.0472912, \
       6146864190777269715907477114280310734534525885000007414638209/\
       793176230249158201524722073600000000000000000000000000000000000000}

これが不正解になるのは，各回の鳴き方は独立ではないためです。

このカエルは確率 で正しく鳴いて，確率 で間違えます。
「P」と鳴いたら素数にいる確率が高く，移動先は高い確率で非素数です。
これは次の鳴き方に影響しますよね。条件つき確率で考えないとダメです。

メモ化再帰

「前の回にどう鳴いたか」≒「前の回にどこにいたか」を考えないといけないので，和のとり方を変えます。

たとえば 5 番からスタートする場合，「P」と鳴いたら 4 番か 6 番に確率 で移動します。
番からスタートして正しく鳴く確率を とおいて， 番で正しく鳴く確率を とおきます。おおまかには次のようになります。

実際にはそれぞれの場所で「PPPPNNPPPNPPNPN」の何番目の鳴き声をあげるかを指定しないといけません。
と の引数を1つ増やしましょう。

この数字が 15 になったら終わり，と言いたいところですが，まだダメです。
を「 番からスタートして～」と定義すると，「開始位置はランダム」という設定のせいで毎回余計な 倍をすることになってしまいます。

の設定を「 番からスタートして正しく鳴く確率の500倍」に変えて，最後に 500 で割るようにします。

\begin{align*}
&p(n,\,k)=q(n,\, k)\cdot \frac{1}{2}\left\{p(n-1,\,k+1)+p(n+1,\, k+1)\right\}\quad (2\leqq n\leqq 499)\\[3pt]
&p(1,\,k)=q(1,\, k)\cdot p(2,\,k+1)\\
&p(500,\,k)=q(500,\, k)\cdot p(499,\, k+1)\\
&p(n,\, 16)=\frac{1}{500}
\end{align*}

実際に書いたコードは次の通りです。*1

    In[]:= AbsoluteTiming[
        croaks = Characters@"PPPPNNPPPNPPNPN";
        q[n_, k_] := q[n, k] =
         If[croaks[[k]] == "P", If[PrimeQ@n, 2/3, 1/3], 
         If[PrimeQ@n, 1/3, 2/3]];
        p[n_, 16] := p[n, 16] = 1/500;
        p[1, k_] := p[1, k] = q[1, k]*p[2, k + 1];
        p[500, k_] := p[500, k] = q[500, k]*p[499, k + 1];
        p[n_, k_] := 
         p[n, k] = q[n, k]*1/2 (p[n + 1, k + 1] + p[n - 1, k + 1]);
        ans = Sum[p[k, 1], {k, 1, 500}];]

    Out[]= {0.0403924, Null}

python で解く

同じことを python でもやりました。計算時間はほぼ同じでした。

    import time
    from sympy import isprime
    from fractions import Fraction
    from functools import lru_cache
    
    @lru_cache(maxsize=None)
    def q(n, k):
        croaks = list("PPPPNNPPPNPPNPN")
        if croaks[k - 1] == "P":
            return Fraction(2, 3) if isprime(n) else Fraction(1, 3)
        else:
            return Fraction(1, 3) if isprime(n) else Fraction(2, 3)
    
    @lru_cache(maxsize=None)
    def p(n, k):
        if k == 16:
            return Fraction(1, 500)
        if n == 1:
            return q(1, k) * p(2, k + 1)
        if n == 500:
            return q(500, k) * p(499, k + 1)
        return q(n, k) * Fraction(1, 2) * (p(n + 1, k + 1) + p(n - 1, k + 1))
    
    def solve():
        return sum([p(k, 1) for k in range(1, 501)])
    
    if __name__ == '__main__':
        start = time.time()
        print(solve())
        elapsed_time = time.time() - start
        print("elapsed_time: {0} [sec]".format(elapsed_time))

*1:101問以降なので答えの数値は省略。