fix: avoid overflow and underflow for scales far from 1 in `stats/base/dists/rayleigh/logpdf` by Abhist17 · Pull Request #15923 · stdlib-js/stdlib

Abhist17 · 2026-10-05T19:19:15Z

Description

What is the purpose of this pull request?

This pull request:

avoids overflow and underflow for scale parameters far from 1 in stats/base/dists/rayleigh/logpdf (JS, factory and C).

For sigma outside about [1.5e-154, 1.3e154], sigma^2 overflows or underflows (and x^2 does for large x), so the log-PDF comes out NaN or a spurious -Infinity where it is finite:

var logpdf = require( '@stdlib/stats/base/dists/rayleigh/logpdf' );

logpdf( 1.0e200, 1.0e200 );
// develop: NaN
// this PR: -461.01701859880916   (= ln(1e-200) - 1/2)

logpdf( 1.0e-200, 1.0e-200 );
// develop: NaN
// this PR: 460.01701859880916

logpdf( 1.0e150, 1.0e200 );
// develop: NaN
// this PR: -575.6462732485114

This is the same class of problem @0PrashantYadav0 flagged in rayleigh/cdf (#15772) and rayleigh/logcdf (#15798), in the sibling function.

Why it's a hybrid rather than a plain switch to x/sigma. I first rewrote it entirely in terms of the ratio. That fixed the extremes, but it regressed one point in the existing small_scale fixture (x = 0.536, sigma = 0.260). There logpdf ≈ -0.057 is a small difference of two terms near 2.1. Both formulas cancel at that point; the old one just happened to round inside the tolerance. So the existing expression is kept verbatim wherever it is safe — sigma^2 normal and finite, x^2 finite, x/sigma^2 normal — and the function only falls through to the ratio form otherwise. On every one of the 3,000 points in the existing fixtures the result is bit-for-bit identical to develop.

In the ratio branch, ln( x/sigma^2 ) is taken directly while that ratio is normal, falling back to ln( r ) - ln( sigma ) only when it isn't: subtracting two large logarithms cancels when x is close to sigma^2. In factory.js, ln( sigma ) and the fast-path test are computed once per sigma.

Measured against mpmath at 80 digits on 5,876 points, with sigma log-uniform in [1e-300, 1e300] and x/sigma in [1e-3, 10] (plus 15% scaled down by 1e-100):

	NaN or wrong `±Infinity`	> 3 ULP	max ULP (finite)
`develop`	2,784	5	71
this PR	0	5	71

The five points over 3 ULP and the 71-ULP maximum are the same points with the same values in both rows: the old expression's own cancellation in the normal range, which this PR doesn't touch.

Tests. gen in runner.jl now evaluates the textbook formula in 2048-bit BigFloat, matching the pattern from #15774 and #15772. Distributions.jl's logpdf works in Float64 and overflows exactly where the bug is, so it can't produce these fixtures. Two new fixtures, both with x/sigma in [1e-3, 10]:

tiny_scale.json: sigma in [1e-200, 1e-155], where sigma^2 underflows.
huge_scale.json: sigma in [1e155, 1e200], where sigma^2 overflows.

They're restricted to those ranges on purpose. An earlier draft sampled the whole [1e-200, 1e200] and picked up two ordinary-scale points from the unchanged branch that sit at 3 and 5 ULP on develop too. test.logpdf.js, test.factory.js and test.native.js each gain a block per fixture, at the existing tolerances (3*EPS in test.logpdf.js, 2 ULP via isAlmostSameValue in the other two). On develop, all 2,000 new assertions fail in both test.logpdf.js and test.factory.js. With this PR, all three files pass 5,014 assertions. The other three fixtures weren't regenerated.

Related Issues

Does this pull request have any related issues?

None. Related to #15772 (rayleigh/cdf) and #15798 (rayleigh/logcdf).

Questions

Any questions for reviewers of this pull request?

The pdf sibling has the same overflow (pdf( 1e200, 1e200 ) is NaN). I left it out of this PR to keep one package per PR, and will follow up if this approach looks right to you.

Other

Any other information relevant to this pull request? This may include screenshots, references, and/or implementation notes.

No.

Checklist

Please ensure the following tasks are completed before submitting this pull request.

Read, understood, and followed the contributing guidelines.

AI Assistance

When authoring the changes proposed in this PR, did you use any kind of AI assistance?

Yes
No

If you answered "yes" above, how did you use AI assistance?

Code generation (e.g., when writing an implementation or fixing a bug)
Test/benchmark generation
Documentation (including examples)
Research and understanding

Disclosure

If you answered "yes" to using AI assistance, please provide a short disclosure indicating how you used AI assistance. This helps reviewers determine how much scrutiny to apply when reviewing your contribution. Example disclosures: "This PR was written primarily by Claude Code." or "I consulted ChatGPT to understand the codebase, but the proposed changes were fully authored manually by myself.".

This PR was written primarily by Claude Code. I found the bug by looking for siblings of the sigma^2 issue in the review of #15772.

@stdlib-js/reviewers

0PrashantYadav0

The three failing checks were cancelled before a runner picked them up, so CI hasn't actually run yet and needs a re-run.

0PrashantYadav0 · 2026-10-07T11:41:51Z

 	s2 = pow( sigma, 2.0 );
 	s2i = 1.0 / s2;
-	return ln( s2i * x ) - ( pow( x, 2.0 ) / ( 2.0 * s2 ) );
+	x2 = pow( x, 2.0 );
+	q = s2i * x;
+	if (
+		s2 >= FLOAT64_SMALLEST_NORMAL &&
+		s2 < PINF &&
+		x2 < PINF &&
+		q >= FLOAT64_SMALLEST_NORMAL
+	) {
+		return ln( q ) - ( x2 / ( 2.0 * s2 ) );
+	}
+	// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
+	r = x / sigma;
+	q = r / sigma;
+	if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {
+		lq = ln( q );
+	} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {
+		lq = ln( r ) - ln( sigma );
+	} else {
+		lq = ( ln( x ) - ln( sigma ) ) - ln( sigma );
+	}
+	return lq - ( 0.5 * r * r );


Suggested change

s2 = pow( sigma, 2.0 );

s2i = 1.0 / s2;

return ln( s2i * x ) - ( pow( x, 2.0 ) / ( 2.0 * s2 ) );

x2 = pow( x, 2.0 );

q = s2i * x;

if (

s2 >= FLOAT64_SMALLEST_NORMAL &&

s2 < PINF &&

x2 < PINF &&

q >= FLOAT64_SMALLEST_NORMAL

) {

return ln( q ) - ( x2 / ( 2.0 * s2 ) );

}

// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...

r = x / sigma;

q = r / sigma;

if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {

lq = ln( q );

} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {

lq = ln( r ) - ln( sigma );

} else {

lq = ( ln( x ) - ln( sigma ) ) - ln( sigma );

}

return lq - ( 0.5 * r * r );

r = x / sigma;

if ( r === PINF ) {

return NINF;

}

// Use `ln( x/sigma^2 )` while it is normal, as `ln( r ) - ln( sigma )` cancels when `x` is close to `sigma^2`:

q = r / sigma;

if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {

lq = ln( q );

} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {

lq = ln( r ) - ln( sigma );

} else {

lq = ( ln( x ) - ln( sigma ) ) - ln( sigma );

}

return lq - ( 0.5 * r * r );

2.0*s2 overflows near sigma = 1e154: logpdf( 1e154, 1e154 ) returns -354.598 (exact -355.098), and the fast path isn't more accurate (the old small_scale value it keeps is itself 43 ULP off). The r === PINF guard fixes logpdf( 1e300, 1e-10 ) returning NaN. Same x/sigma idea as #15798 (comment).

Same in lib/factory.js:88-103 and src/main.c:58-80; then drop s2, s2i, x2, normal and pow (also from manifest.json).

0PrashantYadav0 · 2026-10-07T11:41:52Z

+	z = setprecision( BigFloat, 2048 ) do
+		X = BigFloat.( x );
+		S = BigFloat.( sigma );
+		Float64.( log.( X ./ ( S .^ 2 ) ) .- ( ( X .^ 2 ) ./ ( 2 .* ( S .^ 2 ) ) ) );


small_scale.json, medium_scale.json and large_scale.json still hold the old Distributions.jl values, so they don't match this runner.jl. small_scale is up to 43 ULP off. Can you regenerate all five? With regenerated fixtures and the change above, every existing tolerance passes. Distributions can also be removed from REQUIRE.

…e/dists/rayleigh/logpdf` For `sigma` outside about [1.5e-154, 1.3e154], `sigma^2` overflows or underflows, and so does `x^2` for large `x`, so the log-PDF came out NaN or a spurious -Infinity where it is finite: `logpdf( 1e200, 1e200 )` is about -461.017 and returned NaN. Keep the existing expression wherever it is safe — it is unchanged, bit for bit, on every point of the existing fixtures — and otherwise evaluate via the ratio `x/sigma`, taking `ln( x/sigma^2 )` directly while that ratio is normal so the two logarithms don't cancel. New `tiny_scale` and `huge_scale` fixtures cover both sides; develop fails all 2,000 of their points.

Abhist17 · 2026-10-07T21:16:34Z

@0PrashantYadav0 thanks — I can't re-run upstream workflows from a fork, so I rebased onto current develop (8251feb) to trigger a fresh run. No code changes; test.logpdf.js and test.factory.js still pass 5,014 each locally.

stdlib-bot · 2026-10-07T21:20:09Z

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/rayleigh/logpdf	$\\color{red}338/346$ $\\color{green}+0.00\\%$	$\\color{red}43/45$ $\\color{green}+0.00\\%$	$\\color{green}4/4$ $\\color{green}+0.00\\%$	$\\color{red}338/346$ $\\color{green}+0.00\\%$

The above coverage report was generated for the changes in this PR.

0PrashantYadav0 · 2026-10-08T02:42:06Z

Thanks for the rebase, CI is green now. The two inline comments are still open, though: logpdf( 1e154, 1e154 ) still returns -354.598 (exact -355.098), and small_scale.json, medium_scale.json and large_scale.json still need to be regenerated. Could you take a look?

…se/dists/rayleigh/logpdf` At `sigma = 1e154`, `sigma^2` is finite but `2*sigma^2` overflows, so the `x^2` term vanished: `logpdf( 1e154, 1e154 )` gave -354.598 (exact -355.098). Also regenerates the small, medium and large scale fixtures with the extended-precision `gen`.

Abhist17 · 2026-10-08T22:23:34Z

@0PrashantYadav0 sorry, I'd missed those two. Both are in 2a7b5aa:

logpdf( 1e154, 1e154 ): sigma^2 is finite there but 2*sigma^2 overflows, so the x^2 term vanished on the direct path. The guard now requires 2.0 * s2 < PINF (JS, factory and C), and it returns -355.09810432108304 in both main and factory.
small_scale.json, medium_scale.json and large_scale.json are regenerated with the extended-precision gen.

test.logpdf.js / test.factory.js pass 5,014 each.

0PrashantYadav0 · 2026-10-09T05:18:52Z

+		return ln( q ) - ( x2 / ( 2.0 * s2 ) );
+	}
+	// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
+	r = x / sigma;


Suggested change

r = x / sigma;

r = x / sigma;

if ( r === PINF ) {

return NINF;

}

Thanks, logpdf( 1e154, 1e154 ) is right now. The r === PINF guard from the earlier suggestion is still missing, though: logpdf( 1e300, 1e-10 ) returns NaN (Infinity - Infinity) instead of -Infinity. Same in lib/factory.js:94 and src/main.c:71. Could you add a test for it too? Also, Distributions can be dropped from test/fixtures/julia/REQUIRE now.

Abhist17 requested a review from a team October 5, 2026 19:19

stdlib-bot added Statistics Issue or pull request related to statistical functionality. Needs Review A pull request which needs code review. labels Oct 5, 2026

0PrashantYadav0 self-requested a review October 7, 2026 11:24

0PrashantYadav0 suggested changes Oct 7, 2026

View reviewed changes

Abhist17 force-pushed the fix/rayleigh-logpdf-scale branch from 593df45 to 8251feb Compare October 7, 2026 21:16

Abhist17 mentioned this pull request Oct 8, 2026

fix: avoid overflow and underflow for scales far from 1 in stats/base/dists/rayleigh/pdf #15994

Open

0PrashantYadav0 reviewed Oct 9, 2026

View reviewed changes

Navigation Menu

Repository navigation

fix: avoid overflow and underflow for scales far from 1 in `stats/base/dists/rayleigh/logpdf` - #15923

fix: avoid overflow and underflow for scales far from 1 in `stats/base/dists/rayleigh/logpdf`#15923
Abhist17 wants to merge 2 commits into
stdlib-js:developfrom
Abhist17:fix/rayleigh-logpdf-scale

Abhist17 commented Oct 5, 2026

Uh oh!

0PrashantYadav0 left a comment

Uh oh!

0PrashantYadav0 Oct 7, 2026

Uh oh!

0PrashantYadav0 Oct 7, 2026

Uh oh!

Abhist17 commented Oct 7, 2026

Uh oh!

stdlib-bot commented Oct 7, 2026 •

edited

Loading

Uh oh!

0PrashantYadav0 commented Oct 8, 2026

Uh oh!

Abhist17 commented Oct 8, 2026

Uh oh!

0PrashantYadav0 Oct 9, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Conversation

Abhist17 commented Oct 5, 2026

Description

Related Issues

Questions

Other

Checklist

AI Assistance

Disclosure

Uh oh!

0PrashantYadav0 left a comment

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 7, 2026

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 7, 2026

Choose a reason for hiding this comment

Uh oh!

Abhist17 commented Oct 7, 2026

Uh oh!

stdlib-bot commented Oct 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Coverage Report

Uh oh!

0PrashantYadav0 commented Oct 8, 2026

Uh oh!

Abhist17 commented Oct 8, 2026

Uh oh!

0PrashantYadav0 Oct 9, 2026

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

stdlib-bot commented Oct 7, 2026 •

edited

Loading