Project Euler 339

发表于 更新于 分类于 阅读次数： 本文字数： 4.3k 阅读时长 ≈ 4 分钟

Project Euler 339

题目

Peredur fab Efrawg

“And he came towards a valley, through which ran a river; and the borders of the valley were wooded, and on each side of the river were level meadows. And on one side of the river he saw a flock of white sheep, and on the other a flock of black sheep. And whenever one of the white sheep bleated, one of the black sheep would cross over and become white; and when one of the black sheep bleated, one of the white sheep would cross over and become black.”

en.wikisource.org

Initially each flock consists of \(n\) sheep. Each sheep (regardless of colour) is equally likely to be the next sheep to bleat. After a sheep has bleated and a sheep from the other flock has crossed over, Peredur may remove a number of white sheep in order to maximize the expected final number of black sheep. Let \(E(n)\) be the expected final number of black sheep if Peredur uses an optimal strategy.

You are given that \(E(5) = 6.871346\) rounded to \(6\) places behind the decimal point.

Find \(E(10 000)\) and give your answer rounded to \(6\) places behind the decimal point.

解决方案

把“河这边白羊、河那边黑羊”的信息抽象为一个二元状态：

  • 当前黑羊数为 \(b\)
  • 当前白羊数为 \(w\)

每次“咩叫”发生时，从全部 \(b+w\) 只羊里等概率选出下一只咩叫的羊：

  • 以概率 \(\dfrac{b}{b+w}\)：黑羊咩叫，则有一只白羊渡河变黑，\((b,w)\to(b+1,w-1).\)
  • 以概率 \(\dfrac{w}{b+w}\)：白羊咩叫，则有一只黑羊渡河变白，\((b,w)\to(b-1,w+1).\)

每次渡河变色之后，Peredur 可以立刻移走若干只白羊（只能减少 \(w\)），以最大化最终黑羊数的期望。过程在某一侧羊群为空时结束；终局黑羊数就是当时的 \(b\)。

令\(V(b,w)\)表示在状态 \((b,w)\)、且 下一次咩叫尚未发生之前，Peredur 采用最优策略能获得的“最终黑羊数期望”。边界显然为\(V(b,0)=b, V(0,w)=0.\)

在不移走白羊的前提下，从 \((b,w)\) 经过一次咩叫与变色后，黑羊数的期望变化为 \[ \mathbb E[b'] = \frac{b}{b+w}(b+1)+\frac{w}{b+w}(b-1) = b+\frac{b-w}{b+w}. \]

因此：

  • 若 \(w\ge b\)，则 \(\mathbb E[b']\le b\)，黑羊“没有正向漂移”，局面对黑羊不利。
  • 若 \(w\le b-1\)，则 \(\mathbb E[b']\ge b+\dfrac{1}{2b-1}\)，黑羊有正向漂移。

Peredur 唯一能做的就是减少 \(w\)。为了让局面具备正向漂移，至少应确保\(w\le b-1.\)另一方面，在满足 \(w\le b-1\) 的前提下，保留更多的白羊不会降低黑羊的当步优势（仍然是 \(w<b\)），并且保留的白羊未来仍有机会通过“黑羊咩叫”变成黑羊，因而不应额外移走。

这导出一个非常稳定的最优策略结构：若 \(w\ge b\)，最优是立刻移走白羊使 \(w=b-1\)；若 \(w<b\)，最优是不移走任何白羊。 下面的计算都建立在这一策略上；在此策略下，除了起始的 \((n,n)\) 之外，系统会被“钳制”在区域 \(w\le b-1\) 中运行，并且只在“刚好触碰到 \(w=b\)（或更大）”时触发一次性移走到 \(w=b-1\)。

在最优策略下，一个非常重要的族是\(g(k)=V(k,k-1),(k\ge 1),\)即“黑羊比白羊多 \(1\)”的临界状态。

当处于 \((k,k-1)\) 时：

  • 以概率 \(\dfrac{k}{2k-1}\) 黑羊咩叫：\((k,k-1)\to(k+1,k-2).\)此时仍满足 \(w<b\)，不会移走白羊。
  • 以概率 \(\dfrac{k-1}{2k-1}\) 白羊咩叫：\((k,k-1)\to(k-1,k).\)这时 \(w\ge b\)，按策略立刻移走白羊到 \(w=b-1=k-2\)，得到\((k-1,k)\to(k-1,k-2),\)其值为 \(g(k-1)\)。

因此有

\[ g(k)=\frac{k}{2k-1}V(k+1,k-2)+\frac{k-1}{2k-1}g(k-1). \tag{1} \]

接下来要把 \(V(k+1,k-2)\) 用 \(g(k)\) 表示出来。

考虑总羊数固定为 \(2k-1\) 的那条线上的状态\((k+d,k-1-d), d=0,1,\dots,k-1,\)其中 \(d=k-1\) 对应吸收态 \((2k-1,0)\)，其价值为\(V(2k-1,0)=2k-1.\)

对 \(d\ge 1\) 的这些状态，下一步不会触发“移走白羊”（因为 \(w<b\) 仍保持），所以在这条线内部它是一个单纯的生灭链。特别地，从\((k+1,k-2)(d=1)\)出发，要么先回到 \(d=0\)（即回到 \((k,k-1)\)，价值为 \(g(k)\)），要么先到达吸收端 \(d=k-1\)（即 \((2k-1,0)\)，价值为 \(2k-1\)）。

接下来只关心从状态 \(d=1\) 出发时，这条生灭链最终会先走到哪一端。把“在到达 \(d=0\) 之前，先到达右端吸收态 \(d=k-1\)”的概率记为 \(p_k\)。换句话说，\(p_k\) 就是从 \((k+1,k-2)\) 起步，最终先抵达 \((2k-1,0)\) 而不是回到 \((k,k-1)\) 的可能性。

那么由“调和函数=边界值的命中概率线性组合”得到

\[ V(k+1,k-2)=(1-p_k)\cdot g(k)+p_k\cdot (2k-1). \tag{2} \]

现在只需计算 \(p_k\)。 在 \(d=1,2,\dots,k-2\) 处，向上（\(d\to d+1\)）与向下（\(d\to d-1\)）的转移概率分别是

  • 向上概率：\(\alpha_d=\dfrac{k+d}{2k-1},\)
  • 向下概率：\(\beta_d=\dfrac{k-1-d}{2k-1}.\)

令 \(r_d\) 为“从 \(d\) 出发先到 \(k-1\) 而非 \(0\)”的概率，则满足\(r_0=0, r_{k-1}=1,r_d=\alpha_d r_{d+1}+\beta_d r_{d-1}\ (1\le d\le k-2).\)

令差分 \(s_d=r_d-r_{d-1}\)，代入可得经典化简：\(\displaystyle{s_{d+1}=\frac{\beta_d}{\alpha_d}s_d=\frac{k-1-d}{k+d}s_d.}\)因此\(\displaystyle{s_{d} = s_1\prod_{t=1}^{d-1}\frac{k-1-t}{k+t}.}\)并且因为\(\displaystyle{1=r_{k-1}-r_0=\sum_{d=1}^{k-1}s_d=s_1\left(1+\sum_{j=1}^{k-2}\prod_{t=1}^{j}\frac{k-1-t}{k+t}\right).}\) 所以

\[ r_1=p_k =\frac{1}{ 1+\sum_{j=1}^{k-2}\prod_{t=1}^{j}\frac{k-1-t}{k+t} }. \tag{3} \]

把式 \((2)\) 代回 \((1)\)：\(g(k)=\dfrac{k}{2k-1}((1-p_k)g(k)+p_k(2k-1))+\dfrac{k-1}{2k-1}g(k-1).\)移项解出 \(g(k)\)，可整理为

\[ g(k)=\frac{(2k-1)kp_k+(k-1)g(k-1)}{(k-1)+kp_k}. \tag{4} \]

初值：\(g(1)=V(1,0)=1.\)由 \((3)(4)\) 就能自底向上算出 \(g(1),g(2),\dots,g(n)\)。

题目要求从初始 \((n,n)\) 出发。第一次咩叫发生前不能移走白羊：

  • 以概率 \(\dfrac12\) 黑羊咩叫：\((n,n)\to(n+1,n-1),\)其期望为 \(V(n+1,n-1)\)。
  • 以概率 \(\dfrac12\) 白羊咩叫：\((n,n)\to(n-1,n+1),\)随后触发策略移走白羊到 \(w=b-1\)，即移走 \(3\) 只白羊得到\((n-1,n+1)\to(n-1,n-2),\)其期望为 \(g(n-1)\)。

因此 \[ E(n)=\frac12V(n+1,n-1)+\frac12g(n-1). \tag{5} \]

剩下只需计算 \(V(n+1,n-1)\)。在总数固定为 \(2n\) 的线上定义\(u_d=V(n+d,n-d), d=1,2,\dots,n,\)其中 \(u_n=V(2n,0)=2n\)。

对 \(2\le d\le n-1\)，一步递推（不触发移走）给出三项关系

\[ u_d=\frac{n+d}{2n}u_{d+1}+\frac{n-d}{2n}u_{d-1}. \tag{6} \]

而对 \(d=1\) 的状态 \((n+1,n-1)\)，若白羊咩叫则到 \((n,n)\) 并立刻移走到 \((n,n-1)\)，其值为 \(g(n)\)；若黑羊咩叫则到 \((n+2,n-2)\) 即 \(u_2\)，于是 \[ u_1=\frac{n+1}{2n}u_2+\frac{n-1}{2n}g(n). \tag{7} \]

联立 \((6)(7)\) 与边界 \(u_n=2n\)，得到一个三对角线性方程组，可用 Thomas 算法解三对角线性方程组 在 \(O(n)\)时间内解出 \(u_1\)，再由 \((5)\) 得到 \(E(n)\)。

代码

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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#include <bits/stdc++.h>
using namespace std;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const int n = 10000;

    // g[k] = V(k, k-1) for k>=1
    vector<double> g(n + 1, 0.0);
    g[1] = 1.0;

    for (int k = 2; k <= n; k++) {
        // p_k = 1 / (1 + sum_{j=1}^{k-2} prod_{t=1}^j (k-1-t)/(k+t))
        double S = 1.0;       // includes j=0 term
        double ratio = 1.0;   // running product
        for (int t = 1; t <= k - 2; t++) {
            ratio *= (double)(k - 1 - t) / (double)(k + t);
            if (ratio < 1e-320) break; // negligible tail
            S += ratio;
        }
        double pk = 1.0 / S;

        // g(k) recurrence
        double num = (double)(2 * k - 1) * (double)k * pk + (double)(k - 1) * g[k - 1];
        double den = (double)(k - 1) + (double)k * pk;
        g[k] = num / den;
    }

    // Solve even-total line for u1 = V(n+1, n-1)
    // Unknowns: u_d for d=1..n-1, with boundary u_n = 2n
    int m = n - 1;
    vector<double> a(m + 1, 0.0), b(m + 1, 1.0), c(m + 1, 0.0), d(m + 1, 0.0);

    double twon = 2.0 * n;

    // d=1 equation: u1 - (n+1)/(2n) u2 = (n-1)/(2n) g(n)
    c[1] = -(double)(n + 1) / twon;
    d[1] = (double)(n - 1) / twon * g[n];

    // 2 <= d <= n-2
    for (int i = 2; i <= m - 1; i++) {
        a[i] = -(double)(n - i) / twon;
        c[i] = -(double)(n + i) / twon;
        d[i] = 0.0;
    }

    // d = n-1 equation uses u_n = 2n:
    // u_{n-1} - (n+(n-1))/(2n) u_n - (n-(n-1))/(2n) u_{n-2} = 0
    // => a[m] = -1/(2n), rhs = (2n-1)/(2n)*(2n) = 2n-1
    a[m] = -1.0 / twon;
    d[m] = 2.0 * n - 1.0;

    // Thomas forward elimination
    for (int i = 2; i <= m; i++) {
        double w = a[i] / b[i - 1];
        b[i] -= w * c[i - 1];
        d[i] -= w * d[i - 1];
    }

    // Back substitution
    vector<double> x(m + 1, 0.0);
    x[m] = d[m] / b[m];
    for (int i = m - 1; i >= 1; i--) {
        x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
    }

    double u1 = x[1];

    // Final answer
    double ans = 0.5 * u1 + 0.5 * g[n - 1];

    cout.setf(std::ios::fixed);
    cout << setprecision(6) << ans << "\n";
    return 0;
}