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

docs/reference/cost-model.md

@@ -1473,9 +1473,9 @@ per-element counts below even for a float32 input (e.g. `stats.uniform.pdf` on a
 | `stats.laplace.pdf` | 22 | DERIVED: \|x−loc\|(3) + exp(−z)(17) + /(2·scale)(2); weight 1.0 |
 | `stats.laplace.cdf` | 40 | DERIVED composite: two eager exp branches + arithmetic/select; weight 1.0 |
 | `stats.laplace.ppf` | 51 | DERIVED composite: two eager log branches + edge selects; weight 1.0 |
-| `stats.truncnorm.pdf` | 28 | DERIVED composite: norm.pdf + cdf normalization; weight 1.0 |
-| `stats.truncnorm.cdf` | 51 | DERIVED composite: affine + norm.cdf + boundary selects; weight 1.0 |
-| `stats.truncnorm.ppf` | 81 | DERIVED composite: affine + rational + Newton with erf+exp; weight 1.0 |
+| `stats.truncnorm.pdf` | 315 | DERIVED upper bound: domain/masks 39 + max(narrow 211, log-mass density 276); weight 1.0 |
+| `stats.truncnorm.cdf` | 844 | DERIVED upper bound: domain/masks 44 + max(narrow 374, three log masses and eager tail selection 800); weight 1.0 |
+| `stats.truncnorm.ppf` | 1392 | DERIVED upper bound: domain/masks 49 + max(narrow fixed Newton 1037, log-tail inverse 1343); weight 1.0 |
 | `stats.lognorm.pdf` | 62 | DERIVED composite: log + exp + arithmetic per element; weight 1.0 |
 | `stats.lognorm.cdf` | 70 | DERIVED composite: log + erf rational approx + arithmetic; weight 1.0 |
 | `stats.lognorm.ppf` | 106 | DERIVED composite: ndtri + exp; weight 1.0 |
@@ -1484,6 +1484,12 @@ per-element counts below even for a float32 input (e.g. `stats.uniform.pdf` on a
 
 Source: `src/flopscope/stats/`.
 
+Truncated-normal counts are fixed analytical numerical bounds, including the
+fixed four-step inverse and eight-node narrow-interval quadrature. They are
+not average timings or hardware-counter calibrations. See the
+[complete branch derivation and numerical limits](truncnorm-cost.md).
+The existing empty-output minimum of one billed element is retained.
+
 ---
 
 ### Window

docs/reference/truncnorm-cost.md

@@ -0,0 +1,59 @@
+# Truncated-normal numerical cost bound
+
+This is an analytical **upper bound for numerical work**, not a hardware calibration or a literal count of NumPy memory operations. It applies to `src/flopscope/stats/_truncnorm_kernels.py`, with four fixed Newton steps, eight fixed quadrature nodes, and import-time scalar thresholds.
+
+Conventions come from `docs/reference/cost-model.md`: FMA = 2; arithmetic, square root, comparison and conditional selection = 1; exp/expm1/log/log1p/logaddexp = 16; a fixed square = one multiplication; each output assignment/fill = 1. A zero/uninitialized allocation and indexing/view mechanics are not numerical operations. Input coercion and Python dispatch are excluded, as in the stats composite convention. All comparison masks are counted even when a branch is empty. Each `any` reduction is bounded by one operation per input element. This is deliberately conservative; mutually exclusive arithmetic branches use their maximum, and eager expressions (including both `where` values) are both counted.
+
+All constants below are per element entering the named helper, before the public float64 dtype multiplier. Write costs in common paths may conservatively cover incompatible endpoint cases simultaneously.
+
+## Helper bounds
+
+| Helper | Derivation | Bound |
+|---|---|---:|
+| `_erf`, restricted to abs(x) < 1.25 | masks/reductions/sign/fill 19 + max(first-region arithmetic/write 23, second-region arithmetic/write 29) | 48 |
+| `_erf`, unrestricted | same common 19 + max tail arithmetic/write 57 | 76 |
+| `_ndtri`, unrestricted | masks, endpoint writes, reductions 25 + max rational seed 40 + Newton correction 104 | 169 |
+| `_narrow` | two finite predicates, two logical conjunctions, subtraction, two absolutes, two maxima, multiplication, comparison | 11 |
+| `_local_integral` | 41 multiplications + one negation + eight subtractions + eight exp at 16 + seven reduction additions | 185 |
+| `log_ndtr`, nonpositive input | common masks/normalization/writes 20 + max(moderate branch 68, first rational tail 60, second rational tail 56) | 88 |
+| `log_ndtr`, arbitrary input | nonpositive bound 88 + optional exp/negation/log1p/write 34 | 122 |
+| `log_mass` | common predicates/fills 26 + max(narrow 207, reflected same-tail 218, central 228) | 254 |
+| `ndtri_log_lower` | common 10 + max(middle seed 186, tail seed 24) + 4 Newton steps of 146 | 780 |
+
+`log_mass` only calls `log_ndtr` on nonpositive arguments: its same-tail branch reflects positive bounds and its central branch uses `a < 0` and `-b < 0`. The restricted erf bound applies inside `log_ndtr`'s moderate branch because it explicitly filters abs(x/sqrt(2)) < 1.25. The general inverse's Newton step conservatively uses the unrestricted `log_ndtr` bound, even though ordinary iterates are nonpositive.
+
+The inverse-normal Newton sub-bound 104 is: CDF 79 (= division 1 + unrestricted erf 76 + add/multiply 2), PDF 19, safety comparison 1, eager subtract/divide/select 3, final subtraction/write 2. The log-domain Newton sub-bound 146 is: log CDF 122, derivative arithmetic/exp 20, correction subtract/divide/subtract/write 4.
+
+## Public numerical kernels
+
+| Method | Common | Narrow branch | Regular branch | Bound |
+|---|---:|---:|---:|---:|
+| PDF | 39 | 211 | 276 | **315** |
+| CDF | 44 | 374 | 800 | **844** |
+| PPF | 49 | 1037 | 1343 | **1392** |
+
+PDF's regular branch is one log mass (254), square/multiply/two subtractions (4), exp (16), scale division (1), and write (1). Its narrow branch uses a 185-cost integral, 21-cost exponential density factor, width/offset subtractions (2), two divisions and write.
+
+CDF's regular branch is three interval log masses (762), two subtractions, and 36 for the final comparison, both eager exp/expm1 values, negation, selection and write. Its narrow branch is two integrals (370), two offset subtractions, division and write.
+
+PPF's regular branch is log(q)/log1p(-q) (33), four log CDFs (488), two logaddexp combinations including their additions (36), side comparison (1), the one selected inverse branch including optional sign and write (782), two full-subset inversions of the right-side mask (2), and output write (1). Its narrow branch has width and initial-point arithmetic (2), one integral (185), four steps of 212, then final offset addition/write (2). A narrow Newton step has integral (185), multiply/subtract (2), exponential derivative (21), division/subtraction (2), and clip (2).
+
+The common costs cover input-domain predicates (12), method masks and endpoint writes, and applicable normalization/clipping/affine operations. Empty inputs retain the existing wrapper minimum of one billed element. The total charges are therefore 315, 844, and 1392 times `max(numel(broadcast(x,a,b,loc,scale)),1)`, subsequently multiplied by the configured float64 dtype rate and stats weight. With packaged defaults those factors are 2 and 1.
+
+## Numerical scope
+
+The implementation retains normal tail probabilities in logarithmic form and
+uses a factored local integral for nearly coincident bounds. Regression checks
+cover positive and negative 9- and 40-sigma intervals, semi-infinite intervals,
+central intervals, and adjacent representable bounds. The log-CDF helper is
+checked on a grid through absolute standardized arguments of 100. Narrow-interval
+checks use independent adaptive quadrature because subtracting two rounded CDF
+values is not an adequate reference there.
+
+This is not an arbitrary-finite-double accuracy guarantee. At extremely large
+standardized bounds, squaring and subtracting log probabilities can still lose
+precision or overflow; an interval near 1e10 can still produce a NaN quantile.
+For very narrow intervals, the quantile's representable spacing can itself
+prevent exact inversion. Correct domain handling returns NaN for invalid
+probabilities, nonpositive scales and invalid bounds. The numerical cost stays
+fixed across these cases, so input values cannot select a cheaper billed path.

flopscope-server/tests/test_dispatch_truncnorm_tails.py

@@ -0,0 +1,41 @@
+"""The existing stats dispatch must retain stable tail results and metadata."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+from flopscope_server._request_handler import RequestHandler
+from flopscope_server._session import Session
+from scipy import stats
+
+
+@pytest.mark.parametrize("method", ["pdf", "cdf", "ppf"])
+def test_tail_stats_dispatch(method):
+    values = (
+        np.array([0.1, 0.5, 0.9]) if method == "ppf" else np.array([40.0, 40.01, 40.1])
+    )
+    session = Session(flop_budget=10**9)
+    try:
+        handler = RequestHandler(session)
+        handle = session.store_array(values)
+        before = session.budget_remaining
+        response = handler.handle(
+            {
+                "op": "stats.truncnorm." + method,
+                "args": [handle, 40.0, 41.0],
+                "kwargs": {},
+            }
+        )
+        assert response["status"] == "ok", response
+        metadata = response["result"]
+        result = np.asarray(session.get_array(metadata["id"]))
+        assert result.shape == (3,) and result.dtype == np.float64
+        assert session.budget_remaining < before
+        np.testing.assert_allclose(
+            result,
+            getattr(stats.truncnorm, method)(values, 40, 41),
+            rtol=5e-13,
+            atol=5e-13,
+        )
+    finally:
+        session.close()

scripts/generate_api_docs.py

@@ -321,9 +321,14 @@ def load_registry() -> dict[str, dict]:
             | `logistic` | $8n$ | $5n$ | $5n$ |
             | `laplace` | $5n$ | $5n$ | $5n$ |
             | `lognorm` | $15n$ | $25n$ | $45n$ |
-            | `truncnorm` | $30n$ | $30n$ | $50n$ |
+            | `truncnorm` | $315n$ | $844n$ | $1392n$ |
 
             where $n$ = `numel(x)` (or `numel(q)` for ppf).
+            For `truncnorm`, $n$ is `max(numel(broadcast(x, a, b, loc, scale)), 1)`;
+            the displayed constants are numerical upper bounds at weight 1.
+            The packaged float64 dtype rate doubles these base counts. The
+            kernel uses log-domain tails and fixed narrow-interval quadrature;
+            arbitrary extreme finite bounds are not an accuracy guarantee.
 
             ## Examples
 
@@ -633,9 +638,9 @@ def generate_api_page(page_path: str, page_info: dict) -> None:
     "stats.lognorm.pdf": ("15n", r"$15n$"),
     "stats.lognorm.cdf": ("25n", r"$25n$"),
     "stats.lognorm.ppf": ("45n", r"$45n$"),
-    "stats.truncnorm.pdf": ("30n", r"$30n$"),
-    "stats.truncnorm.cdf": ("30n", r"$30n$"),
-    "stats.truncnorm.ppf": ("50n", r"$50n$"),
+    "stats.truncnorm.pdf": ("315n", r"$315n$"),
+    "stats.truncnorm.cdf": ("844n", r"$844n$"),
+    "stats.truncnorm.ppf": ("1392n", r"$1392n$"),
 }
 
 CATEGORY_COST_LATEX: dict[str, tuple[str, str]] = {

src/flopscope/stats/_truncnorm.py

@@ -2,24 +2,15 @@
 
 from __future__ import annotations
 
-import numpy as _np
-
+from flopscope.stats import _truncnorm_kernels
 from flopscope.stats._base import ContinuousDistribution
-from flopscope.stats._erf import _erf
-from flopscope.stats._ndtri import _ndtri
-
-_SQRT2 = _np.sqrt(2.0)
-_INV_SQRT_2PI = 1.0 / _np.sqrt(2.0 * _np.pi)
-
-
-def _std_norm_cdf(x):
-    """Standard normal CDF (no budget deduction)."""
-    return 0.5 * (1.0 + _erf(x / _SQRT2))
-
 
-def _std_norm_pdf(x):
-    """Standard normal PDF (no budget deduction)."""
-    return _INV_SQRT_2PI * _np.exp(-0.5 * x * x)
+# Fixed analytical bounds for the numerical kernels, before dtype pricing.
+# See docs/reference/cost-model.md. Four Newton steps and eight quadrature nodes
+# are part of this derivation; changing them requires revisiting these bounds.
+_PDF_COST = 315
+_CDF_COST = 844
+_PPF_COST = 1392
 
 
 class TruncnormDistribution(ContinuousDistribution):
@@ -40,14 +31,17 @@ class TruncnormDistribution(ContinuousDistribution):
     -----
     ``a`` and ``b`` are standardized lower and upper bounds. The truncated
     support is ``[a * scale + loc, b * scale + loc]``, and both bounds appear
-    before ``loc`` and ``scale`` to match SciPy's signature. pdf deducts
-    ``28 * numel(broadcast(input, a, b, loc, scale))`` FLOPs (composite:
-    z(2)+std_norm_pdf(20)+div(1)+bounds(5), FMA=2, weight 1.0; calibrated
-    alpha 28.0). cdf deducts ``51 * numel(broadcast(input, a, b, loc, scale))``
-    FLOPs (composite: z(2)+std_norm_cdf(46)+result(3)+2 where(4), FMA=2,
-    weight 1.0; calibrated alpha 50.6). ppf deducts
-    ``81 * numel(broadcast(input, a, b, loc, scale))`` FLOPs (composite:
-    erf + ndtri rational approx + arithmetic, weight 1.0; audit-2 verified).
+    before ``loc`` and ``scale`` to match SciPy's signature. Invalid bounds
+    (``a >= b``), nonpositive scales, and NaN inputs return NaN; ppf also
+    returns NaN for probabilities outside ``[0, 1]``.
+
+    Base composite costs are ``315*n`` for pdf, ``844*n`` for cdf, and
+    ``1392*n`` for ppf, where
+    ``n = max(numel(broadcast(input, a, b, loc, scale)), 1)``. These are
+    analytical numerical upper bounds (FMA=2, stats weight 1.0), not hardware
+    calibrations. Configured dtype pricing applies to the float64 output.
+    The kernels use stable log-tail probabilities, eight fixed quadrature
+    nodes for narrow intervals, and four fixed inverse-refinement steps.
     """
 
     def __init__(self):
@@ -78,9 +72,9 @@ def pdf(self, x, a, b, loc=0, scale=1):
         Notes
         -----
         Equivalent to ``scipy.stats.truncnorm.pdf(x, a, b, loc, scale)``.
-        FLOP cost: ``28 * numel(broadcast(x, a, b, loc, scale))`` (composite:
-        z(2)+std_norm_pdf(20)+div(1)+bounds(5), FMA=2, weight 1.0; calibrated
-        alpha 28.0).
+        Base FLOP cost: ``315 * max(numel(broadcast(x, a, b, loc, scale)), 1)``.
+        Analytical numerical upper bound, FMA=2, weight 1.0, before the
+        configured float64 dtype multiplier.
 
         Examples
         --------
@@ -90,7 +84,7 @@ def pdf(self, x, a, b, loc=0, scale=1):
         >>> np.round(flops.stats.truncnorm.pdf(x, a=-1.0, b=1.0), 3)
         array([0.516, 0.584, 0.516])
         """
-        return self._deduct_and_call("pdf", 28, x, a, b, loc=loc, scale=scale)
+        return self._deduct_and_call("pdf", _PDF_COST, x, a, b, loc=loc, scale=scale)
 
     def cdf(self, x, a, b, loc=0, scale=1):
         """Evaluate the cumulative distribution function.
@@ -117,9 +111,9 @@ def cdf(self, x, a, b, loc=0, scale=1):
         Notes
         -----
         Equivalent to ``scipy.stats.truncnorm.cdf(x, a, b, loc, scale)``.
-        FLOP cost: ``51 * numel(broadcast(x, a, b, loc, scale))`` (composite:
-        z(2)+std_norm_cdf(46)+result(3)+2 where(4), FMA=2, weight 1.0;
-        calibrated alpha 50.6).
+        Base FLOP cost: ``844 * max(numel(broadcast(x, a, b, loc, scale)), 1)``.
+        Analytical numerical upper bound, FMA=2, weight 1.0, before the
+        configured float64 dtype multiplier.
 
         Examples
         --------
@@ -129,7 +123,7 @@ def cdf(self, x, a, b, loc=0, scale=1):
         >>> np.round(flops.stats.truncnorm.cdf(x, a=-1.0, b=1.0), 3)
         array([0.22, 0.5 , 0.78])
         """
-        return self._deduct_and_call("cdf", 51, x, a, b, loc=loc, scale=scale)
+        return self._deduct_and_call("cdf", _CDF_COST, x, a, b, loc=loc, scale=scale)
 
     def ppf(self, q, a, b, loc=0, scale=1):
         """Evaluate the percent-point function.
@@ -156,8 +150,9 @@ def ppf(self, q, a, b, loc=0, scale=1):
         Notes
         -----
         Equivalent to ``scipy.stats.truncnorm.ppf(q, a, b, loc, scale)``.
-        FLOP cost: ``81 * numel(broadcast(q, a, b, loc, scale))`` (composite:
-        erf + ndtri rational approx + arithmetic, weight 1.0).
+        Base FLOP cost: ``1392 * max(numel(broadcast(q, a, b, loc, scale)), 1)``.
+        Analytical numerical upper bound, FMA=2, weight 1.0, before the
+        configured float64 dtype multiplier.
 
         Examples
         --------
@@ -167,31 +162,16 @@ def ppf(self, q, a, b, loc=0, scale=1):
         >>> np.round(flops.stats.truncnorm.ppf(q, a=-1.0, b=1.0), 3)
         array([-0.442,  0.   ,  0.442])
         """
-        return self._deduct_and_call("ppf", 81, q, a, b, loc=loc, scale=scale)
+        return self._deduct_and_call("ppf", _PPF_COST, q, a, b, loc=loc, scale=scale)
 
     def _compute_pdf(self, x, a, b, loc=0, scale=1):
-        z = (x - loc) / scale
-        phi_a = _std_norm_cdf(a)
-        phi_b = _std_norm_cdf(b)
-        denom = scale * (phi_b - phi_a)
-        result = _std_norm_pdf(z) / denom
-        return _np.where((z >= a) & (z <= b), result, 0.0)
+        return _truncnorm_kernels.pdf(x, a, b, loc=loc, scale=scale)
 
     def _compute_cdf(self, x, a, b, loc=0, scale=1):
-        z = (x - loc) / scale
-        phi_a = _std_norm_cdf(a)
-        phi_b = _std_norm_cdf(b)
-        result = (_std_norm_cdf(z) - phi_a) / (phi_b - phi_a)
-        result = _np.where(z < a, 0.0, result)
-        result = _np.where(z > b, 1.0, result)
-        return result
+        return _truncnorm_kernels.cdf(x, a, b, loc=loc, scale=scale)
 
     def _compute_ppf(self, q, a, b, loc=0, scale=1):
-        phi_a = _std_norm_cdf(a)
-        phi_b = _std_norm_cdf(b)
-        inner = phi_a + q * (phi_b - phi_a)
-        z = _ndtri(inner)
-        return loc + scale * z
+        return _truncnorm_kernels.ppf(q, a, b, loc=loc, scale=scale)
 
 
 truncnorm = TruncnormDistribution()