Computing RD(sqrt(x)) with a FPU in RU mode
Intervals of floating-point bounds can be used to over-approximate sets of reals, as long as the upper bound of any result interval is computed in round-upwards and the lower bound in round-downwards.
One recommended trick is to actually compute the negation of the lower bound. This allows to keep the FPU in round-upwards at all times (for instance, “Handbook of Floating-Point Arithmetic”, 2.9.2).
This works well for addition and for multiplication. The square root operation, on the other hand, is not symmetrical in the ways addition and multiplication are.
It occurs to me that in order to compute sqrtRD, for the lower bound, the following idiom, despite its complications, might be faster on an ordinary platform with IEEE 754 double-precision and FLT_EVAL_METHOD defined to 0 than changing the rounding mode twice:
#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;
}
I am wondering whether this is better, and whether if it is the fastest yet. As one possible alternative, but still not necessarily the fastest, it seems to me that FMARU(candidate, candidate, -l) might perhaps not be always exact (because of the directed rounding) but might be accurate enough around 0 for the following to work:
/* 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;
}
What other inexpensive ways are there to detect that sqrt was inexact? What combination of floating-point operations leads to the fastest computation of sqrt_rd on a modern FPU set to round upwards?
I think you should be able to use:
/* 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;
}
The justification here being that if u is inexact, then it will be strictly greater than √ l , which in turn implies that w >= u 2 > l (since w is also calculated in RU mode). And if u is exact, then so is w (since we know it must be representable as a double).
fma calculates the result with infinite precision, then applies the rounding mode.
If your candidate is too large, then the infinitely precise result is greater than 0, and since you are rounding up, it will be rounded up. Even if it is only a tiny little bit larger than zero. To verify this, first try l = 1 + 2eps, where (1 + eps) = sqrt (1 + 2eps + eps^2) is just a tiny bit too large; then scale l down by a negative power of 4 so that the eps^2 is way beyond the resolution of denormalised numbers, and check that as well.
链接地址: http://www.djcxy.com/p/85614.html上一篇: 正确舍入的自然记录