使用浮点数计算极坐标中两点之间的距离
Computing the distance between two points in polar coordinates using floating-point
当平面中两点的坐标以极坐标形式提供为r1, a1和r2, a2时，其中r1, r2, a1, a2是浮点数，目标是计算两点之间的距离作为浮点数，使用下面的数学公式可能很诱人：
D = sqrt(r1*r1 + r2*r2 - 2*r1*r2*cos(a1-a2))
这个公式可以在this and several other answers on Whosebug. The link to the source provided in the linked answer is dead, but you find this formula in a lot of mathematical resources, for instance this website中找到。
作为浮点公式，此公式有许多不良特性，按严重程度递减列出如下。简而言之，我的问题是：在上述条件下计算距离 D 的更好公式是什么？
A/ sqrt 当 r1 = r2 和 r1*r1 次正规时应用于负数
当r1*r1为次正规且r1 = r2时，r1*r1的舍入与2*r1*r2的舍入不同。在这个反例范围的一个极端，r1*r1 为零而 2*r1*r2 为非零。在另一个极端，r1*r1 低于正态并且向下四舍五入有些残酷，而 2*r1*r2 是正常的并且四舍五入不那么残酷。无论哪种方式，sqrt 中的表达式最终都是负数，D 的浮点公式的结果是 NaN.
双精度示例值：
double r1 = sqrt(DBL_MIN * DBL_EPSILON) * sqrt(0.45);
double r2 = r1;
B/sqrt应用于负数时r1接近r2而不同
当r1和r2很接近时，数学乘积r1*r1和r1*r2之间的距离与数学乘积[=]之间的距离很近38=] 和 r2*r2。当这个公共距离对应一个小的、奇数个 half-ULPs 时，可能会出现浮点乘法 r1*r1 向下舍入，r1*r2 向上舍入，r2*r2 向下舍入的情况再次。同样，在这些条件下，常见的公式取负数的平方根。
双精度示例值：
double r1 = 0x1.1c71c71c71d3fp0;
double r2 = 0x1.1c71c71c71dcbp0;
C/ 当r1接近r2时，减法放大了乘法
中发生的近似值
这实际上是问题 A/ 和 B/ 的相同根本原因的更良性症状。当 r1 非常接近 r2 时，就会出现称为 catastrophic cancellation 的现象。减法的相对精度很糟糕（例如，a1 = a2 和 r1 接近 r2，乘法期间的舍入可以使减法的结果为 0.0，尽管最佳答案是 fabs(r1 - r2)，相对而言是可表示的并且无限准确。请注意，结果的绝对准确度是 r1 和 r2 的 ULP 的数量级，这可能仍然没问题。
D/当r1或r2的量级过大时溢出
如果只有 r1*r1 或 r2*r2 溢出，结果计算为 +inf，这可能不是数学距离的最佳表示近似值。
如果r1*r2溢出，则减法结果为NaN，因此距离计算为NaN.
问题 A/ 和 B/ 使本来不必如此的结果变得毫无意义，可以通过计算 dq > 0 ? sqrt(dq) : 0 而不是 dq 来解决。对于导致它们的输入，此更改改为 0.0 的答案。这个结果有无限的相对误差，其他输入的结果也是如此，因为这不能解决问题 C/。
如果程序员希望计算将在触发它的条件下使用，则问题 D/ 可以通过缩放来解决。就此而言，问题 A/ 也可以通过缩放来解决，但这不能解决问题 B/。
解决所有问题 A-D 的完整解决方案可能会涉及更多计算。可能存在几个最佳点，它们只解决了一些问题，或者或多或少彻底地解决了问题 C/，计算的距离精确到 10 ULP，或 3，或 1。任何比起点改进的解决方案是值得回答。
Guillaume Melquiond 已经在场外指出，下面的公式在数学方面等同于原始公式，并且它显然避开了问题 A/ 和 B/，因为平方根的参数是非负项的总和：
D = sqrt((r1-r2)*(r1-r2) + 2*r1*r2*(1 - cos(a2-a1)))
在此解决方案中，灾难性取消发生在 1 - cos(a2-a1)，因此问题 C/ 的某些方面仍然存在（尽管使用此公式的计算对于 a1=a2 和 r1 关闭是最佳的到 r2，因为那时 r1-r2 和 cos(a2-a1) 是准确的）。关于问题 D/ 的情况有所改善，但仍然存在结果可表示为有限值且公式计算 +inf.
的情况
我不是这类数值分析方面的专家，但想指出的是，现在有一些计算机程序试图给出关于浮点数的更多 "stable" 版本的公式问题。更好的系统之一是所谓的 Herbie 系统：https://herbie.uwplse.org/
他们确实有一个网络演示，当我插入你的方程式时，它会得出：
上面的 link 提供了大量关于它推荐的翻译的信息。如果 web-link 消失了，这里是它提出的最终公式的屏幕截图：
您也可以在 LaTeX 或 C 中获得相同的输出。它声称错误已从 31.5 减少到 18.5；尽管我不确定这些数字的确切含义。他们确实有一个简短的教程来帮助您入门：https://herbie.uwplse.org/doc/latest/tutorial.html
希望对您有所帮助！
Let b = (a1-a2)/2
then using
cos( a1-a2) = 1 - 2*sin(b)*sin(b)
D = sqrt( (r1-r2)*(r1-r2) + 4*r1*r2*sin(b)*sin(b))
这至少消除了负数的平方根，但仍然会出现溢出问题。
也许
x = 2*sqrt(r1)*sqrt(r2)*sin(b)
D = hypot( r1-r2, x)
能解决吗？
如果我们旋转整个图形，距离保持不变。
使用 pol2cart 的想法，我们可以有一些变化，如：
x1=r1 cos((a2-a1)/2),y1=-r1 sin((a2-a1)/2)
x2=r2 cos((a2-a1)/2),y2= r2 sin((a2-a1)/2)
dist = hypot((r2-r1) cos((a2-a1)/2),(r2+r1) sin((a2-a1)/2)))
如果角度没有限制而不是适当减小，则有必要进一步开发 sin/cos 公式，使该公式不太有用...
r2+r1 也可能溢出，在这种情况下我们可以应用简单的缩放，如 2*hypot((r1/2+r2/2)...)
当平面中两点的坐标以极坐标形式提供为r1, a1和r2, a2时，其中r1, r2, a1, a2是浮点数，目标是计算两点之间的距离作为浮点数，使用下面的数学公式可能很诱人：
D = sqrt(r1*r1 + r2*r2 - 2*r1*r2*cos(a1-a2))
这个公式可以在this and several other answers on Whosebug. The link to the source provided in the linked answer is dead, but you find this formula in a lot of mathematical resources, for instance this website中找到。
作为浮点公式，此公式有许多不良特性，按严重程度递减列出如下。简而言之，我的问题是：在上述条件下计算距离 D 的更好公式是什么？
A/ sqrt 当 r1 = r2 和 r1*r1 次正规时应用于负数
当r1*r1为次正规且r1 = r2时，r1*r1的舍入与2*r1*r2的舍入不同。在这个反例范围的一个极端，r1*r1 为零而 2*r1*r2 为非零。在另一个极端，r1*r1 低于正态并且向下四舍五入有些残酷，而 2*r1*r2 是正常的并且四舍五入不那么残酷。无论哪种方式，sqrt 中的表达式最终都是负数，D 的浮点公式的结果是 NaN.
双精度示例值：
double r1 = sqrt(DBL_MIN * DBL_EPSILON) * sqrt(0.45);
double r2 = r1;
B/sqrt应用于负数时r1接近r2而不同
当r1和r2很接近时，数学乘积r1*r1和r1*r2之间的距离与数学乘积[=]之间的距离很近38=] 和 r2*r2。当这个公共距离对应一个小的、奇数个 half-ULPs 时，可能会出现浮点乘法 r1*r1 向下舍入，r1*r2 向上舍入，r2*r2 向下舍入的情况再次。同样，在这些条件下，常见的公式取负数的平方根。
双精度示例值：
double r1 = 0x1.1c71c71c71d3fp0;
double r2 = 0x1.1c71c71c71dcbp0;
C/ 当r1接近r2时，减法放大了乘法
中发生的近似值
这实际上是问题 A/ 和 B/ 的相同根本原因的更良性症状。当 r1 非常接近 r2 时，就会出现称为 catastrophic cancellation 的现象。减法的相对精度很糟糕（例如，a1 = a2 和 r1 接近 r2，乘法期间的舍入可以使减法的结果为 0.0，尽管最佳答案是 fabs(r1 - r2)，相对而言是可表示的并且无限准确。请注意，结果的绝对准确度是 r1 和 r2 的 ULP 的数量级，这可能仍然没问题。
D/当r1或r2的量级过大时溢出
如果只有 r1*r1 或 r2*r2 溢出，结果计算为 +inf，这可能不是数学距离的最佳表示近似值。
如果r1*r2溢出，则减法结果为NaN，因此距离计算为NaN.
问题 A/ 和 B/ 使本来不必如此的结果变得毫无意义，可以通过计算 dq > 0 ? sqrt(dq) : 0 而不是 dq 来解决。对于导致它们的输入，此更改改为 0.0 的答案。这个结果有无限的相对误差，其他输入的结果也是如此，因为这不能解决问题 C/。
如果程序员希望计算将在触发它的条件下使用，则问题 D/ 可以通过缩放来解决。就此而言，问题 A/ 也可以通过缩放来解决，但这不能解决问题 B/。
解决所有问题 A-D 的完整解决方案可能会涉及更多计算。可能存在几个最佳点，它们只解决了一些问题，或者或多或少彻底地解决了问题 C/，计算的距离精确到 10 ULP，或 3，或 1。任何比起点改进的解决方案是值得回答。
Guillaume Melquiond 已经在场外指出，下面的公式在数学方面等同于原始公式，并且它显然避开了问题 A/ 和 B/，因为平方根的参数是非负项的总和：
D = sqrt((r1-r2)*(r1-r2) + 2*r1*r2*(1 - cos(a2-a1)))
在此解决方案中，灾难性取消发生在 1 - cos(a2-a1)，因此问题 C/ 的某些方面仍然存在（尽管使用此公式的计算对于 a1=a2 和 r1 关闭是最佳的到 r2，因为那时 r1-r2 和 cos(a2-a1) 是准确的）。关于问题 D/ 的情况有所改善，但仍然存在结果可表示为有限值且公式计算 +inf.
我不是这类数值分析方面的专家，但想指出的是，现在有一些计算机程序试图给出关于浮点数的更多 "stable" 版本的公式问题。更好的系统之一是所谓的 Herbie 系统：https://herbie.uwplse.org/
他们确实有一个网络演示，当我插入你的方程式时，它会得出：
上面的 link 提供了大量关于它推荐的翻译的信息。如果 web-link 消失了，这里是它提出的最终公式的屏幕截图：
您也可以在 LaTeX 或 C 中获得相同的输出。它声称错误已从 31.5 减少到 18.5；尽管我不确定这些数字的确切含义。他们确实有一个简短的教程来帮助您入门：https://herbie.uwplse.org/doc/latest/tutorial.html
希望对您有所帮助！
Let b = (a1-a2)/2
then using
cos( a1-a2) = 1 - 2*sin(b)*sin(b)
D = sqrt( (r1-r2)*(r1-r2) + 4*r1*r2*sin(b)*sin(b))
这至少消除了负数的平方根，但仍然会出现溢出问题。
也许
x = 2*sqrt(r1)*sqrt(r2)*sin(b)
D = hypot( r1-r2, x)
能解决吗？
如果我们旋转整个图形，距离保持不变。
使用 pol2cart 的想法，我们可以有一些变化，如：
x1=r1 cos((a2-a1)/2),y1=-r1 sin((a2-a1)/2)
x2=r2 cos((a2-a1)/2),y2= r2 sin((a2-a1)/2)
dist = hypot((r2-r1) cos((a2-a1)/2),(r2+r1) sin((a2-a1)/2)))
如果角度没有限制而不是适当减小，则有必要进一步开发 sin/cos 公式，使该公式不太有用...
r2+r1 也可能溢出，在这种情况下我们可以应用简单的缩放，如 2*hypot((r1/2+r2/2)...)