fix: avoid underflow for small `x` in `stats/base/dists/rayleigh/logcdf` by Abhist17 · Pull Request #15798 · stdlib-js/stdlib

Abhist17 · 2026-10-02T08:15:52Z

Description

What is the purpose of this pull request?

This pull request:

avoids the underflow of x^2 for small x in stats/base/dists/rayleigh/logcdf (JS, factory and C).

For x below about 1e-154, x^2 underflows to 0 and the function returns -Infinity, although the log-CDF is finite there:

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

logcdf( 1.0e-300, 1.0 );
// develop: -Infinity
// this PR: -1382.2442029769873 (= ln(1e-600 / 2))

With t = x^2/(2*sigma^2), ln(1 - e^{-t}) = ln(t) + ln(1 - t/2 + ...). Once ln(t) = 2*ln(x) - ln(2*sigma^2) is below -36, the correction is less than half an ULP of ln(t). So the function returns ln(t) there, which never underflows. Above it, the existing log1p/expm1 paths are unchanged, and x = 0 still gives -Infinity. This is the same change as #15797 for weibull/logcdf.

Tests: runner.jl gains a gen_small function. It evaluates the textbook formula in 2048-bit BigFloat (written with expm1, so it stays exact at that precision) and writes small_x.json, with x log-spaced over [1e-300, 0.01]. The new tests pass at 1 ULP in JS and native. On develop, 487 of the 4,011 assertions in test.logcdf.js fail, most of them returning -Infinity. The existing fixtures never reach the new branch and still pass at 0 ULP.

Related Issues

Does this pull request have any related issues?

None. Related to #15772 (rayleigh/cdf).

Questions

Any questions for reviewers of this pull request?

No.

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 evaluating every stats/base/dists logcdf just above the lower end of its support.

@stdlib-js/reviewers

`x^2` underflows to `0` for small `x`, so the function returned `-Infinity` although the log-CDF is `ln(x^2/(2*sigma^2))` to working precision there. Return `2*ln(x) - ln(2*sigma^2)` once it is below `-36`.

stdlib-bot · 2026-10-02T08:18:39Z

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/rayleigh/logcdf	$\\color{red}328/330$ $\\color{red}-0.61\\%$	$\\color{red}37/39$ $\\color{red}-5.13\\%$	$\\color{green}4/4$ $\\color{green}+0.00\\%$	$\\color{red}328/330$ $\\color{red}-0.61\\%$

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

0PrashantYadav0

Question: #15779 and #15797 regenerate the existing fixtures with BigFloat, while this PR keeps the Distributions ones. Is the difference intentional?

0PrashantYadav0 · 2026-10-02T11:25:36Z

 		return NINF;
 	}
 	s2 = pow( sigma, 2.0 );
+	lt = ( 2.0*ln( x ) ) - ln( 2.0*s2 );


Suggested change

lt = ( 2.0*ln( x ) ) - ln( 2.0*s2 );

r = x / sigma;

if ( r < FLOAT64_SMALLEST_NORMAL ) {

lt = ( 2.0*( ln( x ) - ln( sigma ) ) ) - LN2;

} else {

lt = ( 2.0*ln( r ) ) - LN2;

}

2*ln(x) - ln(2*s2) cancels when sigma is far from 1. logcdf( 1e-108, 1e-100 ) is now 5 ULP off (correctly rounded on develop). logcdf( 1e-170, 1e-160 ) is wrong because sigma^2 is subnormal, and logcdf( 1e-300, 1e200 ) still returns -Infinity. With x/sigma, as in #15797, all of these are within 1 ULP of mpmath for sigma in [1e-300, 1e300].

Same in lib/factory.js:82. Needs var r;, @stdlib/constants/float64/ln-two and @stdlib/constants/float64/smallest-normal.

0PrashantYadav0 · 2026-10-02T11:25:36Z

 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
 	s2 = stdlib_base_pow( sigma, 2.0 );
+	lt = ( 2.0 * stdlib_base_ln( x ) ) - stdlib_base_ln( 2.0 * s2 );


Suggested change

lt = ( 2.0 * stdlib_base_ln( x ) ) - stdlib_base_ln( 2.0 * s2 );

r = x / sigma;

if ( r < STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {

lt = ( 2.0 * ( stdlib_base_ln( x ) - stdlib_base_ln( sigma ) ) ) - STDLIB_CONSTANT_FLOAT64_LN2;

} else {

lt = ( 2.0 * stdlib_base_ln( r ) ) - STDLIB_CONSTANT_FLOAT64_LN2;

}

Same as lib/main.js:85. Also needs double r;, the ln_two.h/smallest_normal.h includes, and both deps in manifest.json.

0PrashantYadav0 · 2026-10-02T11:25:36Z

+x = exp10.( range( -300.0, stop = -2.0, length = 1000 ) );
+sigma = ( rand( 1000 ) .* 10.0 ) .+ 0.1;


Suggested change

x = exp10.( range( -300.0, stop = -2.0, length = 1000 ) );

sigma = ( rand( 1000 ) .* 10.0 ) .+ 0.1;

sigma = exp10.( ( rand( 1000 ) .* 200.0 ) .- 100.0 );

x = sigma .* exp10.( range( -150.0, stop = -4.0, length = 1000 ) );

With sigma in [0.1, 10.1], this fixture can't catch the cancellation above. With this change, the current code fails 49/1000 points at 1 ULP and the x/sigma version passes all of them.

Co-authored-by: 0PrashantYadav0 <0PrashantYadav0@users.noreply.github.com>

Abhist17 · 2026-10-03T18:30:48Z

@0PrashantYadav0 thanks, applied in 0dbba04:

r = x/sigma with the ln(x) - ln(sigma) fallback for a subnormal ratio, in JS, factory and C, plus the ln-two and smallest-normal dependencies.
small_x.json now uses your wide sigma range. The previous commit fails 54 of its points, and this one passes all of them at 1 ULP in JS and native. logcdf( 1e-108, 1e-100 ) is -37.53450866846468 and logcdf( 1e-300, 1e200 ) is -2303.278240174606.

On the fixtures: yes, that was intentional. The old Distributions fixtures here never reach the new branch (x is at least about 1e-3 there), and they still pass at 0 ULP, so I left them alone to keep the diff to the fix. In #15779 and #15797 the change does affect values the existing fixtures cover, and those fixtures had the same double-precision expression, so they had to be regenerated. I'm happy to regenerate these too for consistency if you prefer.

0PrashantYadav0 · 2026-10-04T04:11:29Z

@Abhist17 Consistency is preferred but we have to discuss this with others first ( to find out which one is required according to our needs )

fix: avoid underflow for small x in stats/base/dists/rayleigh/logcdf

59187c7

`x^2` underflows to `0` for small `x`, so the function returned `-Infinity` although the log-CDF is `ln(x^2/(2*sigma^2))` to working precision there. Return `2*ln(x) - ln(2*sigma^2)` once it is below `-36`.

Abhist17 requested a review from a team October 2, 2026 08:15

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

0PrashantYadav0 suggested changes Oct 2, 2026

View reviewed changes

fix: evaluate ln(t) from x/sigma

0dbba04

Co-authored-by: 0PrashantYadav0 <0PrashantYadav0@users.noreply.github.com>

Navigation Menu

Uh oh!

fix: avoid underflow for small `x` in `stats/base/dists/rayleigh/logcdf` - #15798

fix: avoid underflow for small `x` in `stats/base/dists/rayleigh/logcdf`#15798
Abhist17 wants to merge 2 commits into
stdlib-js:developfrom
Abhist17:fix/rayleigh-logcdf-small-x

Abhist17 commented Oct 2, 2026

Uh oh!

stdlib-bot commented Oct 2, 2026 •

edited

Loading

Uh oh!

0PrashantYadav0 left a comment

Uh oh!

0PrashantYadav0 Oct 2, 2026

Uh oh!

0PrashantYadav0 Oct 2, 2026

Uh oh!

0PrashantYadav0 Oct 2, 2026

Uh oh!

Abhist17 commented Oct 3, 2026

Uh oh!

0PrashantYadav0 commented Oct 4, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

-	lt = ( 2.0*ln( x ) ) - ln( 2.0*s2 );
+	r = x / sigma;
+	if ( r < FLOAT64_SMALLEST_NORMAL ) {
+		lt = ( 2.0*( ln( x ) - ln( sigma ) ) ) - LN2;
+	} else {
+		lt = ( 2.0*ln( r ) ) - LN2;
+	}

		x = exp10.( range( -300.0, stop = -2.0, length = 1000 ) );
		sigma = ( rand( 1000 ) .* 10.0 ) .+ 0.1;

Uh oh!

Conversation

Abhist17 commented Oct 2, 2026

Description

Related Issues

Questions

Other

Checklist

AI Assistance

Disclosure

Uh oh!

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

Uh oh!

Coverage Report

Uh oh!

0PrashantYadav0 left a comment

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 2, 2026

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 2, 2026

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 2, 2026

Choose a reason for hiding this comment

Uh oh!

Abhist17 commented Oct 3, 2026

Uh oh!

0PrashantYadav0 commented Oct 4, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

stdlib-bot commented Oct 2, 2026 •

edited

Loading