On 2012-02-13 15:00:54 -0500, Geert Bosch wrote:
> Properties:
>
> [ ] Conforms to C99 for exceptional values
> (accepting/producing NaNs, infinities)
>
> [ ] Handles non-default rounding modes,
> trapping math, errno, etc.
>
> [ ] Requires IEEE compliant binary64 arithmetic
> (no implicit extended range or precision)
>
> [ ] Requires IEEE compliant binary80 arithmetic
> (I know, not official format, but YKWIM)
Please do not use the term binary80, as it is confusing (and
there is a difference between this format and the formats of
the IEEE binary{k} class concerning the implicit bit). You'd
rather say: Intel (or x86/x87) extended precision. As some
platforms use binary128 instead of the Intel extended precision,
you may want to add:
[ ] Requires IEEE compliant binary128 arithmetic
> Accuracy level:
>
> 0 - Correctly rounded
>
> 1 - Faithfully rounded, preserving symmetry and monotonicity
Between 1 and 2, I would add:
Faithfully rounded, preserving symmetry
The reason is that ensuring monotonicity can be difficult for some
functions, where the scaled derivative can be very small (IIRC, if
you don't control monotonicity in the algorithm and only rely on
the error bound, it can be as difficult as correct rounding). That's
why, in the IEEE 754R discussions, when I proposed it as a fallback
when correct rounding wasn't supported, it was rejected.
However there are functions for which statically-proved correct
rounding is very difficult but monotonicity can easily be proved:
sin, cos, tan.
Other other hand, symmetry is always trivial to support and generally
is a direct consequence of the first simplifications in an algorithm.
> 2 - Tightly approximated, meeting prescribed relative error
> bounds. Conforming to OpenCL and Ada Annex G "strict mode"
> numerical bounds.
>
> 3 - Unspecified error bounds
>
> Note that currently of all different operating systems we (AdaCore)
> support for the GNAT compiler, I don't know of any where we can rely
> on the system math library to meet level 2 requirements for all
> functions and over all of the range! Much of this is due to us
> supporting OS versions as long as there is vendor support, so while
> SPARC Solaris 9 and later are fine, Solaris 8 had some issues.
> GNU Linux is quite good, but has issues with the "pow" function for
> large exponents, even in current versions,
pow(x,y) where x is close to 1 and y is large:
http://sourceware.org/bugzilla/show_bug.cgi?id=706
(waiting for a patch).
> and even though Ada
> allows a relative error of (4.0 + |Right · log(Left)| / 32.0)
> for this function, or up to 310 ulps at the end of the range.
> Similarly, for trigonometric functions, the relative error for
> arguments larger than some implementation-defined angle threshold
> is not specified, though the angle threshold needs to be at
> least radix**(mantissa bits / 2), so 2.0**12 for binary32 or 2.0**32
> for binary80. OpenCL doesn't specify an angle threshold, but I
> doubt they intend to require sufficient accurate reduction over
> the entire range to have a final error of 4 ulps: that doesn't fit
> with the rest of the requirements.
P1788 (interval arithmetic standard) will also have definitions for
accuracy levels, so that sin([x,x]) for x large can still be very
fast (except in the tighest level, which basically corresponds to
correct rounding). There will not be a threshold, but the accuracy
will be defined by considering something like [prev(x),succ(x)].
In any case, the spec should be sufficiently clear so that one can
do a proved error analysis.
> The Ada test suite (ACATS) already has quite extensive tests for (2),
> which are automatically parameterized for any type (including those
> with non-binary radix). That, and the fact that this level apparently
> is still a challenge for most system math libraries seems like a good
> reason for aiming for this level as a base level.
>
> We should consider any libm function not meeting these minimal
> criteria to be buggy (or: only usable with -ffast-math, if actually fast),
> and use a replacement. We could have -fstrict-math select the most
> accurate functions we have, with the actual level being more for
> documentation purposes than anything else.
IEEE 754 recommends correct rounding (which should not be much slower
than a function accurate to a few ulp's, in average) in the full range.
I think this should be the default. The best compromise between speed
and accuracy depends on the application, and the compiler can't guess
anyway.
--
Vincent Lefèvre <vinc...@vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)
