Set errno in tgamma() per C99/POSIX (#325)

tgamma() previously left errno unset on overflow, at the pole x == 0, and
for negative-integer arguments, violating C99/POSIX and diverging from
glibc. Wrap the numeric implementation (now static __tgamma) in a thin
public tgamma() that classifies the (argument, result) pair and sets:

  - ERANGE when the result is +-Inf for a finite argument
    (overflow, and the pole at x == 0);
  - EDOM when the result is NaN for a finite argument
    (negative-integer arguments).

In-range results and NaN/Inf inputs leave errno untouched. The internal
recursive calls in neg_gam() use __tgamma directly so intermediate
computations cannot perturb errno. Values returned are unchanged.

Behaviour now matches system glibc for x in {0, -0.0, -2, 200, 5,
subnormal}. Adds test/regression/test-325.c covering the ERANGE, EDOM,
and errno-untouched cases; it fails pre-fix and passes post-fix.

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

Author: ViralBShah

bsdsrc/b_tgamma.c

@@ -38,6 +38,7 @@
  * acknowledged.
  */
 
+#include <errno.h>
 #include <float.h>
 #include <openlibm_math.h>
 
@@ -77,6 +78,7 @@
  *	Maximum observed error < 4ulp in 1,000,000 trials.
  */
 
+static double __tgamma(double);
 static double neg_gam(double);
 static double small_gam(double);
 static double smaller_gam(double);
@@ -123,9 +125,35 @@ static struct Double ratfun_gam(double, double);
 
 static const double zero = 0., one = 1.0, tiny = 1e-300;
 
+/*
+ * Public entry point.  The numeric work is done by __tgamma(); this wrapper
+ * inspects the (argument, result) pair and sets errno as required by C99/POSIX,
+ * matching the classification used by glibc:
+ *	- pole at x == 0 (result +-Inf):	ERANGE
+ *	- overflow (finite x, result +-Inf):	ERANGE
+ *	- negative integer (result NaN):	EDOM
+ * Underflow to +-0 is left alone (consistent with glibc, which does not set
+ * errno there), and normal in-range results never touch errno.
+ */
 OLM_DLLEXPORT double
 tgamma(x)
 	double x;
+{
+	double r = __tgamma(x);
+
+	if (isinf(r) && isfinite(x)) {
+		/* Pole error at x == 0, or overflow for finite x. */
+		errno = ERANGE;
+	} else if (isnan(r) && isfinite(x)) {
+		/* Domain error: negative integer argument. */
+		errno = EDOM;
+	}
+	return (r);
+}
+
+static double
+__tgamma(x)
+	double x;
 {
 	struct Double u;
 
@@ -307,9 +335,9 @@ neg_gam(x)
 	}
 	y = one-x;
 	if (one-y == x)
-		y = tgamma(y);
+		y = __tgamma(y);	/* internal: do not perturb errno */
 	else		/* 1-x is inexact */
-		y = -x*tgamma(-x);
+		y = -x*__tgamma(-x);
 	if (sgn < 0) y = -y;
 	return (M_PI / (y*z));
 }

test/regression/test-325.c

@@ -0,0 +1,55 @@
+/*
+ * Regression test for issue #325: tgamma() does not set errno.
+ * https://github.com/JuliaMath/openlibm/issues/325
+ *
+ * C99/POSIX require tgamma to report:
+ *	- a range error (errno == ERANGE) on overflow and at the pole x == 0,
+ *	  where the result is +-HUGE_VAL;
+ *	- a domain error (errno == EDOM) for negative-integer arguments, where
+ *	  the result is NaN.
+ * In-range results must leave errno untouched.
+ *
+ * This pins that behaviour so the missing-errno regression cannot return.
+ */
+#include <errno.h>
+#include <math.h>
+
+#include "regress-util.h"
+
+int
+main(void)
+{
+	double r;
+
+	/* Pole error at x == 0: +Inf, ERANGE. */
+	errno = 0;
+	r = tgamma(0.0);
+	CHECK(isinf(r) && r > 0.0);
+	CHECK(errno == ERANGE);
+
+	/* Pole error at x == -0.0: -Inf, ERANGE. */
+	errno = 0;
+	r = tgamma(-0.0);
+	CHECK(isinf(r) && r < 0.0);
+	CHECK(errno == ERANGE);
+
+	/* Domain error for a negative integer: NaN, EDOM. */
+	errno = 0;
+	r = tgamma(-2.0);
+	CHECK(isnan(r));
+	CHECK(errno == EDOM);
+
+	/* Overflow for large x: +Inf, ERANGE. */
+	errno = 0;
+	r = tgamma(200.0);
+	CHECK(isinf(r) && r > 0.0);
+	CHECK(errno == ERANGE);
+
+	/* Normal in-range value: exact result, errno untouched. */
+	errno = 0;
+	r = tgamma(5.0);
+	CHECK(r == 24.0);
+	CHECK(errno == 0);
+
+	return 0;
+}