Fix powl() returning NaN on extreme underflow (#334)

In the ld80 (x86 80-bit) powl(), when |y*log2(x)| is enormous the
extended-precision reduction reducl(G) computes ldexpl(G, LNXT), which
overflows to +-Inf.  The following Inf-Inf / Inf+(-Inf) arithmetic
(Gb = G - Ga, Ha = reducl(Fb + Gb), Ga + Ha) yields NaN, so the
w > MEXP / w < MNEXP over/underflow tests both fail and powl() returns
NaN instead of Inf (overflow) or +0 (underflow).

Reproducer:
  powl(0x1.98p-3072l, 0xa.ab43b6dba9a6383p+16364l)
returned NaN; glibc correctly returns +0 (underflow).

This is a latent bug shared with upstream FreeBSD msun / OpenBSD; there
is no upstream fix to import.  G already holds y*log2(x) in 1/NXT units,
so the final exponent is ldexpl(G, LNXT) to within O(1).  Test G for
over/underflow before reducl(G) can overflow.  A generous 2*MEXP/NXT
margin keeps the early guard far from the genuine over/underflow
boundary, so borderline values still fall through to the original
precise test on w = ldexpl(Ga+Ha, LNXT); only the pathological huge-|G|
case is short-circuited.

Verified the reproducer now returns +0 and a battery of normal,
overflow, underflow, negative-base, non-integer-exponent and boundary
cases all match glibc.  Full suite (test-double, test-float,
regression) passes.

Adds test/regression/test-334.c (skips with status 77 when
LDBL_MANT_DIG < 64, since the ld80 path is unused there); it fails on
the pre-fix code and passes after.

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

Author: ViralBShah

ld80/e_powl.c

@@ -404,6 +404,27 @@ Fa = reducl(F);
 Fb = F - Fa;
 
 G = Fa + w * ya;
+
+/*
+ * Guard against intermediate overflow in reducl(G).  When |y*log2(x)|
+ * is enormous, reducl(G) computes ldexpl(G, LNXT) which overflows to
+ * +-Inf; the subsequent Inf-Inf / Inf+(-Inf) arithmetic then yields
+ * NaN, which defeats the w > MEXP / w < MNEXP tests below and makes
+ * powl() return NaN instead of Inf (overflow) or +0 (underflow).
+ * G holds y*log2(x) in 1/NXT units, so the final exponent is
+ * ldexpl(G, LNXT) to within O(1).  Genuine over/underflow occurs at
+ * |G| ~ MEXP/NXT ~ 16448, vastly below the |G| ~ LDBL_MAX/NXT at
+ * which reducl(G) overflows.  A generous margin (2*MEXP/NXT) keeps
+ * this guard well clear of the real boundary, so borderline cases
+ * still fall through to the precise test on w = ldexpl(Ga+Ha, LNXT)
+ * below, while the pathological huge-|G| case is caught here before
+ * reducl(G) can overflow.
+ */
+if( G > (2.0L * MEXP / NXT) )
+	return (huge * huge);		/* overflow */
+if( G < (2.0L * MNEXP / NXT) )
+	return (twom10000 * twom10000);	/* underflow */
+
 Ga = reducl(G);
 Gb = G - Ga;
 

test/regression/test-334.c

@@ -0,0 +1,47 @@
+/*
+ * Regression test for issue #334: powl() returns NaN instead of +0 on
+ * extreme underflow.
+ * https://github.com/JuliaMath/openlibm/issues/334
+ *
+ * For very large |y * log2(x)|, the extended-precision reduction step in
+ * ld80/e_powl.c (reducl(G), which computes ldexpl(G, LNXT)) overflowed to
+ * +-Inf.  The following Inf-Inf / Inf+(-Inf) arithmetic produced NaN, which
+ * defeated the w > MEXP / w < MNEXP over/underflow tests and made powl()
+ * return NaN.  The mathematically correct result here is +0 (underflow),
+ * which is what glibc returns.
+ *
+ * This manifests only when long double has more precision than double (the
+ * x86 80-bit / ld80 path); otherwise powl == pow and the bug cannot occur.
+ */
+#include <float.h>
+#include <math.h>
+
+#include "regress-util.h"
+
+int
+main(void)
+{
+	/*
+	 * ld80-specific (the fix is in ld80/e_powl.c).  openlibm provides long
+	 * double routines only on x86 and (linux) aarch64, so restrict to the
+	 * 80-bit format; other arches compile to a bare skip with no powl
+	 * reference (powl may be absent, which would break the link).
+	 */
+#if LDBL_MANT_DIG != 64
+	return REGRESS_SKIP;
+#else
+	long double x = 0x1.98p-3072l;
+	long double y = 0xa.ab43b6dba9a6383p+16364l;
+	long double z = powl(x, y);
+
+	/* Must underflow to +0, not return NaN. */
+	CHECK(!isnan(z));
+	CHECK(z == 0.0L);
+	CHECK(!signbit(z));   /* result is +0, not -0 */
+
+	/* Exercise the symmetric early-overflow guard: a result far beyond the
+	   range (here 10^1e6) must return +Inf, matching glibc. */
+	CHECK(isinf(powl(10.0L, 1000000.0L)) == 1);
+	return 0;
+#endif
+}