fix: avoid cancellation in stats/base/dists/kumaraswamy/logcdf

[Abhist17]

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

[stdlib-bot]

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/kumaraswamy/logcdf	$\\color{green}365/365$ $\\color{green}+100.00\\%$	$\\color{green}35/35$ $\\color{green}+100.00\\%$	$\\color{green}4/4$ $\\color{green}+100.00\\%$	$\\color{green}365/365$ $\\color{green}+100.00\\%$

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