fix: use expm1 in stats/base/dists/lognormal/skewness

[Abhist17]

Description

What is the purpose of this pull request?

This pull request:

• evaluates stats/base/dists/lognormal/skewness as ( exp(sigma^2) + 2 ) * sqrt( expm1(sigma^2) ) instead of ( exp(sigma^2) + 2 ) * sqrt( exp(sigma^2) - 1 ) (JS and C).
• returns 3*sigma for sigma < 1e-8, where it is within 1 ULP of the series 3*sigma*(1 + 7*sigma^2/12 + ...); below sigma ≈ 1.5e-154 sigma^2 underflows and the formula above decays to 0.

exp(sigma^2) - 1 cancels when sigma is small, and is 0 for sigma^2 < 2^-53:

var skewness = require( '@stdlib/stats/base/dists/lognormal/skewness' );

skewness( 0.0, 1.0e-9 );
// develop: 0.0
// this PR: 3.0000000000000004e-9

skewness( 0.0, 1.0e-5 );
// develop: 0.000030000001242105537
// this PR: 0.00003000000000175

Fixtures and tests:

• runner.jl now evaluates the textbook formula in 2048-bit BigFloat instead of Distributions.skewness, which uses the same double-precision expression. Distributions is dropped from REQUIRE, and data.json is regenerated with the same input ranges.
• There's a new small_sigma.json with sigma log-spaced over [1e-150, 0.1]. It passes at 1 ULP in JS and native.
• On data.json, the max error goes from 4,598 ULP to 22 ULP. The remaining 22 comes at sigma ≈ 4.6, where the rounding of sigma^2 is multiplied by the size of exp(sigma^2). The tolerance goes from 0 to 22. The old 0 only held because the fixture used the same double-precision formula.
• On develop, 1,005 of the 1,110 assertions in test.js fail. With this change they all pass.

Related Issues

Does this pull request have any related issues?

None. Same cancellation as #15774 (lognormal/variance).

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 exp(...) - 1.

---

@stdlib-js/reviewers

[stdlib-bot]

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/lognormal/skewness	$\\color{green}184/184$ $\\color{green}+0.00\\%$	$\\color{green}11/11$ $\\color{green}+0.00\\%$	$\\color{green}2/2$ $\\color{green}+0.00\\%$	$\\color{green}184/184$ $\\color{green}+0.00\\%$

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

[0PrashantYadav0]

Suggested change

	z = setprecision( BigFloat, 2048 ) do
	es2 = exp.( BigFloat.( sigma ) .^ 2 );
	Float64.( ( es2 .+ 2 ) .* sqrt.( es2 .- 1 ) )
	# The skewness is evaluated in extended precision using the textbook formula `(exp(sigma^2) + 2) * sqrt(exp(sigma^2) - 1)`, as subtracting `1` from `exp(sigma^2)` in double precision cancels for small `sigma`.
	z = setprecision( BigFloat, 2048 ) do
	es2 = exp.( BigFloat.( sigma ) .^ 2 );
	Float64.( ( es2 .+ 2 ) .* sqrt.( es2 .- 1 ) );

Also drop the docstring paragraph at lines 26-27. This matches the docs: move comment and style: add semicolon commits on #15774.

[Abhist17]

Applied in cffe135: the comment moved into gen, the docstring paragraph dropped, and the semicolon added. The precision is 4096 rather than 2048 because the fixtures now go down to 1e-300; see the other thread.

[0PrashantYadav0]

Suggested change

	s2 = sigma * sigma;
	return ( exp( s2 ) + 2.0 ) * sqrt( expm1( s2 ) );
	if ( sigma < 1.0e-8 ) {
	// For small `sigma`, skewness = 3*sigma*(1 + 7*sigma^2/12 + ...):
	return 3.0 * sigma;
	}
	s2 = sigma * sigma;
	return ( exp( s2 ) + 2.0 ) * sqrt( expm1( s2 ) );

On the question in the description: the branch is small enough to include here. Right now skewness( 0.0, 1.0e-163 ) returns 0, and skewness( 0.0, 1.0e-160 ) returns 2.999983300727547e-160 (exact 3.0e-160). 3.0*sigma is within 1 ULP of the exact value for all sigma < 1e-8. The same change applies to src/main.c:45-46, and small_sigma.json could then go down to 1e-300. Okay to fold in, or would you rather do it as a follow-up?

[Abhist17]

Folded in (cffe135): 3.0 * sigma for sigma < 1.0e-8 in main.js and src/main.c, and small_sigma.json now runs from 1e-300 to 1e-1, so both branches are under the 1-ULP test.

One catch when extending the fixtures: at sigma = 1e-300, exp(sigma^2) - 1 cancels even at 2048 bits (about 616 digits against 1 + 1e-600), and 7 of the 1000 fixtures came out 1 ULP off against an 8192-bit evaluation. The runner now uses 4096 bits, which matches 16384 bits on all 1000.