8

我目前正在收紧浮点数字以估计值。（这是：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)不是先决条件。这只是我使用的一种“可能足够好”的方法，但没有为它提供坚实的基础。我可以说分子和分母不会是特殊值：NaN、Infs 等，分母也不会是零。

4

1 回答 1

10

第一：我知道您不想设置舍入模式，但确实应该说，就精度而言，正如其他人所指出的那样，设置舍入模式将产生尽可能好的答案。具体来说，假设x和y都是肯定的（似乎是这种情况，但在您的问题中没有明确说明），以下是具有预期效果的标准 C 代码片段[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/yreal 总是大于或等于无限精确的 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.
于 2012-05-10T14:47:59.377 回答