0353 危险的月球
* * * *
拉格朗日计划
* * * *
危险的月球

把月球视为以$(0,0,0)$为中心、半径为r的球体$C(r)$。

月球上的站点设立在C(r)表面坐标为整数的点处。位于坐标$(0,0,r)$处的站点称为北极站,位于坐标$(0,0,-r)$处的站点称为南极站。

站点两两之间都通过以球心为圆心的圆上的圆弧路径(取弧长较短的一支)相连。站点之间的旅行是危险的,若两站点之间路径的长度为d,则在此两站点之间旅行的危险程度为$(\frac{d}{\pi r})^2$。

若一趟旅程包括了多个站点,则其危险程度为其中每段路径的危险程度之和。

直接从北极站到南极站的旅行,长度为$\pi r$,因此危险程度为1。

从北极站先到$(0,r,0)$站点再到南极站的旅行,尽管总长度不变,但危险程度降低为 $$\left(\frac{\frac{\pi r}{2}}{\pi r}\right)^2+\left(\frac{\frac{\pi r}{2}}{\pi r}\right)^2=\frac{1}{2}.$$ 把从北极站到南极站的所有路径中最低的总危险程度记为$M(r)$,可以验证$M(7)\approx 0.1784943998$(四舍五入到小数点后10位)。

求$\sum_{n=1}^{15}M(2^n-1)$,结果四舍五入到小数点后10位。

本题难度:



解答

不难想到需要先生成所有的整点,再用Dijkstra法求出最短路径(最小风险)。

但本题的上界达到$r=2^{15}-1=32767$,若直接计算,那么整点数量可以达到$O(10^5)$,再用Dijkstra时复杂度就达到$O(10^10)$,因此需要一些优化。

由对称性很容易看出不妨只考虑$x\ge 0, y\ge 0$的情况,进一步地可以验证在$z < 0$的情况也无需计算,这是因为最优路径必定关于赤道对称,只需找出北极点到北半球各点的最短路径再沿赤道作镜面反射即可(此处有最优路径上存在赤道上的整点以及最优路径上不存在赤道上的整点两种情况,因此需要检验北半球所有的整点及其反射点)。

最后,同样由对称性可以只考察$x\le y$的情况,最终结果是1.2759860331。

注:因用到Walrus算符(即:=)以及fstring等新特性,以下代码为Python 3。此外,代码中还打印了进度信息。

import math

risk=lambda p,q,r:(math.acos((p[0]*q[0]+p[1]*q[1]+p[2]*q[2])/r/r)/math.pi)**2

s=0
for i in range(1,16):
    r=(1 << i)-1
    pts=[(x,y,z) for x in range(r+1) for y in range(x,math.isqrt(r*r-x*x)+1) if x*x+y*y+(z:=math.isqrt(r*r-x*x-y*y))*z==r*r]
    n=len(pts)
    tick=n//100
    print(r,n)
    dist=[risk((0,0,r),q,r) for q in pts]
    visited=[0]*n
    unvisited=n
    while unvisited:
        d,i=min((d,i) for i,d in enumerate(dist) if not visited[i])
        visited[i]=1
        unvisited-=1
        for j in range(n):
            if not visited[j]: 
                dist[j]=min(dist[j],d+risk(pts[i],pts[j],r))
        if n>=10000 and unvisited%tick==0:
            print(f"r={r}, with {unvisited//tick} percent left")
    d=min(2*d+risk(pts[i],(pts[i][0],pts[i][1],-pts[i][2]),r) for i,d in enumerate(dist))
    print(r,n,f"{d:.10f}")
    s+=d

print(f"{s:.10f}")