issue with tgammaf

issue with tgammaf

[Paul Zimmermann]

       Hi,

in https://sourceware.org/pipermail/newlib/2020/017920.html I reported an
issue with tgammaf(-0).

With Newlib 3.3.0, there seems to be the same issue for negative subnormal
values of the input, where tgammaf gives +inf instead of -inf.

Best regards,
Paul Zimmermann

Re: issue with tgammaf

[Keith Packard]

[-- Attachment #1.1: Type: text/plain, Size: 436 bytes --]

Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:

>        Hi,
>
> in https://sourceware.org/pipermail/newlib/2020/017920.html I reported an
> issue with tgammaf(-0).
>
> With Newlib 3.3.0, there seems to be the same issue for negative subnormal
> values of the input, where tgammaf gives +inf instead of -inf.

We've been fixing a bunch of gamma problems; this is one that got missed
when the library is configured to report errno.


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1.2: 0001-libm-Make-tgamma-small-INFINITY.patch --]
[-- Type: text/x-diff, Size: 876 bytes --]

From 8cefdba40a4382de6b3f44a5d35880026727f1dd Mon Sep 17 00:00:00 2001
From: Keith Packard <keithp@keithp.com>
Date: Tue, 15 Sep 2020 10:10:07 -0700
Subject: [PATCH] libm: Make tgamma(-small) = -INFINITY

Need to copy the argument sign to the output for tgamma(finite)
overflow case.

Signed-off-by: Keith Packard <keithp@keithp.com>
---
 newlib/libm/math/k_standard.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/newlib/libm/math/k_standard.c b/newlib/libm/math/k_standard.c
index 906412ba7..beab51e16 100644
--- a/newlib/libm/math/k_standard.c
+++ b/newlib/libm/math/k_standard.c
@@ -331,7 +331,7 @@ static double zero = 0.0;	/* used as const */
 	    case 40:
 	    case 140:
 		/* gamma(finite) overflow */
-		retval = HUGE_VAL;
+		retval = copysign(HUGE_VAL, x);
 		errno = ERANGE;
 		break;
 	    case 41:
-- 
2.28.0


[-- Attachment #1.3: Type: text/plain, Size: 15 bytes --]


-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: issue with tgammaf

[Paul Zimmermann]

thank you Keith. For your information I have updated my comparison with
Newlib 3.3.0:

https://members.loria.fr/PZimmermann/papers/accuracy.pdf

The differences with the previous (unknown) version of Newlib I used are
that j1 is less accurate (largest error goes from 2.66e6 ulps to 1.68e7 ulps),
and y1 is more accurate (largest error goes from 4.65e8 ulps to 6.18e6 ulps).
(However, those huge errors simply tell us that j1f is completely off, in the
first case for x=0x6.cfd78p+100 j1f returns 0 instead of -0x1p-52.)

I am still unable to compile a program using lgammaf:

$ cat f.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>

#ifdef NEWLIB
/* without this, newlib says: undefined reference to `__errno' */
int errno;
int* __errno () { return &errno; }
#endif

int
main (int argc, char *argv[])
{
  float x = 0x1p-149f;
  float y = lgammaf (x);
  printf ("x=%.9e y=%.9e\n", x, y);
  return 0;
}

$ gcc -DNEWLIB -no-pie f.c /localdisk/zimmerma/newlib-3.3.0/libm.a
/usr/bin/ld: /localdisk/zimmerma/newlib-3.3.0/libm.a(lib_a-wf_lgamma.o): in function `lgammaf':
wf_lgamma.c:(.text+0x7): undefined reference to `_impure_ptr'
collect2: error: ld returned 1 exit status

If someone manages to compile this program on x86_64/Linux, please tell me!

Paul

Re: issue with tgammaf

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 984 bytes --]

Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:

> thank you Keith. For your information I have updated my comparison with
> Newlib 3.3.0:
>
> https://members.loria.fr/PZimmermann/papers/accuracy.pdf
>
> The differences with the previous (unknown) version of Newlib I used are
> that j1 is less accurate (largest error goes from 2.66e6 ulps to 1.68e7 ulps),
> and y1 is more accurate (largest error goes from 4.65e8 ulps to 6.18e6 ulps).
> (However, those huge errors simply tell us that j1f is completely off, in the
> first case for x=0x6.cfd78p+100 j1f returns 0 instead of -0x1p-52.)
>
> I am still unable to compile a program using lgammaf:

You might have more luck with lgammaf_r as that avoids use of
_impure_ptr in newlib.

Alternatively, you can use the picolibc fork of newlib which uses native
thread local storage instead of a giant struct full of library global
data. That has built-in support for compiling on x86_64 for testing.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: issue with tgammaf

[Paul Zimmermann]

> You might have more luck with lgammaf_r as that avoids use of
> _impure_ptr in newlib.

thank you very much Keith, that made it. The results for lgammaf_r are
the same as those with OpenLibm and Musl:

https://members.loria.fr/PZimmermann/papers/accuracy.pdf

Best regards,
Paul

Re: issue with tgammaf

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 1034 bytes --]

Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:

>> You might have more luck with lgammaf_r as that avoids use of
>> _impure_ptr in newlib.
>
> thank you very much Keith, that made it. The results for lgammaf_r are
> the same as those with OpenLibm and Musl:

That isn't surprising; so many libraries math support can be traced back
to the original SunPro code.

Thanks much for doing this analysis; it's nice to see careful
measurements of how accurate each of these C libraries is. Have you
published the test framework so we can replicate these results and
monitor their change over time?

Do you know that newlib has two different sets of 32-bit math functions?
Arm supplied a new set of many common math functions that use 64-bit
doubles for the intermediate computations? I'd be very interested in
measuring the effect of that code on accuracy. The new paths are
selected by compiling newlib with -D__OBSOLETE_MATH_DEFAULT=0 while the
old can be selected with -D__OBSOLETE_MATH_DEFAULT=1

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: issue with tgammaf

[Paul Zimmermann]

       Dear Keith,

> Do you know that newlib has two different sets of 32-bit math functions?
> Arm supplied a new set of many common math functions that use 64-bit
> doubles for the intermediate computations? I'd be very interested in
> measuring the effect of that code on accuracy. The new paths are
> selected by compiling newlib with -D__OBSOLETE_MATH_DEFAULT=0 while the
> old can be selected with -D__OBSOLETE_MATH_DEFAULT=1

here are the differences. For each function, the first line is newlib-4.1.0
with default configuration, the second line is newlib-4.1.0 with
-D__OBSOLETE_MATH_DEFAULT=1. Functions that do not appear (like acos) give
the same results in both cases (thus probably there is no new code).

For each line, "errors" is the number of binary32 inputs that are incorrectly
rounded, "errors2" is the number of inputs with an error >= 2 ulps, and "maxerr"
is the maximal error.

In summary, the old code (__OBSOLETE_MATH_DEFAULT=1) is always better in terms
of accuracy, with huge differences for some functions (cos, exp, exp10, exp2,
log, log2, sin, pow).

For bivariate functions (atan2, hypot, pow), the tests are not exhaustive.

Best regards,
Paul

acosh:
Total: errors=244658623 (5.72%) errors2=2698 maxerr=2.01e+00 ulp(s)
Total: errors=243413455 (5.69%) errors2=2698 maxerr=2.01e+00 ulp(s)

asinh:
Total: errors=542122908 (12.67%) errors2=2748 maxerr=1.78e+00 ulp(s)
Total: errors=539621602 (12.61%) errors2=2748 maxerr=1.78e+00 ulp(s)

cos:
Total: errors=209833072 (4.90%) errors2=6 maxerr=2.91e+00 ulp(s)
Total: errors=28209648 (0.66%) errors2=0 maxerr=7.76e-01 ulp(s)

cosh:
Total: errors=23905668 (0.56%) errors2=7706 maxerr=2.51e+00 ulp(s)
Total: errors=17868390 (0.42%) errors2=3558 maxerr=1.89e+00 ulp(s)

erf:
Total: errors=126741900 (2.96%) errors2=0 maxerr=9.68e-01 ulp(s)
Total: errors=126734730 (2.96%) errors2=0 maxerr=9.68e-01 ulp(s)

erfc:
Total: errors=21247299 (0.50%) errors2=1131209 maxerr=6.39e+01 ulp(s)
Total: errors=20982747 (0.49%) errors2=1067774 maxerr=6.39e+01 ulp(s)

exp:
Total: errors=17982847 (0.42%) errors2=0 maxerr=9.11e-01 ulp(s)
Total: errors=170646 (0.00%) errors2=0 maxerr=5.02e-01 ulp(s)

exp10:
Total: errors=18423203 (0.43%) errors2=0 maxerr=1.06e+00 ulp(s)
Total: errors=172143 (0.00%) errors2=0 maxerr=5.03e-01 ulp(s)

exp2:
Total: errors=18401203 (0.43%) errors2=0 maxerr=1.02e+00 ulp(s)
Total: errors=168362 (0.00%) errors2=0 maxerr=5.02e-01 ulp(s)

j0:
Total: errors=1338235574 (31.28%) errors2=279528826 maxerr=6.18e+06 ulp(s)
Total: errors=1334176546 (31.19%) errors2=269351612 maxerr=6.18e+06 ulp(s)

j1:
Total: errors=1818091384 (42.50%) errors2=1376362116 maxerr=1.68e+07 ulp(s)
Total: errors=1816660580 (42.46%) errors2=1373004084 maxerr=1.68e+07 ulp(s)

lgamma:
Total: errors=510903809 (11.94%) errors2=13277834 maxerr=7.50e+06 ulp(s)
Total: errors=506865167 (11.85%) errors2=13400241 maxerr=7.50e+06 ulp(s)

log:
Total: errors=13363494 (0.31%) errors2=0 maxerr=8.88e-01 ulp(s)
Total: errors=416908 (0.01%) errors2=0 maxerr=8.18e-01 ulp(s)

log10:
Total: errors=30061115 (0.70%) errors2=91958 maxerr=2.10e+00 ulp(s)
Total: errors=29787060 (0.70%) errors2=62225 maxerr=2.07e+00 ulp(s)

log2:
Total: errors=602745869 (14.09%) errors2=258 maxerr=1.65e+00 ulp(s)
Total: errors=313550 (0.01%) errors2=0 maxerr=7.52e-01 ulp(s)

sin:
Total: errors=206155238 (4.82%) errors2=0 maxerr=1.37e+00 ulp(s)
Total: errors=29362804 (0.69%) errors2=0 maxerr=5.61e-01 ulp(s)

sinh:
Total: errors=74587762 (1.74%) errors2=38924 maxerr=2.51e+00 ulp(s)
Total: errors=71328304 (1.67%) errors2=34776 maxerr=1.89e+00 ulp(s)

tgamma:
Total: errors=2028164922 (47.41%) errors2=1833526367 maxerr=2.39e+02 ulp(s)
Total: errors=2026865970 (47.38%) errors2=1832940352 maxerr=2.39e+02 ulp(s)

y0:
Total: errors=1306144386 (30.53%) errors2=191859954 maxerr=4.84e+06 ulp(s)
Total: errors=1304302538 (30.49%) errors2=187031173 maxerr=4.84e+06 ulp(s)

y1:
Total: errors=1201178797 (28.08%) errors2=153321647 maxerr=6.18e+06 ulp(s)
Total: errors=1199229206 (28.03%) errors2=148430893 maxerr=6.18e+06 ulp(s)

pow:
pow 0 -1 0x1.d55902p-1,-0x1.fe037ep+9 [1.69e+02] 168.527 168.5268728137016
pow 0 -1 0x1.025736p+0,0x1.309f94p+13 [0.817] 0.816736 0.816736102104187

Re: issue with tgammaf

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 1309 bytes --]

Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:

> In summary, the old code (__OBSOLETE_MATH_DEFAULT=1) is always better in terms
> of accuracy, with huge differences for some functions (cos, exp, exp10, exp2,
> log, log2, sin, pow).

Thanks so much for your evaluation. It's very good to know that by
selecting the original SunPro math functions we're not losing anything
in precision while we target machines with less capable hardware.

Are you able to share the code which performs these tests? I'd love to
be able to run these on a regular basis to ensure that no regressions
occur as we maintain the code.

This also provides a potential starting point for someone interested in
improving the accuracy of this library -- we could clearly use better
bessel and gamma functions.

> For bivariate functions (atan2, hypot, pow), the tests are not
> exhaustive.

Yeah, the search space is a bit large :-)

> tgamma:
> Total: errors=2028164922 (47.41%) errors2=1833526367 maxerr=2.39e+02 ulp(s)
> Total: errors=2026865970 (47.38%) errors2=1832940352 maxerr=2.39e+02 ulp(s)

Given that around half of the possible input values (> 35, or near
negative integers) generate an overflow, it seems like tgamma
essentially *never* gives us an accurate finite result...

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: issue with tgammaf

[Paul Zimmermann]

       Dear Keith,

> Are you able to share the code which performs these tests? I'd love to
> be able to run these on a regular basis to ensure that no regressions
> occur as we maintain the code.

I will send you privately the code, since it is not yet ready for a wide
distribution.

> > tgamma:
> > Total: errors=2028164922 (47.41%) errors2=1833526367 maxerr=2.39e+02 ulp(s)
> > Total: errors=2026865970 (47.38%) errors2=1832940352 maxerr=2.39e+02 ulp(s)
> 
> Given that around half of the possible input values (> 35, or near
> negative integers) generate an overflow, it seems like tgamma
> essentially *never* gives us an accurate finite result...

on a sample of 1/1000 of all values, it appears most errors come from exponents
(ilogbf) between -64 and -7. For example with exponent -64 I get no incorrect
rounding with glibc-2.32 out of 16776 samples, but 16410 with newlib.

For example for x=0x1.00e8p-64 newlib gives 0x1.fe31d6p+63 whereas the correct
rounding is 0x1.fe31a4p+63 according to MPFR.

Best regards,
Paul

Re: issue with tgammaf

[Jeff Johnston]

Patch applied.  Thanks Keith.

-- Jeff J.

On Tue, Sep 15, 2020 at 1:28 PM Keith Packard via Newlib <
newlib@sourceware.org> wrote:

> Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:
>
> >        Hi,
> >
> > in https://sourceware.org/pipermail/newlib/2020/017920.html I reported
> an
> > issue with tgammaf(-0).
> >
> > With Newlib 3.3.0, there seems to be the same issue for negative
> subnormal
> > values of the input, where tgammaf gives +inf instead of -inf.
>
> We've been fixing a bunch of gamma problems; this is one that got missed
> when the library is configured to report errno.
>
>
> --
> -keith
>