fix: avoid cancellation in `stats/base/dists/kumaraswamy/logcdf` by Abhist17 · Pull Request #15768 · stdlib-js/stdlib

Abhist17 · 2026-10-01T23:31:24Z

Description

What is the purpose of this pull request?

This pull request:

fixes two cancellations in stats/base/dists/kumaraswamy/logcdf (JS, factory and C),
adds the fixtures the tests had a TODO for.

The function computed ln( 1 - pow( 1 - pow( x, a ), b ) ):

For small x^a, 1 - x^a rounds away its low-order bits, and 1 - (1 - x^a)^b cancels. For x^a < 2^-53 the result is -Infinity.
For x close to 1, 1 - x^a cancels (x^a was already rounded), and ln of a value close to 1 loses the digits of a small result. It returns 0 once (1 - x^a)^b < 2^-53.

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

logcdf( 1.0e-10, 2.0, 3.0 );
// develop: -Infinity
// this PR: -44.95308957121281 (= ln(3e-20))

logcdf( 0.9985187737256925, 1.668922808849776, 5.306636551597105 );
// develop: -1.4654943925052174e-14
// exact:   -1.4613769592956135e-14  (this PR)

logcdf( 0.9999999999910474, 2.3437963327045774, 1.5267961353982975 );
// develop: 0.0
// this PR: -4.9736148447726656e-17

The function now computes v = ln(1 - x^a) as log1p( -x^a ) when x^a < 0.5, and as ln( -expm1( a*ln(x) ) ) otherwise, so 1 - x^a is never formed. It then returns log1mexp( b*v ).

Tests:

There's a new test/fixtures/julia/runner.jl that evaluates the textbook formula ln(1 - (1 - x^a)^b) in 2048-bit BigFloat. It writes three fixtures:
- data.json: random x, with a and b in [0.5, 5.5],
- small_x.json: x log-spaced over [1e-100, 0.1], with a <= 3 so that x^a doesn't underflow,
- large_x.json: 1 - x log-spaced over [1e-15, 0.1].
The "evaluates the logcdf" tests only checked that the result is a number (// TODO: Add fixtures). They are replaced with tests against the fixtures.
Max error, JS and native, which sets the tolerances:

fixture develop this PR

data ~1.3e13 ULP 25 ULP

small_x -Infinity 1 ULP

large_x 0 148 ULP

The remaining error comes from rounding in b * v. A 1-ULP change in b moves the exact result by |b ln(1 - x^a)| ULP, which reaches ~150 at the far end of large_x. So this is about as close as double precision gets without extended arithmetic. On develop, 2,140 of the 3,027 assertions in test.logcdf.js fail. With this change they all pass.

Related Issues

Does this pull request have any related issues?

None. Companion to #15767 (kumaraswamy/cdf).

Questions

Any questions for reviewers of this pull request?

When x^a underflows (e.g. x = 1e-100, a = 5), the result is still -Infinity, as on develop. It could return ln(b) + a*ln(x) there. I left that out to keep this PR to the cancellation, but I'm happy to add it.

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 reading the stats/base/dists functions for 1 - u cancellation.

@stdlib-js/reviewers

`ln( 1 - (1 - x^a)^b )` cancels twice: `1 - x^a` rounds away the low-order bits of a small `x^a`, and `1 - (1 - x^a)^b` loses the digits of a small CDF, as `ln` of a value close to `1` does for a CDF close to `1`. Evaluate `ln(1-x^a)` with `log1p` or `expm1`, depending on the size of `x^a`, and finish with `log1mexp`. Add fixtures generated from the textbook formula in extended precision.

stdlib-bot · 2026-10-01T23:33:54Z

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/kumaraswamy/logcdf	$\\color{red}381/385$ $\\color{red}-1.04\\%$	$\\color{red}41/43$ $\\color{red}-4.65\\%$	$\\color{green}4/4$ $\\color{green}+0.00\\%$	$\\color{red}381/385$ $\\color{red}-1.04\\%$

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

0PrashantYadav0 · 2026-10-08T11:35:56Z

+	u = pow( x, a );
+	if ( u < 0.5 ) {
+		v = log1p( -u );
+	} else {
+		v = ln( -expm1( a * ln( x ) ) );
+	}
+	return log1mexp( b * v );


Suggested change

u = pow( x, a );

if ( u < 0.5 ) {

v = log1p( -u );

} else {

v = ln( -expm1( a * ln( x ) ) );

}

return log1mexp( b * v );

u = pow( x, a );

if ( u < FLOAT64_SMALLEST_NORMAL ) {

// `x^a` is subnormal or underflows, so `ln(1-(1-x^a)^b) = ln(b) + a*ln(x)`:

return ln( b ) + ( a * ln( x ) );

}

if ( u < 0.5 ) {

v = log1p( -u );

} else {

v = ln( -expm1( a * ln( x ) ) );

}

t = b * v;

if ( t > -FLOAT64_SMALLEST_NORMAL ) {

return ln( b ) + ln( -v );

}

return log1mexp( t );

On your question: yes, let's add it here, plus the b*v underflow (same as #15797 (comment)). logcdf( 1.0e-200, 3.0, 2.0 ), logcdf( 0.9, 1.0e5, 2.0 ) and logcdf( 1.0e-100, 2.0, 1.0e-200 ) return -Infinity (exact ≈ -1380.858, -10535.358, -921.034); with this they're within 2 ULP. Same in lib/factory.js:86-92 and src/main.c:59-65; needs var t; and FLOAT64_SMALLEST_NORMAL (plus the manifest.json dep).

0PrashantYadav0 · 2026-10-08T11:35:56Z

+# Small `x` (with `a` small enough that `x^a` does not underflow):
+x = exp10.( range( -100.0, stop = -1.0, length = 1000 ) );
+a = ( rand( 1000 ) .* 2.5 ) .+ 0.5;


Suggested change

# Small `x` (with `a` small enough that `x^a` does not underflow):

x = exp10.( range( -100.0, stop = -1.0, length = 1000 ) );

a = ( rand( 1000 ) .* 2.5 ) .+ 0.5;

# Small `x`:

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

a = ( rand( 1000 ) .* 5.0 ) .+ 0.5;

So the fixture covers the underflow above. Needs setprecision( BigFloat, 8192 ) on line 36 (at 2048 bits 1 - x^a rounds to 1 for x^a < 2^-2048), and the small_x tolerance goes to 2 ULP.

…wamy/logcdf` Per review: when `x^a` is subnormal the log-CDF is `ln(b) + a*ln(x)`, and when `b*v` underflows it is `ln(b) + ln(-v)`. `logcdf( 1e-200, 3, 2 )`, `logcdf( 0.9, 1e5, 2 )` and `logcdf( 1e-100, 2, 1e-200 )` returned -Infinity. The small-`x` fixture now reaches `1e-300` with `a` up to 5.5, evaluated at 8192 bits, at 2 ULP.

Abhist17 · 2026-10-08T22:26:48Z

@0PrashantYadav0 thanks — applied in c83cde0:

Your branch structure in lib/main.js, lib/factory.js and src/main.c: the subnormal-x^a return and the b*v underflow branch, with var t; / double t;, the smallest-normal require and include, and the dep in all three manifest.json lists. Your three cases now give -1380.8579086158677, -10535.358418602069 and -921.0340371976183 (same in factory).
runner.jl: setprecision( BigFloat, 8192 ), and the small-x block over [1e-300, 0.1] with a in [0.5, 5.5]. Only small_x.json is regenerated; the other two come from rand, so I left them.
small_x tolerance is 2 ULP in all three test files; the fixture's max is 2.

test.logcdf.js / test.factory.js / test.native.js pass 3,027 / 3,031 / 3,027.

0PrashantYadav0

LGTM

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

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

Abhist17 mentioned this pull request Oct 1, 2026

fix: avoid cancellation for small p in stats/base/dists/kumaraswamy/quantile #15769

Open

5 of 7 tasks

0PrashantYadav0 suggested changes Oct 8, 2026

View reviewed changes

0PrashantYadav0 mentioned this pull request Oct 8, 2026

fix: avoid cancellation for small x in stats/base/dists/kumaraswamy/cdf #15767

Open

5 of 7 tasks

0PrashantYadav0 approved these changes Oct 9, 2026

View reviewed changes

Navigation Menu

Repository navigation

fix: avoid cancellation in `stats/base/dists/kumaraswamy/logcdf` - #15768

fix: avoid cancellation in `stats/base/dists/kumaraswamy/logcdf`#15768
Abhist17 wants to merge 2 commits into
stdlib-js:developfrom
Abhist17:fix/kumaraswamy-logcdf-small-x

Abhist17 commented Oct 1, 2026

Uh oh!

stdlib-bot commented Oct 1, 2026 •

edited

Loading

Uh oh!

0PrashantYadav0 Oct 8, 2026

Uh oh!

0PrashantYadav0 Oct 8, 2026

Uh oh!

Abhist17 commented Oct 8, 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

fixture	`develop`	this PR
`data`	~1.3e13 ULP	25 ULP
`small_x`	`-Infinity`	1 ULP
`large_x`	`0`	148 ULP

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 8, 2026

Choose a reason for hiding this comment

Uh oh!

0PrashantYadav0 Oct 8, 2026

Choose a reason for hiding this comment

Uh oh!

Abhist17 commented Oct 8, 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