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

Question

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

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

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

我突然想到，为了计算下限的 sqrt_RD，尽管它很复杂，但在具有 IEEE 754 double 和 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

我认为你应该能够使用:

/* 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 也是精确的(因为我们知道它必须可以表示为 double )。