symplex 0.11.1

//! Core node types for the symplex symbolic expression graph.
//!
//! Every expression in symplex is represented as a node in a directed acyclic
//! graph (DAG).  Nodes are identified by [`ExprId`] handles that index into an
//! arena held by the expression context.  The actual payload of each node is
//! described by the [`ExprNode`] enum.
//!
//! Numeric literals and symbolic names are stored in separate side‐tables and
//! referenced via [`NumId`] and [`SymbolId`] respectively, keeping the core
//! node type small and cheap to clone.

use smallvec::{SmallVec, smallvec};
use std::fmt;

// ---------------------------------------------------------------------------
// Id newtypes
// ---------------------------------------------------------------------------

/// Opaque handle that identifies an expression node inside an expression
/// context.
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct ExprId(pub u32);

impl fmt::Debug for ExprId {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "e{}", self.0)
    }
}

/// Index into the numeric literal side‐table of an expression context.
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct NumId(pub u32);

impl fmt::Debug for NumId {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "n{}", self.0)
    }
}

/// Index into the symbol name table of an expression context.
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct SymbolId(pub u32);

impl fmt::Debug for SymbolId {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "s{}", self.0)
    }
}

/// Identifies an expression context (arena).
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct CtxId(pub u32);

impl fmt::Debug for CtxId {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "ctx{}", self.0)
    }
}

// ---------------------------------------------------------------------------
// ExprNode
// ---------------------------------------------------------------------------

/// The payload of a single node in the expression DAG.
///
/// Variants fall into the following categories:
///
/// * **Atoms** – leaves that carry no child expressions: `Num`, `Symbol`,
///   the mathematical constants (`Pi`, `E`, `ImaginaryUnit`, `EulerGamma`,
///   `Catalan`, `GoldenRatio`, `PhysicalConstant`), the special values
///   (`Infinity`, `NegInfinity`, `ComplexInfinity`, `NaN`), the boolean
///   atoms (`BoolTrue`, `BoolFalse`) and the set atoms (`EmptySet`,
///   `UniversalSet`).
/// * **N-ary arithmetic** – `Add`, `Mul`, `Min`, `Max`.
/// * **Binary arithmetic** – `Pow`.
/// * **Elementary functions** – `Neg`, `Floor`, `Ceiling`, trigonometric,
///   hyperbolic and their inverses, `Exp`, `Ln`, `Abs`, `Sign`,
///   `Heaviside`, `DiracDelta`.
/// * **Complex analysis** – `Re`, `Im`, `Conjugate`, `Arg`.
/// * **Special functions** – `Gamma`, `LogGamma`, `Digamma`, `Polygamma`,
///   `Erf`, `Erfc`, `LambertW`, `Beta`, `Si`, `Ci`, `Ei`, `Li`, `Zeta`,
///   `KroneckerDelta`.
/// * **Combinatorial** – `Factorial`, `Binomial`.
/// * **Boolean / relational / logical** – `Gt`, `Ge`, `Eq_`, `Ne`, `And`,
///   `Or`, `Not`, `Piecewise`.
/// * **Function application** – `Apply` (user-defined or library functions
///   identified by name).
/// * **Calculus / formal (unevaluated) forms** – `Derivative`, `Integral`,
///   `Sum`, `Product_`, `Limit`, `Series`, `LaplaceTransform`,
///   `InverseLaplaceTransform`, `Residue`, `DSolve`, `ConditionSet`.
/// * **Algebraic answers** – `RootOf`, `RootSum` (complete closed-form
///   descriptions of polynomial roots; *not* unevaluated).
/// * **Sets** – `Interval`, `FiniteSet`, `SetUnion`, `SetIntersection`,
///   `SetComplement`.
#[derive(Clone, PartialEq, Eq, Hash)]
pub enum ExprNode {
    // -- atoms ---------------------------------------------------------------
    /// A numeric literal.  The actual value lives in the number side‐table
    /// and is looked up via the contained `NumId`.
    Num(NumId),

    /// A symbolic name (variable, constant, or function name).  The actual
    /// string is stored in the symbol table and looked up via `SymbolId`.
    Symbol(SymbolId),

    /// The mathematical constant π ≈ 3.14159…
    Pi,

    /// Euler's number *e* ≈ 2.71828…
    E,

    /// The imaginary unit *i*, satisfying *i*² = −1.
    ImaginaryUnit,

    /// The Euler–Mascheroni constant γ ≈ 0.57721 56649…
    ///
    /// Whether γ is rational or irrational is an open problem, so the
    /// assumption system deliberately leaves those properties unknown.
    EulerGamma,

    /// Catalan's constant G = Σ (−1)ⁿ/(2n+1)² ≈ 0.91596 55941…
    ///
    /// Irrationality of G is unproven; only positivity, realness and
    /// finiteness are asserted.
    Catalan,

    /// The golden ratio φ = (1 + √5)/2 ≈ 1.61803 39887…
    ///
    /// φ is algebraic (root of x² − x − 1) and irrational.
    GoldenRatio,

    /// A named physical constant with a known exact value.
    /// Displays as name, evaluates to value. E.g., speed of light, Planck's constant.
    PhysicalConstant(SymbolId, ExprId),

    /// Positive infinity (+∞).
    Infinity,

    /// Negative infinity (−∞).
    NegInfinity,

    /// Complex infinity (∞̃) — magnitude is infinite, direction is undefined.
    ComplexInfinity,

    /// Not‐a‐number, representing an indeterminate or undefined result.
    NaN,

    // -- n‐ary operators -----------------------------------------------------
    /// An n‐ary sum: `a + b + c + …`
    Add(SmallVec<[ExprId; 6]>),

    /// An n‐ary product: `a · b · c · …`
    Mul(SmallVec<[ExprId; 6]>),

    // -- n-ary min/max -------------------------------------------------------
    /// N-ary minimum: min(a, b, c, ...).
    Min(SmallVec<[ExprId; 4]>),

    /// N-ary maximum: max(a, b, c, ...).
    Max(SmallVec<[ExprId; 4]>),

    // -- binary operators ----------------------------------------------------
    /// Exponentiation: `base ^ exponent`.
    Pow(ExprId, ExprId),

    // -- unary operators -----------------------------------------------------
    /// Unary arithmetic negation: `−x`.
    Neg(ExprId),

    /// Floor function: ⌊x⌋ (greatest integer ≤ x).
    Floor(ExprId),

    /// Ceiling function: ⌈x⌉ (least integer ≥ x).
    Ceiling(ExprId),

    /// Sine function: `sin(x)`.
    Sin(ExprId),

    /// Cosine function: `cos(x)`.
    Cos(ExprId),

    /// Tangent function: `tan(x)`.
    Tan(ExprId),

    /// Natural exponential function: `eˣ`.
    Exp(ExprId),

    /// Natural logarithm: `ln(x)`.
    Ln(ExprId),

    /// Absolute value (or complex modulus): `|x|`.
    Abs(ExprId),

    /// Inverse sine: `asin(x)` (arcsin).
    Asin(ExprId),

    /// Inverse cosine: `acos(x)` (arccos).
    Acos(ExprId),

    /// Inverse tangent: `atan(x)` (arctan).
    Atan(ExprId),

    /// Two-argument arctangent: `atan2(y, x)` — gives the angle in (-π, π].
    /// Correctly handles all four quadrants, unlike `atan(y/x)`.
    Atan2(ExprId, ExprId),

    /// Hyperbolic sine: `sinh(x)`.
    Sinh(ExprId),

    /// Hyperbolic cosine: `cosh(x)`.
    Cosh(ExprId),

    /// Hyperbolic tangent: `tanh(x)`.
    Tanh(ExprId),

    /// Inverse hyperbolic sine: `asinh(x)`.
    Asinh(ExprId),

    /// Inverse hyperbolic cosine: `acosh(x)`.
    Acosh(ExprId),

    /// Inverse hyperbolic tangent: `atanh(x)`.
    Atanh(ExprId),

    /// Sign function: `sign(x)` = 1 if x > 0, -1 if x < 0, 0 if x = 0.
    Sign(ExprId),

    /// Heaviside step function: H(x) = 0 for x<0, 1/2 for x=0, 1 for x>0.
    Heaviside(ExprId),

    /// Dirac delta distribution: δ(x) = 0 for x≠0, symbolic at x=0.
    DiracDelta(ExprId),

    // -- complex analysis ----------------------------------------------------
    /// Real part of a complex quantity: `re(z)`.
    ///
    /// Always real-valued. Only constructed when the real part cannot be
    /// determined from the structure of `z` and the assumption system.
    Re(ExprId),

    /// Imaginary part of a complex quantity: `im(z)`, a real number such
    /// that `z = re(z) + i·im(z)`.
    Im(ExprId),

    /// Complex conjugate: `conjugate(z)`.
    ///
    /// Distributes over `Add`/`Mul`/integer `Pow` and commutes with
    /// real-analytic functions at construction; the node itself is only
    /// produced for arguments whose realness is unknown.
    Conjugate(ExprId),

    /// Principal complex argument `arg(z) ∈ (−π, π]`.
    Arg(ExprId),

    // -- special functions ---------------------------------------------------
    /// Gamma function: Γ(x) = ∫₀^∞ t^(x-1) e^(-t) dt.
    Gamma(ExprId),

    /// Log-gamma function: ln(Γ(x)).
    LogGamma(ExprId),

    /// Digamma function: ψ(x) = Γ'(x)/Γ(x).
    Digamma(ExprId),

    /// Error function: erf(x) = 2/√π ∫₀ˣ e^(-t²) dt.
    Erf(ExprId),

    /// Complementary error function: erfc(x) = 1 - erf(x).
    Erfc(ExprId),

    /// Lambert W function (principal branch): W(x)·exp(W(x)) = x.
    LambertW(ExprId),

    /// Beta function: B(a,b) = Γ(a)Γ(b)/Γ(a+b).
    Beta(ExprId, ExprId),

    /// Sine integral: Si(x) = ∫₀ˣ sin(t)/t dt.
    Si(ExprId),

    /// Cosine integral: Ci(x) = γ + ln(x) + ∫₀ˣ (cos(t) − 1)/t dt.
    Ci(ExprId),

    /// Exponential integral (Cauchy principal value):
    /// Ei(x) = −∫_{−x}^{∞} e^{−t}/t dt = γ + ln|x| + Σ_{k≥1} xᵏ/(k·k!).
    Ei(ExprId),

    /// Logarithmic integral: li(x) = ∫₀ˣ dt/ln(t) = Ei(ln x).
    Li(ExprId),

    /// Riemann zeta function ζ(s).
    Zeta(ExprId),

    /// Polygamma function ψ⁽ⁿ⁾(x) = dⁿ⁺¹/dxⁿ⁺¹ ln Γ(x): `Polygamma(n, x)`.
    ///
    /// `Polygamma(0, x)` is canonicalised to [`Digamma`](ExprNode::Digamma).
    Polygamma(ExprId, ExprId),

    /// Kronecker delta δᵢⱼ = 1 if i = j, else 0: `KroneckerDelta(i, j)`.
    ///
    /// Arguments are stored in canonical (sorted) order since δ is symmetric.
    KroneckerDelta(ExprId, ExprId),

    // -- combinatorial -------------------------------------------------------
    /// Factorial: `n!`
    Factorial(ExprId),

    /// Binomial coefficient: `C(n, k)` = n! / (k! * (n-k)!)
    Binomial(ExprId, ExprId),

    // -- boolean atoms -------------------------------------------------------
    /// Boolean true.
    BoolTrue,
    /// Boolean false.
    BoolFalse,

    // -- relational operators ------------------------------------------------
    /// Greater than: `a > b`.
    Gt(ExprId, ExprId),
    /// Greater than or equal: `a >= b`.
    Ge(ExprId, ExprId),
    /// Mathematical equality test: `a == b` (boolean-valued).
    Eq_(ExprId, ExprId),
    /// Not equal: `a != b`.
    Ne(ExprId, ExprId),

    // -- logical connectives -------------------------------------------------
    /// Logical conjunction (n-ary): `a && b && c`.
    And(SmallVec<[ExprId; 6]>),
    /// Logical disjunction (n-ary): `a || b || c`.
    Or(SmallVec<[ExprId; 6]>),
    /// Logical negation: `!a`.
    Not(ExprId),

    // -- piecewise -----------------------------------------------------------
    /// Piecewise function: each pair is (value, condition).
    Piecewise(SmallVec<[(ExprId, ExprId); 3]>),

    // -- composite forms -----------------------------------------------------
    /// Application of a user‐defined or library function identified by
    /// `SymbolId` to a list of argument expressions.
    Apply(SymbolId, SmallVec<[ExprId; 2]>),

    /// Formal derivative: d/d(var) of an expression.
    ///
    /// `Derivative(body, var)` represents ∂/∂`var` `body`.
    Derivative(ExprId, ExprId),

    /// Formal indefinite integral: ∫ expr d(var).
    ///
    /// `Integral(body, var)` represents ∫ `body` d`var`.
    Integral(ExprId, ExprId),

    /// Formal definite integral: ∫ₗₒʰⁱ expr d(var).
    ///
    /// `DefiniteIntegral(body, var, lo, hi)` represents `∫_lo^hi body dvar`.
    /// It is the unevaluated form returned by definite integration when no
    /// closed form can be established, so that the bounds are preserved for
    /// display, differentiation (Leibniz rule) and numeric quadrature.
    /// `var` is bound inside `body` only; `lo` and `hi` are in the outer
    /// scope.
    DefiniteIntegral(ExprId, ExprId, ExprId, ExprId),

    /// Symbolic summation: Sum(body, var, lower, upper).
    Sum(ExprId, ExprId, ExprId, ExprId),

    /// Symbolic product: Product(body, var, lower, upper).
    Product_(ExprId, ExprId, ExprId, ExprId),

    /// Formal limit: lim_{var → point} body.
    Limit(ExprId, ExprId, ExprId),

    /// Formal series expansion: Series(body, var, point, order).
    Series(ExprId, ExprId, ExprId, ExprId),

    /// Formal Laplace transform: ℒ{body}(t → s).
    LaplaceTransform(ExprId, ExprId, ExprId),

    /// Formal inverse Laplace transform: ℒ⁻¹{body}(s → t).
    InverseLaplaceTransform(ExprId, ExprId, ExprId),

    /// Formal residue: Res_{var=point} body.
    Residue(ExprId, ExprId, ExprId),

    /// Root of a polynomial: the index-th root of poly.
    RootOf(ExprId, ExprId),

    /// Sum over roots of a polynomial: Σ_{α: poly(α)=0} body(α, x).
    ///
    /// Used by Rothstein-Trager when the resultant has irreducible factors
    /// of degree ≥ 5 whose roots cannot be expressed in radicals.
    ///
    /// `RootSum(poly, body, sumvar)` where:
    /// - `poly` is the polynomial expression whose roots are summed over
    /// - `body` is the expression evaluated at each root (contains `sumvar`)
    /// - `sumvar` is the bound summation variable symbol
    RootSum(ExprId, ExprId, ExprId),

    /// Formal ODE solution: DSolve(expr=0, func, var).
    DSolve(ExprId, ExprId, ExprId),

    /// Condition set: {var | condition}.
    ConditionSet(ExprId, ExprId),

    // -- set atoms --------------------------------------------------------
    /// The empty set ∅.
    EmptySet,
    /// The universal set (all values).
    UniversalSet,

    // -- set constructors -------------------------------------------------
    /// A closed/open interval [a, b] with flags encoding open/closed.
    /// Bits: 0x01 = left_open, 0x02 = right_open.
    Interval(ExprId, ExprId, u8),

    /// A finite set of elements {a, b, c, ...}, sorted and deduplicated.
    FiniteSet(SmallVec<[ExprId; 4]>),

    /// Union of sets: A ∪ B ∪ C ∪ ...
    SetUnion(SmallVec<[ExprId; 4]>),

    /// Intersection of sets: A ∩ B ∩ C ∩ ...
    SetIntersection(SmallVec<[ExprId; 4]>),

    /// Set complement (relative): A \ B.
    SetComplement(ExprId, ExprId),
}

/// Interval flag: left endpoint is open (excluded).
pub const INTERVAL_LEFT_OPEN: u8 = 0x01;
/// Interval flag: right endpoint is open (excluded).
pub const INTERVAL_RIGHT_OPEN: u8 = 0x02;
/// Interval flag: both endpoints are open.
pub const INTERVAL_BOTH_OPEN: u8 = 0x03;
/// Interval flag: both endpoints are closed.
pub const INTERVAL_BOTH_CLOSED: u8 = 0x00;

impl ExprNode {
    /// Returns all child `ExprId`s contained in this node.
    ///
    /// For atoms (numbers, symbols, constants, and special values) the
    /// returned collection is empty.  For compound nodes the children are
    /// returned in the order they appear in the variant.
    pub fn children(&self) -> SmallVec<[ExprId; 6]> {
        match self {
            // atoms — no children
            ExprNode::Num(_)
            | ExprNode::Symbol(_)
            | ExprNode::Pi
            | ExprNode::E
            | ExprNode::ImaginaryUnit
            | ExprNode::EulerGamma
            | ExprNode::Catalan
            | ExprNode::GoldenRatio
            | ExprNode::PhysicalConstant(_, _)
            | ExprNode::Infinity
            | ExprNode::NegInfinity
            | ExprNode::ComplexInfinity
            | ExprNode::NaN
            | ExprNode::BoolTrue
            | ExprNode::BoolFalse
            | ExprNode::EmptySet
            | ExprNode::UniversalSet => smallvec![],

            // n‐ary
            ExprNode::Add(ids) | ExprNode::Mul(ids) | ExprNode::And(ids) | ExprNode::Or(ids) => {
                ids.clone()
            }

            // n-ary min/max (SmallVec<[ExprId; 4]> → SmallVec<[ExprId; 6]>)
            ExprNode::Min(ids) | ExprNode::Max(ids) => ids.iter().copied().collect(),

            // n-ary set ops (SmallVec<[ExprId; 4]> → SmallVec<[ExprId; 6]>)
            ExprNode::FiniteSet(ids) | ExprNode::SetUnion(ids) | ExprNode::SetIntersection(ids) => {
                ids.iter().copied().collect()
            }

            // piecewise — flatten pairs into children list
            ExprNode::Piecewise(pairs) => {
                let mut result = SmallVec::new();
                for &(val, cond) in pairs {
                    result.push(val);
                    result.push(cond);
                }
                result
            }

            // binary
            ExprNode::Pow(a, b)
            | ExprNode::Atan2(a, b)
            | ExprNode::Binomial(a, b)
            | ExprNode::Beta(a, b)
            | ExprNode::Polygamma(a, b)
            | ExprNode::KroneckerDelta(a, b)
            | ExprNode::Gt(a, b)
            | ExprNode::Ge(a, b)
            | ExprNode::Eq_(a, b)
            | ExprNode::Ne(a, b)
            | ExprNode::Derivative(a, b)
            | ExprNode::Integral(a, b)
            | ExprNode::SetComplement(a, b)
            | ExprNode::RootOf(a, b)
            | ExprNode::ConditionSet(a, b) => {
                smallvec![*a, *b]
            }

            // ternary-ish: Interval(a, b, _flags)
            ExprNode::Interval(a, b, _) => {
                smallvec![*a, *b]
            }

            // ternary
            ExprNode::Limit(a, b, c)
            | ExprNode::LaplaceTransform(a, b, c)
            | ExprNode::InverseLaplaceTransform(a, b, c)
            | ExprNode::Residue(a, b, c)
            | ExprNode::DSolve(a, b, c)
            | ExprNode::RootSum(a, b, c) => smallvec![*a, *b, *c],

            // 4-ary: Sum, Product_, Series, DefiniteIntegral
            ExprNode::Sum(a, b, c, d)
            | ExprNode::Product_(a, b, c, d)
            | ExprNode::Series(a, b, c, d)
            | ExprNode::DefiniteIntegral(a, b, c, d) => {
                smallvec![*a, *b, *c, *d]
            }

            // unary
            ExprNode::Neg(x)
            | ExprNode::Floor(x)
            | ExprNode::Ceiling(x)
            | ExprNode::Sin(x)
            | ExprNode::Cos(x)
            | ExprNode::Tan(x)
            | ExprNode::Exp(x)
            | ExprNode::Ln(x)
            | ExprNode::Abs(x)
            | ExprNode::Asin(x)
            | ExprNode::Acos(x)
            | ExprNode::Atan(x)
            | ExprNode::Sinh(x)
            | ExprNode::Cosh(x)
            | ExprNode::Tanh(x)
            | ExprNode::Asinh(x)
            | ExprNode::Acosh(x)
            | ExprNode::Atanh(x)
            | ExprNode::Sign(x)
            | ExprNode::Heaviside(x)
            | ExprNode::DiracDelta(x)
            | ExprNode::Re(x)
            | ExprNode::Im(x)
            | ExprNode::Conjugate(x)
            | ExprNode::Arg(x)
            | ExprNode::Gamma(x)
            | ExprNode::LogGamma(x)
            | ExprNode::Digamma(x)
            | ExprNode::Erf(x)
            | ExprNode::Erfc(x)
            | ExprNode::LambertW(x)
            | ExprNode::Si(x)
            | ExprNode::Ci(x)
            | ExprNode::Ei(x)
            | ExprNode::Li(x)
            | ExprNode::Zeta(x)
            | ExprNode::Factorial(x)
            | ExprNode::Not(x) => smallvec![*x],

            // function application
            ExprNode::Apply(_, args) => {
                // Promote the SmallVec<[ExprId; 2]> into SmallVec<[ExprId; 6]>
                args.iter().copied().collect()
            }
        }
    }

    /// Calls `f` for each child `ExprId` without allocating.
    ///
    /// This is the zero-allocation alternative to [`children()`](Self::children).
    /// Prefer this in hot paths (tree walks, display) where the allocation
    /// from `children()` would be significant.
    #[inline]
    pub fn for_each_child(&self, mut f: impl FnMut(ExprId)) {
        match self {
            // atoms — no children
            ExprNode::Num(_)
            | ExprNode::Symbol(_)
            | ExprNode::Pi
            | ExprNode::E
            | ExprNode::ImaginaryUnit
            | ExprNode::EulerGamma
            | ExprNode::Catalan
            | ExprNode::GoldenRatio
            | ExprNode::PhysicalConstant(_, _)
            | ExprNode::Infinity
            | ExprNode::NegInfinity
            | ExprNode::ComplexInfinity
            | ExprNode::NaN
            | ExprNode::BoolTrue
            | ExprNode::BoolFalse
            | ExprNode::EmptySet
            | ExprNode::UniversalSet => {}

            // n-ary
            ExprNode::Add(ids) | ExprNode::Mul(ids) | ExprNode::And(ids) | ExprNode::Or(ids) => {
                for &id in ids {
                    f(id);
                }
            }

            // n-ary min/max
            ExprNode::Min(ids) | ExprNode::Max(ids) => {
                for &id in ids {
                    f(id);
                }
            }

            // n-ary set ops
            ExprNode::FiniteSet(ids) | ExprNode::SetUnion(ids) | ExprNode::SetIntersection(ids) => {
                for &id in ids {
                    f(id);
                }
            }

            // piecewise — flatten pairs
            ExprNode::Piecewise(pairs) => {
                for &(val, cond) in pairs {
                    f(val);
                    f(cond);
                }
            }

            // binary
            ExprNode::Pow(a, b)
            | ExprNode::Atan2(a, b)
            | ExprNode::Binomial(a, b)
            | ExprNode::Beta(a, b)
            | ExprNode::Polygamma(a, b)
            | ExprNode::KroneckerDelta(a, b)
            | ExprNode::Gt(a, b)
            | ExprNode::Ge(a, b)
            | ExprNode::Eq_(a, b)
            | ExprNode::Ne(a, b)
            | ExprNode::Derivative(a, b)
            | ExprNode::Integral(a, b)
            | ExprNode::SetComplement(a, b)
            | ExprNode::RootOf(a, b)
            | ExprNode::ConditionSet(a, b) => {
                f(*a);
                f(*b);
            }

            ExprNode::Interval(a, b, _) => {
                f(*a);
                f(*b);
            }

            // ternary
            ExprNode::Limit(a, b, c)
            | ExprNode::LaplaceTransform(a, b, c)
            | ExprNode::InverseLaplaceTransform(a, b, c)
            | ExprNode::Residue(a, b, c)
            | ExprNode::DSolve(a, b, c)
            | ExprNode::RootSum(a, b, c) => {
                f(*a);
                f(*b);
                f(*c);
            }

            // 4-ary
            ExprNode::Sum(a, b, c, d)
            | ExprNode::Product_(a, b, c, d)
            | ExprNode::Series(a, b, c, d)
            | ExprNode::DefiniteIntegral(a, b, c, d) => {
                f(*a);
                f(*b);
                f(*c);
                f(*d);
            }

            // unary
            ExprNode::Neg(x)
            | ExprNode::Floor(x)
            | ExprNode::Ceiling(x)
            | ExprNode::Sin(x)
            | ExprNode::Cos(x)
            | ExprNode::Tan(x)
            | ExprNode::Exp(x)
            | ExprNode::Ln(x)
            | ExprNode::Abs(x)
            | ExprNode::Asin(x)
            | ExprNode::Acos(x)
            | ExprNode::Atan(x)
            | ExprNode::Sinh(x)
            | ExprNode::Cosh(x)
            | ExprNode::Tanh(x)
            | ExprNode::Asinh(x)
            | ExprNode::Acosh(x)
            | ExprNode::Atanh(x)
            | ExprNode::Sign(x)
            | ExprNode::Heaviside(x)
            | ExprNode::DiracDelta(x)
            | ExprNode::Re(x)
            | ExprNode::Im(x)
            | ExprNode::Conjugate(x)
            | ExprNode::Arg(x)
            | ExprNode::Gamma(x)
            | ExprNode::LogGamma(x)
            | ExprNode::Digamma(x)
            | ExprNode::Erf(x)
            | ExprNode::Erfc(x)
            | ExprNode::LambertW(x)
            | ExprNode::Si(x)
            | ExprNode::Ci(x)
            | ExprNode::Ei(x)
            | ExprNode::Li(x)
            | ExprNode::Zeta(x)
            | ExprNode::Factorial(x)
            | ExprNode::Not(x) => f(*x),

            // function application
            ExprNode::Apply(_, args) => {
                for &id in args {
                    f(id);
                }
            }
        }
    }

    /// Returns the number of children without allocating.
    #[inline]
    pub fn child_count(&self) -> usize {
        match self {
            ExprNode::Num(_)
            | ExprNode::Symbol(_)
            | ExprNode::Pi
            | ExprNode::E
            | ExprNode::ImaginaryUnit
            | ExprNode::EulerGamma
            | ExprNode::Catalan
            | ExprNode::GoldenRatio
            | ExprNode::PhysicalConstant(_, _)
            | ExprNode::Infinity
            | ExprNode::NegInfinity
            | ExprNode::ComplexInfinity
            | ExprNode::NaN
            | ExprNode::BoolTrue
            | ExprNode::BoolFalse
            | ExprNode::EmptySet
            | ExprNode::UniversalSet => 0,
            ExprNode::Add(ids) | ExprNode::Mul(ids) | ExprNode::And(ids) | ExprNode::Or(ids) => {
                ids.len()
            }
            ExprNode::Min(ids) | ExprNode::Max(ids) => ids.len(),
            ExprNode::FiniteSet(ids) | ExprNode::SetUnion(ids) | ExprNode::SetIntersection(ids) => {
                ids.len()
            }
            ExprNode::Piecewise(pairs) => pairs.len() * 2,
            ExprNode::Pow(..)
            | ExprNode::Atan2(..)
            | ExprNode::Binomial(..)
            | ExprNode::Beta(..)
            | ExprNode::Polygamma(..)
            | ExprNode::KroneckerDelta(..)
            | ExprNode::Gt(..)
            | ExprNode::Ge(..)
            | ExprNode::Eq_(..)
            | ExprNode::Ne(..)
            | ExprNode::Derivative(..)
            | ExprNode::Integral(..)
            | ExprNode::SetComplement(..)
            | ExprNode::Interval(..)
            | ExprNode::RootOf(..)
            | ExprNode::ConditionSet(..) => 2,
            ExprNode::Limit(..)
            | ExprNode::LaplaceTransform(..)
            | ExprNode::InverseLaplaceTransform(..)
            | ExprNode::Residue(..)
            | ExprNode::DSolve(..)
            | ExprNode::RootSum(..) => 3,
            ExprNode::Sum(..)
            | ExprNode::Product_(..)
            | ExprNode::Series(..)
            | ExprNode::DefiniteIntegral(..) => 4,
            ExprNode::Neg(_)
            | ExprNode::Floor(_)
            | ExprNode::Ceiling(_)
            | ExprNode::Sin(_)
            | ExprNode::Cos(_)
            | ExprNode::Tan(_)
            | ExprNode::Exp(_)
            | ExprNode::Ln(_)
            | ExprNode::Abs(_)
            | ExprNode::Asin(_)
            | ExprNode::Acos(_)
            | ExprNode::Atan(_)
            | ExprNode::Sinh(_)
            | ExprNode::Cosh(_)
            | ExprNode::Tanh(_)
            | ExprNode::Asinh(_)
            | ExprNode::Acosh(_)
            | ExprNode::Atanh(_)
            | ExprNode::Sign(_)
            | ExprNode::Heaviside(_)
            | ExprNode::DiracDelta(_)
            | ExprNode::Re(_)
            | ExprNode::Im(_)
            | ExprNode::Conjugate(_)
            | ExprNode::Arg(_)
            | ExprNode::Gamma(_)
            | ExprNode::LogGamma(_)
            | ExprNode::Digamma(_)
            | ExprNode::Erf(_)
            | ExprNode::Erfc(_)
            | ExprNode::LambertW(_)
            | ExprNode::Si(_)
            | ExprNode::Ci(_)
            | ExprNode::Ei(_)
            | ExprNode::Li(_)
            | ExprNode::Zeta(_)
            | ExprNode::Factorial(_)
            | ExprNode::Not(_) => 1,
            ExprNode::Apply(_, args) => args.len(),
        }
    }

    /// Returns `true` if this node is an atom (a leaf with no child
    /// expressions).
    ///
    /// Atoms are: [`Num`](ExprNode::Num), [`Symbol`](ExprNode::Symbol),
    /// [`Pi`](ExprNode::Pi), [`E`](ExprNode::E),
    /// [`ImaginaryUnit`](ExprNode::ImaginaryUnit),
    /// [`EulerGamma`](ExprNode::EulerGamma), [`Catalan`](ExprNode::Catalan),
    /// [`GoldenRatio`](ExprNode::GoldenRatio),
    /// [`PhysicalConstant`](ExprNode::PhysicalConstant),
    /// [`Infinity`](ExprNode::Infinity),
    /// [`NegInfinity`](ExprNode::NegInfinity),
    /// [`ComplexInfinity`](ExprNode::ComplexInfinity),
    /// [`NaN`](ExprNode::NaN), the boolean atoms and the set atoms.
    pub fn is_atom(&self) -> bool {
        matches!(
            self,
            ExprNode::Num(_)
                | ExprNode::Symbol(_)
                | ExprNode::Pi
                | ExprNode::E
                | ExprNode::ImaginaryUnit
                | ExprNode::EulerGamma
                | ExprNode::Catalan
                | ExprNode::GoldenRatio
                | ExprNode::PhysicalConstant(_, _)
                | ExprNode::Infinity
                | ExprNode::NegInfinity
                | ExprNode::ComplexInfinity
                | ExprNode::NaN
                | ExprNode::BoolTrue
                | ExprNode::BoolFalse
                | ExprNode::EmptySet
                | ExprNode::UniversalSet
        )
    }

    /// Returns `true` if this node is a set-valued expression.
    pub fn is_set_node(&self) -> bool {
        matches!(
            self,
            ExprNode::EmptySet
                | ExprNode::UniversalSet
                | ExprNode::Interval(..)
                | ExprNode::FiniteSet(_)
                | ExprNode::SetUnion(_)
                | ExprNode::SetIntersection(_)
                | ExprNode::SetComplement(..)
        )
    }
}

impl fmt::Debug for ExprNode {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            ExprNode::Num(id) => write!(f, "Num({id:?})"),
            ExprNode::Symbol(id) => write!(f, "Symbol({id:?})"),
            ExprNode::Pi => write!(f, "Pi"),
            ExprNode::E => write!(f, "E"),
            ExprNode::ImaginaryUnit => write!(f, "ImaginaryUnit"),
            ExprNode::EulerGamma => write!(f, "EulerGamma"),
            ExprNode::Catalan => write!(f, "Catalan"),
            ExprNode::GoldenRatio => write!(f, "GoldenRatio"),
            ExprNode::PhysicalConstant(name, val) => {
                write!(f, "PhysicalConstant({name:?}, {val:?})")
            }
            ExprNode::Infinity => write!(f, "Infinity"),
            ExprNode::NegInfinity => write!(f, "NegInfinity"),
            ExprNode::ComplexInfinity => write!(f, "ComplexInfinity"),
            ExprNode::NaN => write!(f, "NaN"),
            ExprNode::Add(ids) => f.debug_tuple("Add").field(ids).finish(),
            ExprNode::Mul(ids) => f.debug_tuple("Mul").field(ids).finish(),
            ExprNode::Pow(base, exp) => f.debug_tuple("Pow").field(base).field(exp).finish(),
            ExprNode::Neg(x) => f.debug_tuple("Neg").field(x).finish(),
            ExprNode::Sin(x) => f.debug_tuple("Sin").field(x).finish(),
            ExprNode::Cos(x) => f.debug_tuple("Cos").field(x).finish(),
            ExprNode::Tan(x) => f.debug_tuple("Tan").field(x).finish(),
            ExprNode::Exp(x) => f.debug_tuple("Exp").field(x).finish(),
            ExprNode::Ln(x) => f.debug_tuple("Ln").field(x).finish(),
            ExprNode::Abs(x) => f.debug_tuple("Abs").field(x).finish(),
            ExprNode::Asin(x) => f.debug_tuple("Asin").field(x).finish(),
            ExprNode::Acos(x) => f.debug_tuple("Acos").field(x).finish(),
            ExprNode::Atan(x) => f.debug_tuple("Atan").field(x).finish(),
            ExprNode::Atan2(y, x) => f.debug_tuple("Atan2").field(y).field(x).finish(),
            ExprNode::Sinh(x) => f.debug_tuple("Sinh").field(x).finish(),
            ExprNode::Cosh(x) => f.debug_tuple("Cosh").field(x).finish(),
            ExprNode::Tanh(x) => f.debug_tuple("Tanh").field(x).finish(),
            ExprNode::Asinh(x) => f.debug_tuple("Asinh").field(x).finish(),
            ExprNode::Acosh(x) => f.debug_tuple("Acosh").field(x).finish(),
            ExprNode::Atanh(x) => f.debug_tuple("Atanh").field(x).finish(),
            ExprNode::Floor(x) => f.debug_tuple("Floor").field(x).finish(),
            ExprNode::Ceiling(x) => f.debug_tuple("Ceiling").field(x).finish(),
            ExprNode::Min(ids) => f.debug_tuple("Min").field(ids).finish(),
            ExprNode::Max(ids) => f.debug_tuple("Max").field(ids).finish(),
            ExprNode::Sign(id) => write!(f, "Sign({id:?})"),
            ExprNode::Heaviside(x) => f.debug_tuple("Heaviside").field(x).finish(),
            ExprNode::DiracDelta(x) => f.debug_tuple("DiracDelta").field(x).finish(),
            ExprNode::Re(x) => f.debug_tuple("Re").field(x).finish(),
            ExprNode::Im(x) => f.debug_tuple("Im").field(x).finish(),
            ExprNode::Conjugate(x) => f.debug_tuple("Conjugate").field(x).finish(),
            ExprNode::Arg(x) => f.debug_tuple("Arg").field(x).finish(),
            ExprNode::Gamma(x) => f.debug_tuple("Gamma").field(x).finish(),
            ExprNode::LogGamma(x) => f.debug_tuple("LogGamma").field(x).finish(),
            ExprNode::Digamma(x) => f.debug_tuple("Digamma").field(x).finish(),
            ExprNode::Erf(x) => f.debug_tuple("Erf").field(x).finish(),
            ExprNode::Erfc(x) => f.debug_tuple("Erfc").field(x).finish(),
            ExprNode::LambertW(x) => f.debug_tuple("LambertW").field(x).finish(),
            ExprNode::Beta(a, b) => f.debug_tuple("Beta").field(a).field(b).finish(),
            ExprNode::Si(x) => f.debug_tuple("Si").field(x).finish(),
            ExprNode::Ci(x) => f.debug_tuple("Ci").field(x).finish(),
            ExprNode::Ei(x) => f.debug_tuple("Ei").field(x).finish(),
            ExprNode::Li(x) => f.debug_tuple("Li").field(x).finish(),
            ExprNode::Zeta(x) => f.debug_tuple("Zeta").field(x).finish(),
            ExprNode::Polygamma(n, x) => f.debug_tuple("Polygamma").field(n).field(x).finish(),
            ExprNode::KroneckerDelta(i, j) => {
                f.debug_tuple("KroneckerDelta").field(i).field(j).finish()
            }
            ExprNode::Factorial(id) => write!(f, "Factorial({id:?})"),
            ExprNode::Binomial(n, k) => write!(f, "Binomial({n:?}, {k:?})"),
            ExprNode::BoolTrue => write!(f, "BoolTrue"),
            ExprNode::BoolFalse => write!(f, "BoolFalse"),
            ExprNode::Gt(a, b) => f.debug_tuple("Gt").field(a).field(b).finish(),
            ExprNode::Ge(a, b) => f.debug_tuple("Ge").field(a).field(b).finish(),
            ExprNode::Eq_(a, b) => f.debug_tuple("Eq_").field(a).field(b).finish(),
            ExprNode::Ne(a, b) => f.debug_tuple("Ne").field(a).field(b).finish(),
            ExprNode::And(ids) => f.debug_tuple("And").field(ids).finish(),
            ExprNode::Or(ids) => f.debug_tuple("Or").field(ids).finish(),
            ExprNode::Not(x) => f.debug_tuple("Not").field(x).finish(),
            ExprNode::Piecewise(pairs) => {
                let mut d = f.debug_tuple("Piecewise");
                for pair in pairs {
                    d.field(pair);
                }
                d.finish()
            }
            ExprNode::Apply(sym, args) => f.debug_tuple("Apply").field(sym).field(args).finish(),
            ExprNode::Derivative(body, var) => {
                f.debug_tuple("Derivative").field(body).field(var).finish()
            }
            ExprNode::Integral(body, var) => {
                f.debug_tuple("Integral").field(body).field(var).finish()
            }
            ExprNode::DefiniteIntegral(body, var, lo, hi) => f
                .debug_tuple("DefiniteIntegral")
                .field(body)
                .field(var)
                .field(lo)
                .field(hi)
                .finish(),
            ExprNode::Sum(body, var, lo, hi) => f
                .debug_tuple("Sum")
                .field(body)
                .field(var)
                .field(lo)
                .field(hi)
                .finish(),
            ExprNode::Product_(body, var, lo, hi) => f
                .debug_tuple("Product_")
                .field(body)
                .field(var)
                .field(lo)
                .field(hi)
                .finish(),
            ExprNode::EmptySet => write!(f, "EmptySet"),
            ExprNode::UniversalSet => write!(f, "UniversalSet"),
            ExprNode::Interval(a, b, flags) => f
                .debug_tuple("Interval")
                .field(a)
                .field(b)
                .field(flags)
                .finish(),
            ExprNode::FiniteSet(ids) => f.debug_tuple("FiniteSet").field(ids).finish(),
            ExprNode::SetUnion(ids) => f.debug_tuple("SetUnion").field(ids).finish(),
            ExprNode::SetIntersection(ids) => f.debug_tuple("SetIntersection").field(ids).finish(),
            ExprNode::SetComplement(a, b) => {
                f.debug_tuple("SetComplement").field(a).field(b).finish()
            }
            ExprNode::Limit(body, var, point) => f
                .debug_tuple("Limit")
                .field(body)
                .field(var)
                .field(point)
                .finish(),
            ExprNode::Series(body, var, point, order) => f
                .debug_tuple("Series")
                .field(body)
                .field(var)
                .field(point)
                .field(order)
                .finish(),
            ExprNode::LaplaceTransform(body, t, s) => f
                .debug_tuple("LaplaceTransform")
                .field(body)
                .field(t)
                .field(s)
                .finish(),
            ExprNode::InverseLaplaceTransform(body, s, t) => f
                .debug_tuple("InverseLaplaceTransform")
                .field(body)
                .field(s)
                .field(t)
                .finish(),
            ExprNode::Residue(body, var, point) => f
                .debug_tuple("Residue")
                .field(body)
                .field(var)
                .field(point)
                .finish(),
            ExprNode::RootOf(poly, idx) => f.debug_tuple("RootOf").field(poly).field(idx).finish(),
            ExprNode::DSolve(expr, func, var) => f
                .debug_tuple("DSolve")
                .field(expr)
                .field(func)
                .field(var)
                .finish(),
            ExprNode::RootSum(poly, body, sumvar) => f
                .debug_tuple("RootSum")
                .field(poly)
                .field(body)
                .field(sumvar)
                .finish(),
            ExprNode::ConditionSet(var, cond) => f
                .debug_tuple("ConditionSet")
                .field(var)
                .field(cond)
                .finish(),
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn debug_expr_id() {
        let id = ExprId(42);
        assert_eq!(format!("{id:?}"), "e42");
    }

    #[test]
    fn debug_num_id() {
        let id = NumId(7);
        assert_eq!(format!("{id:?}"), "n7");
    }

    #[test]
    fn debug_symbol_id() {
        let id = SymbolId(3);
        assert_eq!(format!("{id:?}"), "s3");
    }

    #[test]
    fn debug_ctx_id() {
        let id = CtxId(0);
        assert_eq!(format!("{id:?}"), "ctx0");
    }

    #[test]
    fn atom_has_no_children() {
        let node = ExprNode::Pi;
        assert!(node.is_atom());
        assert!(node.children().is_empty());
    }

    #[test]
    fn num_is_atom() {
        let node = ExprNode::Num(NumId(0));
        assert!(node.is_atom());
        assert!(node.children().is_empty());
    }

    #[test]
    fn symbol_is_atom() {
        let node = ExprNode::Symbol(SymbolId(1));
        assert!(node.is_atom());
    }

    #[test]
    fn special_values_are_atoms() {
        for node in [
            ExprNode::Infinity,
            ExprNode::NegInfinity,
            ExprNode::ComplexInfinity,
            ExprNode::NaN,
            ExprNode::ImaginaryUnit,
            ExprNode::E,
        ] {
            assert!(node.is_atom(), "{node:?} should be an atom");
            assert!(node.children().is_empty());
        }
    }

    #[test]
    fn add_children() {
        let ids: SmallVec<[ExprId; 6]> = smallvec![ExprId(1), ExprId(2), ExprId(3)];
        let node = ExprNode::Add(ids.clone());
        assert!(!node.is_atom());
        assert_eq!(node.children(), ids);
    }

    #[test]
    fn mul_children() {
        let ids: SmallVec<[ExprId; 6]> = smallvec![ExprId(4), ExprId(5)];
        let node = ExprNode::Mul(ids.clone());
        assert!(!node.is_atom());
        assert_eq!(node.children(), ids);
    }

    #[test]
    fn pow_children() {
        let node = ExprNode::Pow(ExprId(10), ExprId(20));
        assert!(!node.is_atom());
        let kids = node.children();
        assert_eq!(kids.len(), 2);
        assert_eq!(kids[0], ExprId(10));
        assert_eq!(kids[1], ExprId(20));
    }

    #[test]
    fn unary_children() {
        for node in [
            ExprNode::Neg(ExprId(1)),
            ExprNode::Sin(ExprId(1)),
            ExprNode::Cos(ExprId(1)),
            ExprNode::Tan(ExprId(1)),
            ExprNode::Exp(ExprId(1)),
            ExprNode::Ln(ExprId(1)),
            ExprNode::Abs(ExprId(1)),
        ] {
            assert!(!node.is_atom());
            let kids = node.children();
            assert_eq!(kids.len(), 1, "{node:?} should have exactly 1 child");
            assert_eq!(kids[0], ExprId(1));
        }
    }

    #[test]
    fn apply_children() {
        let args: SmallVec<[ExprId; 2]> = smallvec![ExprId(5), ExprId(6)];
        let node = ExprNode::Apply(SymbolId(0), args);
        assert!(!node.is_atom());
        let kids = node.children();
        assert_eq!(kids.len(), 2);
        assert_eq!(kids[0], ExprId(5));
        assert_eq!(kids[1], ExprId(6));
    }

    #[test]
    fn derivative_children() {
        let node = ExprNode::Derivative(ExprId(3), ExprId(7));
        assert!(!node.is_atom());
        let kids = node.children();
        assert_eq!(kids.len(), 2);
        assert_eq!(kids[0], ExprId(3));
        assert_eq!(kids[1], ExprId(7));
    }

    #[test]
    fn integral_children() {
        let node = ExprNode::Integral(ExprId(8), ExprId(9));
        assert!(!node.is_atom());
        let kids = node.children();
        assert_eq!(kids.len(), 2);
        assert_eq!(kids[0], ExprId(8));
        assert_eq!(kids[1], ExprId(9));
    }

    #[test]
    fn named_constants_are_atoms() {
        for node in [
            ExprNode::EulerGamma,
            ExprNode::Catalan,
            ExprNode::GoldenRatio,
        ] {
            assert!(node.is_atom(), "{node:?} should be an atom");
            assert!(node.children().is_empty());
            assert_eq!(node.child_count(), 0);
            let mut count = 0;
            node.for_each_child(|_| count += 1);
            assert_eq!(count, 0);
        }
        assert_eq!(format!("{:?}", ExprNode::EulerGamma), "EulerGamma");
        assert_eq!(format!("{:?}", ExprNode::Catalan), "Catalan");
        assert_eq!(format!("{:?}", ExprNode::GoldenRatio), "GoldenRatio");
    }

    #[test]
    fn complex_and_special_unary_children() {
        for node in [
            ExprNode::Re(ExprId(3)),
            ExprNode::Im(ExprId(3)),
            ExprNode::Conjugate(ExprId(3)),
            ExprNode::Arg(ExprId(3)),
            ExprNode::Si(ExprId(3)),
            ExprNode::Ci(ExprId(3)),
            ExprNode::Ei(ExprId(3)),
            ExprNode::Li(ExprId(3)),
            ExprNode::Zeta(ExprId(3)),
        ] {
            assert!(!node.is_atom());
            assert_eq!(node.child_count(), 1);
            let kids = node.children();
            assert_eq!(kids.len(), 1, "{node:?} should have exactly 1 child");
            assert_eq!(kids[0], ExprId(3));
            let mut seen = Vec::new();
            node.for_each_child(|c| seen.push(c));
            assert_eq!(seen, vec![ExprId(3)]);
        }
    }

    #[test]
    fn polygamma_and_kronecker_children() {
        for node in [
            ExprNode::Polygamma(ExprId(1), ExprId(2)),
            ExprNode::KroneckerDelta(ExprId(1), ExprId(2)),
        ] {
            assert!(!node.is_atom());
            assert_eq!(node.child_count(), 2);
            let kids = node.children();
            assert_eq!(kids.len(), 2);
            assert_eq!(kids[0], ExprId(1));
            assert_eq!(kids[1], ExprId(2));
            let mut seen = Vec::new();
            node.for_each_child(|c| seen.push(c));
            assert_eq!(seen, vec![ExprId(1), ExprId(2)]);
        }
        assert_eq!(
            format!("{:?}", ExprNode::Polygamma(ExprId(1), ExprId(2))),
            "Polygamma(e1, e2)"
        );
        assert_eq!(
            format!("{:?}", ExprNode::KroneckerDelta(ExprId(1), ExprId(2))),
            "KroneckerDelta(e1, e2)"
        );
    }
}