fix: use `log1mexp` in `stats/base/dists/exponential/logcdf` by Abhist17 · Pull Request #15771 · stdlib-js/stdlib

Abhist17 · 2026-10-01T23:36:38Z

Description

What is the purpose of this pull request?

This pull request:

evaluates stats/base/dists/exponential/logcdf as log1mexp( lambda*x ) instead of log1p( -exp( -lambda*x ) ) (JS, factory and C).

For small lambda*x, exp(-lambda*x) is close to 1, and log1p of minus that loses the digits of lambda*x. For lambda*x < 2^-53 the result is -Infinity:

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

logcdf( 1.0e-20, 1.0 );
// develop: -Infinity
// this PR: -46.051701859880914 (= ln(1e-20))

logcdf( 1.0e-10, 1.0 );
// develop: -23.02585084720009
// this PR: -23.025850929990458

log1mexp uses ln( -expm1( -x ) ) for x <= ln(2) and log1p( -exp( -x ) ) above it, so large lambda*x takes the same path as before.

Tests: there's a new test/fixtures/julia/small_x.json (x log-spaced over [1e-300, 0.1], lambda in [0.1, 10.1]). runner.jl evaluates ln(1 - exp(-lambda*x)) in 2048-bit BigFloat. The new tests in test.logcdf.js, test.factory.js and test.native.js pass at 1 ULP. On develop, 996 of those assertions fail.

The R fixture keeps its 511-ULP tolerance. Its worst cases have lambda*x ≈ 267, where the result is -exp(-267), and the half-ULP rounding of the product lambda*x alone is multiplied by lambda*x. That's the conditioning of the function, not this formula.

Related Issues

Does this pull request have any related issues?

None. Companion to #15770 (exponential/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 scanning stats/base/dists for 1 - exp(...).

@stdlib-js/reviewers

`log1p( -exp( -lambda*x ) )` cancels for small `lambda*x` and is `-Infinity` for `lambda*x < 2^-53`.

stdlib-bot · 2026-10-01T23:38:41Z

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/exponential/logcdf	$\\color{green}291/291$ $\\color{green}+0.00\\%$	$\\color{green}24/24$ $\\color{green}+0.00\\%$	$\\color{green}4/4$ $\\color{green}+0.00\\%$	$\\color{green}291/291$ $\\color{green}+0.00\\%$

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

0PrashantYadav0 · 2026-10-05T15:33:35Z

 		return NINF;
 	}
-	return log1p( -exp( -lambda * x ) );
+	return log1mexp( lambda * x );


Suggested change

return log1mexp( lambda * x );

t = lambda * x;

if ( t < FLOAT64_SMALLEST_NORMAL ) {

return ln( x ) + ln( lambda );

}

return log1mexp( t );

lambda * x can underflow: logcdf( 1e-200, 1e-200 ) returns -Infinity (exact ≈ -921.034). Same fix as #15797 (comment). Same in lib/factory.js:71; needs var t;, ln and FLOAT64_SMALLEST_NORMAL.

0PrashantYadav0 · 2026-10-05T15:33:35Z

 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
-	return stdlib_base_log1p( -stdlib_base_exp( -lambda * x ) );
+	return stdlib_base_log1mexp( lambda * x );


Suggested change

return stdlib_base_log1mexp( lambda * x );

t = lambda * x;

if ( t < STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {

return stdlib_base_ln( x ) + stdlib_base_ln( lambda );

}

return stdlib_base_log1mexp( t );

Same as lib/main.js:74; also needs double t;, the ln.h/smallest_normal.h includes and both deps in manifest.json.

0PrashantYadav0 · 2026-10-05T15:33:36Z

+
+# Small `x`:
+x = exp10.( range( -300.0, stop = -1.0, length = 1000 ) );
+lambda = ( rand( 1000 ) .* 10.0 ) .+ 0.1;


Suggested change

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

lambda = exp10.( ( rand( 1000 ) .* 41.0 ) .- 40.0 );

With lambda in [0.1, 10.1], lambda*x never gets subnormal, so the case above isn't tested (with this range the current code fails 33/1000). Same as #15797 (comment).

…ogcdf` `lambda * x` underflows for small arguments, so `logcdf( 1e-200, 1e-200 )` returned -Infinity instead of about -921.034. Below the smallest normal, `ln(1 - exp(-t))` is `ln(t)` to within rounding, so evaluate it as `ln(x) + ln(lambda)` there, as `weibull/logcdf` does. The small-`x` fixture now draws `lambda` from [1e-40, 1e1] so it reaches that branch; the previous code fails 36 of its points.

Abhist17 · 2026-10-05T19:07:35Z

@0PrashantYadav0 thanks — applied in 21b42d6:

The t = lambda * x / ln( x ) + ln( lambda ) branch in lib/main.js, lib/factory.js and src/main.c, with var t; / double t;, the ln and smallest-normal requires and includes, and both deps added to all three manifest.json dependency lists.
runner.jl: lambda = exp10.( ( rand( 1000 ) .* 41.0 ) .- 40.0 );, and small_x.json regenerated.

logcdf( 1e-200, 1e-200 ) is now -921.0340371976183 in JS and native. The boundaries are unchanged: logcdf( 0.0, 1.0 ) and logcdf( 1.0, 0.0 ) still return -Infinity, since t = 0 takes the new branch and ln( 0 ) is -Infinity either way. Against the old code the regenerated fixture fails 36 of 1000 points (close to your 33 — different draw); with the change test.logcdf.js / test.factory.js / test.native.js pass 2013 each.

0PrashantYadav0

LGTM

fix: use log1mexp in stats/base/dists/exponential/logcdf

8d4291f

`log1p( -exp( -lambda*x ) )` cancels for small `lambda*x` and is `-Infinity` for `lambda*x < 2^-53`.

Abhist17 requested a review from a team October 1, 2026 23:36

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

0PrashantYadav0 suggested changes Oct 5, 2026

View reviewed changes

0PrashantYadav0 approved these changes Oct 7, 2026

View reviewed changes

Navigation Menu

Repository navigation

fix: use `log1mexp` in `stats/base/dists/exponential/logcdf` - #15771

fix: use `log1mexp` in `stats/base/dists/exponential/logcdf`#15771
Abhist17 wants to merge 2 commits into
stdlib-js:developfrom
Abhist17:fix/exponential-logcdf-small-x

Abhist17 commented Oct 1, 2026

Uh oh!

stdlib-bot commented Oct 1, 2026 •

edited

Loading

Uh oh!

0PrashantYadav0 Oct 5, 2026

Uh oh!

0PrashantYadav0 Oct 5, 2026

Uh oh!

0PrashantYadav0 Oct 5, 2026

Uh oh!

Abhist17 commented Oct 5, 2026

Uh oh!

0PrashantYadav0 left a comment

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

-	return log1mexp( lambda * x );
+	t = lambda * x;
+	if ( t < FLOAT64_SMALLEST_NORMAL ) {
+		return ln( x ) + ln( lambda );
+	}
+	return log1mexp( t );

-	return stdlib_base_log1mexp( lambda * x );
+	t = lambda * x;
+	if ( t < STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
+		return stdlib_base_ln( x ) + stdlib_base_ln( lambda );
+	}
+	return stdlib_base_log1mexp( t );

	lambda = ( rand( 1000 ) .* 10.0 ) .+ 0.1;
	lambda = exp10.( ( rand( 1000 ) .* 41.0 ) .- 40.0 );

Conversation

Abhist17 commented Oct 1, 2026

Description

Related Issues

Questions

Other

Checklist

AI Assistance

Disclosure

Uh oh!

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

Uh oh!

Coverage Report

Uh oh!

0PrashantYadav0 Oct 5, 2026

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 5, 2026

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 5, 2026

Choose a reason for hiding this comment

Uh oh!

Abhist17 commented Oct 5, 2026

Uh oh!

0PrashantYadav0 left a comment

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 1, 2026 •

edited

Loading