在 RU 模式下使用 FPU 计算 RD(sqrt(x))

Question

浮点区间的区间可用于过度逼近实数集，只要任何结果区间的上限是向上舍入计算的，下限是向下舍入的。

一个推荐的技巧是实际计算下界的否定。这允许 FPU 始终保持向上舍入（例如，“Handbook of Floating-Point Arithmetic”，2.9.2）。

这适用于加法和乘法。另一方面，平方根运算不像加法和乘法那样对称。

我突然想到，为了计算 sqrt_RD 的下界，以下惯用语尽管很复杂，但在 IEEE 754 的普通平台上可能更快双精度和 FLT_EVAL_METHOD 定义为 0 比更改舍入模式两次：

#include <fenv.h>
#include <math.h>
#pragma STDC FENV_ACCESS ON
…
/* assumes round-upwards */
double sqrt_rd(double l) { 
  feclearexcept(FE_INEXACT);
  double candidate = sqrt(l);
  if (fetestexcept(FE_INEXACT))
    return nextafter(candidate, 0);
  return candidate;
}

我想知道这是否更好，是否是最快的。作为一种可能的替代方案，但不一定是最快的，在我看来 FMA_RU(candidate, candidate, -l) 可能并不总是准确的（因为有向舍入）但可能在 0 左右足够准确，以便以下工作：

/* assumes round-upwards */
double sqrt_rd(double l) { 
  double candidate = sqrt(l);
  if (fma(candidate, candidate, -l) != 0.0)
    return nextafter(candidate, 0);
  return candidate;
}

还有哪些廉价的方法可以检测到 sqrt 不准确？什么样的浮点运算组合可以在设置为向上舍入的现代 FPU 上实现最快的 sqrt_rd 计算？

Answer 1

fma 以无限精度计算结果，然后应用舍入模式。

如果你的candidate太大，那么无限精度的结果就大于0，既然你是四舍五入，那就四舍五入。即使它只比零大一点点。要验证这一点，首先尝试 l = 1 + 2eps，其中 (1 + eps) = sqrt (1 + 2eps + eps^2) 只是有点太大了；然后按 4 的负幂缩小 l ，使 eps^2 远远超出非规范化数字的分辨率，并检查它。

Answer 2

我认为你应该可以使用：

/* assumes round-upwards */
double sqrt_rd(double l) { 
  double u = sqrt(l);
  double w = u*u;
  if (w != l)
    return nextafter(u, 0);
  return u;
}

这里的理由是，如果 u 不精确，那么它将严格大于 √l，这反过来意味着 w >= u ² > l（因为 w 也在 RU 模式下计算）。如果 u 是精确的，那么 w 也是精确的（因为我们知道它必须可以表示为双精度数）。

在 RU 模式下使用 FPU 计算 RD(sqrt(x))

Computing RD(sqrt(x)) with a FPU in RU mode

c

floating-point

c99

ieee-754