时间戳

首页 > 解决方案 > 使用浮点计算极坐标中两点之间的距离

floating-point - 使用浮点计算极坐标中两点之间的距离

问题描述

当平面中两点的坐标以极坐标形式提供时，其中r1, a1和是浮点数，并且目标是将两点之间的距离计算为浮点数，使用下面的数学公式：r2, a2r1, r2, a1, a2

D = sqrt(r1*r1 + r2*r2 - 2*r1*r2*cos(a1-a2))

这个公式可以在这个和 StackOverflow 上的其他几个答案中找到。链接答案中提供的源链接已失效，但您可以在很多数学资源中找到这个公式，例如这个网站。

作为一个浮点公式，这个公式有许多不受欢迎的性质，下面按严重程度递减列出。简而言之，我的问题是：D在上述条件下计算距离的更好公式是什么？

A/在和次正规sqrt时应用于负数r1 = r2r1*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应用于负数r1r2

当r1和r2非常接近时，数学积r1*r1和之间的距离非常接近数学积和r1*r2之间的距离。当这个公共距离对应于一个小的、奇数个半 ULP 时，可能会出现浮点乘法向下舍入、向上舍入和再次向下舍入的情况。在这些条件下，常见的公式再次取负数的平方根。r1*r2r2*r2r1*r1r1*r2r2*r2

双精度的示例值：

  double r1 = 0x1.1c71c71c71d3fp0;
  double r2 = 0x1.1c71c71c71dcbp0;

在编译器资源管理器上运行它。

C/ 当r1接近 时r2，减法放大了乘法中发生的近似值

这实际上是问题 A/ 和 B/ 的相同根本原因的更良性症状。当r1非常接近 时r2，就会发生称为灾难性抵消的现象。减法的相对精度很糟糕（例如，在乘法过程中舍入a1 = a2和r1接近r2可以使减法的结果为 0.0，尽管最佳答案fabs(r1 - r2)是可表示的，并且相对而言更准确。请注意，绝对精度的结果是和的 ULP 的顺序r1，r2这可能仍然没问题。

r1D/或r2大小过大时溢出

如果只有r1*r1或r2*r2溢出，则计算结果+inf可能不是数学距离的最佳可表示近似。

如果r1*r2溢出，则减法的结果为NaN，因此距离计算为NaN。

问题 A/ 和 B/ 使原本不必存在的结果变得毫无意义，并且可以通过计算dq > 0 ? sqrt(dq) : 0而不是dq. 对于导致它们的输入，此更改0.0代替了答案。这个结果有一个无限的相对误差，其他输入的结果也是如此，因为这不能解决问题 C/。

如果程序员期望计算将在触发它的条件下使用，则可以通过缩放来解决问题 D/。就此而言，问题 A/ 也可以通过缩放来解决，但这不会解决问题 B/。

解决所有 AD 问题的完整解决方案可能会涉及更多计算。可能存在几个仅解决部分问题的最佳点，或者或多或少彻底解决问题 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。

标签： floating-pointlanguage-agnosticieee-754numerical-analysis

解决方案

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)

会解决吗？

推荐阅读