ld80: fix spurious tgammal overflow for large negative arguments (#223)

tgammal(-0xd.b6e8f5c28f5c29p+7L) (approx -1755.455) returned +Inf, but
the true value is the tiny subnormal 0x0.01dbd551da54538p-16385L. The
ld80 implementation made two mistakes near the gamma-overflow threshold
MAXGAML:

  1. The `q > MAXGAML` early-out was applied to *both* signs of x.
     MAXGAML is the threshold at which the *positive* gamma overflows;
     for large negative x, |Gamma(x)| is instead tiny, so this wrongly
     forced negative arguments down the +Inf overflow path.

  2. Just below MAXGAML the reflection code computed
     `fabsl(z) * stirf(q)` first. For q near 1755, tgamma(q) is enormous
     and the product overflowed to Inf, after which `PIL/Inf` collapsed
     to 0 instead of the correct tiny subnormal.

Fix: restrict the MAXGAML overflow check to positive x (the only side
that overflows), and evaluate the reflection formula as
`(PIL / stirf(q)) / z` so the huge tgamma(q) divides pi down to a tiny
value before the multiply, avoiding the intermediate overflow. When
stirf(q) itself overflows to Inf the quotient correctly underflows
toward zero.

With this change tgammal(-0xd.b6e8f5c28f5c29p+7L) returns
0x0.01dbd551da54538p-16385L, matching glibc bit-for-bit, and nearby
large negative non-integers that previously flushed to 0 now yield the
correct subnormals. Positive overflow behavior is unchanged.

Adds test/regression/test-223.c (skipped when LDBL_MANT_DIG < 64).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: ViralBShah

ld80/e_tgammal.c

@@ -233,8 +233,6 @@ q = fabsl(x);
 if( q > 13.0L )
 	{
 	int sign = 1;
-	if( q > MAXGAML )
-		goto goverf;
 	if( x < 0.0L )
 		{
 		p = floorl(q);
@@ -250,16 +248,22 @@ if( q > 13.0L )
 			z = q - p;
 			}
 		z = q * sinl( PIL * z );
-		z = fabsl(z) * stirf(q);
-		if( z <= PIL/LDBL_MAX )
-			{
-goverf:
-			return( sign * INFINITY);
-			}
-		z = PIL/z;
+		z = fabsl(z);
+		/*
+		 * Reflection formula: |tgamma(x)| = pi / (z * tgamma(q)).
+		 * For large negative x, tgamma(q) is enormous, so |tgamma(x)|
+		 * is tiny.  Divide pi by tgamma(q) *before* dividing by z to
+		 * avoid an intermediate overflow of z * tgamma(q) that would
+		 * spuriously yield Inf (issue #223).  When tgamma(q) is so
+		 * large that stirf(q) overflows to Inf, the quotient correctly
+		 * underflows toward zero.
+		 */
+		z = (PIL / stirf(q)) / z;
 		}
 	else
 		{
+		if( q > MAXGAML )
+			return( INFINITY );
 		z = stirf(x);
 		}
 	return( sign * z );

test/regression/test-223.c

@@ -0,0 +1,59 @@
+/*
+ * Regression test for issue #223: tgammal spuriously overflows for
+ * large-ish negative non-integer arguments.
+ * https://github.com/JuliaMath/openlibm/issues/223
+ *
+ * For a large negative non-integer x, |tgamma(x)| is tiny (the reflection
+ * formula divides pi by the enormous tgamma(1-x)).  The ld80 implementation
+ * applied its positive-overflow threshold (MAXGAML) to negative arguments
+ * too and, just below that threshold, let z * tgamma(q) overflow to Inf
+ * before the final division -- both paths wrongly produced +Inf (or a
+ * flush-to-zero) instead of the correct tiny subnormal.
+ *
+ * glibc on x86-64 (the correctness oracle) returns the correctly-rounded
+ * subnormal 0x0.01dbd551da54538p-16385L for the argument below; pin that.
+ */
+#include <float.h>
+
+#include "regress-util.h"
+
+/*
+ * ld80-specific (the fix is in ld80/e_tgammal.c, and the exact hex below holds
+ * only for the 80-bit format).  openlibm provides long double routines only on
+ * x86 and (linux) aarch64, so every other format must compile to a bare skip
+ * with no reference to tgammal (it may be absent, which would break the link).
+ */
+#if LDBL_MANT_DIG != 64
+
+int
+main(void)
+{
+	return REGRESS_SKIP;
+}
+
+#else
+
+#include <math.h>
+
+int
+main(void)
+{
+	long double x = -0xd.b6e8f5c28f5c29p+7L;	/* approx -1755.455 */
+	long double y = tgammal(x);
+
+	/* Must be a finite, positive, tiny value -- not Inf and not zero. */
+	CHECK(isfinite(y));
+	CHECK(y > 0.0L);
+	CHECK(y < 1.0L);
+
+	/* Correctly-rounded result observed from glibc tgammal on x86-64. */
+	CHECK(y == 0x0.01dbd551da54538p-16385L);
+
+	/* The positive-overflow path (moved out of the shared code) must still
+	   return +Inf for large positive x. */
+	CHECK(isinf(tgammal(2000.0L)) == 1);
+
+	return 0;
+}
+
+#endif