On Dec 4, 2006, at 20:19, Howard Hinnant wrote:
If that is the question, I'm afraid your answer is not accurate.
In the example I showed the difference is 2 ulp. The difference
appears to grow with the magnitude of the argument. On my systems,
when the argument is DBL_MAX, the difference is 75 ulp.
pow(DBL_MAX, 1./3.) = 0x1.428a2f98d7240p+341
cbrt(DBL_MAX) = 0x1.428a2f98d728bp+341
And yes, I agree with you about the C99 standard. It allows the
vendor to compute pretty much any answer it wants from either pow
or cbrt. Accuracy is not mandated. And I'm not trying to mandate
accuracy for Gfortran either. I just had a knee jerk reaction when
I read that pow(x, 1./3.) could be optimized to cbrt(x) (and on re-
reading, perhaps I inferred too much right there). This isn't just
an optimization. It is also an approximation. Perhaps that is
acceptable. I'm only highlighting the fact in case it might be
important but not recognized.
This is really a very similar case as using a fused multiply-add
instead of two separate instructions with a separate rounding.
In computing cbrt(x) instead of pow(x, 1.0/3.0), we omit the rounding
of 1.0/3.0 and only round after computing the cube root.
As you should for IEEE double, for cbrt, the relative difference is
still quite small. For fused multiply-add, cancellation can make
the relative difference as arbitrarily large.
Still this is not nearly as bad as the re-association of summands,
which throws away any hope of being able to determine accuracy of
computations due to catastrophic cancellation. This is the
one of the biggest issues with -funsafe-math-optimizations.
Maybe we should have an option between IEEE math and
-funsafe--math-optimizations that says for the standard
arithmetic operations *, / and sqrt, the result can either
be correctly rounded floating-point value or the exact result?
This would allow useful optimizations such as
X * 4.0 / 2.0 -> X + X (otherwise invalid in case of overflow)
X / 2.0 * 4.0 -> X + X (otherwise invalid in case of overflow)
X * Y + Z -> fma(X, Y, Z)
X + Y - Y -> X
X / (1.0 / Y) -> X * Y
pow((sqrt (X), 2.0) -> X
sqrt (X * X) -> X
pow (x, 1.0 / 3.0) -> cbrt (x)
In addition we should have a built-in identity function that
only explicitly evaluate its argument to a correctly rounded
floating-point number. That would allow users (or front ends)
to selectively avoid contractions where necessary. Currently
the only brute force workaround is using volatile variables.
It would be great to keep the value in registers instead,
but still prevent other clever optimizations from wreaking havoc.
-Geert
