极坐标下两点浮点距离计算的最优公式问询

嘿，处理极坐标下两点距离的浮点计算时，经典公式 D = sqrt(r1*r1 + r2*r2 - 2*r1*r2*cos(a1-a2)) 确实藏着不少容易踩坑的浮点精度问题，我把这些问题按严重程度从高到低梳理出来：

• A/ 亚正规数场景下的NaN错误：当r1=r2且r1*r1是亚正规数时，根号内的参数会变成负数，直接返回NaN；
• B/ 相近不等半径的NaN错误：当r1和r2非常接近但不完全相等时，同样可能导致根号内参数为负，结果返回NaN；
• C/ 灾难性抵消导致精度丢失：当r1与r2数值接近时，公式中的减法操作会引发灾难性抵消，最终结果的相对精度极差；
• D/ 大数值溢出问题：当r1或r2的幅度过大时，计算过程中会出现溢出，结果变成+inf或者NaN。

这些问题靠临时补丁（比如给根号内参数加个极小值兜底）没法彻底解决，现有的改进公式也或多或少存在局限性。这里分享几个更鲁棒的计算方案：

方案1：半角公式转换 #

把原公式通过三角恒等式转换为基于半角的形式：

D = 2 * sqrt(r1*r2) * sin(|a1 - a2| / 2)

这个公式能完美避开r1和r2接近时的灾难性抵消，同时大幅降低溢出风险——因为sqrt(r1*r2)的增长速度远慢于r1*r1或r2*r2。不过要注意两个细节：

• 当角度差为π时，sin(π/2)=1，此时结果为2*sqrt(r1*r2)，和原公式的r1+r2一致，精度没问题；
• 单独处理r1或r2为0的边界情况：此时直接返回非零的那个半径即可，避免不必要的计算。

方案2：转直角坐标计算 #

先把极坐标转成直角坐标，再用直角坐标距离公式：

x1 = r1 * cos(a1), y1 = r1 * sin(a1)
x2 = r2 * cos(a2), y2 = r2 * sin(a2)
D = sqrt((x1 - x2)^2 + (y1 - y2)^2)

这个方法逻辑直观，能避开原公式的大部分缺陷。虽然极坐标转直角坐标时可能存在精度损耗，但现代编译器的数学库对三角函数的实现已经足够鲁棒，只要不是极端场景，精度都能满足需求。另外，计算(x1-x2)^2 + (y1-y2)^2时，可以用一些优化技巧（比如先计算较大的差值再平方）进一步提升精度。

方案3：混合策略计算 #

结合前两种方案的优势，针对不同场景选择最优计算方式：

1. 先判断r1和r2的大小，假设r1 >= r2；
2. 如果r2远小于r1（比如r2 < r1 * 1e-10），直接返回r1（误差可忽略）；
3. 若角度差|a1-a2|很小（比如小于π/10），用半角公式；角度差接近π（比如大于9π/10），用直角坐标转换；中间区间两种方法都可以；
4. 单独处理边界情况：角度差为0时返回|r1-r2|，角度差为π时返回r1+r2，任一半径为0时返回另一半径。

这种混合策略能最大化不同场景下的计算精度和稳定性，是工业级代码中常用的方案。

内容的提问来源于stack exchange，提问作者Pascal Cuoq