fix(stats): stabilize truncated-normal tails and narrow intervals

[AminMohamed-3]

Valid tail intervals such as [9, 10] make truncnorm.pdf and ppf return infinities and can make the scalar CDF raise ZeroDivisionError: subtracting two saturated normal CDFs gives a zero normalization mass. This implements a fix for #259.

Keep interval probabilities in logarithmic form, use a factored eight-node quadrature for nearly coincident bounds, and invert log probabilities with four fixed Newton steps. The pure-NumPy kernels preserve the existing public wrapper and server dispatch path. Invalid probabilities, bounds and nonpositive scales now produce SciPy-style NaNs; broadcast, dtype-promotion, source-refusal and budget checks remain covered.

The numerical work is more expensive, so this proposal replaces the old base costs 28/51/81 with 315/844/1392 for PDF/CDF/PPF per broadcast output element. These are conservative analytical numerical bounds at weight 1, before configured dtype pricing, not hardware calibrations or measured averages. The complete branch derivation is in docs/reference/truncnorm-cost.md; source docstrings, generated-documentation formulas and regenerated ops.json agree. Eight quadrature nodes and four inverse iterations are part of that derivation.

Validation:

• 641 selected tests passed with the repository's frozen dependencies on Windows/Python 3.12, covering statistics, broadcast/dtype billing, cost conformance and server dispatch. A separate NumPy 2.4.6 run passed 201 stats/dispatch tests plus 21 relevant cost tests.
• Reported positive and negative tail cases, semi-infinite intervals, narrow/adjacent bounds, quantile monotonicity/inversion, joint broadcasting and invalid inputs are checked against SciPy and independent adaptive quadrature. SciPy is used only by tests.
• Full pyright src/flopscope tests passed with zero errors and four existing export warnings. Repository Ruff lint/format, diff check, generated-doc verification and the ops.json synchronization check passed. The full Linux CI matrix has not run locally.

Numerical limits: this addresses the reported 9- and 40-sigma cases and tested narrow intervals; it is not an arbitrary-finite-double accuracy guarantee. At extremely large standardized bounds (for example near 1e10), cancellation or overflow can still produce failures, and quantile spacing limits inversion on very narrow intervals.

Developed, independently reviewed and tested with AI assistance. The fdlibm-derived coefficients retain their original license notice.