alkahest-cas 3.12.0

//! Probability distributions as symbolic objects, and expectations over them.
//!
//! A [`Distribution`] is a *name plus symbolic parameters*. It carries its own
//! density, its support, and the parameter constraints that make it a
//! distribution at all (`σ > 0`, `a < b`, `0 ≤ p ≤ 1`). Everything derived from
//! it — [`Distribution::mean`], [`Distribution::variance`],
//! [`Distribution::moment`], [`Distribution::cdf`], [`Distribution::quantile`],
//! and the general [`expectation`] — is **checked against a numerical
//! quadrature of its own defining integral before it is returned**. A closed
//! form the checker cannot confirm is withheld ([`ProbError::Unverified`],
//! `E-PROB-005`), never returned with a caveat.
//!
//! # Why a separate module rather than a wider `integrate`
//!
//! An expectation is not an integral the user wrote down; it is an integral
//! this module *constructs* from a distribution, and the construction is where
//! the mistakes live — the wrong Jacobian, the wrong support, the wrong
//! normalising constant, a parameter outside its constraint. Those are exactly
//! the mistakes a symbolic integrator cannot catch, because by the time it sees
//! the integrand the error is already baked in and it will integrate the wrong
//! thing perfectly. So the construction is checked *end to end*, against the
//! definition, numerically, every time.
//!
//! # How an expectation is computed
//!
//! `E[f(X)] = ∫_S f(x) p(x) dx` is real but rarely tractable in `x`. Each
//! continuous distribution therefore declares a **reduction**: an increasing
//! change of variables `x = ξ(z)` under which the density becomes a standard
//! kernel the existing integrator has a chance with
//! (`e^{-z²/2}` for `Normal`/`LogNormal`, `e^{-u}` for `Exponential`/`Gamma`,
//! the identity for `Uniform`/`Beta`). The integrand is pushed through the
//! reduction, expanded into a sum, each term's `z`-free factor is pulled out,
//! and what is left is handed to [`crate::integrate::integrate_definite`].
//!
//! That pipeline can go wrong in four places, so it is **verified against the
//! untransformed definition**: the numeric gate quadratures `f(x) p(x)` over
//! the original support in `x`, not over `z`. A wrong Jacobian, a wrong
//! inverse, a dropped constant or a mis-mapped bound all show up as a
//! disagreement and the answer is refused. See the `verify` module (private;
//! `--document-private-items` to read it) and the [`Evidence`] it returns.
//!
//! # What refuses, and why that is right
//!
//! | input | refusal |
//! |---|---|
//! | `Normal(μ, 0)`, `Uniform(1, 0)`, `Beta(-1, 2)` | `E-PROB-001`: the parameter is a *number* and it violates the constraint. There is no distribution to ask about |
//! | `Gamma(k, θ).cdf(x)` for symbolic `k` | `E-PROB-004`: the lower incomplete gamma function is not in this library's primitive set, so there is no closed form to return. Integer `k` (Erlang) does close and is returned |
//! | `Normal(μ, σ).quantile(p)` | `E-PROB-004`: needs `erf⁻¹`, which is not a registered primitive |
//! | `E[e^{X²}]` under `Normal(0, 1)` | `E-PROB-006`: the reduction integral `∫ e^{z²} e^{-z²/2} dz` **diverges**. The gate sees the quadrature fail to converge and says so rather than returning the integrator's formal answer |
//! | an integrand the integrator declines | `E-PROB-003`, naming the exact integral that did not close. Not an unevaluated `Integral` object dressed as an answer |
//! | a closed form the quadrature contradicts | `E-PROB-005`. This is the prime-directive case: something *was* computed and is being withheld |
//! | `Uniform(a, b).quantile(2)` | `E-PROB-001`: `p` is a probability. The closed form returns `a + 2(b-a)` there — a point outside the support, offered as a quantile |
//!
//! # What this module deliberately does not do
//!
//! There is no joint distribution, no conditioning, and no dependence
//! structure. [`expectation_affine`] applies linearity — `E[aX + bY] = aE[X] +
//! bE[Y]`, which needs no independence at all — and
//! [`variance_affine_independent`] applies the variance rule under an
//! independence assumption the caller states by choosing that function's name.
//! Anything beyond those two would require a joint-law design, and guessing at
//! one produces wrong covariances silently.

use std::fmt;

use crate::deriv::SideCondition;
use crate::errors::AlkahestError;
use crate::kernel::{Domain, ExprId, ExprPool};

mod charfun;
mod cplx;
mod dists;
mod entropy;
mod expect;
mod genfun;
mod moments;
mod quad;
#[cfg(test)]
mod tests;
mod verify;

pub use charfun::characteristic_function;
pub use entropy::{
    cross_entropy, entropy, kl_divergence, mutual_information_independent, MAX_BINOMIAL_ENTROPY_N,
};
pub use expect::{expectation, expectation_affine, variance_affine_independent};
pub use genfun::{
    cumulant, cumulant_generating_function, excess_kurtosis, factorial_moment,
    moment_generating_function, probability_generating_function, skewness, MAX_CUMULANT_ORDER,
};
pub use verify::Evidence;

pub(crate) use quad::numeric_ball;

// ---------------------------------------------------------------------------
// Undischarged hypotheses
// ---------------------------------------------------------------------------

std::thread_local! {
    static PROB_CONDITIONS: std::cell::RefCell<Vec<SideCondition>> =
        const { std::cell::RefCell::new(Vec::new()) };
}

/// The hypotheses the most recent [`Distribution::cdf`],
/// [`Distribution::quantile`] or [`expectation`] call on this thread had to
/// **assume** in order to return the answer it did.
///
/// The same out-of-band channel as
/// [`crate::transform::take_transform_side_conditions`], and for the same
/// reason: the value is an `ExprId`, so there is nowhere in band to hang a
/// hypothesis, and a caller who cannot see one cannot tell a theorem from a
/// conditional.
///
/// Non-empty in exactly the cases the argument could not be *placed*:
///
/// * `cdf(x)` for a symbolic `x` — the closed forms are the in-support branch,
///   and `P(X ≤ x)` is `0`/`1` outside it. A decidable argument is answered
///   exactly instead and records nothing; a `Normal`, whose support is the
///   whole line, never records anything at all.
/// * `quantile(p)` for a symbolic `p` — `F⁻¹` is defined on `[0, 1]`, and off
///   it the closed forms return a number rather than failing.
/// * `expectation` of a payoff whose kink position cannot be placed inside the
///   support — `E[max(S - K, 0)]` with a symbolic strike needs `K > 0`.
///
/// Consuming, so one call's hypotheses cannot be read as a later call's, and
/// **reset by every one of those three entry points** — including on the
/// refusal path, so a stale list is never attributed to a call that failed.
/// An empty list means every branch taken was forced by the input, not that
/// none was taken.
pub fn take_prob_side_conditions() -> Vec<SideCondition> {
    PROB_CONDITIONS.with(|c| std::mem::take(&mut *c.borrow_mut()))
}

/// Replace the channel with `conds`. Called at the *entry* of every route that
/// can record one, so a refusal clears rather than leaks.
pub(crate) fn stash_prob_side_conditions(conds: Vec<SideCondition>) {
    PROB_CONDITIONS.with(|c| *c.borrow_mut() = conds);
}

/// Collect every side condition a derivation recorded, in order, without
/// repeats.
///
/// A split payoff attaches the *same* hypothesis to each of its pieces —
/// `E[max(S - K, 0)]` cuts the support once and both sides carry `K > 0` — so
/// without this the caller reads `["K > 0", "K > 0"]` and has to work out
/// whether that is one hypothesis or two. `ExprId` is interned, so equality is
/// a pointer comparison.
pub(crate) fn conditions_of(log: &crate::deriv::DerivationLog) -> Vec<SideCondition> {
    let mut out: Vec<SideCondition> = Vec::new();
    for c in log.0.iter().flat_map(|s| s.side_conditions.iter()) {
        if !out.contains(c) {
            out.push(c.clone());
        }
    }
    // `e > 0` subsumes `e ≥ 0`, and a route can record both: the entropy of a
    // `Bernoulli(p)` needs the strict `p > 0` (the closed form is `NaN` at the
    // endpoint) on top of the constructor's `p ≥ 0`. Publishing the weaker one
    // alongside the stronger is not wrong, it is noise — and a caller reading
    // four hypotheses where there are two stops reading them.
    let strict: Vec<ExprId> = out
        .iter()
        .filter_map(|c| match c {
            SideCondition::Positive(e) => Some(*e),
            _ => None,
        })
        .collect();
    out.retain(|c| {
        !matches!(c,
        SideCondition::InDomain(e, Domain::NonNegative) if strict.contains(e))
    });
    out
}

// ---------------------------------------------------------------------------
// Errors
// ---------------------------------------------------------------------------

/// Why a probability query did not return a value.
///
/// `#[non_exhaustive]` from birth: a future refusal reason has to be an
/// additive change rather than a semver break.
#[derive(Clone, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum ProbError {
    /// A parameter evaluates to a number that violates the distribution's own
    /// constraint — `σ ≤ 0`, `a ≥ b`, `p ∉ [0, 1]`, a non-integer `n` for
    /// `Binomial`. `E-PROB-001`.
    ///
    /// A *symbolic* parameter is never reported here: it cannot be decided, so
    /// it is carried as a side condition on
    /// [`Distribution::constraints`] instead.
    InvalidParameter {
        /// Which parameter, by name (`"sigma"`, `"a"`, `"p"`, …).
        parameter: &'static str,
        /// The constraint that failed, spelled out.
        requirement: &'static str,
    },
    /// The query is outside this module's supported shape — a general `f` over
    /// a distribution with infinite discrete support, a non-affine argument to
    /// [`expectation_affine`], a moment order this table does not carry.
    /// `E-PROB-002`.
    Unsupported(String),
    /// The reduction integral was built and [`crate::integrate::integrate`]
    /// declined it. The message names the integral. `E-PROB-003`.
    ///
    /// Distinct from [`ProbError::NoClosedForm`] on purpose: this one says
    /// *this* integrator could not do it, which may change; that one says the
    /// closed form needs a function this library does not have.
    IntegralDidNotClose(String),
    /// The quantity has no closed form within the primitive set — the `Gamma`
    /// CDF needs `γ(k, x)`, the `Beta` CDF needs `I_x(α, β)`, the `Normal`
    /// quantile needs `erf⁻¹`. `E-PROB-004`.
    NoClosedForm {
        /// The quantity asked for.
        quantity: &'static str,
        /// The function that would be needed to express it.
        missing: &'static str,
    },
    /// A closed form **was** computed and then **withheld**, because the
    /// numeric gate could not confirm it. `E-PROB-005`.
    ///
    /// Never downgraded to a warning: a caller cannot tell a checked value
    /// from an unchecked one once it is in their hands.
    Unverified(UnverifiedReason),
    /// The defining integral (or series) does not converge, so the quantity
    /// does not exist. `E-PROB-006`.
    ///
    /// The formal symbolic answer for a divergent expectation is typically a
    /// clean, plausible, wrong number — which is the whole reason the gate
    /// runs before the answer is returned rather than after.
    Divergent(String),
}

/// Why the numeric gate declined to confirm a computed closed form.
#[derive(Clone, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum UnverifiedReason {
    /// The closed form and the quadrature of the defining integral disagree by
    /// more than the quadrature's own error estimate allows. The strings are
    /// the two values and the parameter point at which they parted.
    Disagreement(String),
    /// The closed form could not be evaluated to a number at any admissible
    /// parameter point — an unregistered head, or a branch that is not real
    /// there. An unevaluable "answer" reports success and then turns into
    /// `NaN`, which is the failure mode this library exists to avoid.
    ClaimNotEvaluable,
    /// The quadrature did not converge to the precision the comparison needs,
    /// so it cannot *contradict* the closed form either. Silence is not
    /// evidence: the form is withheld rather than passed on a check that could
    /// not fail.
    QuadratureInconclusive,
    /// No admissible parameter point could be constructed — the constraints
    /// involve symbols this module cannot sample.
    NoAdmissiblePoint,
}

impl fmt::Display for ProbError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            ProbError::InvalidParameter {
                parameter,
                requirement,
            } => write!(
                f,
                "prob: parameter `{parameter}` violates the constraint {requirement}"
            ),
            ProbError::Unsupported(msg) => write!(f, "prob: unsupported: {msg}"),
            ProbError::IntegralDidNotClose(msg) => {
                write!(f, "prob: the reduction integral did not close: {msg}")
            }
            ProbError::NoClosedForm { quantity, missing } => write!(
                f,
                "prob: no closed form for {quantity} — it requires {missing}, \
                 which is not a registered primitive"
            ),
            ProbError::Unverified(reason) => match reason {
                UnverifiedReason::Disagreement(detail) => write!(
                    f,
                    "prob: a closed form was computed and withheld — it disagrees \
                     with quadrature of its own defining integral ({detail})"
                ),
                UnverifiedReason::ClaimNotEvaluable => write!(
                    f,
                    "prob: a closed form was computed and withheld — it could not be \
                     evaluated to a number at any admissible parameter point"
                ),
                UnverifiedReason::QuadratureInconclusive => write!(
                    f,
                    "prob: a closed form was computed and withheld — the verifying \
                     quadrature did not converge, so the check could not have failed"
                ),
                UnverifiedReason::NoAdmissiblePoint => write!(
                    f,
                    "prob: a closed form was computed and withheld — no admissible \
                     parameter point could be constructed to check it at"
                ),
            },
            ProbError::Divergent(msg) => write!(f, "prob: the quantity does not exist: {msg}"),
        }
    }
}

impl std::error::Error for ProbError {}

impl AlkahestError for ProbError {
    fn code(&self) -> &'static str {
        match self {
            ProbError::InvalidParameter { .. } => "E-PROB-001",
            ProbError::Unsupported(_) => "E-PROB-002",
            ProbError::IntegralDidNotClose(_) => "E-PROB-003",
            ProbError::NoClosedForm { .. } => "E-PROB-004",
            ProbError::Unverified(_) => "E-PROB-005",
            ProbError::Divergent(_) => "E-PROB-006",
        }
    }

    fn remediation(&self) -> Option<&'static str> {
        match self {
            ProbError::InvalidParameter { .. } => Some(
                "pass a parameter inside the distribution's constraint, or leave it symbolic \
                 (a symbolic parameter is carried as a side condition on `constraints` \
                 rather than decided)",
            ),
            ProbError::Unsupported(_) => Some(
                "reduce the query to a supported shape: a polynomial or exponential `f` for \
                 `expectation`, an affine combination for `expectation_affine`",
            ),
            ProbError::IntegralDidNotClose(_) => Some(
                "the named integral has no antiderivative this integrator found; integrate it \
                 numerically, or ask for a moment/CDF that is in the closed-form table",
            ),
            ProbError::NoClosedForm { .. } => Some(
                "the closed form needs a special function this library does not implement; \
                 use an integer shape parameter where that collapses the special function to \
                 a finite sum (Erlang, integer-parameter Beta), or integrate numerically",
            ),
            ProbError::Unverified(_) => Some(
                "an unconfirmed closed form is withheld rather than returned; report the \
                 distribution and query as a minimal failing example",
            ),
            ProbError::Divergent(_) => Some(
                "the integral defining this expectation does not converge — there is no value \
                 to return; restrict `f` so that `f(x)·p(x)` is integrable over the support",
            ),
        }
    }
}

// ---------------------------------------------------------------------------
// Support
// ---------------------------------------------------------------------------

/// Where a distribution puts its mass.
#[derive(Clone, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum Support {
    /// The whole real line.
    Real,
    /// `(0, ∞)`.
    Positive,
    /// The closed interval `[lo, hi]`, endpoints symbolic.
    Interval(ExprId, ExprId),
    /// `{0, 1, …, n}`.
    IntegersUpTo(u32),
    /// `{0, 1, 2, …}` — countably infinite.
    NonNegativeIntegers,
}

impl Support {
    /// Is this a discrete support?
    pub fn is_discrete(&self) -> bool {
        matches!(
            self,
            Support::IntegersUpTo(_) | Support::NonNegativeIntegers
        )
    }
}

// ---------------------------------------------------------------------------
// Distribution
// ---------------------------------------------------------------------------

/// Which distribution. Parameters live on [`Distribution`], not here.
///
/// `#[non_exhaustive]`: adding a distribution must not be a semver break.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
#[non_exhaustive]
pub enum DistKind {
    /// `Normal(μ, σ)`, `σ > 0`, density `e^{-(x-μ)²/(2σ²)} / (σ√(2π))` on `ℝ`.
    Normal,
    /// `LogNormal(μ, σ)`, `σ > 0`: `log X ~ Normal(μ, σ)`. Density
    /// `e^{-(log x - μ)²/(2σ²)} / (x σ√(2π))` on `(0, ∞)`.
    LogNormal,
    /// `Uniform(a, b)`, `a < b`, density `1/(b-a)` on `[a, b]`.
    Uniform,
    /// `Exponential(λ)`, `λ > 0`, density `λ e^{-λx}` on `(0, ∞)`.
    Exponential,
    /// `Gamma(k, θ)`, `k > 0`, `θ > 0` — **shape–scale**, density
    /// `x^{k-1} e^{-x/θ} / (Γ(k) θ^k)` on `(0, ∞)`.
    Gamma,
    /// `Beta(α, β)`, `α > 0`, `β > 0`, density
    /// `x^{α-1}(1-x)^{β-1} / B(α, β)` on `[0, 1]`.
    Beta,
    /// `Bernoulli(p)`, `0 ≤ p ≤ 1`, on `{0, 1}`.
    Bernoulli,
    /// `Binomial(n, p)`, `n` a non-negative integer, `0 ≤ p ≤ 1`, on `{0…n}`.
    Binomial,
    /// `Poisson(λ)`, `λ > 0`, on `{0, 1, 2, …}`.
    Poisson,
}

impl DistKind {
    /// The distribution's name, as it is printed and as it appears in
    /// derivation logs.
    pub fn name(self) -> &'static str {
        match self {
            DistKind::Normal => "Normal",
            DistKind::LogNormal => "LogNormal",
            DistKind::Uniform => "Uniform",
            DistKind::Exponential => "Exponential",
            DistKind::Gamma => "Gamma",
            DistKind::Beta => "Beta",
            DistKind::Bernoulli => "Bernoulli",
            DistKind::Binomial => "Binomial",
            DistKind::Poisson => "Poisson",
        }
    }
}

/// A distribution with symbolic parameters.
///
/// Constructed through the named constructors ([`Distribution::normal`],
/// [`Distribution::gamma`], …), which is what makes it impossible to hold one
/// whose *numeric* parameters violate its own constraints. Symbolic parameters
/// are not decidable and are carried on [`Distribution::constraints`] instead.
///
/// Fields are private: a publicly-constructible struct cannot gain a field
/// without a semver break, and this one will gain fields.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Distribution {
    kind: DistKind,
    params: Vec<ExprId>,
}

impl Distribution {
    /// Which distribution this is.
    pub fn kind(&self) -> DistKind {
        self.kind
    }

    /// The parameters, in the order of the constructor.
    pub fn params(&self) -> &[ExprId] {
        &self.params
    }

    /// `Normal(μ, σ)` — requires `σ > 0`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `σ` is a number that is not
    /// positive. A symbolic `σ` is accepted and carried on
    /// [`Distribution::constraints`].
    pub fn normal(mu: ExprId, sigma: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_positive(sigma, "sigma", pool)?;
        Ok(Distribution {
            kind: DistKind::Normal,
            params: vec![mu, sigma],
        })
    }

    /// `LogNormal(μ, σ)` — `log X ~ Normal(μ, σ)`, requires `σ > 0`.
    ///
    /// Note the convention: `μ` and `σ` are the parameters of the *underlying
    /// normal*, not the mean and standard deviation of `X`. `E[X] =
    /// e^{μ+σ²/2}`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `σ` is a non-positive number.
    pub fn log_normal(mu: ExprId, sigma: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_positive(sigma, "sigma", pool)?;
        Ok(Distribution {
            kind: DistKind::LogNormal,
            params: vec![mu, sigma],
        })
    }

    /// `Uniform(a, b)` — requires `a < b`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `a` and `b` are both numbers and
    /// `a ≥ b`.
    pub fn uniform(a: ExprId, b: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_less(a, b, "b", "a < b", pool)?;
        Ok(Distribution {
            kind: DistKind::Uniform,
            params: vec![a, b],
        })
    }

    /// `Exponential(λ)` — rate parametrisation, mean `1/λ`, requires `λ > 0`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `λ` is a non-positive number.
    pub fn exponential(lambda: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_positive(lambda, "lambda", pool)?;
        Ok(Distribution {
            kind: DistKind::Exponential,
            params: vec![lambda],
        })
    }

    /// `Gamma(k, θ)` — **shape–scale**, mean `kθ`, requires `k > 0`, `θ > 0`.
    ///
    /// The other common convention is shape–*rate*; pass `1/rate` as `θ`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when either parameter is a non-positive
    /// number.
    pub fn gamma(k: ExprId, theta: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_positive(k, "k", pool)?;
        dists::require_positive(theta, "theta", pool)?;
        Ok(Distribution {
            kind: DistKind::Gamma,
            params: vec![k, theta],
        })
    }

    /// `Beta(α, β)` — requires `α > 0`, `β > 0`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when either parameter is a non-positive
    /// number.
    pub fn beta(alpha: ExprId, beta: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_positive(alpha, "alpha", pool)?;
        dists::require_positive(beta, "beta", pool)?;
        Ok(Distribution {
            kind: DistKind::Beta,
            params: vec![alpha, beta],
        })
    }

    /// `Bernoulli(p)` — requires `0 ≤ p ≤ 1`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `p` is a number outside `[0, 1]`.
    pub fn bernoulli(p: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_probability(p, "p", pool)?;
        Ok(Distribution {
            kind: DistKind::Bernoulli,
            params: vec![p],
        })
    }

    /// `Binomial(n, p)` — `n` trials, requires `0 ≤ p ≤ 1`.
    ///
    /// `n` must be a literal non-negative integer: the support `{0…n}` is
    /// enumerated by every route in this module, so a symbolic `n` is not a
    /// distribution this module can compute with and is refused rather than
    /// half-supported.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `n` is not a literal non-negative
    /// integer, or `p` is a number outside `[0, 1]`.
    pub fn binomial(n: ExprId, p: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_count(n, "n", pool)?;
        dists::require_probability(p, "p", pool)?;
        Ok(Distribution {
            kind: DistKind::Binomial,
            params: vec![n, p],
        })
    }

    /// `Poisson(λ)` — requires `λ > 0`.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] when `λ` is a non-positive number.
    pub fn poisson(lambda: ExprId, pool: &ExprPool) -> Result<Self, ProbError> {
        dists::require_positive(lambda, "lambda", pool)?;
        Ok(Distribution {
            kind: DistKind::Poisson,
            params: vec![lambda],
        })
    }

    /// Where the mass lives.
    pub fn support(&self, pool: &ExprPool) -> Support {
        dists::support(self, pool)
    }

    /// The parameter constraints, as predicate expressions.
    ///
    /// A constraint whose parameters are numeric has already been *checked* by
    /// the constructor, so it appears here as a record rather than an
    /// obligation. A constraint over a symbolic parameter has **not** been
    /// decided and is the caller's to discharge — that is the "carried as a
    /// side condition" half of the contract, and it is why this returns
    /// predicates rather than a bool.
    pub fn constraints(&self, pool: &ExprPool) -> Vec<ExprId> {
        dists::constraints(self, pool)
    }

    /// The density at `x` (continuous) or the mass at `k` (discrete).
    ///
    /// For a continuous distribution this is the density *on the support only*
    /// — it is not extended by zero outside it, because a `Piecewise` density
    /// defeats every downstream symbolic route for no gain. Use
    /// [`Distribution::support`] to know where it is valid.
    pub fn pdf(&self, x: ExprId, pool: &ExprPool) -> ExprId {
        dists::pdf(self, x, pool)
    }

    /// `E[X]`, in closed form, verified.
    ///
    /// # Errors
    ///
    /// See [`ProbError`]; in particular [`ProbError::Unverified`] if the
    /// numeric gate could not confirm the table value against quadrature.
    pub fn mean(&self, pool: &ExprPool) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        moments::mean(self, pool)
    }

    /// `Var[X] = E[X²] - E[X]²`, in closed form, verified.
    ///
    /// # Errors
    ///
    /// See [`ProbError`].
    pub fn variance(
        &self,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        moments::variance(self, pool)
    }

    /// The raw moment `E[Xⁿ]`, in closed form, verified.
    ///
    /// # Errors
    ///
    /// See [`ProbError`]. `n` past [`MAX_MOMENT_ORDER`] is
    /// [`ProbError::Unsupported`]: the closed forms stay exact but the numeric
    /// gate loses its ability to *fail*, and an unfailable check is decoration.
    pub fn moment(
        &self,
        n: u32,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        moments::moment(self, n, pool)
    }

    /// `P(X ≤ x)`, in closed form, verified.
    ///
    /// # Outside the support
    ///
    /// Unlike [`Distribution::pdf`], this is **not** "the formula on the
    /// support, and you are on your own elsewhere". `P(X ≤ x)` is `0` below the
    /// support and `1` above it, and the in-support closed form does not merely
    /// fail there — it produces a clean wrong number offered as a probability
    /// (the `Erlang(3, 4/5)` form at `x = -5` is `-7396.87`). So:
    ///
    /// * an argument that can be *decided* to sit outside the support returns
    ///   the exact `0` or `1`;
    /// * an argument that cannot be decided returns the in-support branch with
    ///   the restriction recorded as a [`crate::deriv::SideCondition`] on the
    ///   `prob_cdf_argument_inside_support` step, rather than silently assumed.
    ///
    /// The second case is a weaker guarantee than the first, and it is the one
    /// a symbolic argument gets.
    ///
    /// # Errors
    ///
    /// [`ProbError::NoClosedForm`] where the CDF needs a special function this
    /// library does not have — the incomplete gamma for a non-integer `Gamma`
    /// shape, the incomplete beta for non-integer `Beta` parameters.
    pub fn cdf(
        &self,
        x: ExprId,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        moments::cdf(self, x, pool)
    }

    /// The characteristic function `φ_X(t) = E[e^{itX}]`, in closed form,
    /// verified — against the defining integral **and**, where
    /// [`crate::transform::fourier_transform`] has a rule for the density,
    /// against `F{p}(-t/2π)`.
    ///
    /// `φ` is a table entry keyed on the distribution, not a transform of
    /// [`Distribution::pdf`]: a law can have a closed-form `φ` and no
    /// closed-form density at all, and a design that derives one from the other
    /// cannot express that. See [`crate::prob::characteristic_function`] — in
    /// particular its note on **inversion**, of which there is none here.
    ///
    /// # Errors
    ///
    /// [`ProbError::NoClosedForm`] for `LogNormal` (whose `φ` has no closed
    /// form at all) and `Beta` (which needs `₁F₁`).
    pub fn characteristic_function(
        &self,
        t: ExprId,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        charfun::characteristic_function(self, t, pool)
    }

    /// The entropy — Shannon `H = -Σ p log p` on a discrete support,
    /// differential `h = -∫ f log f` on a continuous one — in nats unless
    /// `base` says otherwise, verified against its own defining integral.
    ///
    /// **The two are different quantities**, and this method returns whichever
    /// the support calls for rather than pretending they are one. `h` is not
    /// the limit of `H`, is not non-negative (`h(Uniform(0, ½)) = log ½`), and
    /// is not invariant under a change of variables. See
    /// [`crate::prob::entropy`] before treating a value from here as "the
    /// information content" of anything.
    ///
    /// # Errors
    ///
    /// [`ProbError::NoClosedForm`] for `Poisson`, whose entropy needs a closed
    /// form for `e^{-λ}Σ λ^k log(k!)/k!` that does not exist; see
    /// [`crate::prob::entropy`] for the rest.
    pub fn entropy(
        &self,
        base: Option<ExprId>,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        entropy::entropy(self, base, pool)
    }

    /// The quantile `F⁻¹(p)`, in closed form, verified.
    ///
    /// `p` is a probability and is checked as one: a numeric `p` outside
    /// `[0, 1]` is [`ProbError::InvalidParameter`], because the closed forms
    /// return a number there rather than failing — `Uniform(-2, 3)` at `p = 2`
    /// evaluates to `8`, outside the support it claims to be a point of. A
    /// symbolic `p` cannot be decided and carries `0 ≤ p ≤ 1` as a
    /// [`crate::deriv::SideCondition`] instead.
    ///
    /// # Errors
    ///
    /// [`ProbError::InvalidParameter`] for a numeric `p` outside `[0, 1]`, and
    /// [`ProbError::NoClosedForm`] for every distribution whose quantile needs
    /// `erf⁻¹` or a numerical inversion — which is most of them. `Uniform` and
    /// `Exponential` close; nothing else here does.
    pub fn quantile(
        &self,
        p: ExprId,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        moments::quantile(self, p, pool)
    }

    /// The moment generating function `M_X(t) = E[e^{tX}]`, in closed form,
    /// verified — **with its region of convergence reported rather than
    /// assumed**.
    ///
    /// This is the route where a table lookup and a correct answer part
    /// company. `M_X(t) = λ/(λ - t)` for an `Exponential(λ)` holds only on
    /// `t < λ`; outside that strip the same expression is finite, clean and
    /// not the value of any integral. So an argument that can be *decided* to
    /// sit outside the strip is [`ProbError::Divergent`], an argument that
    /// cannot carries `λ - t > 0` on
    /// [`take_prob_side_conditions`], and a law whose MGF is entire records
    /// that fact in the derivation log — an empty condition list means
    /// "checked", not "not looked at".
    ///
    /// # Errors
    ///
    /// [`ProbError::Divergent`] for `LogNormal` at a positive or undecidable
    /// `t` (`E[e^{tX}] = ∞` for every `t > 0`) and for an argument outside an
    /// `Exponential`/`Gamma` strip; [`ProbError::NoClosedForm`] for `Beta`,
    /// whose `M` is `₁F₁(α; α+β; t)`.
    pub fn moment_generating_function(
        &self,
        t: ExprId,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::moment_generating_function(self, t, pool)
    }

    /// The cumulant generating function `K_X(t) = log M_X(t)`, in closed form,
    /// verified as `e^{K(t)} = E[e^{tX}]`.
    ///
    /// Carries the same convergence strip as
    /// [`Distribution::moment_generating_function`], through the same channel:
    /// the logarithm of a divergent expectation is not a cumulant generating
    /// function.
    ///
    /// # Errors
    ///
    /// As [`Distribution::moment_generating_function`]. `Beta` refuses here
    /// while [`Distribution::cumulant`] does not — `K` is analytic at the
    /// origin for a Beta, it simply has no name in this library.
    pub fn cumulant_generating_function(
        &self,
        t: ExprId,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::cumulant_generating_function(self, t, pool)
    }

    /// The probability generating function `G_X(z) = E[z^X] = Σ_k z^k P(X=k)`,
    /// in closed form, verified against that sum.
    ///
    /// # Errors
    ///
    /// [`ProbError::Unsupported`] for every law not supported on the
    /// non-negative integers. `E[z^X]` for a `Normal` is a category error, not
    /// a harder integral: the formal rewrite `E[e^{X log z}]` returns a number,
    /// and that number is the MGF at `log z` and says nothing about any
    /// `P(X = k)`.
    pub fn probability_generating_function(
        &self,
        z: ExprId,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::probability_generating_function(self, z, pool)
    }

    /// The factorial moment `E[X(X-1)⋯(X-n+1)] = G_X^{(n)}(1)`, in closed
    /// form, verified against the defining sum.
    ///
    /// # Errors
    ///
    /// [`ProbError::Unsupported`] for a continuous law — there is no `G` to
    /// differentiate — and for `n` past [`MAX_MOMENT_ORDER`].
    pub fn factorial_moment(
        &self,
        n: u32,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::factorial_moment(self, n, pool)
    }

    /// The `n`-th cumulant `κ_n = K_X^{(n)}(0)`, in closed form, verified.
    ///
    /// `κ₁ = E[X]`, `κ₂ = Var[X]`, and for a `Normal` every `κ_n` with `n ≥ 3`
    /// is exactly `0` — the sharpest available test that the moment–cumulant
    /// recursion has no off-by-one in it, and one this module runs.
    ///
    /// # Errors
    ///
    /// [`ProbError::Divergent`] for `LogNormal`, whose `K` exists on no
    /// neighbourhood of the origin: the recursion runs, and what it computes
    /// is the coefficient of a divergent series. Use
    /// [`Distribution::skewness`] and [`Distribution::excess_kurtosis`], which
    /// are defined from central moments and do exist there.
    /// [`ProbError::Unsupported`] for `n = 0` and for `n` past
    /// [`MAX_CUMULANT_ORDER`].
    pub fn cumulant(
        &self,
        n: u32,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::cumulant(self, n, pool)
    }

    /// `γ₁ = E[((X - μ)/σ)³] = κ₃/σ³`, in closed form, verified against that
    /// expectation.
    ///
    /// # Errors
    ///
    /// See [`ProbError`].
    pub fn skewness(
        &self,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::skewness(self, pool)
    }

    /// `γ₂ = E[((X - μ)/σ)⁴] - 3 = κ₄/σ⁴`, in closed form, verified.
    ///
    /// **Excess** kurtosis — `0` for a normal, not `3`.
    ///
    /// # Errors
    ///
    /// See [`ProbError`].
    pub fn excess_kurtosis(
        &self,
        pool: &ExprPool,
    ) -> Result<crate::deriv::DerivedExpr<ExprId>, ProbError> {
        genfun::excess_kurtosis(self, pool)
    }
}

/// Highest moment order the numeric gate can still *fail* on, and therefore
/// the highest [`Distribution::moment`] will return.
///
/// Not a capability limit — a verification limit. `E[X¹⁶]` under a unit normal
/// is `2027025`; the quadrature that checks it integrates `z¹⁶ e^{-z²/2}`,
/// whose mass sits where `z¹⁶` spans twenty orders of magnitude, and the
/// double-exponential rule's relative error there stops being small enough to
/// distinguish the right answer from a neighbouring wrong one. Past this the
/// check would pass unconditionally, so the moment is refused instead of
/// returned unchecked.
pub const MAX_MOMENT_ORDER: u32 = 12;

// ---------------------------------------------------------------------------
// Small shared helpers
// ---------------------------------------------------------------------------

/// `π`, the way the rest of the crate spells it.
pub(crate) fn pi(pool: &ExprPool) -> ExprId {
    pool.symbol("pi", Domain::Real)
}

/// `a - b`.
pub(crate) fn sub(a: ExprId, b: ExprId, pool: &ExprPool) -> ExprId {
    let neg = pool.integer(-1);
    pool.add(vec![a, pool.mul(vec![neg, b])])
}

/// `a / b`.
pub(crate) fn div(a: ExprId, b: ExprId, pool: &ExprPool) -> ExprId {
    let inv = pool.pow(b, pool.integer(-1));
    pool.mul(vec![a, inv])
}