Project Euler 314

发表于 2026-01-07 分类于 Project Euler 阅读次数：本文字数： 4.1k 阅读时长 ≈ 4 分钟

Project Euler 314

题目

The Mouse on the Moon

The moon has been opened up, and land can be obtained for free, but there is a catch. You have to build a wall around the land that you stake out, and building a wall on the moon is expensive. Every country has been allotted a $500 \text{m}$ by $500 \text{m}$ square area, but they will possess only that area which they wall in. $251001$ posts have been placed in a rectangular grid with $1$ meter spacing. The wall must be a closed series of straight lines, each line running from post to post.

The bigger countries of course have built a $2000 \text{m}$ wall enclosing the entire $250 000 \text{m}^2$ area. The Duchy of Grand Fenwick, has a tighter budget, and has asked you (their Royal Programmer) to compute what shape would get best maximum enclosed-area/wall-length ratio.

You have done some preliminary calculations on a sheet of paper.

For a $2000$ meter wall enclosing the $250 000 \text{m}^2$ area the enclosed-area/wall-length ratio is $125$.

Although not allowed , but to get an idea if this is anything better: if you place a circle inside the square area touching the four sides the area will be equal to $\pi\cdot 250^2 \text{m}^2$ and the perimeter will be $\pi\cdot 500 \text{m}$, so the enclosed-area/wall-length ratio will also be $125$.

However, if you cut off from the square four triangles with sides $75 \text{m}$, $75 \text{m}$ and $75\sqrt{2} \text{m}$ the total area becomes $238750 \text{m}^2$ and the perimeter becomes $1400+300\sqrt{2} \text{m}$. So this gives an enclosed-area/wall-length ratio of $130.87$, which is significantly better.

Find the maximum enclosed-area/wall-length ratio.

Give your answer rounded to $8$ places behind the decimal point in the form $abc.defghijk$.

解决方案

Dinkelbach 迭代

Dinkelbach 迭代用于求解经典的分式规划问题（最大值或最小值均可，不失一般性，这里使用最大值讲解），例如

\[ \max_{x\in\mathcal X}\frac{f(x)}{g(x)}, g(x)>0. \]

对任意实数 $r$，定义

\[ \phi(r)=\max_{x\in\mathcal X}(f(x)-r,g(x)). \]

这个函数有关键性质：

设最优比值为 $r^*$，则$\phi(r^*)=0.$
若 $r<r^{\ast}$，则存在某个 $x$ 使 $\dfrac{f(x)}{g(x)}>r$，于是 $f(x)-r g(x)>0$，从而$\phi(r)>0.$
若 $r>r^{\ast}$，则对所有 $x$ 有 $\dfrac{f(x)}{g(x)}\le r$，于是 $f(x)-r g(x)\le 0$，从而$\phi(r)\le 0.$

因此，最大化 $\dfrac{f(x)}{g(x)}$ 等价于寻找使 $\phi(r)=0$ 的 $r$。

迭代步骤：

从一个初值 $r_0$ 开始，重复：

解子问题$x_k\in\arg\max_{x\in\mathcal X}\bigl(f(x)-r_k g(x)\bigr).$
用该解更新比值$r_{k+1}=\dfrac{f(x_k)}{g(x_k)}.$

常用停止条件是子问题的最优值接近 0：$f(x_k)-r_k g(x_k)\approx 0$；或相邻两次比值变化很小：$|r_{k+1}-r_k|<\varepsilon$。

回到本题，令$R=250$，可见，所构造的图形一定是一个凸多边形，这是因为若边界存在凹角，把该凹角的两条边用其端点的弦替换：

面积不减（凹进去的部分被“填平”）；
周长不增（由三角不等式，折线长度 $\ge$ 端点直线段长度）。

反复消除凹角可得到凸格点多边形，且 $A/L$ 不下降。因此最优解可限制在凸的格点多边形上。

此处我们做出一个假设：最优形状具有 $D_4$ 对称性（关于两坐标轴与两条对角线对称）。

由于路径必定包含格点$(0,R),(R,0)$，因此菱形$(0,R)\rightarrow(R,0)\rightarrow(0,-R)\rightarrow(-R,0)\rightarrow(0,R)$里面的所有点都肯定不是轨迹的一部分。

因此在第一象限中，我们只需要枚举满足 $y\ge x$ 且 $x+y\ge R$（即 $y\ge \max(x,R-x)$）的格点。原因是：在假设边界具有八重对称时，第一象限内的边界从 $(0,R)$ 连接到 $(R,0)$，并且关于对角线 $y=x$ 镜像对称。

于是我们固定起点为 $S=(0,R)$，在上述允许点集中寻找一条从 $S$ 到某个点 $P=(x,y)$ 的单调折线链；得到这段链后，再将其关于 $y=x$ 反射并拼接，从而补全第一象限的整段边界。

设原点为 $O=(0,0)$。令链为$V_0=S,V_1,\dots,V_k=P,$那么从$(0,R)$到$(R,0)$的路径为$V_0,V_1,\dots,V_k,V_k',\dots,V_1',V_0'$

其中每个 $V_i=(x_i,y_i)$ 都是格点，$V_i'=(y_i,x_i)$是$V_i$关于$x=y$的对称点。

对于折线链 $V_0,V_1,\dots,V_k$，考虑相邻两点 $V=(x_2,y_2)$ 与 $U=(x_1,y_1)$。

面积增量：原点 $O$ 与线段 $VU$ 组成的有向三角形面积为$\Delta A=\dfrac{x_1y_2-x_2y_1}{2}.$
长度增量：该段边的长度为$\Delta L=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}.$

因此，链的总面积$A_1$与总长度$L_1$都可以通过对这些增量分别求和得到。

到达 $P=(x,y)$ 后，连接到对称点 $P'=(y,x)$：

这段长度为$ L_{2}==(y-x).$
对应三角形的面积（$y\ge x$）为$A_{2}=\dfrac{y^2-x^2}{2}.$

由于关于对角线对称，第一象限的总面积与总周长为$A_0=2A_1+A_2,L=2L_1+L_2$

整体四象限旋转复制后，比例不变，因此最终目标等价于最大化$\dfrac{A}{L}=\dfrac{A_0}{L_0}.$

直接最大化 $\frac{A}{L}$ 不易做动态规划。这时考虑使用Dinkelbach 迭代对任意猜测比值 $r$，定义$F_r = A_{0} - rL_{0}.$

设$g(r)=\max F_r.$

则有：

若存在形状使 $\dfrac{A}{L}>r$，则该形状满足 $F_r>0$，从而 $g(r)>0$；
若所有形状都满足 $\frac{A}{L}\le r$，则对所有形状 $F_r\le0$，从而 $g(r)\le0$；
最优比值 $r^*$ 满足 $g(r^*)=0$。

可见，$F_r=2(A_1-rL_1)+(A_2-rL_2).$

因此 DP 先求每个点 $(x,y)$ 的最优$S_{1}(x,y)=A_{1}-rL_{1},$并同时记录对应的 $A_{1}$ 与 $L_{1}$ 以便最后恢复 $A_{2},L_{2}$。

因此，为了保持凸性，$x$ 单调增，$y$ 单调不增。

以下还有两个转移剪枝，对从前点 $(x_2,y_2)$ 到当前点 $(x_1,y_1)$，要求：

方向不陡过 $-45^\circ$，也就是$y_2 \le x_1+y_1-x_2$，这是因为把第一象限的整条路径恢复后，必定违反凸性。
前点位于起点 $S=(0,R)$ 到当前点 $(x_1,y_1)$ 的连线之上：直线 $S\rightarrow(x_1,y_1)$ 在 $x=x_2$ 处的值为$y_{0} = R + (y_1-R)\dfrac{x_2}{x_1},$那么要求$y_2 > y_{0}.$这是避免在枚举路径中违反凸性。

代码

import math

R = 250
SQRT2 = math.sqrt(2.0)

# 预计算所有整数位移的距离：dist[dx][dy] = sqrt(dx^2 + dy^2)
dist = [[0.0] * (R + 1) for _ in range(R + 1)]
for dx in range(R + 1):
    for dy in range(R + 1):
        dist[dx][dy] = math.hypot(dx, dy)

# 1/16 区域下界：y >= max(x, R-x)
y_low = [max(x, R - x) for x in range(R + 1)]


def best_for_ratio(ratio: float):
    """
    固定 ratio=r，最大化 F_r = A_quad - r * L_quad
    返回 (best_S, best_A, best_L) 对应第一象限(1/4)的 (F_r, A_quad, L_quad)
    """
    inf = float('inf')
    dpS = [[-inf] * (R + 1) for _ in range(R + 1)]
    dpA = [[0.0] * (R + 1) for _ in range(R + 1)]
    dpL = [[0.0] * (R + 1) for _ in range(R + 1)]

    dpS[0][R] = 0.0  # 起点 (0,R)

    for x1 in range(1, R + 1):
        y1_min = y_low[x1]
        for y1 in range(R, y1_min - 1, -1):
            # 基础：从 (0,R) 直接连到 (x1,y1)
            bestA = x1 * R * 0.5
            bestL = dist[x1][abs(y1 - R)]
            bestS = bestA - ratio * bestL


            # 从所有可能前驱 (x2,y2) 转移
            for x2 in range(1, x1):
                # 约束1：斜率 >= -1
                y2_hi = x1 + y1 - x2 - 1
                if y2_hi > R:
                    y2_hi = R

                # 约束2：y2 > 线性插值 y_line = R + (y1-R)*x2/x1
                y2_lo = R + (y1 - R) * x2 // x1  # floor

                if y2_hi <= y2_lo:
                    continue

                # 还要满足 y2 在 1/16 区域内
                y2_reg = y_low[x2]
                if y2_hi < y2_reg:
                    continue
                if y2_lo < y2_reg - 1:
                    y2_lo = y2_reg - 1

                dx = x1 - x2
                for y2 in range(y2_hi, y2_lo, -1):
                    prevS = dpS[x2][y2]
                    dA = (x1 * y2 - x2 * y1) * 0.5
                    dL = dist[dx][y2 - y1]  # y2 >= y1
                    newS = prevS + dA - ratio * dL
                    if newS > bestS:
                        bestS = newS
                        bestA = dpA[x2][y2] + dA
                        bestL = dpL[x2][y2] + dL

            dpS[x1][y1] = bestS
            dpA[x1][y1] = bestA
            dpL[x1][y1] = bestL

    # 用对角线反射补全到第一象限(1/4)
    bestTS = -inf
    bestTA = 0.0
    bestTL = 1.0

    for x in range(1, R + 1):
        y_min = y_low[x]
        for y in range(R, y_min - 1, -1):
            partS = dpS[x][y]
            if partS <= -inf / 2:
                continue
            connA = (y * y - x * x) * 0.5
            connL = (y - x) * SQRT2
            totalS = 2.0 * partS + connA - ratio * connL
            if totalS > bestTS:
                bestTS = totalS
                bestTA = 2.0 * dpA[x][y] + connA
                bestTL = 2.0 * dpL[x][y] + connL

    return bestTS, bestTA, bestTL


def solve():
    # Dinkelbach 迭代
    ratio = 125
    for _ in range(20):
        _, A, L = best_for_ratio(ratio)
        new_ratio = A / L
        if abs(new_ratio - ratio) < 1e-10:
            return new_ratio
        ratio = new_ratio
    return ratio


if __name__ == "__main__":
    ans = solve()
    print(f"{ans:.8f}")