Project Euler 584

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

Project Euler 584

题目

Birthday Problem Revisited

A long long time ago in a galaxy far far away, the Wimwians, inhabitants of planet WimWi, discovered an unmanned drone that had landed on their planet. On examining the drone, they uncovered a device that sought the answer for the so called “Birthday Problem”. The description of the problem was as follows:

If people on your planet were to enter a very large room one by one, what will be the expected number of people in the room when you first find \(\mathbf{3}\) people with Birthdays within \(\mathbf{1}\) day from each other.

The description further instructed them to enter the answer into the device and send the drone into space again. Startled by this turn of events, the Wimwians consulted their best mathematicians. Each year on Wimwi has \(10\) days and the mathematicians assumed equally likely birthdays and ignored leap years (leap years in Wimwi have \(11\) days), and found \(5.78688636\) to be the required answer. As such, the Wimwians entered this answer and sent the drone back into space.

After traveling light years away, the drone then landed on planet Joka. The same events ensued except this time, the numbers in the device had changed due to some unknown technical issues. The description read:

If people on your planet were to enter a very large room one by one, what will be the expected number of people in the room when you first find \(\mathbf{3}\) people with Birthdays within \(\mathbf{7}\) days from each other.

With a \(100\)-day year on the planet, the Jokars (inhabitants of Joka) found the answer to be \(8.48967364\) (rounded to \(8\) decimal places because the device allowed only \(8\) places after the decimal point) assuming equally likely birthdays. They too entered the answer into the device and launched the drone into space again.

This time the drone landed on planet Earth. As before the numbers in the problem description had changed. It read:

If people on your planet were to enter a very large room one by one, what will be the expected number of people in the room when you first find \(\mathbf{4}\) people with Birthdays within \(\mathbf{7}\) days from each other.

What would be the answer (rounded to eight places after the decimal point) the people of Earth have to enter into the device for a year with \(365\) days? Ignore leap years. Also assume that all birthdays are equally likely and independent of each other.

解决方案

简而言之，题目是说一年共有 \(Y=365\) 天，生日均匀独立分布在 \({0,1,\dots,Y-1}\) 上，并按圆环距离定义相隔不超过 \(w=7\) 天。

题目要求的触发条件是：存在 \(K=4\) 个人，他们的生日都落在某个长度为 \(w+1\) 的连续区间内（按圆环取模）。

记随机变量 \(N\) 为第一个触发条件出现时，房间里的人数，求 \(\mathbb E[N]\)。

设事件 \(B_n\) 表示“前 \(n\) 个人中尚未出现触发条件”，则

  • \(\mathbb P(N>n)=\mathbb P(B_n)\)
  • \(\mathbb P(N=n)=\mathbb P(B_{n-1})-\mathbb P(B_n)\)

于是\(\displaystyle{\mathbb E[N]=\sum_{n\ge1} n\big(\mathbb P(B_{n-1})-\mathbb P(B_n)\big)= \sum_{n\ge0}\mathbb P(B_n).}\)因此问题化为：计算所有 \(n\) 的 \(\mathbb P(B_n)\) 并求和。

令\(\displaystyle{x_0,x_1,\dots,x_{Y-1}\ge 0, \sum_{i=0}^{Y-1}x_i=n,}\)其中 \(x_i\) 表示第 \(i\) 天生日的人数。在给定计数向量 \(x\) 时，其概率为多项分布：

\[ \mathbb P(x)=\frac{n!}{Y^n\prod_{i=0}^{Y-1} x_i!}. \]

因此，触发条件：存在某个长度 \(w+1\) 的连续区间内至少有 \(K\) 人可以转换成另一种表达：对所有 \(j\)（按模 \(Y\)），都有\(\displaystyle{\sum_{t=0}^w x_{(j+t)\%Y}\le K-1}\)。

因此

\[ \mathbb P(B_n)= \sum_{\substack{x_0+\cdots+x_{Y-1}=n\\ \forall j:\sum_{t=0}^{w}x_{(j+t)\%Y}\le K-1}} \frac{n!}{Y^n\prod_{i=0}^{Y-1} x_i!}. \]

代入 \(\displaystyle{\mathbb E[N]=\sum_{n\ge0}\mathbb P(B_n)}\) 得

\[ \mathbb E[N]= \sum_{\substack{x_0,\dots,x_{Y-1}\ge 0\\ \forall j:\sum_{t=0}^{w}x_{(j+t)\%Y}\le K-1}} \frac{\big(\sum_{i=0}^{Y-1}x_i\big)!}{Y^{\sum_{i=0}^{Y-1}x_i}\prod_{i=0}^{Y-1} x_i!}. \]

为了简便计算，我们使用伽马函数把阶乘去除。可见，有\(\displaystyle{n! = \int_{0}^{\infty} e^{-u}u^ndu,}\)则每个配置的贡献\(\displaystyle{\frac{n!}{Y^n\prod_{i=0}^{Y-1} x_i!}=\int_{0}^{\infty} e^{-u}\frac{(u/Y)^n}{\prod_{i=0}^{Y-1} x_i!}du.}\)

交换求和与积分，得到\(\displaystyle{\mathbb E[N]=\int_{0}^{\infty} e^{-u}H\left(\frac{u}{Y}\right)du}\)，其中\(\displaystyle{H(z)=\sum_{\substack{x_0,\dots,x_{Y-1}\ge 0\\\forall j:\sum_{t=0}^{w}x_{j+t}\le K-1}}\frac{z^{\sum_{i=0}^{Y-1} x_i}}{\prod_{i=0}^{Y-1} x_i!}.}\)

所以核心变成：对给定 \(z\) 高效计算 \(H(z)\)。

令\(L=w+1, B=K-1.\)约束是“任意长度 \(L\) 的连续窗口总和 \(\le B\)”。定义状态为最近 \(L-1\) 天的计数：\(\displaystyle{s=(a_1,a_2,\dots,a_{L-1}), a_i\ge 0, \sum_{i=1}^{L-1}a_i\le B.}\)

如果下一天放入 \(t\) 个人，则新状态为\(s'=\text{shift}(s,t)=(a_2,a_3,\dots,a_{L-1},t),\)并且必须满足窗口约束\(\displaystyle{\sum_{i=1}^{L-1}a_i+t\le B.}\)

每次添加 \(t\) 贡献对应权重：\(\dfrac{z^t}{t!}.\)因此形成一个转移矩阵 \(T(z)\)（大小为状态数 \(S\)）：

\[ T(z)_{s',s}= \begin{cases} \dfrac{z^t}{t!} &\text{if} \quad s'=\text{shift}(s,t) \land \displaystyle{\sum_{i=1}^{L-1} a_i}+t\le B\\ 0,&\text{else}. \end{cases} \]

然而，上面的转移只保证了线性序列内部的所有长度 \(L\) 窗口都合法。

但一年是圆环，还需要检查跨越边界的窗口，即：

把开头的 \(L-1\) 天记作起始状态 \(p\)，末尾的 \(L-1\) 天记作结束状态 \(e\)。

对于每个 \(j=1,2,\dots,L-1\)，跨边界窗口包含：

  • 结束状态 \(e\) 的最后 \(j\) 项
  • 起始状态 \(p\) 的最前 \(L-j\) 项

因此必须满足\(\displaystyle{\sum_{t=L-j}^{L-1} e_t+\sum_{t=1}^{L-j} p_t \le B.}\)记兼容矩阵\(C_{p,e}=1\)当且仅当\((p,e)\)满足所有跨边界窗口约束，否则\(C_{p,e}=0\)。

固定起始状态 \(p\)，其自身对应前 \(L-1\) 天的权重为\(\displaystyle{W_p(z)=\frac{z^{\sum_{t=1}^{L-1} p_t}}{\prod_{t=1}^{L-1} p_t!}.}\)

剩余的 \(Y-(L-1)\) 天用矩阵推进。设\(P(z)=T(z)^{Y-(L-1)}.\)则从 \(p\) 走到 \(e\) 的总权重为 \(P(z)_{e,p}\)。

最后乘上首尾兼容性 \(C_{p,e}\) 并对所有 \(p,e\) 求和：

\[ H(z)=\sum_{p}\sum_{e} W_p(z)P(z)_{e,p}C_{p,e}. \]

最终，我们需要 \[ \mathbb E[N]=\int_{0}^{\infty} e^{-u}H\left(\frac{u}{Y}\right)du. \]

Gauss–Laguerre 公式直接针对权重 \(e^{-u}\)：\(\displaystyle{\int_0^\infty e^{-u} f(u),du \approx \sum_{i=1}^{m} w_i f(x_i),}\)所以只需在若干节点 \(x_i\) 处计算 \(H(x_i/Y)\) 并加权求和即可。

代码

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
from tools import fac
import numpy as np


def build_states(m: int, bound: int):
    """All m-tuples (a0..a_{m-1}) with ai in [0,bound] and sum <= bound."""
    states = []
    for s in range(bound + 1):
        def rec(pos, rem, cur):
            if pos == m:
                if rem == 0:
                    states.append(tuple(cur))
                return
            for x in range(rem + 1):
                cur.append(x)
                rec(pos + 1, rem - x, cur)
                cur.pop()

        rec(0, s, [])
    return states


def prepare(K: int, w: int):
    """
    K: required people in window (K=4)
    w: within w days => window length L=w+1 (w=7 => L=8)
    Y: days in year (365)
    """
    bound = K - 1
    L = w + 1
    m = L - 1  # state length, number of "previous days" we remember

    states = build_states(m, bound)
    idx = {s: i for i, s in enumerate(states)}
    S = len(states)

    # edges: from state j to state i by appending digit x
    edges = []
    for j, s in enumerate(states):
        sm = sum(s)
        for x in range(bound - sm + 1):
            ns = s[1:] + (x,)
            i = idx[ns]
            edges.append((i, j, x))

    # compatibility matrix for wrap-around windows across end and start
    # C[P,E]=1 if for every j=1..m:
    #   sum(last j of E) + sum(first L-j of P) <= bound
    compat = np.zeros((S, S), dtype=np.uint8)
    states_arr = [np.array(s, dtype=np.int16) for s in states]

    for p_idx, P in enumerate(states_arr):
        P_prefix = np.concatenate([[0], np.cumsum(P)])  # P_prefix[t]=sum(P[:t])
        for e_idx, E in enumerate(states_arr):
            E_suffix = np.concatenate([[0], np.cumsum(E[::-1])])  # suffix sums by length
            ok = True
            for j in range(1, m + 1):
                # last j from E + first (L-j) from P
                if E_suffix[j] + P_prefix[L - j] > bound:
                    ok = False
                    break
            compat[p_idx, e_idx] = 1 if ok else 0

    # start weight: z^{sum(P)} / prod(P_i!)
    inv_fact_prod = np.zeros(S, dtype=np.float64)
    sums = np.zeros(S, dtype=np.int16)
    for i, s in enumerate(states):
        sums[i] = sum(s)
        denom = 1
        for x in s:
            denom *= fac(x)
        inv_fact_prod[i] = 1.0 / denom

    inv_fact = [1.0 / fac(x) for x in range(bound + 1)]
    return states, edges, compat, sums, inv_fact_prod, inv_fact, m


def H_of_z(z: float, edges, compat, sums, inv_fact_prod, inv_fact, steps: int):
    """
    Compute H(z) = sum_{valid (a_t)} z^{sum a_t} / prod a_t!
    via transfer matrix and wrap compatibility.
    """
    S = len(sums)
    T = np.zeros((S, S), dtype=np.float64)

    z2 = z * z
    z3 = z2 * z
    zp = (1.0, z, z2, z3)

    # T[next, cur] = z^x / x!
    for i, j, x in edges:
        T[i, j] += zp[x] * inv_fact[x]

    # Tpow[end, start]
    Tpow = np.linalg.matrix_power(T, steps)

    # For each start p:
    # sum over end e: Tpow[e,p] * compat[p,e]
    col_sum = (Tpow * compat.T).sum(axis=0)

    # multiply by starting-day weight z^{sum(P)}/prod(P_i!)
    w = inv_fact_prod * (z ** sums)

    return float(np.dot(w, col_sum))


def expected_value(K: int, w: int, Y: int, quad_n: int = 40):
    """
    E[N] = ∫_0^∞ e^{-u} H(u/Y) du
    Use Gauss-Laguerre quadrature with quad_n points.
    """
    states, edges, compat, sums, inv_fact_prod, inv_fact, m = prepare(K, w)
    steps = Y - m

    nodes, weights = np.polynomial.laguerre.laggauss(quad_n)

    ans = 0.0
    for x, wgt in zip(nodes, weights):
        z = x / Y
        ans += wgt * H_of_z(z, edges, compat, sums, inv_fact_prod, inv_fact, steps)

    return ans


if __name__ == "__main__":
    # Earth case: 4 people within 7 days, 365-day year
    val = expected_value(K=4, w=7, Y=365, quad_n=40)
    print(f"{val:.8f}")