From: Szabolcs Nagy <nsz@port70.net>
Date: Sun, 16 Dec 2012 19:22:17 +0000 (+0100)
Subject: math: tgammal.c fixes
X-Git-Tag: v0.9.9~50
X-Git-Url: https://git.librecmc.org/?a=commitdiff_plain;h=d8a7619e371ff0f226200f6316abb46dd1192f3d;p=oweals%2Fmusl.git

math: tgammal.c fixes

this is not a full rewrite just fixes to the special case logic:
+-0 and non-integer x<INT_MIN inputs incorrectly raised invalid
exception and for +-0 the return value was wrong

so integer test and odd/even test for negative inputs are changed
and a useless overflow test was removed
---

diff --git a/src/math/tgammal.c b/src/math/tgammal.c
index 86782a96..5c1a02a6 100644
--- a/src/math/tgammal.c
+++ b/src/math/tgammal.c
@@ -205,42 +205,36 @@ static long double stirf(long double x)
 long double tgammal(long double x)
 {
 	long double p, q, z;
-	int i, sign;
 
-	if (isnan(x))
-		return NAN;
-	if (x == INFINITY)
-		return INFINITY;
-	if (x == -INFINITY)
-		return x - x;
+	if (!isfinite(x))
+		return x + INFINITY;
+
 	q = fabsl(x);
 	if (q > 13.0) {
-		sign = 1;
-		if (q > MAXGAML)
-			goto goverf;
 		if (x < 0.0) {
 			p = floorl(q);
-			if (p == q)
-				return (x - x) / (x - x);
-			i = p;
-			if ((i & 1) == 0)
-				sign = -1;
 			z = q - p;
-			if (z > 0.5L) {
-				p += 1.0;
-				z = q - p;
-			}
-			z = q * sinl(PIL * z);
-			z = fabsl(z) * stirf(q);
-			if (z <= PIL/LDBL_MAX) {
-goverf:
-				return sign * INFINITY;
+			if (z == 0)
+				return 0 / z;
+			if (q > MAXGAML) {
+				z = 0;
+			} else {
+				if (z > 0.5) {
+					p += 1.0;
+					z = q - p;
+				}
+				z = q * sinl(PIL * z);
+				z = fabsl(z) * stirf(q);
+				z = PIL/z;
 			}
-			z = PIL/z;
+			if (0.5 * p == floorl(q * 0.5))
+				z = -z;
+		} else if (x > MAXGAML) {
+			z = x * 0x1p16383L;
 		} else {
 			z = stirf(x);
 		}
-		return sign * z;
+		return z;
 	}
 
 	z = 1.0;
@@ -268,8 +262,9 @@ goverf:
 	return z;
 
 small:
-	if (x == 0.0)
-		return (x - x) / (x - x);
+	/* z==1 if x was originally +-0 */
+	if (x == 0 && z != 1)
+		return x / x;
 	if (x < 0.0) {
 		x = -x;
 		q = z / (x * __polevll(x, SN, 8));