Project Euler 856

题目

Waiting for a Pair

A standard \(52\)-card deck comprises \(13\) ranks in four suits. A pair is a set of two cards of the same rank.

Cards are drawn, without replacement, from a well shuffled \(52\)-card deck waiting for consecutive cards that form a pair. For example, the probability of finding a pair in the first two draws is \(\dfrac{1}{17}\).

Cards are drawn until either such a pair is found or the pack is exhausted waiting for one. In the latter case we say that all \(52\) cards were drawn.

Find the expected number of cards that were drawn. Give your answer rounded to eight places after the decimal point.

解决方案

令\(N=13\)。可见，本题使用动态规划可以解决。不过，直接记录每个点数出现了几次会有巨大状态数。利用点数对称性，可以把抽牌过程的状态压缩为：\((r_0,r_1,r_2,r_3,r_4,l)\)

其含义为，\(r_i\)表示当前已经抽出的序列中，恰好出现 \(i\) 次的点数有多少个（\(i=0,1,2,3,4\)）；\(l\)则表示当前序列最后一张牌的点数，在序列中总共出现了多少次（取值 1~4）。注意，初始时还没有最后一张，用哨兵 \(l=4\) 表示（只做编码，不表示真的抽到过 4 次）。

不变量：\(\displaystyle r_i\ge 0, \sum_{i=0}^4 r_i=N.\)已抽牌数与剩余牌数：\(\displaystyle u=\sum_{i=0}^4 i\cdot r_i, r=4N-u.\)

我们定义\(F(r_0,r_1,r_2,r_3,r_4,l)\)为从该状态出发，直到停止，还要再抽多少张牌”的期望。

考虑下一张牌来自当前出现 \(i\) 次的某个点数（\(i=0,1,2,3\)），那么记其后继状态将会便化成\(M(r_0,r_1,r_2,r_3,r_4,l,i)\)。这类点数每个还剩 \(4-i\) 张可抽。 若 \(i=l\)，其中恰有 1 个点数是上一张的点数，抽到它会立刻形成相邻对子并停止，贡献概率 \((4-i)/r\)，本步只抽 1 张；除它之外，其余 \(r_i-1\) 个点数可继续游戏。若 \(i\ne l\)，这组内所有 \(r_i\) 个点数都可继续游戏。

设\(t_i=r_i-\mathbf 1\{i=l\}\)，则继续分支的可抽牌数是\((4-i)t_i\)。若 \(t_i>0\)，抽到后状态变为：\(r_i\leftarrow r_i-1, r_{i+1}\leftarrow r_{i+1}+1, l\leftarrow i+1.\)

最终，若 \(r=0\)，显然 \(F=0\)。否则：

\[ F(r_0,r_1,r_2,r_3,r_4,l)=\mathbf 1\{0\le l\le 3\}\frac{4-l}{r}+\sum_{i=0}^{3}\frac{(4-i)t_i}{r} \left(1+F_{M(r_0,r_1,r_2,r_3,r_4,l,i)}\right) \]

其中第一项就表示下一张立刻停止；对第二项，抽 1 张只对 \(t_i>0\) 的项求和。

目标即：\(F(N,0,0,0,0,4).\)

代码

from fractions import Fraction
from functools import lru_cache

N = 13
@lru_cache(None)
def dp(r0, r1, r2, r3, r4, last):
    used = r1 + 2 * r2 + 3 * r3 + 4 * r4
    rem = 52 - used
    if rem == 0:
        return Fraction(0, 1)
    ans = Fraction(0, 1)
    arr = [r0, r1, r2, r3, r4]
    for i in range(4):
        t = arr[i]
        if last == i:
            if t == 0:
                continue
            ans += Fraction(4 - i, rem)
            t -= 1
        if t > 0:
            ways = (4 - i) * t
            nxt = arr[:]
            nxt[i] -= 1
            nxt[i + 1] += 1
            ans += Fraction(ways, rem) * (1 + dp(nxt[0], nxt[1], nxt[2], nxt[3], nxt[4], i + 1))
    return ans

exact = dp(N, 0, 0, 0, 0, 4)
print(f"{exact:.8f}")