From 4d5dbeceac8f2a65c8d4174dacf0c19b4165e9f1 Mon Sep 17 00:00:00 2001
From: Andreas Jaeger <aj@suse.de>
Date: Wed, 21 Aug 2002 08:00:55 +0000
Subject: Merge with mainline.

---
 sysdeps/ieee754/ldbl-96/e_lgammal_r.c | 42 ++++++++++++++++++++++-------------
 1 file changed, 26 insertions(+), 16 deletions(-)

diff --git a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
index f39ef355a1..24c5a4a450 100644
--- a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
+++ b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
@@ -18,9 +18,9 @@
  *
  * Method:
  *   1. Argument Reduction for 0 < x <= 8
- * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
- * 	reduce x to a number in [1.5,2.5] by
- * 		lgamma(1+s) = log(s) + lgamma(s)
+ *	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ *	reduce x to a number in [1.5,2.5] by
+ *		lgamma(1+s) = log(s) + lgamma(s)
  *	for example,
  *		lgamma(7.3) = log(6.3) + lgamma(6.3)
  *			    = log(6.3*5.3) + lgamma(5.3)
@@ -50,13 +50,13 @@
  *	Let z = 1/x, then we approximation
  *		f(z) = lgamma(x) - (x-0.5)(log(x)-1)
  *	by
- *	  			    3       5             11
+ *				    3       5             11
  *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
  *
  *   4. For negative x, since (G is gamma function)
  *		-x*G(-x)*G(x) = pi/sin(pi*x),
- * 	we have
- * 		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
+ *	we have
+ *		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
  *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
  *	Hence, for x<0, signgam = sign(sin(pi*x)) and
  *		lgamma(x) = log(|Gamma(x)|)
@@ -69,7 +69,7 @@
  *		lgamma(1)=lgamma(2)=0
  *		lgamma(x) ~ -log(x) for tiny x
  *		lgamma(0) = lgamma(inf) = inf
- *	 	lgamma(-integer) = +-inf
+ *		lgamma(-integer) = +-inf
  *
  */
 
@@ -84,6 +84,7 @@ static long double
   half = 0.5L,
   one = 1.0L,
   pi = 3.14159265358979323846264L,
+  two63 = 9.223372036854775808e18L,
 
   /* lgam(1+x) = 0.5 x + x a(x)/b(x)
      -0.268402099609375 <= x <= 0
@@ -206,8 +207,7 @@ sin_pi (x)
 
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
-
-  i1 = (ix << 16) | (i0 >> 16);
+  ix = (ix << 16) | (i0 >> 16);
   if (ix < 0x3ffd8000) /* 0.25 */
     return __sinl (pi * x);
   y = -x;			/* x is assume negative */
@@ -219,13 +219,25 @@ sin_pi (x)
   z = __floorl (y);
   if (z != y)
     {				/* inexact anyway */
-      y *= half;
-      y = 2.0 * (y - __floorl (y));	/* y = |x| mod 2.0 */
-      n = (int) (y * 4.0);
+      y  *= 0.5;
+      y = 2.0*(y - __floorl(y));		/* y = |x| mod 2.0 */
+      n = (int) (y*4.0);
     }
   else
     {
-      return (zero + zero);
+      if (ix >= 0x403f8000)  /* 2^64 */
+	{
+	  y = zero; n = 0;                 /* y must be even */
+	}
+      else
+	{
+	if (ix < 0x403e8000)  /* 2^63 */
+	  z = y + two63;	/* exact */
+	GET_LDOUBLE_WORDS (se, i0, i1, z);
+	n = i1 & 1;
+	y  = n;
+	n <<= 2;
+      }
     }
 
   switch (n)
@@ -267,6 +279,7 @@ __ieee754_lgammal_r (x, signgamp)
   int i, ix;
   u_int32_t se, i0, i1;
 
+  *signgamp = 1;
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
 
@@ -276,7 +289,6 @@ __ieee754_lgammal_r (x, signgamp)
   ix = (ix << 16) | (i0 >> 16);
 
   /* purge off +-inf, NaN, +-0, and negative arguments */
-  *signgamp = 1;
   if (ix >= 0x7fff0000)
     return x * x;
 
@@ -292,8 +304,6 @@ __ieee754_lgammal_r (x, signgamp)
     }
   if (se & 0x8000)
     {
-      if (x == __floorl(x))
-	return x / zero;
       t = sin_pi (x);
       if (t == zero)
 	return one / fabsl (t);	/* -integer */
-- 
cgit 1.4.1