Fix powl() returning NaN instead of +0 on extreme underflow (#334)#351
Merged
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #351 +/- ##
==========================================
- Coverage 72.12% 72.11% -0.02%
==========================================
Files 233 233
Lines 6135 6139 +4
Branches 1607 1609 +2
==========================================
+ Hits 4425 4427 +2
- Misses 1417 1419 +2
Partials 293 293 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
ViralBShah
force-pushed
the
fix-334-powl-overflow
branch
2 times, most recently
from
9b11950 to
d437edc
Compare
ViralBShah
force-pushed
the
fix-334-powl-overflow
branch
from
d437edc to
62d3ef1
Compare
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>
ViralBShah
force-pushed
the
fix-334-powl-overflow
branch
from
62d3ef1 to
0654492
Compare
|
good work! I'll check that with the next release |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #334.
powl(0x1.98p-3072L, 0xa.ab43b6dba9a6383p+16364L)returnedNaN; glibc returns+0(underflow).Root cause
In
ld80/e_powl.c,powlformsG = y*log2(x)in 1/NXT units (hereG ≈ -2^16379). The reduction helperreducl(G)computesldexpl(G, LNXT) = G*32 ≈ -2^16384, which overflows to -Inf. That cascades (Gb = G - (-Inf) = +Inf;w = ldexpl(-Inf + +Inf, …) = NaN), so the over/underflow testsw > MEXP/w < MNEXPare both false and the routine falls through returning NaN. The true result underflows to +0. The bug is shared verbatim with current FreeBSD/OpenBSD msun, so there is no upstream patch to import.Fix
Short-circuit over/underflow on
G(which already holds the final exponent to within O(1)) beforereducl(G)can overflow, with a generous2*MEXP/NXTmargin so genuine borderline values still reach the original precise test:Verification
0x0p+0, matching glibc.1^y,x^0, and four over/underflow boundary cases at y=16383.5/16384.5/−16445/−16500) all match glibc, including exact boundary classification.make testpasses (double + float + exception-flag suites).test/regression/test-334.casserts the result is+0and not NaN; verified it fails before this fix and passes after (skips wherelong double == double).Stacked on the regression-harness PR #350 (provides
test/regression/).