fix: avoid overflow in stats/base/dists/poisson/entropy for large lambda

[Abhist17]

Description

What is the purpose of this pull request?

This pull request:

• fixes stats/base/dists/poisson/entropy (JS and C) returning Infinity for lambda >= ~80 and NaN for lambda >= ~745.

var entropy = require( '@stdlib/stats/base/dists/poisson/entropy' );

entropy( 50.0 );   // 3.373266261170244  (ok)
entropy( 80.0 );   // Infinity           (expected ~3.6089)
entropy( 1000.0 ); // NaN                (expected ~4.8727)

The implementation evaluates lambda*(1-ln(lambda)) + exp(-lambda) * sum_k lambda^k ln(k!)/k!. Both lambda^k and k! overflow long before the series converges once lambda is more than a few dozen, so the terms become Infinity and then Infinity/Infinity. The existing fixtures only cover lambda in [0, 20], which is why the tests never saw it. (Julia's entropy( Poisson( lambda ) ), which generated them, switches to an asymptotic expansion at lambda >= 50 for the same reason.)

The fix evaluates H = -sum_k p_k ln(p_k) directly, with weights relative to the mode, w_k = p_k/p_m:

• ln(w_k) comes from the recurrence w_{k+1} = w_k * lambda/(k+1) (and downward, w_{k-1} = w_k * k/lambda), so nothing overflows and there is no cancellation between k*ln(lambda), lambda and ln(k!), each of which is far larger than the entropy.
• Since the p_k sum to one, p_m = 1/W with W = sum w_k, so H = ln(W) - sum w_k ln(w_k) / W. No gammaln is needed.
• Each ln(lambda/(k+1)) near the mode goes through log1p, where the subtraction is exact (the two operands are within a factor of two).
• For lambda >= 5000, it uses the asymptotic expansion ln(2*pi*e*lambda)/2 - 1/(12 lambda) - 1/(24 lambda^2) - 19/(360 lambda^3), whose truncation error there is below 1e-16 relative. That keeps the loop to at most about 1,300 iterations.

Accuracy against a 256-bit reference: at most 3 ULP across lambda in [0.005, 5e5], compared with up to 76 ULP for the old series on [0, 20], where it did work.

Related Issues

Does this pull request have any related issues?

No.

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.

• Fixtures: runner.jl now computes -sum p_k ln(p_k) in BigFloat (256 bits) over lambda ± 40 sd, and no longer uses Distributions.entropy. Its double-precision series was itself up to 76 ULP off on the existing data.json, and it becomes an asymptotic approximation at lambda >= 50. data.json is regenerated, and medium.json (lambda in [20, 5000]) and large.json (lambda in [5e3, 5e5]) are new. The tolerance is tightened from 20 to 4 ULP, and the worst case is 3.
• With the old implementation, 560 of the 806 assertions in test/test.js fail. With this change, all pass in both test.js and test.native.js, and the poisson/ctor tests also pass.
• entropy( Infinity ) now returns Infinity (the expansion's limit) where it used to return NaN.

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. The bug was found by comparing every stats/base/dists function against SciPy and mpmath on random parameters.

---

@stdlib-js/reviewers

[stdlib-bot]

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/poisson/entropy	$\\color{green}257/257$ $\\color{green}+0.00\\%$	$\\color{green}22/22$ $\\color{green}+0.00\\%$	$\\color{green}3/3$ $\\color{green}+0.00\\%$	$\\color{green}257/257$ $\\color{green}+0.00\\%$

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

[kgryte]

@Planeshifter Mind taking a look at this PR?