Project Euler 737

题目

Coin Loops

A game is played with many identical, round coins on a flat table.

Consider a line perpendicular to the table.

The first coin is placed on the table touching the line.

Then, one by one, the coins are placed horizontally on top of the previous coin and touching the line.

The complete stack of coins must be balanced after every placement.

The diagram below [not to scale] shows a possible placement of \(8\) coins where point \(P\) represents the line.

It is found that a minimum of \(31\) coins are needed to form a coin loop around the line, i.e. if in the projection of the coins on the table the centre of the \(n\text{th}\) coin is rotated \(\theta_n\), about the line, from the centre of the \((n-1)\text{th}\) coin then the sum of \(\displaystyle\sum_{k=2}^n \theta_k\) is first larger than \(360^\circ\) when \(n=31\). In general, to loop \(k\) times, \(n\) is the smallest number for which the sum is greater than \(360^\circ k\).

Also, \(154\) coins are needed to loop two times around the line, and \(6947\) coins to loop ten times.

Calculate the number of coins needed to loop \(2020\) times around the line.

解决方案

把桌面投影平面中的那条竖线接触点记为点 \(P\)，并作为原点。设所有硬币半径都为 \(1\)（缩放不影响角度与圈数）。第 \(k\) 枚硬币圆心记为 \(O_k\)。由于每枚硬币都与竖线相切，所以投影中有\(|O_k|=1.\)

设其极角为 \(\alpha_k\)，题目中的转角就是相邻圆心极角差，累计转角为\(\displaystyle F(n)=\sum_{k=1}^{n}(\alpha_{k+1}-\alpha_k)=\alpha_{n+1}-\alpha_1.\)若要绕竖线转 \(K\) 圈，则需找到最小的 \(n\) 使得\(F(n)>2\pi K,\)最终答案是 \(n+1\)。

每次放第 \(k+1\) 枚硬币时，前 \(k\) 枚硬币整体质心必须落在第 \(k+1\) 枚硬币圆盘内，才能保持平衡。

为了让总硬币数最少，需要让每一步新增的转角尽可能大。若质心严格在圆盘内部，那么仍可沿着可行方向继续旋转新硬币一点点，且仍然平衡，从而得到更大转角。这与最优矛盾。

因此最优摆放必须满足 临界平衡：\(|O_{k+1}-W_k|=1,\)其中 \(W_k\) 是前 \(k\) 枚硬币的质心。

把平面视为复平面，令\(\displaystyle O_k=e^{i\alpha_k},S_k=\sum_{j=1}^k O_j,W_k=\frac{S_k}{k}.\)初始取 \(\alpha_1=0\)，即 \(O_1=1\)。由 \(|O_{k+1}|=1\) 且 \(|O_{k+1}-W_k|=1\)，展开平方：\(|O_{k+1}|^2-2\langle O_{k+1},W_k\rangle+|W_k|^2=1,\)得到\(\langle O_{k+1},W_k\rangle=\dfrac{|W_k|^2}{2}.\)记\(s_{k+1}=\langle O_{k+1},S_k\rangle,\)则\(s_{k+1}=k\langle O_{k+1},W_k\rangle=\dfrac{|S_k|^2}{2k}.\)

令 \(\displaystyle H_k=\sum_{j=1}^k\frac1j\) 为调和数。先看递推：\(S_{k+1}=S_k+O_{k+1},\)因此\(|S_{k+1}|^2=|S_k|^2+|O_{k+1}|^2+2\langle S_k,O_{k+1}\rangle=|S_k|^2+1+2s_{k+1}.\)代入上式 \(s_{k+1}=\dfrac{|S_k|^2}{2k}\)，得

\[ |S_{k+1}|^2=\left(1+\frac1k\right)|S_k|^2+1. \]

又有 \(|S_1|^2=1=1\cdot H_1\)。并且由\(S_{k+1}=S_k+O_{k+1}\)以及\(|S_{k+1}|^2=|S_k|^2+1+2\langle S_k,O_{k+1}\rangle=|S_k|^2+1+\dfrac{|S_k|^2}{k}=\left(1+\dfrac1k\right)|S_k|^2+1\)可归纳得到\(|S_k|^2=kH_k.\)于是质心模长为\(|W_k|=\dfrac{|S_k|}{k}=\sqrt{\dfrac{H_k}{k}}.\)

设第 \(k\) 步增量为\(\Delta_k=\alpha_{k+1}-\alpha_k.\)我们把它拆成两个角：\(A_k=\angle(W_k,O_{k+1}),B_k=\angle(W_k,O_k)\)则取最大转角的那个分支时，有\(\Delta_k=A_k-B_k.\)

首先计算\(A_k.\)由前文内积关系可知，\(\cos A_k=\dfrac{\langle O_{k+1},W_k\rangle}{|O_{k+1}||W_k|}=\dfrac{|W_k|}{2}=\dfrac12\sqrt{\dfrac{H_k}{k}}.\)所以\(A_k=\arccos\left(\dfrac12\sqrt{\dfrac{H_k}{k}}\right).\)

接下来计算\(B_k\).先算 \(\langle O_k,S_k\rangle\)。因为 \(S_k=S_{k-1}+O_k\)，\(\langle O_k,S_k\rangle=1+\langle O_k,S_{k-1}\rangle=1+s_k.\)而对 \(s_k\) 使用前面的公式（把 \(k+1\) 改成 \(k\)）：\(s_k=\dfrac{|S_{k-1}|^2}{2(k-1)}=\dfrac{H_{k-1}}2.\)故\(\langle O_k,S_k\rangle=1+\dfrac{H_{k-1}}2=\dfrac{H_k-\frac1k+2}{2}.\)

于是\(\cos B_k=\dfrac{\langle O_k,W_k\rangle}{|O_k||W_k|}=\dfrac{\langle O_k,S_k\rangle/k}{\sqrt{H_k/k}}=\dfrac{H_k-\frac1k+2}{2\sqrt{kH_k}}.\)因此有\(B_k=\arccos\left(\dfrac{H_k-\frac1k+2}{2\sqrt{kH_k}}\right).\)

综上所述，可得，

\[ \Delta_k= \arccos\left(\frac12\sqrt{\frac{H_k}{k}}\right)- \arccos\left(\frac{H_k-\frac1k+2}{2\sqrt{kH_k}}\right). \]

从而\(\displaystyle F(n)=\sum_{k=1}^n \Delta_k.\)

直接逐项累加到答案规模在效率上仍不可行，因此需要对总转角函数 \(F(n)\) 做渐近分析。记\(h=\log n+\gamma,\)其中 \(\gamma\) 为欧拉–马歇罗尼常数。具体做法是：先将单步增量 \(\Delta_n\) 中涉及的量按 \(n\to\infty\) 展开

其中用到调和数展开\(\displaystyle H_n=\log n+\gamma+\frac1{2n}-\frac1{12n^2}+\frac1{120n^4}-\frac1{252n^6}+\cdots,\)以及 \(\arccos(x)\) 在 \(x=0\) 附近的级数展开），再对\(\displaystyle F(n)=\sum_{k=1}^n \Delta_k\)应用 Maclaurin 求和公式，将离散求和转化为积分项与若干端点/高阶修正项，从而得到 \(F(n)\) 的渐近展开。

可以得到 \(F(n)\) 的渐近形式。其主导项为

\[ F(n)\sim \sqrt{2\pi}e^{-\gamma/2}\operatorname{erfi}\left(\sqrt{\frac{h}{2}}\right), \]

这来自于单步增量的主阶行为\(\displaystyle \Delta_n\sim \frac1{\sqrt{n(\log n+\gamma)}}.\)对该主项积分时，令\(u=\sqrt{\dfrac{\log x+\gamma}{2}},\)就会自然出现 \(\operatorname{erfi}(u)\)。进一步保留到较高精度，写成\(F(n)=\Phi(n)+C+\widetilde O(n^{-3/2}),\)

其中（这里把端点修正也并入 \(\Phi\)）

\[ \Phi(n)= \sqrt{2\pi}e^{-\gamma/2}\operatorname{erfi}!\left(\sqrt{\frac h2}\right) +\frac{\sqrt{2\pi}}{24}e^{\gamma/2}\operatorname{erfc}!\left(\sqrt{\frac h2}\right) +\frac{2-3h}{12\sqrt{hn}} +\frac12,\delta_{\mathrm{end}}(n), \]

而

\[ \delta_{\mathrm{end}}(n)= \arccos\left(\frac{H_n^{(*)}}{2\sqrt{nH_n^{(*)}}}\right)= \arccos\left(\frac{H_n^{(*)}-\frac1n+2}{2\sqrt{nH_n^{(*)}}}\right), \]

这里 \(H_n^{(*)}\) 表示用调和数的渐近展开代替 \(H_n\)（在实现中取到 \(n^{-6}\) 项已足够）。

常数 \(C\) 与 \(n\) 无关。为了避免显式求出该常数，可引入一个参考点（锚点）\(M\)，先用精确公式直接累加得到\(\displaystyle F(M)=\sum_{k=1}^M \Delta_k.\)

由\(F(n)=\Phi(n)+C+\varepsilon(n),F(M)=\Phi(M)+C+\varepsilon(M),\)两式相减可得\(F(n)=F(M)+(\Phi(n)-\Phi(M)\bigr)+\bigl(\varepsilon(n)-\varepsilon(M)).\)因此定义近似函数\(\widehat F(n)=F(M)+\Phi(n)-\Phi(M),\)

则未知常数项 \(C\) 被自动消去；同时在 \(n\ge M\) 时，这种写法通常还能改善数值稳定性。随后只需二分求出最小满足\(\widehat F(n)>2\pi K\)的 \(n\)，答案即为 \(n+1\)。

代码

import mpmath as mp

N = 2020
def h_series(n):
    n = mp.mpf(n)
    return mp.log(n) + mp.euler + mp.mpf(1) / (2 * n) - mp.mpf(1) / (12 * n**2) + mp.mpf(1) / (120 * n**4) - mp.mpf(1) / (252 * n**6)

def delta_series_endpoint(n):
    n = mp.mpf(n)
    h = h_series(n)
    den = 2 * mp.sqrt(h * n)
    return mp.acos(h / den) - mp.acos((h - mp.mpf(1) / n + 2) / den)

def A(n):
    n = mp.mpf(n)
    h = mp.log(n) + mp.euler
    x = mp.sqrt(h / 2)
    c0 = mp.sqrt(2 * mp.pi) * mp.e**(-mp.euler / 2)
    c1 = mp.sqrt(2 * mp.pi) * mp.e**( mp.euler / 2) / 24
    return c0 * mp.erfi(x) + c1 * mp.erfc(x) + (2 - 3 * h) / (12 * mp.sqrt(h * n)) + delta_series_endpoint(n) / 2

def exact_F(m):
    h = mp.mpf(0)
    s = mp.mpf(0)
    k = mp.mpf(0)
    for _ in range(m):
        k += 1
        invk = 1 / k
        h += invk
        den = 2 * mp.sqrt(h * k)
        s += mp.acos(h / den) - mp.acos((h - invk + 2) / den)
    return s

def solve(K):
    target = 2 * mp.pi * K
    M = 10 ** 5
    FM = exact_F(M)
    if FM > target:
        h = mp.mpf(0)
        s = mp.mpf(0)
        k = mp.mpf(0)
        for i in range(1, M + 1):
            k += 1
            invk = 1 / k
            h += invk
            den = 2 * mp.sqrt(h * k)
            s += mp.acos(h / den) - mp.acos((h - invk + 2) / den)
            if s > target:
                return i + 1
    shift = FM - A(M)

    def Fest(n):
        return A(n) + shift

    l = M
    r = int(mp.ceil(3 * target * target * mp.log(target + 10))) + 10
    if r <= l:
        r = l + 1
    while Fest(r) <= target:
        r *= 2
    while l < r:
        mid = (l + r) // 2
        if Fest(mid) > target:
            r = mid
        else:
            l = mid + 1
    return l + 1

print(solve(N))