浮点除法 - 避免结果小于“精确”值的偏差

Question

我目前正在收紧浮点数字以估计值。 （对于那些感兴趣的人来说，它是：p(k,t)。）本质上，该实用程序永远不会产生该值的低估：可能的素数生成的安全性取决于数值稳健的实现。虽然输出结果与公布的值一致，但我使用了DBL_EPSILON 值来确保除法，特别是，产生的结果永远不会小于真实值：

考虑：double x, y; /* assigned some values... */

评估：r = x / y; 经常发生，但这些（有限精度）结果可能会从真实结果中截断有效数字 - 可能是无限精度的理性扩展。我目前尝试通过对分子应用偏差来缓解这种情况，即

r = ((1.0 + DBL_EPSILON) * x) / y;

如果您对这个主题有所了解，p(k,t) 通常比大多数估计值要小得多 - 但它根本不足以解决这个“观察”的问题。我当然可以说：

(((1.0 + DBL_EPSILON) * x) / y) >= (x / y)

当然，我需要确保“偏差”结果大于或等于“精确”值。虽然我确信它与操纵或缩放DBL_EPSILON 有关，但我显然希望“有偏差”的结果至少超过“精确”结果 - 在 IEEE-754 算术假设下可以证明。

是的，我查看了 Goldberg 的论文，并寻找了一个可靠的解决方案。 请不要建议操纵舍入模式。理想情况下，我希望得到对浮点定理非常了解的人的回答，或者知道一个很好的示例。

编辑：澄清一下，(((1.0 + DBL_EPSILON) * x) / y) 或表单(((1.0 + c) * x) / y) 不是先决条件。这只是我使用的一种“可能足够好”的方法，但没有为它提供坚实的基础。我可以声明分子和分母不会是特殊值：NaNs、Infs 等，分母也不会是零。

Answer 1

第一：我知道你不想设置舍入模式，但确实应该说 正如其他人所指出的，就精度而言，设置舍入模式将产生尽可能好的答案。具体来说，假设 x 和 y 都是肯定的（似乎是这种情况，但在您的问题中没有明确说明），以下是具有预期效果的标准 C sn-p [1] ：

#include <math.h>
#pragma STDC FENV_ACCESS on

int OldRoundingMode = fegetround();
fesetround(FE_UPWARD);
r = x/y;
fesetround(OldRoundingMode);

现在，除此之外，有正当的理由不想更改舍入模式（某些平台不支持舍入到正无穷大，在某些平台上更改舍入模式会引入较大的序列化停顿等） )，您不应该这样做的愿望不应该如此随便地置之不理。那么，尊重您的问题，我们还能做些什么呢？

如果您的平台支持融合乘加，那么您可以使用一个非常优雅的解决方案：

#include <math.h>
r = x/y;
if (fma(r,y,-x) < 0) r = nextafter(r, INFINITY);

在支持硬件 fma 的平台上，这是非常有效的。即使 fma( ) 是在软件中实现的，它也是可以接受的。这种方法的优点是它可以提供与更改舍入模式相同的结果；也就是尽可能严格的界限。

如果你平台的 C 库是上古时代的，不提供fma，还是有希望的。您声称的陈述是正确的（至少假设没有异常值-我需要更多地考虑异常值会发生什么）； (1.0+DBL_EPSILON)*x/y 确实总是大于或等于无限精确的 x/y。它有时会比具有此属性的最小值大 1 ulp，但这是一个非常小的并且可能可以接受的余量。这些说法的证明非常繁琐，可能不适合 StackOverflow，但我将简要介绍一下：

1. 忽略非规范化，将我们限制在 [1.0, 2.0) 中的 x, y 就足够了。

2. (1.0 + eps)*x >= x + eps > x。要看到这一点，请观察：

   (1.0 + eps)*x = x + x*eps >= x + eps > x.
   

3. 令 P 为数学上精确的 x/y。我们有：

   (1.0 + eps)*x/y >= (x + eps)/y = x/y + eps/y = P + eps/y
   

   现在，y 以 2 为界，所以这给了我们：

   (1.0 + eps)*x/y > P + eps/2
   

   这足以保证结果四舍五入到一个值 >= P。这也向我们展示了更严格的界限。在许多情况下，我们可以改为使用nextafter(x,INFINITY)/y 来获得所需的效果，并具有更严格的界限。 （nextafter(x,INFINITY) 始终是 x + ulp，而(1.0 + eps)*x 将是 x + 2ulp 一半的时间。如果你想避免调用 nextafter 库函数，你可以使用 (x + (0.75*DBL_EPSILON)*x) 来获得相同的结果，在正正常值的工作假设下）。

---

1. 为了真正学究式地正确，这将变得更加复杂。没有人真正写过这样的代码，但它会遵循以下思路：

   #include <math.h>
   #pragma STDC FENV_ACCESS on
   
   #if defined FE_UPWARD
       int OldRoundingMode = fegetround();
       if (OldRoundingMode < 0) goto Error;
       if (fesetround(FE_UPWARD)) goto Error;
       r = x/y;
       if (fesetround(OldRoundingMode)) goto TrulyHosed;
       return r;
   TrulyHosed:
       // we established the desired rounding mode and did our computation,
       // but now we can't set it back to the original mode.  I have no idea
       // how you handle this gracefully.
   Error:
   #else
       // we can't establish the desired rounding mode, so fall back on
       // something else.