Project Euler 815

发表于 分类于 阅读次数： 本文字数： 2.5k 阅读时长 ≈ 2 分钟

Project Euler 815

题目

Group by Value

A pack of cards contains \(4n\) cards with four identical cards of each value. The pack is shuffled and cards are dealt one at a time and placed in piles of equal value. If the card has the same value as any pile it is placed in that pile. If there is no pile of that value then it begins a new pile. When a pile has four cards of the same value it is removed.

Throughout the process the maximum number of non empty piles is recorded. Let \(E(n)\) be its expected value. You are given \(E(2) = 1.97142857\) rounded to \(8\) decimal places.

Find \(E(60)\). Give your answer rounded to \(8\) digits after the decimal point.

解决方案

对每一种点数而言，在任意时刻它已经出现的张数只能是 \(0,1,2,3,4\) 之一。令\(C_i (i=0,1,2,3,4)\)表示当前“已出现 \(i\) 张”的点数种类数量。显然恒有守恒关系\(C_0+C_1+C_2+C_3+C_4=n.\)

桌面上的“非空堆”恰好对应那些出现次数为 \(1,2,3\) 的点数，因此此刻非空堆数量为\(P=C_1+C_2+C_3.\)

其中 \(C_4\) 表示已经完整出现 \(4\) 次的点数，它们的堆已被移除，不再贡献非空堆数。

设已发出 \(t\) 张牌，则剩余牌数为\(L=4n-t.\)若某点数已经出现了 \(i\) 张，则仍有 \(4-i\) 张留在牌堆中；因此“下一张牌来自已出现 \(i\) 张的点数”的总可能牌数为 \((4-i)C_i\)。

于是对于 \(i=0,1,2,3\) 可以得到下一张来自出现\(i\)张的点数的概率为\(\dfrac{(4-i)C_i}{L}\)，并且有一致性恒等式\(4C_0+3C_1+2C_2+C_3=L,\)这恰好保证上述四种情况概率和为 \(1\)。

抽到不同类别的点数，会使其出现次数加一，具体为：

  • \(0\to 1\)：\((C_0,C_1)\mapsto(C_0-1,C_1+1),\)
  • \(1\to 2\)：\((C_1,C_2)\mapsto(C_1-1,C_2+1),\)
  • \(2\to 3\)：\((C_2,C_3)\mapsto(C_2-1,C_3+1),\)
  • \(3\to 4\)（堆移除）：\((C_3,C_4)\mapsto(C_3-1,C_4+1).\)

从非空堆数 \(P=C_1+C_2+C_3\) 的角度看，变化规律为：

  • \(0\to 1\)：\(P\) 增加 \(1\)；
  • \(1\to 2,2\to 3\)：\(P\) 不变；
  • \(3\to 4\)：\(P\) 减少 \(1\)。

为了得到 \(M=\max \{P\}\) 的期望，必须记录历史最大堆数。令\(\displaystyle{m=\max_{0\le s\le t}\{P(s)\}}\)表示截至目前为止出现过的最大非空堆数。于是定义概率函数\(f_t(c_1,c_2,c_3,m)\)表示在发出 \(t\) 张牌后，系统满足\(C_1=c_1,C_2=c_2,C_3=c_3\)，且历史最大非空堆数为 \(m\)的概率。那么初始条件是\(f_0(0,0,0,0)=1.\)

为了使递推只依赖 \((c_1,c_2,c_3)\)，注意到发出的总牌数满足\(t=1\cdot C_1+2\cdot C_2+3\cdot C_3+4\cdot C_4.\)因此当已知 \((t,c_1,c_2,c_3)\) 时，\(C_4=\dfrac{t-(c_1+2c_2+3c_3)}{4}.\)进而由总数守恒得\(C_0=n-c_1-c_2-c_3-C_4.\)

在所有可达状态中，\(C_4\) 必为非负整数，且 \(C_0\ge 0\)。令当前非空堆数为\(p=c_1+c_2+c_3,\)剩余牌数为\(L=4n-t.\)

6. 概率递推方程

从状态 \((c_1,c_2,c_3,m)\) 出发，下一张牌有四种类型，其概率分别为\(\dfrac{4C_0}{L},\dfrac{3c_1}{L},\dfrac{2c_2}{L},\dfrac{c_3}{L}.\)

对应的状态转移为：

  1. 抽到一个“未出现过”的点数（\(0\to 1\)）：\((c_1,c_2,c_3,m)\mapsto(c_1+1,c_2,c_3,\max(m,p+1)).\)
  2. 抽到一个“已出现 1 张”的点数（\(1\to 2\)）：\((c_1,c_2,c_3,m)\mapsto(c_1-1,c_2+1,c_3,m).\)
  3. 抽到一个“已出现 2 张”的点数（\(2\to 3\)）：\((c_1,c_2,c_3,m)\mapsto(c_1,c_2-1,c_3+1,m).\)
  4. 抽到一个“已出现 3 张”的点数（\(3\to 4\)，堆移除）：\((c_1,c_2,c_3,m)\mapsto(c_1,c_2,c_3-1,m).\)

因此递推关系可以写成对所有合法参数的求和形式：

\[ f_{t+1}(c_1',c_2',c_3',m')=\sum_{(c_1,c_2,c_3,m)}f_t(c_1,c_2,c_3,m)\cdot \Pr\{(c_1,c_2,c_3,m)\to(c_1',c_2',c_3',m')\}. \]

更具体地，若把四类转移分别展开，则由上面的四种映射可直接累加得到 \(f_{t+1}\)。

当 \(t=4n\) 时，牌已全部发完，系统过程终止，历史最大堆数 \(m\) 已成为最终随机变量 \(M\) 的取值。因此期望为

\[ E(n)=\sum_{c_1,c_2,c_3,m} m\cdot f_{4n}(c_1,c_2,c_3,m). \]

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
from collections import defaultdict


N = 60

def expected_max_piles(n: int) -> float:
    inf = float('inf')
    f = defaultdict(float)
    f[(0, 0, 0, 0)] = 1.0
    for t in range(4 * n):
        left = 4 * n - t
        nf = defaultdict(float)
        for (c1, c2, c3, mp), p in f.items():
            x = t - (c1 + 2 * c2 + 3 * c3)
            if x < 0 or x & 3:
                continue
            c4 = x >> 2
            c0 = n - c1 - c2 - c3 - c4
            if c0 < 0:
                continue
            piles = c1 + c2 + c3
            if c0:
                q = (4 * c0) / left
                mp2 = max(mp, piles + 1)
                nf[(c1 + 1, c2, c3, mp2)] += p * q
            if c1:
                q = (3 * c1) / left
                nf[(c1 - 1, c2 + 1, c3, mp)] += p * q
            if c2:
                q = (2 * c2) / left
                nf[(c1, c2 - 1, c3 + 1, mp)] += p * q
            if c3:
                q = c3 / left
                nf[(c1, c2, c3 - 1, mp)] += p * q
        f = nf
    ans = 0.0
    for (_, _, _, mp), p in f.items():
        ans += mp * p
    return ans

print(f"{expected_max_piles(N):.8f}")