libs/newlib-cygwin
libs/newlib-cygwin
4
0
Fork 0
mirror of git://sourceware.org/git/newlib-cygwin.git synced 2025-02-08 18:19:08 +08:00
Code Issues Projects Releases Wiki Activity

libm: Fix 'gamma' and 'gammaf' functions. Clean up other gamma code. [v2]

Browse Source 
The current gamma, gamma_r, gammaf and gammaf_r functions return
|gamma(x)| instead of ln(|gamma(x)|) due to a change made back in 2002
to the __ieee754_gamma_r implementation. This patch fixes that, making
all of these functions map too their lgamma equivalents.

To fix the underlying bug, the __ieee754_gamma functions have been
changed to return gamma(x), removing the _r variants as those are no
longer necessary. Their names have been changed to __ieee754_tgamma to
avoid potential confusion from users.

Now that the __ieee754_tgamma functions return the correctly signed
value, the tgamma functions have been modified to use them.

libm.a now exposes the following gamma functions:

    ln(|gamma(x)|):

	__ieee754_lgamma_r
	__ieee754_lgammaf_r

	lgamma
	lgamma_r
	gamma
	gamma_r

	lgammaf
	lgammaf_r
	gammaf
	gammaf_r

	lgammal	(on machines where long double is double)

    gamma(x):

	__ieee754_tgamma
	__ieee754_tgammaf
	tgamma
	tgammaf
	tgammal (on machines where long double is double)

Additional aliases for any of the above functions can be added if
necessary; in particular, I'm not sure if we need to include
__ieee754_gamma*_r functions (which would return ln(|(gamma(x)|).

Signed-off-by: Keith Packard <keithp@keithp.com>

----

v2:
	Switch commit message to ASCII
This commit is contained in:
Keith Packard via Newlib 2020-09-04 09:25:27 -07:00 committed by Ken Brown
parent 2d9bb96b51
commit 63c71d3ae4
13 changed files with 56 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
newlib/libc/sys/linux/cmath/math_private.h
View File

@ -178,9 +178,7 @@ extern double __ieee754_cosh (double);
extern double __ieee754_fmod (double,double);
extern double __ieee754_pow (double,double);
extern double __ieee754_lgamma_r (double,int *);
extern double __ieee754_gamma_r (double,int *);
extern double __ieee754_lgamma (double);
extern double __ieee754_gamma (double);
extern double __ieee754_tgamma (double);
extern double __ieee754_log10 (double);
extern double __ieee754_log2 (double);
extern double __ieee754_sinh (double);
@ -241,9 +239,7 @@ extern float __ieee754_coshf (float);
extern float __ieee754_fmodf (float,float);
extern float __ieee754_powf (float,float);
extern float __ieee754_lgammaf_r (float,int *);
extern float __ieee754_gammaf_r (float,int *);
extern float __ieee754_lgammaf (float);
extern float __ieee754_gammaf (float);
extern float __ieee754_tgammaf (float);
extern float __ieee754_log10f (float);
extern float __ieee754_log2f (float);
extern float __ieee754_sinhf (float);
@ -281,9 +277,7 @@ extern long double __ieee754_coshl (long double);
extern long double __ieee754_fmodl (long double,long double);
extern long double __ieee754_powl (long double,long double);
extern long double __ieee754_lgammal_r (long double,int *);
extern long double __ieee754_gammal_r (long double,int *);
extern long double __ieee754_lgammal (long double);
extern long double __ieee754_gammal (long double);
extern long double __ieee754_tgammal (long double);
extern long double __ieee754_log10l (long double);
extern long double __ieee754_log2l (long double);
extern long double __ieee754_sinhl (long double);

4
newlib/libm/math/Makefile.am
View File

@ -8,7 +8,7 @@ src = k_standard.c k_rem_pio2.c \
k_cos.c k_sin.c k_tan.c \
e_acos.c e_acosh.c e_asin.c e_atan2.c \
e_atanh.c e_cosh.c e_exp.c e_fmod.c \
er_gamma.c e_hypot.c e_j0.c \
e_tgamma.c e_hypot.c e_j0.c \
e_j1.c e_jn.c er_lgamma.c \
e_log.c e_log10.c e_pow.c e_rem_pio2.c e_remainder.c \
e_scalb.c e_sinh.c e_sqrt.c \
@ -31,7 +31,7 @@ fsrc = kf_rem_pio2.c \
kf_cos.c kf_sin.c kf_tan.c \
ef_acos.c ef_acosh.c ef_asin.c ef_atan2.c \
ef_atanh.c ef_cosh.c ef_exp.c ef_fmod.c \
erf_gamma.c ef_hypot.c ef_j0.c \
ef_tgamma.c ef_hypot.c ef_j0.c \
ef_j1.c ef_jn.c erf_lgamma.c \
ef_log.c ef_log10.c ef_pow.c ef_rem_pio2.c ef_remainder.c \
ef_scalb.c ef_sinh.c ef_sqrt.c \

28
newlib/libm/math/Makefile.in
View File

@ -78,7 +78,7 @@ am__objects_1 = lib_a-k_standard.$(OBJEXT) lib_a-k_rem_pio2.$(OBJEXT) \
lib_a-e_acosh.$(OBJEXT) lib_a-e_asin.$(OBJEXT) \
lib_a-e_atan2.$(OBJEXT) lib_a-e_atanh.$(OBJEXT) \
lib_a-e_cosh.$(OBJEXT) lib_a-e_exp.$(OBJEXT) \
lib_a-e_fmod.$(OBJEXT) lib_a-er_gamma.$(OBJEXT) \
lib_a-e_fmod.$(OBJEXT) lib_a-e_tgamma.$(OBJEXT) \
lib_a-e_hypot.$(OBJEXT) lib_a-e_j0.$(OBJEXT) \
lib_a-e_j1.$(OBJEXT) lib_a-e_jn.$(OBJEXT) \
lib_a-er_lgamma.$(OBJEXT) lib_a-e_log.$(OBJEXT) \
@ -111,7 +111,7 @@ am__objects_2 = lib_a-kf_rem_pio2.$(OBJEXT) lib_a-kf_cos.$(OBJEXT) \
lib_a-ef_asin.$(OBJEXT) lib_a-ef_atan2.$(OBJEXT) \
lib_a-ef_atanh.$(OBJEXT) lib_a-ef_cosh.$(OBJEXT) \
lib_a-ef_exp.$(OBJEXT) lib_a-ef_fmod.$(OBJEXT) \
lib_a-erf_gamma.$(OBJEXT) lib_a-ef_hypot.$(OBJEXT) \
lib_a-ef_tgamma.$(OBJEXT) lib_a-ef_hypot.$(OBJEXT) \
lib_a-ef_j0.$(OBJEXT) lib_a-ef_j1.$(OBJEXT) \
lib_a-ef_jn.$(OBJEXT) lib_a-erf_lgamma.$(OBJEXT) \
lib_a-ef_log.$(OBJEXT) lib_a-ef_log10.$(OBJEXT) \
@ -147,7 +147,7 @@ LTLIBRARIES = $(noinst_LTLIBRARIES)
libmath_la_LIBADD =
am__objects_4 = k_standard.lo k_rem_pio2.lo k_cos.lo k_sin.lo k_tan.lo \
e_acos.lo e_acosh.lo e_asin.lo e_atan2.lo e_atanh.lo e_cosh.lo \
e_exp.lo e_fmod.lo er_gamma.lo e_hypot.lo e_j0.lo e_j1.lo \
e_exp.lo e_fmod.lo e_tgamma.lo e_hypot.lo e_j0.lo e_j1.lo \
e_jn.lo er_lgamma.lo e_log.lo e_log10.lo e_pow.lo \
e_rem_pio2.lo e_remainder.lo e_scalb.lo e_sinh.lo e_sqrt.lo \
w_acos.lo w_acosh.lo w_asin.lo w_atan2.lo w_atanh.lo w_cosh.lo \
@ -159,7 +159,7 @@ am__objects_4 = k_standard.lo k_rem_pio2.lo k_cos.lo k_sin.lo k_tan.lo \
s_signif.lo s_sin.lo s_tan.lo s_tanh.lo w_exp2.lo w_tgamma.lo
am__objects_5 = kf_rem_pio2.lo kf_cos.lo kf_sin.lo kf_tan.lo \
ef_acos.lo ef_acosh.lo ef_asin.lo ef_atan2.lo ef_atanh.lo \
ef_cosh.lo ef_exp.lo ef_fmod.lo erf_gamma.lo ef_hypot.lo \
ef_cosh.lo ef_exp.lo ef_fmod.lo ef_tgamma.lo ef_hypot.lo \
ef_j0.lo ef_j1.lo ef_jn.lo erf_lgamma.lo ef_log.lo ef_log10.lo \
ef_pow.lo ef_rem_pio2.lo ef_remainder.lo ef_scalb.lo \
ef_sinh.lo ef_sqrt.lo wf_acos.lo wf_acosh.lo wf_asin.lo \
@ -333,7 +333,7 @@ src = k_standard.c k_rem_pio2.c \
k_cos.c k_sin.c k_tan.c \
e_acos.c e_acosh.c e_asin.c e_atan2.c \
e_atanh.c e_cosh.c e_exp.c e_fmod.c \
er_gamma.c e_hypot.c e_j0.c \
e_tgamma.c e_hypot.c e_j0.c \
e_j1.c e_jn.c er_lgamma.c \
e_log.c e_log10.c e_pow.c e_rem_pio2.c e_remainder.c \
e_scalb.c e_sinh.c e_sqrt.c \
@ -356,7 +356,7 @@ fsrc = kf_rem_pio2.c \
kf_cos.c kf_sin.c kf_tan.c \
ef_acos.c ef_acosh.c ef_asin.c ef_atan2.c \
ef_atanh.c ef_cosh.c ef_exp.c ef_fmod.c \
erf_gamma.c ef_hypot.c ef_j0.c \
ef_tgamma.c ef_hypot.c ef_j0.c \
ef_j1.c ef_jn.c erf_lgamma.c \
ef_log.c ef_log10.c ef_pow.c ef_rem_pio2.c ef_remainder.c \
ef_scalb.c ef_sinh.c ef_sqrt.c \
@ -552,11 +552,11 @@ lib_a-e_fmod.o: e_fmod.c
lib_a-e_fmod.obj: e_fmod.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_fmod.obj `if test -f 'e_fmod.c'; then $(CYGPATH_W) 'e_fmod.c'; else $(CYGPATH_W) '$(srcdir)/e_fmod.c'; fi`
lib_a-er_gamma.o: er_gamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-er_gamma.o `test -f 'er_gamma.c' || echo '$(srcdir)/'`er_gamma.c
lib_a-e_tgamma.o: e_tgamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_tgamma.o `test -f 'e_tgamma.c' || echo '$(srcdir)/'`e_tgamma.c
lib_a-er_gamma.obj: er_gamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-er_gamma.obj `if test -f 'er_gamma.c'; then $(CYGPATH_W) 'er_gamma.c'; else $(CYGPATH_W) '$(srcdir)/er_gamma.c'; fi`
lib_a-e_tgamma.obj: e_tgamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_tgamma.obj `if test -f 'e_tgamma.c'; then $(CYGPATH_W) 'e_tgamma.c'; else $(CYGPATH_W) '$(srcdir)/e_tgamma.c'; fi`
lib_a-e_hypot.o: e_hypot.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_hypot.o `test -f 'e_hypot.c' || echo '$(srcdir)/'`e_hypot.c
@ -948,11 +948,11 @@ lib_a-ef_fmod.o: ef_fmod.c
lib_a-ef_fmod.obj: ef_fmod.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_fmod.obj `if test -f 'ef_fmod.c'; then $(CYGPATH_W) 'ef_fmod.c'; else $(CYGPATH_W) '$(srcdir)/ef_fmod.c'; fi`
lib_a-erf_gamma.o: erf_gamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-erf_gamma.o `test -f 'erf_gamma.c' || echo '$(srcdir)/'`erf_gamma.c
lib_a-ef_tgamma.o: ef_tgamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_tgamma.o `test -f 'ef_tgamma.c' || echo '$(srcdir)/'`ef_tgamma.c
lib_a-erf_gamma.obj: erf_gamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-erf_gamma.obj `if test -f 'erf_gamma.c'; then $(CYGPATH_W) 'erf_gamma.c'; else $(CYGPATH_W) '$(srcdir)/erf_gamma.c'; fi`
lib_a-ef_tgamma.obj: ef_tgamma.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_tgamma.obj `if test -f 'ef_tgamma.c'; then $(CYGPATH_W) 'ef_tgamma.c'; else $(CYGPATH_W) '$(srcdir)/ef_tgamma.c'; fi`
lib_a-ef_hypot.o: ef_hypot.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_hypot.o `test -f 'ef_hypot.c' || echo '$(srcdir)/'`ef_hypot.c

19
newlib/libm/math/er_gamma.c → newlib/libm/math/e_tgamma.c
View File

@ -1,5 +1,5 @@
/* @(#)er_gamma.c 5.1 93/09/24 */
/* @(#)e_tgamma.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -12,9 +12,8 @@
*
*/
/* __ieee754_gamma_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
/* __ieee754_tgamma(x)
* Gamma function. Returns gamma(x)
*
* Method: See __ieee754_lgamma_r
*/
@ -22,11 +21,15 @@
#include "fdlibm.h"
#ifdef __STDC__
double __ieee754_gamma_r(double x, int *signgamp)
double __ieee754_tgamma(double x)
#else
double __ieee754_gamma_r(x,signgamp)
double x; int *signgamp;
double __ieee754_tgamma(x)
double x;
#endif
{
return __ieee754_exp (__ieee754_lgamma_r(x,signgamp));
int signgam_local;
double y = __ieee754_exp(__ieee754_lgamma_r(x, &signgam_local));
if (signgam_local < 0)
y = -y;
return y;
}

18
newlib/libm/math/erf_gamma.c → newlib/libm/math/ef_tgamma.c
View File

@ -1,4 +1,4 @@
/* erf_gamma.c -- float version of er_gamma.c.
/* ef_tgamma.c -- float version of e_tgamma.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
@ -14,9 +14,8 @@
*
*/
/* __ieee754_gammaf_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
/* __ieee754_tgammaf(x)
* Float version the Gamma function. Returns gamma(x)
*
* Method: See __ieee754_lgammaf_r
*/
@ -24,11 +23,14 @@
#include "fdlibm.h"
#ifdef __STDC__
float __ieee754_gammaf_r(float x, int *signgamp)
float __ieee754_tgammaf(float x)
#else
float __ieee754_gammaf_r(x,signgamp)
float x; int *signgamp;
float __ieee754_tgammaf(x)
#endif
{
return __ieee754_expf (__ieee754_lgammaf_r(x,signgamp));
int signgam_local;
float y = __ieee754_expf(__ieee754_lgammaf_r(x, &signgam_local));
if (signgam_local < 0)
y = -y;
return y;
}

17
newlib/libm/math/w_gamma.c
View File

@ -132,7 +132,7 @@ in terms of the base return values, although the <[signgam]> global for
/* double gamma(double x)
* Return the logarithm of the Gamma function of x.
*
* Method: call gamma_r
* Method: call lgamma
*/
#include "fdlibm.h"
@ -148,18 +148,7 @@ in terms of the base return values, although the <[signgam]> global for
double x;
#endif
{
#ifdef _IEEE_LIBM
return __ieee754_gamma_r(x,&(_REENT_SIGNGAM(_REENT)));
#else
double y;
y = __ieee754_gamma_r(x,&(_REENT_SIGNGAM(_REENT)));
if(_LIB_VERSION == _IEEE_) return y;
if(!finite(y)&&finite(x)) {
/* gamma(finite) overflow */
errno = ERANGE;
}
return y;
#endif
}
return lgamma(x);
}
#endif /* defined(_DOUBLE_IS_32BITS) */

4
newlib/libm/math/w_tgamma.c
View File

@ -27,9 +27,7 @@
#endif
{
double y;
int local_signgam;
y = __ieee754_gamma_r(x,&local_signgam);
if (local_signgam < 0) y = -y;
y = __ieee754_tgamma(x);
#ifdef _IEEE_LIBM
return y;
#else

17
newlib/libm/math/wf_gamma.c
View File

@ -25,19 +25,8 @@
float x;
#endif
{
#ifdef _IEEE_LIBM
return __ieee754_gammaf_r(x,&(_REENT_SIGNGAM(_REENT)));
#else
float y;
y = __ieee754_gammaf_r(x,&(_REENT_SIGNGAM(_REENT)));
if(_LIB_VERSION == _IEEE_) return y;
if(!finitef(y)&&finitef(x)) {
/* gammaf(finite) overflow */
errno = ERANGE;
}
return y;
#endif
}
return lgammaf(x);
}
#ifdef _DOUBLE_IS_32BITS
@ -48,7 +37,7 @@
double x;
#endif
{
return (double) gammaf((float) x);
return (double) lgammaf((float) x);
}
#endif /* defined(_DOUBLE_IS_32BITS) */

4
newlib/libm/math/wf_tgamma.c
View File

@ -24,9 +24,7 @@
#endif
{
float y;
int local_signgam;
y = __ieee754_gammaf_r(x,&local_signgam);
if (local_signgam < 0) y = -y;
y = __ieee754_tgammaf(x);
#ifdef _IEEE_LIBM
return y;
#else

18
newlib/libm/math/wr_gamma.c
View File

@ -27,23 +27,7 @@
double x; int *signgamp;
#endif
{
#ifdef _IEEE_LIBM
return __ieee754_gamma_r(x,signgamp);
#else
double y;
y = __ieee754_gamma_r(x,signgamp);
if(_LIB_VERSION == _IEEE_) return y;
if(!finite(y)&&finite(x)) {
if(floor(x)==x&&x<=0.0)
/* gamma(-integer) or gamma(0) */
errno = EDOM;
else
/* gamma(finite) overflow */
errno = ERANGE;
return HUGE_VALF;
} else
return y;
#endif
return lgamma_r(x, signgamp);
}
#endif /* defined(_DOUBLE_IS_32BITS) */

5
newlib/libm/math/wr_lgamma.c
View File

@ -40,9 +40,8 @@
else
/* lgamma(finite) overflow */
errno = ERANGE;
return HUGE_VAL;
} else
return y;
}
return y;
#endif
}

21
newlib/libm/math/wrf_gamma.c
View File

@ -27,22 +27,5 @@
float x; int *signgamp;
#endif
{
#ifdef _IEEE_LIBM
return __ieee754_gammaf_r(x,signgamp);
#else
float y;
y = __ieee754_gammaf_r(x,signgamp);
if(_LIB_VERSION == _IEEE_) return y;
if(!finitef(y)&&finitef(x)) {
if(floorf(x)==x&&x<=0.0f) {
/* gammaf(-integer) or gamma(0) */
errno = EDOM;
} else {
/* gammaf(finite) overflow */
errno = ERANGE;
}
return HUGE_VALF;
} else
return y;
#endif
}
return lgammaf_r(x, signgamp);
}

7
newlib/libm/math/wrf_lgamma.c
View File

@ -41,8 +41,7 @@
/* lgammaf(finite) overflow */
errno = ERANGE;
}
return HUGE_VALF;
} else
return y;
}
return y;
#endif
}
}