Re: Comparisons of GNU math functions

[Ryan Arnold <ryan dot arnold at gmail dot com>]

On Thu, Mar 26, 2009 at 1:18 PM, JohnT <jrt@worldlinc.net> wrote:
> Regarding the odd math results I reported not too long ago, I was
> wondering what the differences are between the math functions of the GNU
> Scientific Library (libgsl, I think) and those in glibc. Does anyone
> compare results of standard "production" libraries to the results of
> high-precision NIST/ISO results with perhaps 256 significant bits or more?

I don't think anyone has expressed interest (i.e. volunteered) in
making these comparisons.  We don't have any volunteers that are
solely focused on math library precision/results so these problems
tend to take a long time to resolve.  Of course, there's going to be
some difference in the last few digits between 'double precision'
float and 256-bit float due to potentially accrued rounding error.

> One potential source of error that gcc indicated by warning messages
> about == and != is that floating-point comparisons are risky because of
> extended-precision bits that some processors use by default for
> calculations. The first 64 bits might be identical while the last 16
> differ, so a native result would differ from an IEEE 64-bit result. An
> "immediate-mode" constant in a binary surely would have no more than 64
> bits of precision, while it might be the result of a source-code
> operation producing a constant result of 80 internal bits. Dropping the
> least significant 16 bits during compilation in GCC 4.3.x as the
> operation is converted to a constant in the object code could lead to
> errors.
>
> A more suitable way than == or != to do floating-point comparisons in
> glibc test routines might be to evaluate (a - b) < epsilon where epsilon
> is the desired accuracy. ÂAre there any published standards on this
> question?
>
> John T

Comparison against a constant is tricky.  I believe you can create C99
hex floating point constants.  So if you know your current rounding
mode and you know how your operation is going to (accrue error)/round
then you can create an appropriate constant to compare equivalence.

I think in GLIBC we often avoid doing direct floating point
equivalence checks when implementing the C-spec math library
functions.

We're more concerned with the attributes of the inputs and outputs,
i.e. is it a NaN, isfinite, does it exceed the max, and we return the
computed result to user apps.

GLIBC is aware that there's a general acceptable level of imprecision
for floating point functions (as defined by per-architecture ulps file
used by make check).  It is probably wise for a user app to evaluate
using the method you've described if direct comparisons are required
against functions which have a natural and/or unpredictable
imprecision. You can try to define your 'epsilon' by deciphering the
GLIBC ulps files to figure out what the ulps are for each function.

Often 'acceptable' precision is defined by: "I was only able to get it
this precise before I ran out of time and or ideas."

Ryan S. Arnold