Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

Keith Packard via Newlib <[email protected]> Thu, 27 Aug 2020 10:05:59 -0700
Newsgroups	gmane.comp.lib.newlib
Message-ID	<[email protected]>
Corinna Vinschen via Newlib <[email protected]> writes:

> Nevertheless, the symbols have to be exported for backward compatibility.
> So you're saying aliasing gamma_r to lgamma_r and gammaf_r to lgammaf_r
> in the DLLs export table would be sufficient?

It depends if you want the pre-2002 functionality or the post-2002
functionality.

Before 2002, gamma_r and gammaf_r were as BSD originally specified them,
returning the log of gamma (ln(|Γ(x)|)) with the sign of gamma stored
through the second, int *, argument. In 2002, this patch was made to
er_gamma.c:

$ git show 0953fe640f177b565578ed7ecc77169ec1a914fa er_gamma.c

diff --git a/newlib/libm/math/er_gamma.c b/newlib/libm/math/er_gamma.c
index a7183c50f..3c0e241e5 100644
@@ -28,5 +28,5 @@
 	double x; int *signgamp;
 #endif
 {
-	return __ieee754_lgamma_r(x,signgamp);
+	return __ieee754_exp (__ieee754_lgamma_r(x,signgamp));
 }

Before this patch, __ieee754_gamma_r was simply another name for
__ieee754_lgamma_r. After this patch, __ieee754_gamma_r became something
that doesn't exist in any other math library: it returns |Γ(x)| and
stores the sign through the second, int *, argument. Internally, this is
used in w_gamma.c, w_tgamma.c, wr_gamma.c.

w_tgamma.c uses it correctly:

	y = __ieee754_gamma_r(x,&local_signgam);
	if (local_signgam < 0) y = -y;

w_gamma.c, which exports the 'gamma' function uses it *incorrectly* --
it didn't change after the above patch, and so it unexpectedly changed
from returning ln(|Γ(x)|) to returning |Γ(x)|. Applications using
'gamma' instead of 'tgamma' will be getting the wrong answer for some
parameters.

I have to assume this was just an oversight. In 2009, another patch to
w_gamma.c adds a comment explicitly stating that gamma and gamma_r
return ln(|Γ(x)|), when in fact they had been changed to
|Γ(x)| due to the change to __ieee754_gamma_r in the above patch.

Even today, wr_gamma.c says:

	double gamma_r(double x, int *signgamp) /* wrapper lgamma_r */

which is incorrect, as it calls __ieee754_gamma_r.

I see a couple of options here:

 1) Leave things as they are. This makes the 'gamma' family of functions
    incompatible with all other libc implementations and creates a trap
    for applications expecting either lgamma or tgamma values.

 2) Finish the work started in 2002 and make the 'gamma' family
    compatible with their 'tgamma' equivalents. This is what my patch
    does. This makes newlib compatible with 4.4 BSD libc, and leaves
    'gamma' returning the same value anytime that value is non-negative.

 3) Revert the 2002 change so that all of the 'gamma' family return the
    same value as their lgamma equivalents. That would make gamma_r
    useful again, as an alias for lgamma_r. This would make newlib
    compatible with glibc.

 4) Remove the 'gamma' family completely. Given that the meaning of
    'gamma' differs between 4.4 BSD and glibc, it would be safest to
    force applications to select between tgamma and lgamma.

I chose option 2 because that offers the best API compatibility -- the
gamma functions return the same value for most arguments (any positive
arguments, and half of negative arguments).

Any code using the gamma_r functions are likely expecting that to return
lgamma_r (as documented in several places). I can't see how we can leave
that function in the library as-is responsibly.

-- 
-keith
signature.asc (application/pgp-signature, 832 B)