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

Question

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

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

这适用于加法和乘法。另一方面，平方根运算在加法和乘法方面是不对称的。

我突然想到，为了计算 sqrt_RD，对于下限，以下习语尽管很复杂，但在具有 IEEE 754 双精度和 @ 的普通平台上可能会更快987654321@定义为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

我认为你应该可以使用：

/* 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² > @ 987654326@（因为w也是在RU模式下计算的）。如果u 是准确的，那么w 也是准确的（因为我们知道它必须可以表示为双精度）。

Answer 2

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

如果你的候选项太大，那么无限精确的结果就大于0，既然你是向上取整，就会向上取整。即使它只比零大一点点。为了验证这一点，首先尝试 l = 1 + 2eps，其中 (1 + eps) = sqrt (1 + 2eps + eps^2) 有点太大了；然后将 l 缩小 4 的负幂，使 eps^2 超出非规范化数字的分辨率，并检查它。