浮点除法 - 偏向于避免结果小于"精确"值

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片段[1]:

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

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

现在,除此之外,有合理的理由不想改变舍入模式(某些平台不支持round-to-infinity,在某些平台上更改舍入模式会引入大型序列化停顿等),以及你不愿意这样做的愿望不应该随便撇在一边.所以,尊重你的问题,我们还能做些什么？

如果您的平台支持融合乘法加法,那么您可以使用非常优雅的解决方案:

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

在具有硬件fma支持的平台上,这非常有效.即使fma()在软件中实现,也可以接受.这种方法的优点在于它将提供与改变舍入模式相同的结果; 也就是说,最严格的约束可能.

如果您的平台的C库是antediluvian并且没有提供fma,那么仍然有希望.您声明的陈述是正确的(假设没有非正规值,至少 - 我需要更多地考虑非正规的情况); (1.0+DBL_EPSILON)*x/y确实总是大于或等于无限精确的x/y.它有时会比具有此属性的最小值大一个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)*xx + 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.