From b6cd89849e57f0caebcd4e45408b62d2004e6e4e Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 8 Jan 2015 09:49:31 +0100
Subject: [PATCH 1/3] Don't let tgammal() modify signgam. Only lgamma*() should
 modify it.

Letting tgammal() modify signgam has two disadvantages:

- It breaks valid code that assumes that the value of signgam is not
  clobbered by calls to tgammal().
- It makes this function depend on the presence of signgam. signgam is
  an X/Open System Interface. It is not part of the C standard.
---
 ld128/e_tgammal.c | 12 +++---------
 ld80/e_tgammal.c  | 19 ++++++-------------
 2 files changed, 9 insertions(+), 22 deletions(-)

diff --git a/ld128/e_tgammal.c b/ld128/e_tgammal.c
index 2fa8002..fb36209 100644
--- a/ld128/e_tgammal.c
+++ b/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) {
-		signgam = 0;
+	if (((i0&0x7fffffffffffffffLL)|i1) == 0)
 		return (1.0/x);
-	}
 
-	if (i0<0 && (u_int64_t)i0<0xffff000000000000ULL && rintl(x)==x) {
-		signgam = 0;
+	if (i0<0 && (u_int64_t)i0<0xffff000000000000ULL && rintl(x)==x)
 		return (x-x)/(x-x);
-	}
 
-	if (i0==0xffff000000000000ULL && i1==0) {
-		signgam = 0;
+	if (i0==0xffff000000000000ULL && i1==0)
 		return (x-x);
-	}
 
 	return expl(lgammal(x));
 }
diff --git a/ld80/e_tgammal.c b/ld80/e_tgammal.c
index 571a0a9..8dfcd65 100644
--- a/ld80/e_tgammal.c
+++ b/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
- * correctly signed, and the sign (+1 or -1) is also
- * returned in a global (extern) variable named signgam.
- * This variable is also filled in by the logarithmic gamma
+ * Returns gamma function of the argument.  The result is correctly
+ * signed.  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 ));

From 24cec16fccb2a42b5c598eb6893b8c6ca1b86a3f Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 8 Jan 2015 11:07:03 +0100
Subject: [PATCH 2/3] Remove references to _DECLARE_C99_LDBL_MATH.

When building openlibm with Clang, I seem to get a lot of warnings in
ld80/ related to some prototypes for long double functions that are
missing. This seems to be because we don't define _DECLARE_C99_LDBL_MATH
anywhere.

It seems that this definition only existed on FreeBSD, as certain C99
math functions were not present yet. The prototypes were simply there as
placeholders. This flag has been removed upstream (FreeBSD SVN r236148).
---
 src/openlibm.h | 26 --------------------------
 1 file changed, 26 deletions(-)

diff --git a/src/openlibm.h b/src/openlibm.h
index 5635612..315a227 100644
--- a/src/openlibm.h
+++ b/src/openlibm.h
@@ -499,35 +499,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);
@@ -539,20 +527,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 */
@@ -562,27 +544,19 @@ 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 */

From 55ac46280872f93e8c5dadede43ca33aac6dd680 Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 8 Jan 2015 11:23:28 +0100
Subject: [PATCH 3/3] 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.c => e_lgammal_r.c} | 25 ++++++++++++-------------
 ld80/Make.files                      |  2 +-
 ld80/{e_lgammal.c => e_lgammal_r.c}  | 13 ++++++-------
 src/Make.files                       |  2 +-
 src/e_lgammal.c                      | 12 ++++++++++++
 src/openlibm.h                       |  6 +++++-
 6 files changed, 37 insertions(+), 23 deletions(-)
 rename ld128/{e_lgammal.c => e_lgammal_r.c} (98%)
 rename ld80/{e_lgammal.c => e_lgammal_r.c} (98%)
 create mode 100644 src/e_lgammal.c

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