引き返す遷移のある双六の期待回数

問題文など

昨日のくしらっちょさん*1のツイートで気になったので、時間計算量 $O(N)$ で解いてみました。怠惰人間なので余りを取っていないです。

今日の東海オンエアの動画で気になった。

3, 5の2種類の数字が等確率で出るサイコロを使ってNマスの双六をする。また、この双六ではあがるためにはピッタリゴールマスに着く必要があり、通り越したらその分だけ引き返す。

ゴールするためにサイコロを振る回数の期待値をmod998.... で求めよ。
— くしらっちょ (@kusirakusira) 2023年8月9日

解法

実装を簡単にするため、 $N$ が $3$ 以下の時は手計算で解き、 $N$ は $4$ 以上とします。

マス $N-4$ ～ $N$ のいずれかに初めて到達すれば、マス $N-4$ ～ $N$ の双六に帰着できて、そこからは連立方程式で解けます。動的計画法の漸化式など、細かい部分を合わせるのが難しくて嵌ってしまいました。

詳細は画像の通りです。1枚目の冒頭にあるように期待値が収束しないゲームを構成することもできるので、 $N$ について答えが収束することを本当は言わないといけないんですが、多分大丈夫でしょうということで（えぇ.....）

数学お強い方いらっしゃいましたら、収束性をご教示いただけますと（他力本願）

実装

import numpy as np
 
def solve(n: int) -> float:
    if n == 1:
        return 1
    if n == 2:
        return 4
    if n == 3:
        A = np.array([[0, -1/2, 1], [-1/2, 1, -1/2], [1, -1, 0]])
        A_inv = np.linalg.inv(A)
        b = np.array([1, 1, 1]).T
        return (A_inv @ b)[2]
    else:
        p = [0] * (n + 1)
        a = [0] * (n + 1)
        p[0] = 1
        a[0] = 0
        def fp(i: int) -> float:
            if 0 <= i and i <= n - 5:
                return p[i]
            return 0
        
        for i in range(1, n + 1):
            p[i] = 1/2 * fp(i - 3) + 1/2 * fp(i - 5)
            if p[i] == 0:
                a[i] = 0
            else:
                a[i] = (1/2 * fp(i - 3) * a[i - 3] + 1/2 * fp(i - 5) * a[i - 5]) / p[i] + 1
        mat_A = np.array([[-1, 0, 0, 1], [0, -1/2, 1, 0], [-1/2, 1, -1/2, 0], [1, -1/2, 0, -1/2]])
        A_inv = np.linalg.inv(mat_A)
        b = np.array([1, 1, 1, 1]).T
        E = [0] + list(A_inv @ b)
        A = [a[n], a[n - 1], a[n - 2], a[n - 3], a[n - 4]]
        P = [p[n], p[n - 1], p[n - 2], p[n - 3], p[n - 4]]
        ans = 0
        for i in range(5):
            ans += P[i] * (A[i] + E[i])
        return ans

if __name__ == '__main__':
    n = int(input())
    print(solve(n))

モンテカルロ法によるランダム試行

確率漸化式の問題は、とりあえずランダム試行をして当たりをつけると、デバッグなどがしやすいと思います。 $N$ 回試行すると、大体 $O(解の大きさ / \sqrt N)$ くらいの誤差に収束すると思います。（この辺あまり詳しくない）

#include <iostream>
#include <random>
#include <cstdio>
#define rep(i, n) for(i = 0; i < n; i++)
using namespace std;

mt19937 mt(252521);
int simurate(int n) {
    int pos = 0, cnt = 0;
    while (true) {
        cnt++;
        if (mt() % 2 == 0) pos += 3;
        else pos += 5;
        pos %= (2 * n);
        if (pos == n) break;
        else if (pos > n) { pos = 2 * n - pos; }
    }
    return cnt;
}

signed main() {
    int n, T = 1000000;
    cin >> n;

    double ans = 0;
    for (int t = 0; t < T; t++) {
        ans += simurate(n);
    }
    ans /= T;

    printf("%.14f\n", ans);
    return 0;
}

感想

ABC-Fにありそう（主観）。
解いたことがなかったので、面白かったです。（もっと計算量の小さい解法、あるのかなぁ？）

*1:し→ちに修正。ごめん。