Merge pull request #75 from NuxiNL/signgam

Clean up handling of signgam

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();
+ * long double x, y, lgammal_r();
- * extern int signgam;
+ * 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
+ * The sign (+1 or -1) of the gamma function is returned through signgamp.
- * global (extern) variable named signgam.
  *
  * 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;

ld128/e_tgammal.c

@@ -26,20 +26,14 @@ tgammal(long double x)
 	int64_t i0,i1;
 
 	GET_LDOUBLE_WORDS64(i0,i1,x);
-	if (((i0&0x7fffffffffffffffLL)|i1) == 0) {
+	if (((i0&0x7fffffffffffffffLL)|i1) == 0)
-		signgam = 0;
 		return (1.0/x);
-	}
 
-	if (i0<0 && (u_int64_t)i0<0xffff000000000000ULL && rintl(x)==x) {
+	if (i0<0 && (u_int64_t)i0<0xffff000000000000ULL && rintl(x)==x)
-		signgam = 0;
 		return (x-x)/(x-x);
-	}
 
-	if (i0==0xffff000000000000ULL && i1==0) {
+	if (i0==0xffff000000000000ULL && i1==0)
-		signgam = 0;
 		return (x-x);
-	}
 
 	return expl(lgammal(x));
 }

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

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 <openlibm.h>
 
 #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;
     }
 

ld80/e_tgammal.c

@@ -25,7 +25,6 @@
  * SYNOPSIS:
  *
  * long double x, y, tgammal();
- * extern int signgam;
  *
  * y = tgammal( x );
  *
@@ -33,10 +32,8 @@
  *
  * DESCRIPTION:
  *
- * Returns gamma function of the argument.  The result is
+ * Returns gamma function of the argument.  The result is correctly
- * correctly signed, and the sign (+1 or -1) is also
+ * signed.  This variable is also filled in by the logarithmic gamma
- * returned in a global (extern) variable named signgam.
- * This variable is also filled in by the logarithmic gamma
  * function lgamma().
  *
  * Arguments |x| <= 13 are reduced by recurrence and the function
@@ -61,7 +58,6 @@
 #include <openlibm.h>
 
 #include "math_private.h"
-extern int signgam;
 
 /*
 tgamma(x+2)  = tgamma(x+2) P(x)/Q(x)
@@ -224,7 +220,6 @@ tgammal(long double x)
 long double p, q, z;
 int i;
 
-signgam = 1;
 if( isnan(x) )
 	return(NAN);
 if(x == INFINITY)
@@ -237,6 +232,7 @@ q = fabsl(x);
 
 if( q > 13.0L )
 	{
+	int sign = 1;
 	if( q > MAXGAML )
 		goto goverf;
 	if( x < 0.0L )
@@ -246,7 +242,7 @@ if( q > 13.0L )
 			return (x - x) / (x - x);
 		i = p;
 		if( (i & 1) == 0 )
-			signgam = -1;
+			sign = -1;
 		z = q - p;
 		if( z > 0.5L )
 			{
@@ -258,7 +254,7 @@ if( q > 13.0L )
 		if( z <= PIL/LDBL_MAX )
 			{
 goverf:
-			return( signgam * INFINITY);
+			return( sign * INFINITY);
 			}
 		z = PIL/z;
 		}
@@ -266,7 +262,7 @@ goverf:
 		{
 		z = stirf(x);
 		}
-	return( signgam * z );
+	return( sign * z );
 	}
 
 z = 1.0L;
@@ -298,8 +294,6 @@ x -= 2.0L;
 p = __polevll( x, P, 7 );
 q = __polevll( x, Q, 8 );
 z = z * p / q;
-if( z < 0 )
-	signgam = -1;
 return z;
 
 small:
@@ -311,7 +305,6 @@ else
 		{
 		x = -x;
 		q = z / (x * __polevll( x, SN, 8 ));
-		signgam = -1;
 		}
 	else
 		q = z / (x * __polevll( x, S, 8 ));

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 \

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));
+}

src/openlibm.h

@@ -515,35 +515,23 @@ float	significandf(float);
  * long double versions of ISO/POSIX math functions
  */
 #if __ISO_C_VISIBLE >= 1999
-#if _DECLARE_C99_LDBL_MATH
 long double	acoshl(long double);
-#endif
 long double	acosl(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	asinhl(long double);
-#endif
 long double	asinl(long double);
 long double	atan2l(long double, long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	atanhl(long double);
-#endif
 long double	atanl(long double);
 long double	cbrtl(long double);
 long double	ceill(long double);
 long double	copysignl(long double, long double) __pure2;
-#if _DECLARE_C99_LDBL_MATH
 long double	coshl(long double);
-#endif
 long double	cosl(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	erfcl(long double);
 long double	erfl(long double);
-#endif
 long double	exp2l(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	expl(long double);
 long double	expm1l(long double);
-#endif
 long double	fabsl(long double) __pure2;
 long double	fdiml(long double, long double);
 long double	floorl(long double);
@@ -555,20 +543,14 @@ long double	frexpl(long double value, int *); /* fundamentally !__pure2 */
 long double	hypotl(long double, long double);
 int		ilogbl(long double) __pure2;
 long double	ldexpl(long double, int);
-#if _DECLARE_C99_LDBL_MATH
 long double	lgammal(long double);
-#endif
 long long	llrintl(long double);
 long long	llroundl(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	log10l(long double);
 long double	log1pl(long double);
 long double	log2l(long double);
-#endif
 long double	logbl(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	logl(long double);
-#endif
 long		lrintl(long double);
 long		lroundl(long double);
 long double	modfl(long double, long double *); /* fundamentally !__pure2 */
@@ -578,31 +560,27 @@ long double	nextafterl(long double, long double);
 double		nexttoward(double, long double);
 float		nexttowardf(float, long double);
 long double	nexttowardl(long double, long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	powl(long double, long double);
-#endif
 long double	remainderl(long double, long double);
 long double	remquol(long double, long double, int *);
 long double	rintl(long double);
 long double	roundl(long double);
 long double	scalblnl(long double, long);
 long double	scalbnl(long double, int);
-#if _DECLARE_C99_LDBL_MATH
 long double	sinhl(long double);
-#endif
 long double	sinl(long double);
 long double	sqrtl(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	tanhl(long double);
-#endif
 long double	tanl(long double);
-#if _DECLARE_C99_LDBL_MATH
 long double	tgammal(long double);
-#endif
 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)