Repository navigation
fix: avoid premature overflow in math/base/special/expm1 - #16239
Open
Thanatos9404 wants to merge 1 commit into
Open
Thanatos9404 wants to merge 1 commit into
Thanatos9404 wants to merge 1 commit into
Conversation
Previously, `expm1` returned `+infinity` when `x` was exactly equal to the overflow threshold `709.782712893384`, which is the double-precision floating-point number nearest to, and just below, `ln(MAX_VALUE)`. The exact result there is about `1.79769313486227322e+308`, which is finite and about 213 ULPs below `1.7976931348623157e+308`, and `exp` already returns `1.7976931348622732e+308` for the same input. As a result, `expm1` overflowed for an input at which the smaller function `exp` did not. The FDLIBM and FreeBSD implementations only return `+infinity` when `x` is strictly greater than the threshold, and the implementation notes in this package already state that `expm1(x)` overflows for `x > 7.09782712893383973096e+02`. This commit updates both the JavaScript and C implementations to use a strict comparison. For `x` equal to the threshold, argument reduction yields `k = 1024`, which is handled by the existing branch that adds `k` to the exponent of `y` directly, so no intermediate value overflows, and the function now returns the correctly rounded value `1.7976931348622732e+308`. Closes: stdlib-js#12468 Ref: http://www.netlib.org/fdlibm/s_expm1.c
stdlib-bot
added
Math
stdlib-bot
removed
the
First-time Contributor
This branch has not been deployed
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
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.
Resolves #12468.
Description
This pull request:
@stdlib/math/base/special/expm1fromy >= OVERFLOW_THRESHOLDtoy > OVERFLOW_THRESHOLDin both the JavaScript and C implementations, matching FDLIBM and FreeBSD (if(x > o_threshold) return huge*huge;) and the existing implementation notes ("if x > 7.09782712893383973096e+02, then expm1(x) will overflow").expm1( 709.782712893384 )now returns1.7976931348622732e+308(the same asexp( 709.782712893384 )andMath.expm1( 709.782712893384 )) instead ofInfinity.Related Issues
This pull request has the following related issues:
math/base/special/expm1overflows prematurely #12468math/base/special/expm1reloverflows prematurely #12436Questions
No.
Other
@kgryte I noticed this is assigned to you, so feel free to close this if you already have a fix in progress.
For
xequal to the threshold, argument reduction yieldsk = 1024, which goes through the existing branch that addskto the exponent ofydirectly (rather than multiplying by2^k), so no intermediate value overflows.Expected values in the new tests come from
mpmathat 400 bits of precision, rounded to nearest:xexpm1(x)709.78271289338391.79769313486206884...e+3081.7976931348620688e+308709.7827128933841.79769313486227322...e+3081.7976931348622732e+308709.78271289338411.79769313486247759...e+308Infinity(exceeds1.7976931348623157e+308)Both implementations were also compared against
mpmathon 7,004 inputs (the threshold, the 2,000 doubles below it, the three above it, and random inputs in[600, 710]and in thek = 1024band): the maximum error is 1 ULP and every overflow matches the reference. The C implementation (built withgcc -std=c99 -Wall -Wextra -pedanticat-O0and-O3) returns the same bits as the JavaScript implementation on all of them.expm1rel(#12436) is unaffected: forxequal toHUGE_VALUEit already takes theexpm1( x/2 )branch, which returns the same value.Checklist
AI Assistance
If you answered "yes" above, how did you use AI assistance?
Disclosure
Written primarily with Claude Code; the new values were checked against a high-precision reference.
@stdlib-js/reviewers