diff --git a/ld128/e_lgammal.c b/ld128/e_lgammal_r.c similarity index 98% rename from ld128/e_lgammal.c rename to ld128/e_lgammal_r.c index 5af2057..855ba4c 100644 --- a/ld128/e_lgammal.c +++ b/ld128/e_lgammal_r.c @@ -16,7 +16,7 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ -/* lgammal +/* lgammal_r * * Natural logarithm of gamma function * @@ -24,10 +24,10 @@ * * SYNOPSIS: * - * long double x, y, lgammal(); - * extern int signgam; + * long double x, y, lgammal_r(); + * int signgam; * - * y = lgammal(x); + * y = lgammal_r(x, &signgam); * * * @@ -35,8 +35,7 @@ * * Returns the base e (2.718...) logarithm of the absolute * value of the gamma function of the argument. - * The sign (+1 or -1) of the gamma function is returned in a - * global (extern) variable named signgam. + * The sign (+1 or -1) of the gamma function is returned through signgamp. * * The positive domain is partitioned into numerous segments for approximation. * For x > 10, @@ -757,12 +756,12 @@ deval (long double x, const long double *p, int n) long double -lgammal(long double x) +lgammal_r(long double x, int *signgamp) { long double p, q, w, z, nx; int i, nn; - signgam = 1; + *signgamp = 1; if (! finite (x)) return x * x; @@ -770,7 +769,7 @@ lgammal(long double x) if (x == 0.0L) { if (signbit (x)) - signgam = -1; + *signgamp = -1; return one / fabsl (x); } @@ -782,9 +781,9 @@ lgammal(long double x) return (one / (p - p)); i = p; if ((i & 1) == 0) - signgam = -1; + *signgamp = -1; else - signgam = 1; + *signgamp = 1; z = q - p; if (z > 0.5L) { @@ -793,7 +792,7 @@ lgammal(long double x) } z = q * sinl (PIL * z); if (z == 0.0L) - return (signgam * huge * huge); + return (*signgamp * huge * huge); w = lgammal (q); z = logl (PIL / z) - w; return (z); @@ -1025,7 +1024,7 @@ lgammal(long double x) } if (x > MAXLGM) - return (signgam * huge * huge); + return (*signgamp * huge * huge); q = ls2pi - x; q = (x - 0.5L) * logl (x) + q; diff --git a/ld80/Make.files b/ld80/Make.files index 6e31da3..1198cfb 100644 --- a/ld80/Make.files +++ b/ld80/Make.files @@ -1,6 +1,6 @@ $(CUR_SRCS) += invtrig.c \ e_acoshl.c e_powl.c k_tanl.c s_exp2l.c \ - e_atanhl.c e_lgammal.c e_sinhl.c s_asinhl.c s_expm1l.c \ + e_atanhl.c e_lgammal_r.c e_sinhl.c s_asinhl.c s_expm1l.c \ e_coshl.c e_log10l.c e_tgammal.c \ e_expl.c e_log2l.c k_cosl.c s_log1pl.c s_tanhl.c \ e_logl.c k_sinl.c s_erfl.c diff --git a/ld80/e_lgammal.c b/ld80/e_lgammal_r.c similarity index 98% rename from ld80/e_lgammal.c rename to ld80/e_lgammal_r.c index 0f0970d..ae5f5c5 100644 --- a/ld80/e_lgammal.c +++ b/ld80/e_lgammal_r.c @@ -25,7 +25,7 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ -/* lgammal(x) +/* lgammal_r(x, signgamp) * Reentrant version of the logarithm of the Gamma function * with user provide pointer for the sign of Gamma(x). * @@ -89,7 +89,6 @@ #include #include "math_private.h" -extern int signgam; static const long double half = 0.5L, @@ -267,20 +266,20 @@ sin_pi(long double x) long double -lgammal(long double x) +lgammal_r(long double x, int *signgamp) { long double t, y, z, nadj, p, p1, p2, q, r, w; int i, ix; u_int32_t se, i0, i1; - signgam = 1; + *signgamp = 1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; if ((ix | i0 | i1) == 0) { if (se & 0x8000) - signgam = -1; + *signgamp = -1; return one / fabsl (x); } @@ -294,7 +293,7 @@ lgammal(long double x) { /* |x|<2**-63, return -log(|x|) */ if (se & 0x8000) { - signgam = -1; + *signgamp = -1; return -logl (-x); } else @@ -307,7 +306,7 @@ lgammal(long double x) return one / fabsl (t); /* -integer */ nadj = logl (pi / fabsl (t * x)); if (t < zero) - signgam = -1; + *signgamp = -1; x = -x; } diff --git a/src/Make.files b/src/Make.files index e4b28d0..edd64e4 100644 --- a/src/Make.files +++ b/src/Make.files @@ -4,7 +4,7 @@ $(CUR_SRCS) = common.c \ e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \ e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \ e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \ - e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \ + e_lgammal.c e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \ e_pow.c e_powf.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \ e_rem_pio2.c e_rem_pio2f.c \ e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c \ diff --git a/src/e_lgammal.c b/src/e_lgammal.c new file mode 100644 index 0000000..e608eb1 --- /dev/null +++ b/src/e_lgammal.c @@ -0,0 +1,12 @@ +#include "cdefs-compat.h" +#include "openlibm.h" +#include "math_private.h" + +extern int signgam; + +DLLEXPORT long double +lgammal(long double x) +{ + + return (lgammal_r(x, &signgam)); +} diff --git a/src/openlibm.h b/src/openlibm.h index 315a227..911ce50 100644 --- a/src/openlibm.h +++ b/src/openlibm.h @@ -558,9 +558,13 @@ long double tanhl(long double); long double tanl(long double); long double tgammal(long double); long double truncl(long double); - #endif /* __ISO_C_VISIBLE >= 1999 */ +/* Reentrant version of lgammal. */ +#if __BSD_VISIBLE +long double lgammal_r(long double, int *); +#endif /* __BSD_VISIBLE */ + #include "openlibm_complex.h" #if defined(__cplusplus)