Project Euler 262

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

Project Euler 262

题目

Mountain Range

The following equation represents the continuous topography of a mountainous region, giving the elevation h at any point \((x,y)\):

\[h = \left(5000 - \frac{x^2 + y^2 + xy}{200} + \frac{25(x + y)}2\right) \cdot e^{-\left|\frac{x^2 + y^2}{1000000} - \frac{3(x + y)}{2000} + \frac 7 {10}\right|}\]

A mosquito intends to fly from \(A(200,200)\) to \(B(1400,1400)\), without leaving the area given by \(0\le x,y\le1600\).

Because of the intervening mountains, it first rises straight up to a point \(A'\), having elevation \(f\).

Then, while remaining at the same elevation \(f\), it flies around any obstacles until it arrives at a point \(B'\) directly above \(B\).

First, determine \(f_{\min}\) which is the minimum constant elevation allowing such a trip from \(A\) to \(B\), while remaining in the specified area.

Then, find the length of the shortest path between \(A'\) and \(B'\), while flying at that constant elevation \(f_{\min}\).

Give that length as your answer, rounded to three decimal places.

Note: For convenience, the elevation function shown above is repeated below, in a form suitable for most programming languages:

1
h=(5000-0.005*(x*x+y*y+x*y)+12.5*(x+y))*exp(-abs(0.000001*(x*x+y*y)-0.0015*(x+y)+0.7))

解决方案

令\(p(x,y)=5000 - \dfrac{x^2 + y^2 + xy}{200} + \dfrac{25(x + y)}{2},g(x,y)=\dfrac{x^2 + y^2}{1000000} - \dfrac{3(x+y)}{2000} + \dfrac{7}{10}\)，那么有\(h(x, y) = p(x,y)\cdot e^{-|g(x,y)|}.\)

此外还需要注意的一点是，\(h(x,y)=h(y,x).\)

\(h(x,y)\)的图像和等高线图如下图所示。所使用的 Mathematica 代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
p[x_, y_] := 5000 - (x^2 + y^2 + x*y)/200 + (25*(x + y))/2;
g[x_, y_] := (x^2 + y^2)/1000000 - (3*(x + y))/2000 + 7/10;
h[x_, y_] := p[x, y]*Exp[-Abs[g[x, y]]];

(*三维曲面图*)
Plot3D[h[x, y], {x, 0, 1600}, {y, 0, 1600}, PlotRange -> All,
 AxesLabel -> {"x", "y", "h"}, PlotPoints -> 50, Mesh -> None,
 ColorFunction -> "Rainbow", PerformanceGoal -> "Quality"]

(*等高线图*)
ContourPlot[h[x, y], {x, 0, 1600}, {y, 0, 1600},
 PlotLegends -> Automatic, ColorFunction -> "Rainbow",
 PlotRange -> All,
 Epilog -> {{Red, PointSize[0.03], Point[{200, 200}]}, {Blue,
    PointSize[0.03], Point[{1400, 1400}]},
   Text[Style["A(200,200)", 12, Bold, Background -> White], {200,
     200}, Scaled[{0, -0.01}]],
   Text[Style["B(1400,1400)", 12, Bold, Background -> White], {1400,
     1400}, Scaled[{0, 0.01}]]}]

若在恒定高度 \(f\) 飞行，则蚊子在平面上允许经过的点必须满足：\(h(x,y)\le f\).

因此对任意 \(f\)，可把\(\Omega_f=\{(x,y)\in[0,1600]^2:\ h(x,y)>f\}\)视为“山体障碍”。

可飞区域为\(D_f=[0,1600]^2\setminus\Omega_f.\)

随着 \(f\) 增大，障碍 \(\Omega_f={h>f}\) 单调缩小，连通性只会“从不连通变为连通”，并在某个临界值首次连通。

在临界时刻\(f=f_{\min}\)，“从不可连通到可连通”通常发生在一个瓶颈被刚好打开的瞬间。对于连续地形的等高线结构，这个瞬间对应两种典型事件之一：

  1. 障碍集合 \(\Omega_f\) 自身发生拓扑变化（比如出现/消失一个洞、两个分支接触等）；
  2. 障碍集合 \(\Omega_f\) 与边界 \([0,1600]^2\) 的边界发生“相切接触”，从而堵住或打开通道。

本题的地形形态（经典解法所依赖的事实）属于第 2 类：临界时障碍的边界 \({h=f_{\min}}\) 恰好与正方形边界相切，形成“刚好卡住/刚好放开”的狭口，这在我们图中可以观察得出。

因此，我们假设，这个相切点落在边界 \(y=0\) 上，记为\(S=(x_0,0).\)

若 \(\Gamma\) 在边界处“刚好接触”而不穿出边界，则沿边界方向高度必须达到局部极值。对边界 (y=0) 而言，相切点满足\(\displaystyle{x_0=\arg\max_{x\in[0,1600]} h(x,0)}\)，那么就有\(f_{\min}=h(x_0,0)\).

固定高度 \(f_{\min}\) 后，障碍边界就是等高线\(\Gamma={(x,y): h(x,y)=f_{\min}}\).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
p[x_, y_] := 5000 - (x^2 + y^2 + x*y)/200 + (25*(x + y))/2;
g[x_, y_] := (x^2 + y^2)/1000000 - (3*(x + y))/2000 + 7/10;
h[x_, y_] := p[x, y]*Exp[-Abs[g[x, y]]];

ContourPlot[h[x, y], {x, 0, 1600}, {y, 0, 1600},
 Contours -> {10396.462193284104},(*特定等高线*)
 ContourStyle -> Directive[Thick, Black],(*等高线样式*)
 ContourLabels ->
  Function[{x, y, z},
   Text[Style[NumberForm[z, {8, 3}], 12, Bold,
     Background -> White], {x, y}]], PlotLegends -> Automatic,
 ColorFunction -> "Rainbow", PlotRange -> All,
 Epilog -> {{Red, PointSize[0.03], Point[{200, 200}]}, {Blue,
    PointSize[0.03], Point[{1400, 1400}]},
   Text[Style["A(200,200)", 12, Bold, Background -> White], {200,
     200}, {0, -8}],
   Text[Style["B(1400,1400)", 12, Bold, Background -> White], {1400,
     1400}, {0, 8}]}]
alt text

在欧氏平面中，绕开光滑障碍的最短路满足经典几何性质：

  1. 在可行区域内部（未贴着障碍边界）的部分，最短路一定是直线段；
  2. 最短路与障碍边界的连接点必须满足相切条件：若不是相切，可在局部“削去拐角”得到更短路径，矛盾。

因此，最短路径必可表示为：\(A \rightarrow T_A \rightarrow T_B \rightarrow B\)：其中\(A \rightarrow T_A, T_B \rightarrow B\)均为直线段，\(T_A \rightarrow T_B\)，且中间一段\(T_A,T_B\in\Gamma\) 沿着等高线 \(\Gamma\) 绕行，并在连接点 \(T_A\)、\(T_B\) 处与 \(Γ\) 满足相切条件（即直线段在该点与 \(\Gamma\) 的切线方向一致）。

因此，\(T_A,T_B\)满足连线方向应当与法向量正交，即\(\nabla h(T_A)\cdot \overrightarrow{T_AA}=0,\nabla h(T_B)\cdot \overrightarrow{T_BB}=0.\)

为计算在\(\Gamma\)上的固定目标切点，不妨假设 \(G\in\{A,B\}.\)

对等高线上的点 \(X\) 定义\(\Phi(X)=\nabla h(X)\cdot \overrightarrow{GX}.\)当\(\Phi(X)=0\)时，连线 \(GX\) 与等高线在 \(X\) 处满足切线条件（也就是连线方向落在切空间内，等价于与法向量正交）。

因此，从瓶颈点 \(S\) 沿 \(\Gamma\) 走动时，只要监控 \(\Phi\) 的符号变化，就能找到“从尚未切到已经越过”的过零点；再对该过零区间做二分，就能把切点逼近到任意精度。

最终，用很小步长沿某个方向推进一个坐标，再在新竖线或新横线上解方程\(h(\cdot,\cdot)=f_{\min}\)把点“拉回”到等高线。因为步长小，求根区间很短，用二分法非常稳健。相邻点的欧氏距离累加就是弧长近似。

分别得到：

  • 从 \(S\) 到 \(T_A\) 的等高线弧长记为 \(L_A^{\text{arc}}\)，直线段长度记为 \(|AT_A|\)；
  • 从 \(S\) 到 \(T_B\) 的等高线弧长记为 \(L_B^{\text{arc}}\)，直线段长度记为 \(|BT_B|\)。

总最短路长度为\(L=(L_A^{\text{arc}}+|AT_A|)+(L_B^{\text{arc}}+|BT_B|).\)

代码

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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
import math  # 导入数学库：用到 exp / abs / hypot 等基础函数


def h(x: float, y: float) -> float:  # 高度函数：输入 (x,y)，输出海拔 h
    """
    地形高度函数 h(x, y)。

    题面给定：
        h = P(x,y) * exp(-|G(x,y)|)

    其中：
        P(x,y) = 5000 - 0.005*(x^2 + y^2 + x*y) + 12.5*(x+y)
        G(x,y) = 1e-6*(x^2 + y^2) - 0.0015*(x+y) + 0.7

    注意：exp(-|G|) 使得高度在某些区域会被指数衰减。
    """
    p = 5000.0 - 0.005 * (x * x + y * y + x * y) + 12.5 * (x + y)  # 多项式部分 P(x,y)
    g = 0.000001 * (x * x + y * y) - 0.0015 * (x + y) + 0.7  # 指数衰减的“指数内部”G(x,y)
    return p * math.exp(-abs(g))  # 组合得到 h = P * exp(-|G|)


def dh(x: float, y: float) -> float:  # 辅助“导数表达式”：主要用于方向判定与切线判据
    """
    计算一个“沿 x 方向的辅助导数表达式”。

    这里返回的是：
        dh = Px - P * Gx
    其中：
        Px = ∂P/∂x
        Gx = ∂G/∂x

    它的用途是：
    - 在 y=0 的截线上，用 dh(x,0)=0 来定位该截线上的极值点位置（峰顶/鞍点方向上的临界点）。
    - 在沿等高线行走时，用 dh(x,y) 与 dh(y,x) 组合成一个“近似梯度方向”的向量，
      用于判断何时达到“切线”条件、以及决定下一步沿哪条坐标方向更新。

    注意：这里并没有完整展开 h 的真实偏导（真实偏导还会带 exp(-|G|) 和 |G| 的符号项）。
    这种写法相当于提取出对“方向/符号判断”最有用的部分，从而使算法更简洁、数值上也较稳。
    """
    p = 5000.0 - 0.005 * (x * x + y * y + x * y) + 12.5 * (x + y)  # 计算 P(x,y)，供后续复用
    px = -0.005 * (2.0 * x + y) + 12.5  # 计算 Px = ∂P/∂x（对 x 求偏导）
    gx = 0.000002 * x - 0.0015  # 计算 Gx = ∂G/∂x（对 x 求偏导）
    return px - gx * p  # 返回 dh = Px - P*Gx（作为“法向/梯度方向”的简化替身）


def dist(ax: float, ay: float, bx: float, by: float) -> float:  # 两点距离函数：用于累计弧长、以及末段直线距离
    """
    欧氏距离。
    用于累计等高线走过的弧长、以及最终切点到目标点的直线距离。
    """
    return math.hypot(ax - bx, ay - by)  # hypot(u,v)=sqrt(u^2+v^2)，数值更稳健


def sgn(v: float) -> float:  # 符号函数：只返回 ±1（这里不区分 0）
    """
    符号函数：返回 +1 或 -1（不区分 0）。
    在构造“用于夹住根”的区间端点时，用它来决定向哪一侧试探。
    """
    return 1.0 if v > 0.0 else -1.0  # v>0 返回 +1，否则返回 -1（v==0 也走 -1）


def findz(eps: float):  # 在 y=0 上二分定位 dh(x,0)=0 的点，并取其高度作为候选 f_min
    """
    在 y=0 的截线上，用二分法寻找 dh(x,0)=0 的位置 x0，
    并返回该点与其高度 (x0, h(x0,0))。

    直观意义：
    - y=0 是边界的一条线。
    - 我们在这条线上找一个“高度达到峰值附近的关键点”，作为等高线追踪的起点 S=(x0,0)。
    - 然后把 z0 = h(x0,0) 当作候选的“最低可通行高度 f_min”。

    二分规则：
    - 对区间中点 c 计算 dh(c,0)。
    - 若 dh(c,0) > 0，表示在该侧高度仍在上升趋势，于是移动左端点到 c；
      否则移动右端点到 c。
    """
    l, r = 0.0, 1600.0  # 搜索区间：[0,1600]，对应题面允许的边界范围

    while abs(l - r) > eps:  # 直到区间长度收缩到 eps 以内（满足精度要求）
        mid = (l + r) * 0.5  # 取中点
        if dh(mid, 0.0) > 0.0:  # 若 dh>0，按约定认为“往右仍在上升”，根在右侧
            l = mid  # 丢弃左半段：l 向右收缩到 mid
        else:  # 否则认为“已不再上升”，根在左侧
            r = mid  # 丢弃右半段：r 向左收缩到 mid

    x0 = (l + r) * 0.5  # 取最终区间中点作为 x0 的近似值
    return x0, h(x0, 0.0)  # 返回 (x0, z0)，其中 z0=h(x0,0)


def bisect_root_for_level(z0: float, fixed_y: float, a: float, b: float, eps: float) -> float:  # 在一条竖线/横线上用二分把点“拉回”到等高线
    """
    在固定 fixed_y 的情况下，求解方程：
        h(x, fixed_y) = z0
    在给定区间 [a,b] 上做二分（要求端点函数值异号才能保证夹住根）。

    返回值：
    - 若 [a,b] 的确夹住根：返回该根的近似位置。
    - 若没夹住根：返回端点中更接近等高线的那个点（更“贴近” z0）。

    这个函数的目的：
    - 在等高线行走时，我们会先把某个坐标推进一步得到新的 x（或 y），
      然后需要在该新 x（或新 y）上反求对应的 y（或 x），
      让点重新落回到同一条等高线上 h=z0。
    - 由于步长不大，通常 [a,b] 是一个很短的邻域区间，因此二分非常稳定。
    """
    fa = h(a, fixed_y) - z0  # 计算端点 a 处的“等高线残差”：h(a,fixed_y)-z0
    fb = h(b, fixed_y) - z0  # 计算端点 b 处的“等高线残差”：h(b,fixed_y)-z0

    if fa == 0.0:  # 若 a 点恰好在等高线上
        return a  # 直接返回 a
    if fb == 0.0:  # 若 b 点恰好在等高线上
        return b  # 直接返回 b

    # 若 fa 与 fb 同号，说明 [a,b] 未必穿过等高线（可能没夹住根）
    if fa * fb > 0.0:  # 同号：无法保证区间内存在根
        return a if abs(fa) < abs(fb) else b  # 退而求其次：返回更接近等高线的端点

    l, r = a, b  # 初始化二分区间端点
    fl, fr = fa, fb  # 保存端点函数值，避免重复计算

    while abs(r - l) > eps:  # 二分直到区间缩到 eps
        m = (l + r) * 0.5  # 中点
        fm = h(m, fixed_y) - z0  # 中点残差
        if fl * fm <= 0.0:  # 若 [l,m] 异号（或中点为零），则根在左半区间
            r, fr = m, fm  # 收缩右端点到 m
        else:  # 否则根在 [m,r]
            l, fl = m, fm  # 收缩左端点到 m

    return (l + r) * 0.5  # 返回最终区间中点作为根的近似


def get_length(step: float, eps: float, x0: float, z0: float, gx: float, gy: float) -> float:  # 从 S 沿等高线追踪到切点，再直线连到目标点 (gx,gy)
    """
    从起点 S=(x0,0) 出发，沿着等高线 h=z0 进行“折线式追踪”，
    直到到达一个切点 T，使得从 T 到目标点 (gx,gy) 的直线
    与等高线在 T 处“刚好相切/不再穿入障碍”。

    返回值：
        等高线弧长( S -> ... -> T ) + 直线距离( T -> (gx,gy) )

    关键思想：
    - 在高度限制为 z0 时，可飞行区域是 { (x,y) | h(x,y) <= z0 }。
      障碍物边界就是等高线 h=z0。
    - 最短路径会由两段组成：从目标到切点的直线 + 沿边界滑行的弧线（反过来一样）。
      切点处满足“直线方向与边界切线相切”的几何条件（等价于：目标方向与法向的点积为 0）。
    - 本函数用一个简单的“切线判据”来检测何时到了切点。

    追踪方法：
    - 在当前点 (x,y) 处，构造一个近似法向量 (dx,dy)：
        dx = dh(x,y)
        dy = dh(y,x)
    - 若 (dx,dy)·((x,y)-(gx,gy)) >= 0，则认为已经到达切点：
      再沿等高线走只会“更不利于直线连到目标”，于是结束追踪并改走直线。
    - 否则继续沿等高线走一步：
      先推进一个坐标（x 或 y），再在对应的竖线/横线上解 h(·,·)=z0，
      把点“拉回”到等高线。
    """
    x, y = x0, 0.0  # 当前追踪点从 S=(x0,0) 出发（在边界 y=0 上）
    t = 0.0  # 用于累计等高线段的“折线弧长”（相邻点欧氏距离之和）

    # 迭代上限：防止步长/精度不合适导致无止境循环
    for _ in range(5_000_000):  # 最多走这么多步，理论上足够覆盖整条相关等高线段
        dx = dh(x, y)  # 近似法向量的 x 分量（辅助导数表达式）
        dy = dh(y, x)  # 近似法向量的 y 分量（利用对称交换 x/y）

        # 切线判据：当法向量与“从目标指向当前点”的向量点积 >= 0 时，认为到达切点
        if dx * (x - gx) + dy * (y - gy) >= 0.0:  # 等价于：沿等高线继续前进将不再“朝向目标”
            return t + dist(gx, gy, x, y)  # 返回：已累计弧长 t + 最后一段直线距离

        # 决策：优先更新“变化更陡”的方向（|d| 更大），通常更稳，减少数值抖动
        if abs(dx) < abs(dy):  # 若 y 方向更“陡”，则优先更新 x（让“拉回”过程更稳定）
            s = -step if x > gx else step  # x 朝向目标方向推进：若当前 x>gx，就往左走，否则往右走
            nx = x + s  # 先把 x 推进一步，得到新的竖线 x=nx

            # 构造一个 y 的试探端点，目的是在 [y, y_guess] 上尽量跨过等高线以便二分夹根
            y_guess = y - sgn(dy * dx) * s  # 方向由 dy*dx 的符号决定（经验性选择，使区间更可能异号）

            # 在固定 x=nx 的竖线上，解出 y 使 h(nx,y)=z0，从而把点拉回等高线
            ny = bisect_root_for_level(z0, fixed_y=nx, a=y, b=y_guess, eps=eps)  # 二分求根（或取更贴近端点）
            nx2, ny2 = nx, ny  # 得到下一步等高线点 P_next=(nx2,ny2)
        else:  # 否则 x 方向更“陡”，则优先更新 y（对称处理）
            s = -step if y > gy else step  # y 朝向目标方向推进：若当前 y>gy，就往下走，否则往上走
            ny = y + s  # 先把 y 推进一步，得到新的横线 y=ny
            x_guess = x - sgn(dy * dx) * s  # 构造 x 的试探端点，用于夹住根（同理，经验性方向选择）
            nx = bisect_root_for_level(z0, fixed_y=ny, a=x, b=x_guess, eps=eps)  # 在固定 y=ny 上解 h(x,ny)=z0
            nx2, ny2 = nx, ny  # 得到下一步等高线点 P_next=(nx2,ny2)

        t += dist(x, y, nx2, ny2)  # 用相邻点欧氏距离近似弧长微元并累加（折线弧长）
        x, y = nx2, ny2  # 更新当前点：继续从新的等高线点出发追踪

    raise RuntimeError("rawrr: exceeded max steps (did not converge)")  # 若走完上限仍未到切点，说明参数需调整


def solve(step: float = 1.0, eps: float = 1e-7) -> float:  # 总流程：先定 f_min，再计算两侧长度并相加
    """
    主流程：
    1) 在边界 y=0 上找到关键点 S=(x0,0)，并取 z0=h(x0,0) 作为最低通行高度候选。
    2) 分别计算：
       - 从 S 沿等高线走到“对 A(200,200) 的切点”所需的弧长 + 切点到 A 的直线距离
       - 从 S 沿等高线走到“对 B(1400,1400) 的切点”所需的弧长 + 切点到 B 的直线距离
    3) 两段相加，得到从 A' 到 B' 的最短路径长度（在高度 z0 的等高线上绕障 + 两端直线连接）。

    step：
    - 等高线追踪的步长，越小越精细但更慢；越大越快但可能丢精度/不稳。
    eps：
    - 二分与关键点定位的收敛阈值，越小越精确但更慢。
    """
    x0, z0 = findz(eps)  # 在 y=0 上定位关键点 S=(x0,0)，并得到候选最低高度 z0=h(x0,0)
    return (  # 返回两段“从同一 S 出发到两个目标”的长度之和（对称拼成 A'->B' 的最短路径长度）
        get_length(step, eps, x0, z0, 200.0, 200.0)  # 计算：S -> 切点(对A) 的弧长 + 切点 -> A 的直线
        + get_length(step, eps, x0, z0, 1400.0, 1400.0)  # 计算：S -> 切点(对B) 的弧长 + 切点 -> B 的直线
    )


if __name__ == "__main__":  # Python 入口：脚本直接运行时才执行下面内容
    ans = solve(step=1.0, eps=1e-7)  # 调用主流程求解；step/eps 控制速度与精度
    print(f"{ans:.3f}")  # 按题意保留三位小数输出