[PATCH] Constant folding for cabs and fabs.

[Joseph S. Myers]

On Thu, 5 Jun 2003, Roger Sayle wrote:

> The most important issue to remember is that unlike sqrt where the
> relevant standards require 0.5 ulp for float and double, C99 cabs
> and hypot only require 1 ulp precision.  In GCC real.c terms this is

And without defining __STDC_IEC_559__ or __STDC_IEC_559_COMPLEX__ (which
we don't; some glibc headers may, though this is dodgy on x86), almost
everything about floating point is implementation-defined or unspecified;  
this doesn't mean we should use that latitude in the implementation.  For
quality-of-implementation, any floating-point constant folding done
without -funsafe-math-optimizations should be yield 0.5 ulp correctly
rounded values (taking account of the FENV_ACCESS pragma when we have
support for that to disable folding where the rounding mode might make a
difference).

> with is 113 bit mantissae: On these platforms in exceptionally rare
> cases it is possible to incorrectly round the last bit with a probability
> of approximately 2^-54 or 2^-47 respectively, which is less than
> .0000000000000071.  In *all* remaining cases, we get the perfectly

This involves an incorrect assumption that the mantissas are uniformly
distributed.  This isn't the case for much real-world data (though the
nonuniformness in that case is unlikely to affect the probability) and
isn't necessarily the case for data with some mathematical significance
that might just happen to interact badly with the imperfect algorithm.

Also, that probability means there exist a nonzero number of examples
which are incorrectly rounded - which is worse than the probabilistic
methods described in the paper of Ziv referenced in a previous discussion
(ACM Transactions on Mathematical Software, Vol. 17, No. 3, 1991, pp
410-423) (which yield a bounded computation that probably yields correctly
rounded results for all values of a floating point type, and yields an
exception where it is not known to be correctly rounded).

> I'm unaware of a default run-time implementation, either in hardware
> or via software emulation, that provides anywhere near this level of
> accuracy by default for cabs/hypot on any of GCC's supported platforms.

The methods described in Ziv's paper seem to be intended for runtime use,
though I don't know how widely they've been implemented for it.

(For GCC the constraints of bounding resource usage for the computation in
all cases don't apply, so an unlimited number of iterations would be fine.
But for cabs/hypot we can bound the precision needed, it's just that you
might need about 4 times the highest precision of a floating point format
supported for it to be immediately obvious that the results are always
correct.)

-- 
Joseph S. Myers
jsm28@cam.ac.uk