fix: use expm1 in stats/base/dists/rayleigh/cdf

[Abhist17]

Description

What is the purpose of this pull request?

This pull request:

• evaluates stats/base/dists/rayleigh/cdf as -expm1( -x^2 / (2*sigma^2) ) instead of 1 - exp( -x^2 / (2*sigma^2) ) (JS, factory and C).

Subtracting the exponential from 1 cancels whenever the CDF is small, and gives 0 once x^2/(2*sigma^2) < 2^-53:

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

cdf( 1.0e-10, 1.0 );
// develop: 0.0
// this PR: 5.0000000000000005e-21

cdf( 1.0e-5, 1.0 );
// develop: 5.000000413701855e-11
// this PR: 4.999999999875001e-11

rayleigh/logcdf already uses expm1 for the same reason.

Fixtures and tests:

• runner.jl now evaluates the textbook formula in 2048-bit BigFloat instead of Distributions.cdf. Distributions has the same 1 - exp cancellation, which is why the old fixtures passed at 0 ULP. Distributions is dropped from REQUIRE, and the three fixtures are regenerated with the same input ranges.
• Against the exact values, develop is off by up to ~2.9 million ULP on medium_scale.json. With this change all fixtures pass at 2 ULP in JS and native, so the tolerances go from 0 to 2.
• There's a new small_x.json with x log-spaced over [1e-150, 0.1].
• On develop, 1,689 of the 4,011 assertions in test.cdf.js fail. With this change they all pass.

Related Issues

Does this pull request have any related issues?

None.

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

[stdlib-bot]

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/rayleigh/cdf	$\\color{green}284/284$ $\\color{green}+0.00\\%$	$\\color{green}27/27$ $\\color{green}+0.00\\%$	$\\color{green}4/4$ $\\color{green}+0.00\\%$	$\\color{green}284/284$ $\\color{green}+0.00\\%$

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

[0PrashantYadav0]

Same x/sigma issue as #15798 (comment).

[0PrashantYadav0]

Suggested change

	s2 = pow( sigma, 2.0 );
	return 1.0 - exp( -pow( x, 2.0 ) / ( 2.0 * s2 ) );
	return -expm1( -pow( x, 2.0 ) / ( 2.0 * s2 ) );
	return -expm1( -0.5 * pow( x / sigma, 2.0 ) );

Squaring x and sigma separately overflows/underflows: cdf( 1e200, 1e200 ) and cdf( 1e-200, 1e-200 ) return NaN, cdf( 1e150, 1e200 ) returns 0. Taking the ratio first (as in weibull/cdf) stays within 2 ULP of mpmath for sigma in [1e-200, 1e200].

Same in lib/factory.js:75; drop s2 in both.

[0PrashantYadav0]

Suggested change

	s2 = stdlib_base_pow( sigma, 2.0 );
	return 1.0 - stdlib_base_exp( -stdlib_base_pow( x, 2.0 ) / ( 2.0 * s2 ) );
	return -stdlib_base_expm1( -stdlib_base_pow( x, 2.0 ) / ( 2.0 * s2 ) );
	return -stdlib_base_expm1( -0.5 * stdlib_base_pow( x / sigma, 2.0 ) );

Same as lib/main.js:77; also drop double s2;.

[0PrashantYadav0]

Suggested change

	# Small `x`:
	x = exp10.( range( -150.0, stop = -1.0, length = 1000 ) );
	sigma = ( rand( 1000 ) .* 10.0 ) .+ 0.1;
	gen( x, sigma, "small_x.json" );
	# Small `x`:
	sigma = exp10.( ( rand( 1000 ) .* 400.0 ) .- 200.0 );
	x = sigma .* exp10.( range( -150.0, stop = -1.0, length = 1000 ) );
	gen( x, sigma, "small_x.json" );

With sigma in [0.1, 10.1] this fixture can't catch the issue above. With this range the current code fails 415/1000 points at 2 ULP (115 NaN), the ratio version none.

[0PrashantYadav0]

Suggested change

	function gen( x, sigma, name )
	z = Array{Float64}( undef, length(x) );
	for i in eachindex(x)
	z[ i ] = cdf( Rayleigh( sigma[i] ), x[i] );
	z = setprecision( BigFloat, 2048 ) do
	Float64.( 1 .- exp.( .-( BigFloat.( x ) .^ 2 ) ./ ( 2 .* BigFloat.( sigma ) .^ 2 ) ) )
	end
	function gen( x, sigma, name )
	# The CDF is evaluated in extended precision using the textbook formula `1 - exp(-x^2 / (2*sigma^2))`, as subtracting the exponential from `1` in double precision cancels for small `x/sigma`...
	z = setprecision( BigFloat, 2048 ) do
	Float64.( 1 .- exp.( .-( BigFloat.( x ) .^ 2 ) ./ ( 2 .* BigFloat.( sigma ) .^ 2 ) ) );
	end

And drop the docstring paragraph at line 26. Matches kgryte's follow-up commits to runner.jl in #15774.

[Abhist17]

@0PrashantYadav0 thanks — applied in 8b3cf42:

• -expm1( -0.5 * pow( x / sigma, 2.0 ) ) in lib/main.js, lib/factory.js and src/main.c, with s2 dropped from all three.
• runner.jl: your gen (comment moved inside, semicolon) and small-x block, docstring paragraph removed. I regenerated only small_x.json; the other three fixtures come from rand, so I restored them rather than churn them.

Your three cases now give 0.3934693402873666, 0.3934693402873666 and 5e-101, in JS and in the native addon. Against the old code the regenerated fixture fails 401 of its points; with the change, test.cdf.js / test.factory.js / test.native.js pass 4011 / 4013 / 4011.

[0PrashantYadav0]

LGTM