Fix expm1l returning NaN for large finite arguments (#306)

For arguments just below the overflow threshold, expm1l reconstructed the
result as 2^k * qx + (2^k - 1) with k == LDBL_MAX_EXP (16384).  ldexpl(1, k)
then overflowed to +Inf, and the subsequent (+Inf) + (-Inf) (qx is negative
when the reduced remainder is negative) produced NaN, even though the true
result is finite and near LDBL_MAX.

When k is large enough that 2^k overflows, instead form exp(remainder ln 2)
= qx + 1 (which lies near 1) and scale its exponent directly via scalbnl, so
no intermediate value overflows.  This is algebraically identical to the
existing formula; the normal path is unchanged.

Result now matches glibc bit-for-bit on x86-64 (80-bit long double), and
slightly larger arguments still overflow to +Inf as before.  Adds
regression test test-306.

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

Author: ViralBShah

ld80/s_expm1l.c

@@ -57,6 +57,8 @@
  *
  */
 
+#include <float.h>
+
 #include <openlibm_math.h>
 
 static const long double MAXLOGL = 1.1356523406294143949492E4L;
@@ -131,8 +133,17 @@ qx = x + (0.5 * xx + xx * px / qx);
 
 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
    We have qx = exp(remainder ln 2) - 1, so
-   exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.  */
-px = ldexpl(1.0L, k);
-x = px * qx + (px - 1.0);
+   exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.
+
+   For large k, 2^k overflows even though the final result is finite
+   (the argument is just below the overflow threshold).  In that regime
+   form exp(remainder ln 2) = qx + 1 first (which lies near 1) and scale
+   the exponent directly, so no intermediate value overflows.  */
+if (k >= LDBL_MAX_EXP) {
+  x = scalbnl(qx + 1.0L, k) - 1.0L;
+} else {
+  px = ldexpl(1.0L, k);
+  x = px * qx + (px - 1.0);
+}
 return x;
 }

test/regression/test-306.c

@@ -0,0 +1,42 @@
+/*
+ * Regression test for issue #306: expm1l overflows internally and returns NaN
+ * for large arguments that are still below the overflow threshold.
+ * https://github.com/JuliaMath/openlibm/issues/306
+ *
+ * For x just below MAXLOG the true result is finite (near LDBL_MAX), but the
+ * old reconstruction formed 2^k (with k == LDBL_MAX_EXP) which overflowed to
+ * +Inf; the subsequent (+Inf) + (-Inf) yielded NaN.  The fix scales the
+ * exponent of exp(remainder) directly in that regime.
+ *
+ * The expected hex values below are bit-for-bit what the system glibc expm1l
+ * returns on x86-64 (LDBL_MANT_DIG == 64).  Only meaningful for the 80-bit
+ * long double format, so skip elsewhere.
+ */
+#include <float.h>
+#include <math.h>
+
+#include "regress-util.h"
+
+int
+main(void)
+{
+#if LDBL_MANT_DIG < 64
+	return REGRESS_SKIP;
+#else
+	long double x = 0x2.c5c85fdf170c604cp+12l;	/* ~11346.6, just below threshold */
+	long double y = expm1l(x);
+
+	/* Must be a large finite value, not NaN/Inf. */
+	CHECK(!isnanl(y));
+	CHECK(isfinite(y));
+	CHECK(y > 0.0L);
+
+	/* Exact value matching glibc. */
+	CHECK(y == 0xf.ffffcfce79e56d5p+16380l);
+
+	/* A slightly larger argument must still overflow to +Inf. */
+	CHECK(isinf(expm1l(11357.0L)) == 1);
+
+	return 0;
+#endif
+}