Add lgammal_r().

We already provide lgammaf_r() and lgamma_r(). It's not hard to also add
lgammal_r(), for consistency.

I am currently working on porting openlibm to an environment where
global state, and thus signgam, is not available. By adding lgammal_r(),
I can trivially disable support for signgam by just patching up
src/e_lgamma{f,,l}.c. That way there is no need to patch up the actual
algorithms.

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;

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

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

@@ -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)