largest known errors

largest known errors

[Paul Zimmermann]

       Hi,

on https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html it is written:

"The table lists the maximum error for each function which is exposed by one of the existing tests in the test suite." [*]

Function i686 ix86 x86_64
j0f      9    9    9
j0       5    5    3

As a consequence, people conclude that the maximum error is 9 ulps for all
GNU libc functions, see https://dl.acm.org/doi/pdf/10.1145/3624062.3624166
end of page 893, right column.

But we know the largest error for j0 on x86_64 is at least 4.51e14 ulps [1].
See also https://sourceware.org/bugzilla/show_bug.cgi?id=27670.
Moreover the test suite contains many "xfail" entries, for example:

j0 0x1.31ec02p+1 xfail:binary64

If one adds j0 0x1.33d152e971b4p+1 in auto-libm-test-in and reruns make
regen-ulps, one gets as output:

Failure: Test: j0 (0x2.67a2a5d2e368p+0)
Result:
 is:         -5.7480381767140581e-17  -0x1.0914db68ef3e2p-54
 should be:  -6.1087652597367303e-17  -0x1.19b7921f03c8ep-54
 difference:  3.6072708302267220e-18   0x1.0a2b6b6148ac0p-58
 ulp       :  292656578119852.0000
 max.ulp   :  3.0000

but the final result says:

Test suite completed:
  0 test cases plus 0 tests for exception flags and
    0 tests for errno executed.
  All tests passed successfully.

I believe we should amend the sentence [*], and/or replace the bounds that we
know are wrong by some special value, say U for Unbounded.

Paul

[1] https://members.loria.fr/PZimmermann/papers/accuracy.pdf

Re: largest known errors

[Carlos O'Donell]

On 12/13/23 03:49, Paul Zimmermann wrote:
>        Hi,
> 
> on https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html it is written:
> 
> "The table lists the maximum error for each function which is exposed by one of the existing tests in the test suite." [*]
> 
> Function i686 ix86 x86_64
> j0f      9    9    9
> j0       5    5    3
> 
> As a consequence, people conclude that the maximum error is 9 ulps for all
> GNU libc functions, see https://dl.acm.org/doi/pdf/10.1145/3624062.3624166
> end of page 893, right column.

In general they would be correct to assume that given the manual information.

> But we know the largest error for j0 on x86_64 is at least 4.51e14 ulps [1].
> See also https://sourceware.org/bugzilla/show_bug.cgi?id=27670.
> Moreover the test suite contains many "xfail" entries, for example:
> 
> j0 0x1.31ec02p+1 xfail:binary64
> 
> If one adds j0 0x1.33d152e971b4p+1 in auto-libm-test-in and reruns make
> regen-ulps, one gets as output:
> 
> Failure: Test: j0 (0x2.67a2a5d2e368p+0)
> Result:
>  is:         -5.7480381767140581e-17  -0x1.0914db68ef3e2p-54
>  should be:  -6.1087652597367303e-17  -0x1.19b7921f03c8ep-54
>  difference:  3.6072708302267220e-18   0x1.0a2b6b6148ac0p-58
>  ulp       :  292656578119852.0000
>  max.ulp   :  3.0000
> 
> but the final result says:
> 
> Test suite completed:
>   0 test cases plus 0 tests for exception flags and
>     0 tests for errno executed.
>   All tests passed successfully.
> 
> I believe we should amend the sentence [*], and/or replace the bounds that we
> know are wrong by some special value, say U for Unbounded.

There are two issues here and the nuances around them matter to me.

(a) There are known defects where ULPs may reach values that are not useful
    for talking about the library in general.

(b) There is value in being clear about the worst case known ULPs for an
    implementation of a given algorithm.

If a test is marked as XFAIL then it is clearly (a) and listing that worst
case ULPs in the manual may not be useful.

If the test is not marked as XFAIL then it is clearly in (b) and we should
list it in the manual as the worst case known ULPS because that is what
the currently implemented algorithm does.

Lastly, all XFAIL entries should reference bugs in our bug tracker, and
if they don't then we should create them to track and resolve the bug.

In closing...

How might we amend the sentence to better explain (a) and (b)?
 
> Paul
> 
> [1] https://members.loria.fr/PZimmermann/papers/accuracy.pdf
> 

-- 
Cheers,
Carlos.

Re: largest known errors

[Joseph Myers]

On Wed, 13 Dec 2023, Carlos O'Donell wrote:

> There are two issues here and the nuances around them matter to me.
> 
> (a) There are known defects where ULPs may reach values that are not useful
>     for talking about the library in general.
> 
> (b) There is value in being clear about the worst case known ULPs for an
>     implementation of a given algorithm.
> 
> If a test is marked as XFAIL then it is clearly (a) and listing that worst
> case ULPs in the manual may not be useful.
> 
> If the test is not marked as XFAIL then it is clearly in (b) and we should
> list it in the manual as the worst case known ULPS because that is what
> the currently implemented algorithm does.
> 
> Lastly, all XFAIL entries should reference bugs in our bug tracker, and
> if they don't then we should create them to track and resolve the bug.

Also, a huge table by architecture (when many libm-test-ulps files may not 
get reliably updated) may not be the most helpful way of presenting this 
information.

In most case, it might be better for libm-test-ulps data to be by 
floating-point format (rather than type; for the narrowing functions, 
there are two formats involved, but ulps should be zero for those 
everywhere except when narrowing from IBM long double), but shared between 
architectures.  Although there are some architecture-specific 
implementations and variation between architectures for results (because 
of different implementations, variation in whether fma gets contracted, 
etc.), there aren't so many such variations (especially once we remove 
ia64), and listing the maximum ulps expected for a function for a given 
format might be better than trying to track architecture-specific values 
(just as we got rid of ulps entries for individual test inputs a long time 
ago).  If we made that change, maybe the vector function ulps would still 
sensibly be architecture-specific; and initial entries for an 
architecture-independent libm-test-ulps file might better be determined by 
what ulps actually appear on a few architectures than by taking the 
maximum across existing libm-test-ulps files, many of which are not 
well-maintained.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: largest known errors

[Paul Zimmermann]

       Hi Joseph,

> Date: Wed, 13 Dec 2023 21:40:04 +0000
> From: Joseph Myers <joseph@codesourcery.com>
> CC: Paul Zimmermann <Paul.Zimmermann@inria.fr>, <libc-alpha@sourceware.org>
> 
> On Wed, 13 Dec 2023, Carlos O'Donell wrote:
> 
> > There are two issues here and the nuances around them matter to me.
> > 
> > (a) There are known defects where ULPs may reach values that are not useful
> >     for talking about the library in general.
> > 
> > (b) There is value in being clear about the worst case known ULPs for an
> >     implementation of a given algorithm.
> > 
> > If a test is marked as XFAIL then it is clearly (a) and listing that worst
> > case ULPs in the manual may not be useful.
> > 
> > If the test is not marked as XFAIL then it is clearly in (b) and we should
> > list it in the manual as the worst case known ULPS because that is what
> > the currently implemented algorithm does.
> > 
> > Lastly, all XFAIL entries should reference bugs in our bug tracker, and
> > if they don't then we should create them to track and resolve the bug.
> 
> Also, a huge table by architecture (when many libm-test-ulps files may not 
> get reliably updated) may not be the most helpful way of presenting this 
> information.
> 
> In most case, it might be better for libm-test-ulps data to be by 
> floating-point format (rather than type; for the narrowing functions, 
> there are two formats involved, but ulps should be zero for those 
> everywhere except when narrowing from IBM long double), but shared between 
> architectures.  Although there are some architecture-specific 
> implementations and variation between architectures for results (because 
> of different implementations, variation in whether fma gets contracted, 
> etc.), there aren't so many such variations (especially once we remove 
> ia64), and listing the maximum ulps expected for a function for a given 
> format might be better than trying to track architecture-specific values 
> (just as we got rid of ulps entries for individual test inputs a long time 
> ago).  If we made that change, maybe the vector function ulps would still 
> sensibly be architecture-specific; and initial entries for an 
> architecture-independent libm-test-ulps file might better be determined by 
> what ulps actually appear on a few architectures than by taking the 
> maximum across existing libm-test-ulps files, many of which are not 
> well-maintained.

that would indeed simplify the presentation, but we still need to find a way
to give accurate information to the user.

I've filled a bug (PR #31165) so that we can keep track of this issue.

Paul

Re: largest known errors

[Paul Zimmermann]

       Hi Carlos,

> Date: Wed, 13 Dec 2023 12:42:11 -0500
> From: Carlos O'Donell <carlos@redhat.com>
> Organization: Red Hat
> 
> On 12/13/23 03:49, Paul Zimmermann wrote:
> >        Hi,
> > 
> > on https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html it is written:
> > 
> > "The table lists the maximum error for each function which is exposed by one of the existing tests in the test suite." [*]
> > 
> > Function i686 ix86 x86_64
> > j0f      9    9    9
> > j0       5    5    3
> > 
> > As a consequence, people conclude that the maximum error is 9 ulps for all
> > GNU libc functions, see https://dl.acm.org/doi/pdf/10.1145/3624062.3624166
> > end of page 893, right column.
> 
> In general they would be correct to assume that given the manual information.
> 
> > But we know the largest error for j0 on x86_64 is at least 4.51e14 ulps [1].
> > See also https://sourceware.org/bugzilla/show_bug.cgi?id=27670.
> > Moreover the test suite contains many "xfail" entries, for example:
> > 
> > j0 0x1.31ec02p+1 xfail:binary64
> > 
> > If one adds j0 0x1.33d152e971b4p+1 in auto-libm-test-in and reruns make
> > regen-ulps, one gets as output:
> > 
> > Failure: Test: j0 (0x2.67a2a5d2e368p+0)
> > Result:
> >  is:         -5.7480381767140581e-17  -0x1.0914db68ef3e2p-54
> >  should be:  -6.1087652597367303e-17  -0x1.19b7921f03c8ep-54
> >  difference:  3.6072708302267220e-18   0x1.0a2b6b6148ac0p-58
> >  ulp       :  292656578119852.0000
> >  max.ulp   :  3.0000
> > 
> > but the final result says:
> > 
> > Test suite completed:
> >   0 test cases plus 0 tests for exception flags and
> >     0 tests for errno executed.
> >   All tests passed successfully.
> > 
> > I believe we should amend the sentence [*], and/or replace the bounds that we
> > know are wrong by some special value, say U for Unbounded.
> 
> There are two issues here and the nuances around them matter to me.
> 
> (a) There are known defects where ULPs may reach values that are not useful
>     for talking about the library in general.
>
> (b) There is value in being clear about the worst case known ULPs for an
>     implementation of a given algorithm.
> 
> If a test is marked as XFAIL then it is clearly (a) and listing that worst
> case ULPs in the manual may not be useful.
> 
> If the test is not marked as XFAIL then it is clearly in (b) and we should
> list it in the manual as the worst case known ULPS because that is what
> the currently implemented algorithm does.
> 
> Lastly, all XFAIL entries should reference bugs in our bug tracker, and
> if they don't then we should create them to track and resolve the bug.
> 
> In closing...
> 
> How might we amend the sentence to better explain (a) and (b)?

so that the manual gives accurate information, one solution would be to write
"The table lists the maximum known error for each function, apart
known defects for some functions and some architectures", and then reference
the bugs corresponding to the known defects.

If we don't want to list all bugs, another possibility would be to replace
the bound by some special value (say NA) that would mean the maximum known
error exceeds the target value. This would concern only a few functions, for
example j0/j1/y0/y1/lgamma in double precision for x64_64 according to
Table 3 of https://members.loria.fr/PZimmermann/papers/accuracy.pdf.

Paul