不难想到需要先生成所有的整点,再用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}")
|