vj-sim 0.2.0

//! Real-valued constraints (IEEE 1800-2023 18.5.3 and 18.5.13; plan Task C9).
//!
//! # Why a separate solver
//!
//! The BDD solver ([`super::solve`]) is a bit-vector engine; a `real` is a
//! 64-bit IEEE pattern, and bit-blasting a floating-point comparison is both
//! enormous and meaningless for sampling (nearly every bit pattern is a
//! number far from any bound). Real random variables therefore live in their
//! own namespace and are solved here, by interval arithmetic on the float
//! grid, after the integral solve and with the same threaded PRNG seed.
//!
//! # The domain
//!
//! A real variable ranges over the finite values of its type (`real`: f64,
//! `shortreal`: f32; `-0.0` is the same grid point as `+0.0`). A set of
//! values is a [`RealSet`]: sorted disjoint closed intervals of *grid indices*
//! (the order-preserving integer numbering of the floats), so open and closed
//! bounds, `!=` holes and adjacent floats need no special cases. An
//! unconstrained variable is uniform over `[-MAX, MAX]`.
//!
//! # What a constraint may be
//!
//! A conjunct reads at most one random real variable (affinely: a product
//! has at most one side that reads it, a quotient's denominator reads none),
//! or it is a guard (`->`/`if`) whose predicate reads one variable and whose
//! arms constrain other variables one at a time. Each atom `l op r` is turned
//! into the exact set of grid values for which the IEEE evaluation of the
//! atom is true: the sign of the slope is read off the affine form and the
//! boundary is found by binary search on the grid with the real evaluator, so
//! `x / 3.0 < 1.0` includes exactly the floats for which that rounded
//! comparison holds. (The one inexactness: when *both* sides read the
//! variable and rounding makes `l < r` non-monotone within an ulp of the
//! boundary the search lands on one boundary of that zone.)
//!
//! Anything else (two variables in one non-guard conjunct, a real conjunct
//! reading an integral random variable, a non-affine product) is a loud
//! [`RealError`] or a build error, never a silently wrong distribution.
//!
//! # Sampling
//!
//! A variable is drawn *uniformly by length* over its feasible set (a point
//! set, such as `x == 3.0 || x == 5.0`, uniformly over its points; an
//! isolated point inside a set of positive length has probability zero), or
//! by `dist` weights (a `:/` range's weight is spread over its length, so
//! the part of it other constraints exclude loses weight in proportion, the
//! real analogue of 18.5.3's per-value weights). A guard that couples two
//! variables (`a inside [..] -> b == X`) is solved by the staged rule of
//! 18.5.9: it requires `solve a before b`, `a` is drawn from exactly the
//! values for which some `b` completes the solution, and `b` is then drawn
//! from what remains.

use super::guard::{self, GuardValue};
use super::ir::{Assignment, ConstraintEvalError, Expr, StateId, Ty, VarId};
use super::rng;
use num_bigint::{BigInt, BigUint};
use num_traits::{ToPrimitive, Zero};

/// Index of a random real variable (its own numbering, separate from the
/// integral [`super::ir::VarId`]s).
pub type RealVarId = u32;

/// The two floating-point types a real random variable can have.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, serde::Serialize, serde::Deserialize)]
pub enum FloatKind {
    /// `real` / `realtime`: IEEE 754 binary64.
    F64,
    /// `shortreal`: IEEE 754 binary32.
    F32,
}

impl FloatKind {
    /// The grid index of the largest finite value.
    pub fn max_idx(self) -> i64 {
        match self {
            FloatKind::F64 => 0x7fef_ffff_ffff_ffff,
            FloatKind::F32 => 0x7f7f_ffff,
        }
    }

    /// The largest finite value.
    pub fn max_val(self) -> f64 {
        self.val(self.max_idx())
    }

    /// The grid index of the finite or infinite grid value `x` (`-0.0` is 0).
    /// `x` must already be a value of the type ([`Self::round`]).
    pub fn idx(self, x: f64) -> i64 {
        match self {
            FloatKind::F64 => {
                let b = x.to_bits();
                let m = (b & 0x7fff_ffff_ffff_ffff) as i64;
                if b >> 63 == 1 { -m } else { m }
            }
            FloatKind::F32 => {
                let b = (x as f32).to_bits();
                let m = i64::from(b & 0x7fff_ffff);
                if b >> 31 == 1 { -m } else { m }
            }
        }
    }

    /// The value of grid index `i` (never `-0.0`).
    pub fn val(self, i: i64) -> f64 {
        match self {
            FloatKind::F64 => {
                if i >= 0 {
                    f64::from_bits(i as u64)
                } else {
                    f64::from_bits(i.unsigned_abs() | (1 << 63))
                }
            }
            FloatKind::F32 => {
                if i >= 0 {
                    f64::from(f32::from_bits(i as u32))
                } else {
                    f64::from(f32::from_bits(i.unsigned_abs() as u32 | (1 << 31)))
                }
            }
        }
    }

    /// `x` rounded to the type (a no-op for `real`).
    pub fn round(self, x: f64) -> f64 {
        match self {
            FloatKind::F64 => x,
            FloatKind::F32 => f64::from(x as f32),
        }
    }

    /// The smallest grid index whose value is `>= x`, or `None` if there is
    /// none (`x` above the largest finite value, or NaN).
    pub fn ceil_idx(self, x: f64) -> Option<i64> {
        if x.is_nan() || x > self.max_val() {
            return None;
        }
        if x <= -self.max_val() {
            return Some(-self.max_idx());
        }
        let r = self.round(x);
        let i = self.idx(r);
        Some(if r < x { i + 1 } else { i })
    }

    /// The largest grid index whose value is `<= x`, or `None`.
    pub fn floor_idx(self, x: f64) -> Option<i64> {
        if x.is_nan() || x < -self.max_val() {
            return None;
        }
        if x >= self.max_val() {
            return Some(self.max_idx());
        }
        let r = self.round(x);
        let i = self.idx(r);
        Some(if r > x { i - 1 } else { i })
    }

    /// The value held by the low bits of a field (`words` of an object
    /// property: the f64 pattern, or the f32 pattern in the low 32 bits).
    pub fn from_bits(self, w: &BigUint) -> f64 {
        let b = w.to_u64().unwrap_or(0);
        match self {
            FloatKind::F64 => f64::from_bits(b),
            FloatKind::F32 => f64::from(f32::from_bits(b as u32)),
        }
    }

    /// The field bits of the grid value `x`.
    pub fn to_bits(self, x: f64) -> u64 {
        match self {
            FloatKind::F64 => x.to_bits(),
            FloatKind::F32 => u64::from((x as f32).to_bits()),
        }
    }

    /// The field width in bits.
    pub fn width(self) -> u32 {
        match self {
            FloatKind::F64 => 64,
            FloatKind::F32 => 32,
        }
    }
}

/// A set of grid values: sorted, disjoint, non-adjacent closed intervals of
/// grid indices (module docs).
#[derive(Clone, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub struct RealSet {
    kind: FloatKind,
    iv: Vec<(i64, i64)>,
}

/// `[0, 1)` from a draw.
fn unit(seed: &mut u64) -> f64 {
    (rng::draw(seed) >> 11) as f64 / (1u64 << 53) as f64
}

impl RealSet {
    /// The empty set.
    pub fn empty(kind: FloatKind) -> RealSet {
        RealSet {
            kind,
            iv: Vec::new(),
        }
    }

    /// Every finite value.
    pub fn full(kind: FloatKind) -> RealSet {
        RealSet {
            kind,
            iv: vec![(-kind.max_idx(), kind.max_idx())],
        }
    }

    /// The grid indices `lo..=hi` (empty if `lo > hi`), clipped to the
    /// finite values.
    pub fn from_idx(kind: FloatKind, lo: i64, hi: i64) -> RealSet {
        let (lo, hi) = (lo.max(-kind.max_idx()), hi.min(kind.max_idx()));
        RealSet {
            kind,
            iv: if lo <= hi { vec![(lo, hi)] } else { Vec::new() },
        }
    }

    /// The values `v` of the type with `lo <= v <= hi` (compared as exact
    /// numbers); `lo` or `hi` may be infinite. Empty if `lo > hi` or NaN.
    pub fn range(kind: FloatKind, lo: f64, hi: f64) -> RealSet {
        match (kind.ceil_idx(lo), kind.floor_idx(hi)) {
            (Some(a), Some(b)) => RealSet::from_idx(kind, a, b),
            _ => RealSet::empty(kind),
        }
    }

    /// The single value `x` if it is a value of the type, else the empty set
    /// (a `shortreal` never equals a `real` that needs more than 24 bits).
    pub fn point(kind: FloatKind, x: f64) -> RealSet {
        if x.is_nan() || x.abs() > kind.max_val() || kind.round(x) != x {
            return RealSet::empty(kind);
        }
        let i = kind.idx(x);
        RealSet::from_idx(kind, i, i)
    }

    /// The type of the values.
    pub fn kind(&self) -> FloatKind {
        self.kind
    }

    /// Whether no value is in the set.
    pub fn is_empty(&self) -> bool {
        self.iv.is_empty()
    }

    /// Whether `x` (a value of the type) is in the set.
    pub fn contains(&self, x: f64) -> bool {
        if x.is_nan() || x.abs() > self.kind.max_val() {
            return false;
        }
        let i = self.kind.idx(x);
        self.iv.iter().any(|&(a, b)| a <= i && i <= b)
    }

    /// The intervals, as inclusive grid-index pairs.
    pub fn intervals(&self) -> &[(i64, i64)] {
        &self.iv
    }

    fn normalised(kind: FloatKind, mut iv: Vec<(i64, i64)>) -> RealSet {
        iv.sort_unstable();
        let mut out: Vec<(i64, i64)> = Vec::with_capacity(iv.len());
        for (a, b) in iv {
            match out.last_mut() {
                Some(last) if a <= last.1 + 1 => last.1 = last.1.max(b),
                _ => out.push((a, b)),
            }
        }
        RealSet { kind, iv: out }
    }

    /// The union.
    pub fn union(&self, o: &RealSet) -> RealSet {
        let mut iv = self.iv.clone();
        iv.extend_from_slice(&o.iv);
        RealSet::normalised(self.kind, iv)
    }

    /// The intersection.
    pub fn intersect(&self, o: &RealSet) -> RealSet {
        let mut out = Vec::new();
        let (mut i, mut j) = (0, 0);
        while i < self.iv.len() && j < o.iv.len() {
            let (a, b) = (self.iv[i], o.iv[j]);
            let (lo, hi) = (a.0.max(b.0), a.1.min(b.1));
            if lo <= hi {
                out.push((lo, hi));
            }
            if a.1 < b.1 {
                i += 1;
            } else {
                j += 1;
            }
        }
        RealSet {
            kind: self.kind,
            iv: out,
        }
    }

    /// Every value of the type not in the set.
    pub fn complement(&self) -> RealSet {
        let m = self.kind.max_idx();
        let mut out = Vec::new();
        let mut next = -m;
        for &(a, b) in &self.iv {
            if a > next {
                out.push((next, a - 1));
            }
            next = b + 1;
        }
        if next <= m {
            out.push((next, m));
        }
        RealSet {
            kind: self.kind,
            iv: out,
        }
    }

    /// The Lebesgue measure in half-units (so the full range does not
    /// overflow): `0.5 * hi - 0.5 * lo` summed over the intervals. A single
    /// point has measure 0.
    pub fn measure(&self) -> f64 {
        self.iv
            .iter()
            .map(|&(a, b)| 0.5 * self.kind.val(b) - 0.5 * self.kind.val(a))
            .sum()
    }

    /// The number of grid points.
    fn points(&self) -> BigUint {
        self.iv
            .iter()
            .map(|&(a, b)| BigUint::from((i128::from(b) - i128::from(a) + 1) as u128))
            .sum()
    }

    /// One value drawn uniformly by length (module docs); the set must not be
    /// empty.
    pub fn sample(&self, seed: &mut u64) -> f64 {
        debug_assert!(!self.is_empty(), "sampling an empty set");
        let k = self.kind;
        let total = self.measure();
        if !total.is_finite() || total <= 0.0 {
            // Only points (or a measure that underflows): uniform over the
            // grid points.
            let n = self.points();
            let mut r = rng::below(seed, &n);
            for &(a, b) in &self.iv {
                let len = BigUint::from((i128::from(b) - i128::from(a) + 1) as u128);
                if r < len {
                    let off = r.to_i64().unwrap_or(0);
                    return k.val(a + off);
                }
                r -= len;
            }
            return k.val(self.iv[self.iv.len() - 1].1);
        }
        let mut target = unit(seed) * total;
        let mut chosen = self.iv[self.iv.len() - 1];
        for &(a, b) in &self.iv {
            let len = 0.5 * k.val(b) - 0.5 * k.val(a);
            if target < len {
                chosen = (a, b);
                break;
            }
            target -= len;
        }
        let (lo, hi) = (k.val(chosen.0), k.val(chosen.1));
        let (h_lo, h_hi) = (0.5 * lo, 0.5 * hi);
        let x = 2.0 * (h_lo + unit(seed) * (h_hi - h_lo));
        let x = k.round(x.clamp(lo, hi));
        x.clamp(lo, hi)
    }
}

/// A real value with bitwise equality, so IR containing reals can derive
/// `Eq`.
#[derive(Clone, Copy, Debug, serde::Serialize, serde::Deserialize)]
pub struct Real(u64);

impl Real {
    /// Wraps `x`.
    pub fn new(x: f64) -> Real {
        Real(x.to_bits())
    }
    /// The value.
    pub fn get(self) -> f64 {
        f64::from_bits(self.0)
    }
}

impl PartialEq for Real {
    fn eq(&self, o: &Real) -> bool {
        self.0 == o.0
    }
}
impl Eq for Real {}

/// Why a real solve failed.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum RealError {
    /// The constraints have no solution (`randomize()` returns 0).
    Unsatisfiable,
    /// A constraint guard evaluated to `ERROR` (18.5.12).
    GuardError,
    /// A constraint reads a state variable with no value (X/Z bits, or a
    /// null handle on the way to it).
    StateUnreadable(StateId),
    /// An expression produced NaN, or divided by zero.
    NonFinite(&'static str),
    /// A `dist` weight is negative.
    BadWeight,
    /// The problem is not one this solver handles (a malformed program, or
    /// a constraint shape [`RealProblem::check`] would have rejected).
    Shape(String),
}

impl core::fmt::Display for RealError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            RealError::Unsatisfiable => f.write_str("the real constraints have no solution"),
            RealError::GuardError => f.write_str("a real constraint guard evaluated to ERROR"),
            RealError::StateUnreadable(s) => {
                write!(f, "state variable #{s} of a real constraint cannot be read")
            }
            RealError::NonFinite(what) => write!(f, "real constraint evaluation: {what}"),
            RealError::BadWeight => f.write_str("a real `dist` weight is negative"),
            RealError::Shape(what) => write!(f, "unsupported real constraint: {what}"),
        }
    }
}

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

/// A real-valued expression over real random variables, state variables and
/// constants. `shortreal` operands compute in binary32, anything involving a
/// `real` in binary64.
#[derive(Clone, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub enum RealExpr {
    /// A constant of the given type.
    Const(Real, FloatKind),
    /// A random real variable.
    Var(RealVarId),
    /// A state variable holding the bits of a real of the given type.
    State(StateId, FloatKind),
    /// An integral, variable-free expression converted to `real` (signed or
    /// unsigned by its own type).
    FromInt(Box<Expr>),
    /// Unary minus.
    Neg(Box<RealExpr>),
    /// `+`.
    Add(Box<RealExpr>, Box<RealExpr>),
    /// `-`.
    Sub(Box<RealExpr>, Box<RealExpr>),
    /// `*`.
    Mul(Box<RealExpr>, Box<RealExpr>),
    /// `/`.
    Div(Box<RealExpr>, Box<RealExpr>),
    /// Conversion to `shortreal` (rounds to binary32).
    ToShort(Box<RealExpr>),
}

/// A comparison operator.
#[derive(Clone, Copy, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub enum CmpOp {
    /// `<`
    Lt,
    /// `<=`
    Le,
    /// `>`
    Gt,
    /// `>=`
    Ge,
    /// `==`
    Eq,
    /// `!=`
    Ne,
}

/// A boolean condition over real expressions.
#[derive(Clone, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub enum RealCond {
    /// A constant.
    Const(bool),
    /// An integral condition over state variables only (true when nonzero).
    Int(Expr),
    /// A comparison.
    Cmp(CmpOp, RealExpr, RealExpr),
    /// `!`.
    Not(Box<RealCond>),
    /// `&&` of any number of conditions.
    And(Vec<RealCond>),
    /// `||` of any number of conditions.
    Or(Vec<RealCond>),
    /// `if (pred) then else els` (an implication has the constant-true else
    /// arm): a constraint guard (18.5.12).
    Guard(Box<RealCond>, Box<RealCond>, Box<RealCond>),
}

/// The values of one real `dist` item.
#[derive(Clone, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub enum RealSpan {
    /// A single value.
    Point(RealExpr),
    /// `[lo:hi]`; `None` is `$`, the type's own extreme.
    Range {
        /// The low bound.
        lo: Option<RealExpr>,
        /// The high bound.
        hi: Option<RealExpr>,
    },
}

/// One explicit item of a real `dist`.
#[derive(Clone, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub struct RealDistItem {
    /// The item's values.
    pub span: RealSpan,
    /// The weight (integral, over states and constants). A range's weight is
    /// the item's whole weight (`:/`, the only operator a real range may
    /// use, 18.5.3); a point's weight is that of the value.
    pub weight: Expr,
}

/// `var dist { ... }` on a real random variable.
#[derive(Clone, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub struct RealDist {
    /// The distributed variable.
    pub var: RealVarId,
    /// The explicit items.
    pub items: Vec<RealDistItem>,
    /// The weight of `default :/ w`, spread over the values no item names.
    pub default: Option<Expr>,
}

/// The whole real problem of one class template.
#[derive(Clone, Debug, Default, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub struct RealProblem {
    /// The type of each real random variable.
    pub kinds: Vec<FloatKind>,
    /// Conditions that must all hold.
    pub conds: Vec<RealCond>,
    /// The `dist` constraints, at most one per variable.
    pub dists: Vec<RealDist>,
    /// `solve a before b` edges between real variables.
    pub before: Vec<(RealVarId, RealVarId)>,
}

impl RealExpr {
    fn map_int(&self, f: &mut dyn FnMut(&Expr) -> Expr) -> RealExpr {
        let mut b = |x: &RealExpr| Box::new(x.map_int(f));
        match self {
            RealExpr::Const(..) | RealExpr::Var(_) | RealExpr::State(..) => self.clone(),
            RealExpr::FromInt(x) => RealExpr::FromInt(Box::new(f(x))),
            RealExpr::Neg(x) => RealExpr::Neg(b(x)),
            RealExpr::ToShort(x) => RealExpr::ToShort(b(x)),
            RealExpr::Add(x, y) => RealExpr::Add(b(x), b(y)),
            RealExpr::Sub(x, y) => RealExpr::Sub(b(x), b(y)),
            RealExpr::Mul(x, y) => RealExpr::Mul(b(x), b(y)),
            RealExpr::Div(x, y) => RealExpr::Div(b(x), b(y)),
        }
    }

    fn for_each_int(&self, f: &mut dyn FnMut(&Expr)) {
        match self {
            RealExpr::Const(..) | RealExpr::Var(_) | RealExpr::State(..) => {}
            RealExpr::FromInt(x) => f(x),
            RealExpr::Neg(x) | RealExpr::ToShort(x) => x.for_each_int(f),
            RealExpr::Add(x, y)
            | RealExpr::Sub(x, y)
            | RealExpr::Mul(x, y)
            | RealExpr::Div(x, y) => {
                x.for_each_int(f);
                y.for_each_int(f);
            }
        }
    }
}

impl RealCond {
    fn map_int(&self, f: &mut dyn FnMut(&Expr) -> Expr) -> RealCond {
        match self {
            RealCond::Const(_) => self.clone(),
            RealCond::Int(e) => RealCond::Int(f(e)),
            RealCond::Cmp(op, l, r) => RealCond::Cmp(*op, l.map_int(f), r.map_int(f)),
            RealCond::Not(x) => RealCond::Not(Box::new(x.map_int(f))),
            RealCond::And(v) => RealCond::And(v.iter().map(|x| x.map_int(f)).collect()),
            RealCond::Or(v) => RealCond::Or(v.iter().map(|x| x.map_int(f)).collect()),
            RealCond::Guard(p, t, e) => RealCond::Guard(
                Box::new(p.map_int(f)),
                Box::new(t.map_int(f)),
                Box::new(e.map_int(f)),
            ),
        }
    }

    fn for_each_int(&self, f: &mut dyn FnMut(&Expr)) {
        match self {
            RealCond::Const(_) => {}
            RealCond::Int(e) => f(e),
            RealCond::Cmp(_, l, r) => {
                l.for_each_int(f);
                r.for_each_int(f);
            }
            RealCond::Not(x) => x.for_each_int(f),
            RealCond::And(v) | RealCond::Or(v) => v.iter().for_each(|x| x.for_each_int(f)),
            RealCond::Guard(p, t, e) => {
                p.for_each_int(f);
                t.for_each_int(f);
                e.for_each_int(f);
            }
        }
    }
}

impl RealProblem {
    /// Whether any integral sub-expression of a condition or a `dist` bound
    /// reads an integral *random* variable (IEEE 1800-2023 18.5.13: a real
    /// constraint over an integral `rand` variable is solved after it).
    pub fn reads_int_vars(&self) -> bool {
        let mut any = false;
        let mut f = |e: &Expr| e.for_each_var(&mut |_, _| any = true);
        for c in &self.conds {
            c.for_each_int(&mut f);
        }
        for d in &self.dists {
            for it in &d.items {
                match &it.span {
                    RealSpan::Point(x) => x.for_each_int(&mut f),
                    RealSpan::Range { lo, hi } => {
                        lo.iter().chain(hi).for_each(|x| x.for_each_int(&mut f));
                    }
                }
            }
        }
        any
    }

    /// Whether some top-level condition compares an integral random variable
    /// with real values and reads no real random variable: it only filters the
    /// integral draw, which rejection sampling cannot do reliably (a 32-bit
    /// variable against a real bound), so the builder refuses it by name.
    pub fn has_integral_only_filter(&self) -> bool {
        self.conds.iter().any(|c| {
            let (mut int_var, mut real_var) = (false, false);
            c.for_each_int(&mut |e| e.for_each_var(&mut |_, _| int_var = true));
            c.for_each_var(&mut |_| real_var = true);
            int_var && !real_var
        })
    }

    /// This problem with every integral random variable `v` of type `t` that a
    /// condition or `dist` bound reads replaced by `f(v, t)` (a state variable
    /// holding the value the integral solve drew).
    pub fn with_int_vars(&self, mut f: impl FnMut(VarId, Ty) -> Expr) -> RealProblem {
        let mut m = |e: &Expr| e.substitute_vars(&mut f);
        let mut p = self.clone();
        p.conds = self.conds.iter().map(|c| c.map_int(&mut m)).collect();
        for d in &mut p.dists {
            for it in &mut d.items {
                it.span = match &it.span {
                    RealSpan::Point(x) => RealSpan::Point(x.map_int(&mut m)),
                    RealSpan::Range { lo, hi } => RealSpan::Range {
                        lo: lo.as_ref().map(|x| x.map_int(&mut m)),
                        hi: hi.as_ref().map(|x| x.map_int(&mut m)),
                    },
                };
            }
        }
        p
    }
}

impl RealExpr {
    /// The type the expression computes in.
    pub fn kind(&self, kinds: &[FloatKind]) -> FloatKind {
        match self {
            RealExpr::Const(_, k) | RealExpr::State(_, k) => *k,
            RealExpr::Var(v) => kinds.get(*v as usize).copied().unwrap_or(FloatKind::F64),
            RealExpr::FromInt(_) => FloatKind::F64,
            RealExpr::ToShort(_) => FloatKind::F32,
            RealExpr::Neg(x) => x.kind(kinds),
            RealExpr::Add(a, b)
            | RealExpr::Sub(a, b)
            | RealExpr::Mul(a, b)
            | RealExpr::Div(a, b) => {
                if a.kind(kinds) == FloatKind::F32 && b.kind(kinds) == FloatKind::F32 {
                    FloatKind::F32
                } else {
                    FloatKind::F64
                }
            }
        }
    }

    /// Calls `f` on every random variable read.
    pub fn for_each_var(&self, f: &mut impl FnMut(RealVarId)) {
        match self {
            RealExpr::Var(v) => f(*v),
            RealExpr::Const(..) | RealExpr::State(..) | RealExpr::FromInt(_) => {}
            RealExpr::Neg(x) | RealExpr::ToShort(x) => x.for_each_var(f),
            RealExpr::Add(a, b)
            | RealExpr::Sub(a, b)
            | RealExpr::Mul(a, b)
            | RealExpr::Div(a, b) => {
                a.for_each_var(f);
                b.for_each_var(f);
            }
        }
    }

    /// Calls `f` on every state variable read (including those of an
    /// integral sub-expression).
    pub fn for_each_state(&self, f: &mut impl FnMut(StateId)) {
        match self {
            RealExpr::State(s, _) => f(*s),
            RealExpr::FromInt(e) => e.for_each_state(f),
            RealExpr::Const(..) | RealExpr::Var(_) => {}
            RealExpr::Neg(x) | RealExpr::ToShort(x) => x.for_each_state(f),
            RealExpr::Add(a, b)
            | RealExpr::Sub(a, b)
            | RealExpr::Mul(a, b)
            | RealExpr::Div(a, b) => {
                a.for_each_state(f);
                b.for_each_state(f);
            }
        }
    }
}

impl RealCond {
    /// Calls `f` on every random variable read.
    pub fn for_each_var(&self, f: &mut impl FnMut(RealVarId)) {
        match self {
            RealCond::Const(_) | RealCond::Int(_) => {}
            RealCond::Cmp(_, l, r) => {
                l.for_each_var(f);
                r.for_each_var(f);
            }
            RealCond::Not(x) => x.for_each_var(f),
            RealCond::And(xs) | RealCond::Or(xs) => xs.iter().for_each(|x| x.for_each_var(f)),
            RealCond::Guard(p, t, e) => {
                p.for_each_var(f);
                t.for_each_var(f);
                e.for_each_var(f);
            }
        }
    }

    /// The distinct random variables read, ascending.
    pub fn vars(&self) -> Vec<RealVarId> {
        let mut v = Vec::new();
        self.for_each_var(&mut |x| v.push(x));
        v.sort_unstable();
        v.dedup();
        v
    }
}

/// What an evaluation reads: the type of every variable, the values of the
/// variables fixed for this evaluation, and the state values.
struct Ctx<'a> {
    kinds: &'a [FloatKind],
    /// `Some` for a variable with a value (inactive, or the one being probed).
    vars: Vec<Option<f64>>,
    states: &'a [Option<BigUint>],
}

fn int_value(e: &Expr, states: &[Option<BigUint>]) -> Result<BigUint, RealError> {
    let a = Assignment {
        vars: Vec::new(),
        states: states.to_vec(),
    };
    super::ir::eval(e, &a).map_err(|err| match err {
        ConstraintEvalError::StateUnknown(s) => RealError::StateUnreadable(s),
        ConstraintEvalError::DivideByZero => RealError::NonFinite("integer division by zero"),
        ConstraintEvalError::UnboundVar(_)
        | ConstraintEvalError::UnboundState(_)
        | ConstraintEvalError::Unexpanded => {
            RealError::Shape("an integral sub-expression reads an unbound variable".into())
        }
    })
}

/// A two's-complement `v` of `e`'s type as a number.
fn signed_number(v: &BigUint, e: &Expr) -> BigInt {
    let t = e.ty();
    if t.signed && v.bit(u64::from(t.width) - 1) {
        BigInt::from(v.clone()) - (BigInt::from(1) << t.width)
    } else {
        BigInt::from(v.clone())
    }
}

impl Ctx<'_> {
    fn state(&self, s: StateId) -> Result<&BigUint, RealError> {
        match self.states.get(s as usize) {
            Some(Some(v)) => Ok(v),
            _ => Err(RealError::StateUnreadable(s)),
        }
    }

    /// The value of `e` (module docs for the NaN and division rules).
    fn eval(&self, e: &RealExpr) -> Result<f64, RealError> {
        let finite = |x: f64, k: FloatKind| -> Result<f64, RealError> {
            let x = k.round(x);
            if x.is_nan() {
                Err(RealError::NonFinite("a real expression is NaN"))
            } else {
                Ok(x)
            }
        };
        Ok(match e {
            RealExpr::Const(r, _) => r.get(),
            RealExpr::Var(v) => self
                .vars
                .get(*v as usize)
                .copied()
                .flatten()
                .ok_or_else(|| RealError::Shape(format!("real variable #{v} has no value here")))?,
            RealExpr::State(s, k) => k.from_bits(self.state(*s)?),
            RealExpr::FromInt(x) => {
                let v = int_value(x, self.states)?;
                signed_number(&v, x).to_f64().unwrap_or(f64::NAN)
            }
            RealExpr::Neg(x) => -self.eval(x)?,
            RealExpr::ToShort(x) => finite(self.eval(x)?, FloatKind::F32)?,
            RealExpr::Add(a, b)
            | RealExpr::Sub(a, b)
            | RealExpr::Mul(a, b)
            | RealExpr::Div(a, b) => {
                let (x, y) = (self.eval(a)?, self.eval(b)?);
                let k = e.kind(self.kinds);
                let r = match e {
                    RealExpr::Add(..) => x + y,
                    RealExpr::Sub(..) => x - y,
                    RealExpr::Mul(..) => x * y,
                    _ => {
                        if y == 0.0 {
                            return Err(RealError::NonFinite("division by zero"));
                        }
                        x / y
                    }
                };
                finite(r, k)?
            }
        })
    }

    /// The truth of one comparison.
    fn cmp(&self, op: CmpOp, l: &RealExpr, r: &RealExpr) -> Result<bool, RealError> {
        let (x, y) = (self.eval(l)?, self.eval(r)?);
        Ok(match op {
            CmpOp::Lt => x < y,
            CmpOp::Le => x <= y,
            CmpOp::Gt => x > y,
            CmpOp::Ge => x >= y,
            CmpOp::Eq => x == y,
            CmpOp::Ne => x != y,
        })
    }

    fn expr_reads_free(&self, e: &RealExpr) -> bool {
        let mut any = false;
        e.for_each_var(&mut |v| any |= self.vars.get(v as usize).copied().flatten().is_none());
        any
    }

    /// The Figure 18-3 value of the predicate `p` (module docs of
    /// [`super::guard`]): `Random` exactly when it reads a variable with no
    /// value yet and nothing decides it.
    fn pred_value(&self, p: &RealCond) -> GuardValue {
        match p {
            RealCond::Const(b) => {
                if *b {
                    GuardValue::True
                } else {
                    GuardValue::False
                }
            }
            RealCond::Int(e) => guard::value(
                e,
                &Assignment {
                    vars: Vec::new(),
                    states: self.states.to_vec(),
                },
            ),
            RealCond::Cmp(op, l, r) => {
                let mut bad = false;
                let mut note = |s: StateId| bad |= self.state(s).is_err();
                l.for_each_state(&mut note);
                r.for_each_state(&mut note);
                if bad {
                    GuardValue::Error
                } else if self.expr_reads_free(l) || self.expr_reads_free(r) {
                    GuardValue::Random
                } else {
                    match self.cmp(*op, l, r) {
                        Ok(true) => GuardValue::True,
                        Ok(false) => GuardValue::False,
                        Err(_) => GuardValue::Error,
                    }
                }
            }
            RealCond::Not(x) => self.pred_value(x).not(),
            RealCond::And(xs) => xs
                .iter()
                .fold(GuardValue::True, |acc, x| acc.and(self.pred_value(x))),
            RealCond::Or(xs) => xs
                .iter()
                .fold(GuardValue::False, |acc, x| acc.or(self.pred_value(x))),
            RealCond::Guard(c, t, e) => match self.pred_value(c) {
                GuardValue::True => self.pred_value(t),
                GuardValue::False => self.pred_value(e),
                GuardValue::Error => GuardValue::Error,
                GuardValue::Random => GuardValue::Random,
            },
        }
    }

    /// Strict truth of a condition that reads no free variable.
    fn truth(&self, c: &RealCond) -> Result<bool, RealError> {
        Ok(match c {
            RealCond::Const(b) => *b,
            RealCond::Int(e) => !int_value(e, self.states)?.is_zero(),
            RealCond::Cmp(op, l, r) => self.cmp(*op, l, r)?,
            RealCond::Not(x) => !self.truth(x)?,
            RealCond::And(xs) => {
                let mut all = true;
                for x in xs {
                    all &= self.truth(x)?;
                }
                all
            }
            RealCond::Or(xs) => {
                let mut any = false;
                for x in xs {
                    any |= self.truth(x)?;
                }
                any
            }
            RealCond::Guard(p, t, e) => match self.pred_value(p) {
                GuardValue::True => self.truth(t)?,
                GuardValue::False => self.truth(e)?,
                GuardValue::Error => return Err(RealError::GuardError),
                GuardValue::Random => {
                    return Err(RealError::Shape("a guard reads a free variable".into()));
                }
            },
        })
    }
}

/// The values of a monotone predicate over the whole grid: `pred` is true on
/// a prefix (`prefix_true`) or a suffix of the indices, and the boundary is
/// found by bisection.
fn monotone_set(
    kind: FloatKind,
    prefix_true: bool,
    mut pred: impl FnMut(f64) -> Result<bool, RealError>,
) -> Result<RealSet, RealError> {
    let m = kind.max_idx();
    let mut at = |i: i64| pred(kind.val(i));
    let (lo_true, hi_true) = (at(-m)?, at(m)?);
    let (first, last) = if prefix_true {
        (lo_true, hi_true)
    } else {
        (hi_true, lo_true)
    };
    // `first` is the predicate at the end where it is expected to hold.
    if !first {
        return Ok(RealSet::empty(kind));
    }
    if last {
        return Ok(RealSet::full(kind));
    }
    // Invariant: `good` satisfies the predicate, `bad` does not.
    let (mut good, mut bad) = if prefix_true { (-m, m) } else { (m, -m) };
    while good.abs_diff(bad) > 1 {
        let mid = ((i128::from(good) + i128::from(bad)) / 2) as i64;
        if at(mid)? {
            good = mid;
        } else {
            bad = mid;
        }
    }
    Ok(if prefix_true {
        RealSet::from_idx(kind, -m, good)
    } else {
        RealSet::from_idx(kind, good, m)
    })
}

impl<'a> Ctx<'a> {
    /// This context with the free variable `v` bound to `x`.
    fn bind(&self, v: RealVarId, x: f64) -> Ctx<'a> {
        let mut vars = self.vars.clone();
        vars[v as usize] = Some(x);
        Ctx {
            kinds: self.kinds,
            vars,
            states: self.states,
        }
    }

    /// The free random variables `c` reads.
    fn free_vars(&self, c: &RealCond) -> Vec<RealVarId> {
        c.vars()
            .into_iter()
            .filter(|v| self.vars.get(*v as usize).copied().flatten().is_none())
            .collect()
    }

    /// The `(slope, offset)` of `e` as an affine function of the free
    /// variable `v`; only the slope's sign and finiteness are used.
    fn lin(&self, e: &RealExpr, v: RealVarId) -> Result<(f64, f64), RealError> {
        let non_affine =
            || RealError::Shape("a product or quotient is not affine in the variable".into());
        Ok(match e {
            RealExpr::Const(r, _) => (0.0, r.get()),
            RealExpr::Var(x) if *x == v => (1.0, 0.0),
            RealExpr::Var(_) | RealExpr::State(..) | RealExpr::FromInt(_) => (0.0, self.eval(e)?),
            RealExpr::Neg(x) => {
                let (c, d) = self.lin(x, v)?;
                (-c, -d)
            }
            RealExpr::ToShort(x) => self.lin(x, v)?,
            RealExpr::Add(a, b) => {
                let ((c1, d1), (c2, d2)) = (self.lin(a, v)?, self.lin(b, v)?);
                (c1 + c2, d1 + d2)
            }
            RealExpr::Sub(a, b) => {
                let ((c1, d1), (c2, d2)) = (self.lin(a, v)?, self.lin(b, v)?);
                (c1 - c2, d1 - d2)
            }
            RealExpr::Mul(a, b) => {
                let ((c1, d1), (c2, d2)) = (self.lin(a, v)?, self.lin(b, v)?);
                if c1 == 0.0 {
                    (d1 * c2, d1 * d2)
                } else if c2 == 0.0 {
                    (d2 * c1, d2 * d1)
                } else {
                    return Err(non_affine());
                }
            }
            RealExpr::Div(a, b) => {
                let ((c1, d1), (c2, d2)) = (self.lin(a, v)?, self.lin(b, v)?);
                if c2 != 0.0 {
                    return Err(non_affine());
                }
                if d2 == 0.0 {
                    return Err(RealError::NonFinite("division by zero"));
                }
                (c1 / d2, d1 / d2)
            }
        })
    }

    /// The set of values of the free variable `v` for which the comparison
    /// `l op r` is true under IEEE evaluation (module docs).
    fn cmp_set(
        &self,
        op: CmpOp,
        l: &RealExpr,
        r: &RealExpr,
        v: RealVarId,
    ) -> Result<RealSet, RealError> {
        let kind = self.kinds[v as usize];
        let mut others = false;
        for e in [l, r] {
            e.for_each_var(&mut |x| {
                others |= x != v && self.vars.get(x as usize).copied().flatten().is_none();
            });
        }
        if others {
            return Err(RealError::Shape(
                "a comparison reads two random real variables".into(),
            ));
        }
        let (c, _) = {
            let ((c1, d1), (c2, d2)) = (self.lin(l, v)?, self.lin(r, v)?);
            (c1 - c2, d1 - d2)
        };
        if !c.is_finite() {
            return Err(RealError::NonFinite("a real coefficient is not finite"));
        }
        let (full, empty) = (RealSet::full(kind), RealSet::empty(kind));
        if c == 0.0 {
            return Ok(if self.bind(v, 0.0).cmp(op, l, r)? {
                full
            } else {
                empty
            });
        }
        let prefix = c > 0.0;
        let lt = || monotone_set(kind, prefix, |x| self.bind(v, x).cmp(CmpOp::Lt, l, r));
        let le = || monotone_set(kind, prefix, |x| self.bind(v, x).cmp(CmpOp::Le, l, r));
        Ok(match op {
            CmpOp::Lt => lt()?,
            CmpOp::Le => le()?,
            CmpOp::Gt => le()?.complement(),
            CmpOp::Ge => lt()?.complement(),
            CmpOp::Eq => le()?.intersect(&lt()?.complement()),
            CmpOp::Ne => le()?.intersect(&lt()?.complement()).complement(),
        })
    }

    /// The values of `v` that satisfy `c` (which reads no free variable
    /// other than `v`), strictly: an unreadable state is an error, but a
    /// guard predicate is evaluated by Figure 18-3.
    fn set_of(&self, c: &RealCond, v: RealVarId) -> Result<RealSet, RealError> {
        let kind = self.kinds[v as usize];
        let all = |b: bool| {
            if b {
                RealSet::full(kind)
            } else {
                RealSet::empty(kind)
            }
        };
        Ok(match c {
            RealCond::Const(b) => all(*b),
            RealCond::Int(_) => all(self.truth(c)?),
            RealCond::Cmp(op, l, r) => {
                if self.expr_reads_free(l) || self.expr_reads_free(r) {
                    self.cmp_set(*op, l, r, v)?
                } else {
                    all(self.cmp(*op, l, r)?)
                }
            }
            RealCond::Not(x) => self.set_of(x, v)?.complement(),
            RealCond::And(xs) => {
                let mut s = all(true);
                for x in xs {
                    s = s.intersect(&self.set_of(x, v)?);
                }
                s
            }
            RealCond::Or(xs) => {
                let mut s = all(false);
                for x in xs {
                    s = s.union(&self.set_of(x, v)?);
                }
                s
            }
            RealCond::Guard(p, t, e) => match self.pred_value(p) {
                GuardValue::True => self.set_of(t, v)?,
                GuardValue::False => self.set_of(e, v)?,
                GuardValue::Error => return Err(RealError::GuardError),
                GuardValue::Random => {
                    let ps = self.pred_set(p, v)?;
                    ps.intersect(&self.set_of(t, v)?)
                        .union(&ps.complement().intersect(&self.set_of(e, v)?))
                }
            },
        })
    }

    /// The values of `v` for which the guard predicate `p` is true (Figure
    /// 18-3 decides what it can; the rest reads `v`).
    fn pred_set(&self, p: &RealCond, v: RealVarId) -> Result<RealSet, RealError> {
        let kind = self.kinds[v as usize];
        match self.pred_value(p) {
            GuardValue::True => Ok(RealSet::full(kind)),
            GuardValue::False => Ok(RealSet::empty(kind)),
            GuardValue::Error => Err(RealError::GuardError),
            GuardValue::Random => match p {
                RealCond::Not(x) => Ok(self.pred_set(x, v)?.complement()),
                RealCond::And(xs) => {
                    let mut s = RealSet::full(kind);
                    for x in xs {
                        s = s.intersect(&self.pred_set(x, v)?);
                    }
                    Ok(s)
                }
                RealCond::Or(xs) => {
                    let mut s = RealSet::empty(kind);
                    for x in xs {
                        s = s.union(&self.pred_set(x, v)?);
                    }
                    Ok(s)
                }
                RealCond::Cmp(op, l, r) => self.cmp_set(*op, l, r, v),
                RealCond::Guard(..) => self.set_of(p, v),
                RealCond::Const(_) | RealCond::Int(_) => {
                    unreachable!("a state-only predicate is never RANDOM")
                }
            },
        }
    }
}

/// A non-negative weight (module docs of [`super::dist`]).
fn weight_of(e: &Expr, states: &[Option<BigUint>]) -> Result<f64, RealError> {
    let n = signed_number(&int_value(e, states)?, e);
    if n < BigInt::zero() {
        return Err(RealError::BadWeight);
    }
    Ok(n.to_f64().unwrap_or(f64::INFINITY))
}

/// A real `dist` evaluated under the current state values.
struct DistPlan {
    /// Each explicit item's values and weight.
    items: Vec<(RealSet, f64)>,
    /// The default's values (everything no item names) and weight.
    default: Option<(RealSet, f64)>,
    /// The values with a positive weight: where the variable may land.
    support: RealSet,
}

fn dist_plan(ctx: &Ctx<'_>, d: &RealDist, kind: FloatKind) -> Result<DistPlan, RealError> {
    let mut items = Vec::with_capacity(d.items.len());
    let mut named = RealSet::empty(kind);
    let mut support = RealSet::empty(kind);
    for it in &d.items {
        let set = match &it.span {
            RealSpan::Point(e) => RealSet::point(kind, ctx.eval(e)?),
            RealSpan::Range { lo, hi } => {
                let lo = lo.as_ref().map_or(Ok(f64::NEG_INFINITY), |e| ctx.eval(e))?;
                let hi = hi.as_ref().map_or(Ok(f64::INFINITY), |e| ctx.eval(e))?;
                RealSet::range(kind, lo, hi)
            }
        };
        let w = weight_of(&it.weight, ctx.states)?;
        named = named.union(&set);
        if w > 0.0 {
            support = support.union(&set);
        }
        items.push((set, w));
    }
    let default = match &d.default {
        Some(e) => {
            let w = weight_of(e, ctx.states)?;
            let region = named.complement();
            if w > 0.0 {
                support = support.union(&region);
            }
            Some((region, w))
        }
        None => None,
    };
    Ok(DistPlan {
        items,
        default,
        support,
    })
}

/// The share of `item` that `f` leaves (module docs: a `:/` weight is spread
/// over the item's length).
fn share(item: &RealSet, f: &RealSet) -> f64 {
    let inter = item.intersect(f);
    if inter.is_empty() {
        return 0.0;
    }
    let m = item.measure();
    if m > 0.0 {
        (inter.measure() / m).min(1.0)
    } else {
        1.0
    }
}

impl DistPlan {
    /// A value of `f` (a non-empty subset of the support) drawn by the
    /// items' weights.
    fn sample(&self, f: &RealSet, seed: &mut u64) -> Result<f64, RealError> {
        let all: Vec<(&RealSet, f64)> = self
            .items
            .iter()
            .map(|(s, w)| (s, *w))
            .chain(self.default.iter().map(|(s, w)| (s, *w)))
            .collect();
        let pick = |weights: &[f64], seed: &mut u64| -> Option<usize> {
            let total: f64 = weights.iter().sum();
            if !total.is_finite() || total <= 0.0 {
                return None;
            }
            let mut t = unit(seed) * total;
            let mut last = None;
            for (i, w) in weights.iter().enumerate() {
                if *w > 0.0 {
                    last = Some(i);
                    if t < *w {
                        return Some(i);
                    }
                    t -= *w;
                }
            }
            last
        };
        let eff: Vec<f64> = all.iter().map(|(s, w)| w * share(s, f)).collect();
        let i = match pick(&eff, seed) {
            Some(i) => i,
            None => {
                // Every positive-weight item is cut to a point: weigh by the
                // nominal weights of the items that still hold a value.
                let nominal: Vec<f64> = all
                    .iter()
                    .map(|(s, w)| if s.intersect(f).is_empty() { 0.0 } else { *w })
                    .collect();
                pick(&nominal, seed).ok_or(RealError::Unsatisfiable)?
            }
        };
        Ok(all[i].0.intersect(f).sample(seed))
    }
}

/// One guard that couples its predicate variable `a` to other variables.
struct Coupled {
    a: RealVarId,
    pred: RealSet,
    then: Vec<(RealVarId, RealSet)>,
    other: Vec<(RealVarId, RealSet)>,
}

impl Coupled {
    fn arm(&self, then: bool) -> &[(RealVarId, RealSet)] {
        if then { &self.then } else { &self.other }
    }
}

/// The conjuncts of a problem, classified by the variables they involve.
#[derive(Default)]
struct Lowered {
    unary: Vec<(RealVarId, RealSet)>,
    coupled: Vec<Coupled>,
    /// Conjuncts reading two or more free random reals, solved together by
    /// [`joint_draw`].
    joint: Vec<RealCond>,
    infeasible: bool,
}

impl Ctx<'_> {
    /// Classifies `c` into `out` (see [`solve`]).
    fn lower(&self, c: &RealCond, out: &mut Lowered) -> Result<(), RealError> {
        let free = self.free_vars(c);
        match c {
            RealCond::Guard(p, t, e) => match self.pred_value(p) {
                GuardValue::True => self.lower(t, out),
                GuardValue::False => self.lower(e, out),
                GuardValue::Error => Err(RealError::GuardError),
                GuardValue::Random => {
                    if free.len() == 1 {
                        out.unary.push((free[0], self.set_of(c, free[0])?));
                        return Ok(());
                    }
                    let pf = self.free_vars(p);
                    let [a] = pf[..] else {
                        out.joint.push(c.clone());
                        return Ok(());
                    };
                    let pred = self.pred_set(p, a)?;
                    let mut arms: [Vec<(RealVarId, RealSet)>; 2] = [Vec::new(), Vec::new()];
                    for (k, arm) in [t, e].into_iter().enumerate() {
                        let mut l = Lowered::default();
                        self.lower(arm, &mut l)?;
                        if !l.coupled.is_empty() {
                            return Err(RealError::Shape(
                                "a guard nested in a coupling guard couples two more variables"
                                    .into(),
                            ));
                        }
                        // The arm applies where the predicate holds (then) or
                        // fails (else); a requirement on `a` itself or a
                        // contradiction is a restriction of `a` alone.
                        let applies = if k == 0 {
                            pred.clone()
                        } else {
                            pred.complement()
                        };
                        let off = applies.complement();
                        if l.infeasible {
                            out.unary.push((a, off.clone()));
                        }
                        for (v, s) in l.unary {
                            if v == a {
                                out.unary.push((a, off.union(&s)));
                            } else {
                                arms[k].push((v, s));
                            }
                        }
                    }
                    let [then, other] = arms;
                    if !then.is_empty() || !other.is_empty() {
                        out.coupled.push(Coupled {
                            a,
                            pred,
                            then,
                            other,
                        });
                    }
                    Ok(())
                }
            },
            _ => match free.len() {
                0 => {
                    if !self.truth(c)? {
                        out.infeasible = true;
                    }
                    Ok(())
                }
                1 => {
                    out.unary.push((free[0], self.set_of(c, free[0])?));
                    Ok(())
                }
                _ => match c {
                    RealCond::And(xs) => xs.iter().try_for_each(|x| self.lower(x, out)),
                    _ => {
                        out.joint.push(c.clone());
                        Ok(())
                    }
                },
            },
        }
    }
}

/// Draws attempted by [`joint_draw`] before it gives up.
const JOINT_TRIES: usize = 20_000;

/// Solves the conjuncts that read several free random reals together (LRM
/// 18.5: every legal solution is equally likely) by rejection sampling: each
/// variable is drawn from its own feasible set (its unary constraints and
/// `dist`), and the draw is kept only when every joint conjunct holds, so the
/// accepted tuple is uniform over the joint feasible region. Any expression
/// shape is accepted (`f < g`, `f + g < 1.0`, `f * g > 2.0`); a region so thin
/// that no draw lands in it (an equality between two random reals has measure
/// zero) is refused loudly, never reported as a solution.
fn joint_draw(
    ctx: &Ctx<'_>,
    joint: &[RealCond],
    own: &[RealSet],
    plans: &[Option<DistPlan>],
    seed: &mut u64,
) -> Result<Vec<(usize, f64)>, RealError> {
    let mut vars: Vec<usize> = joint
        .iter()
        .flat_map(|c| ctx.free_vars(c))
        .map(|v| v as usize)
        .collect();
    vars.sort_unstable();
    vars.dedup();
    if vars.iter().any(|&v| own[v].is_empty()) {
        return Err(RealError::Unsatisfiable);
    }
    let mut probe = Ctx {
        kinds: ctx.kinds,
        vars: ctx.vars.clone(),
        states: ctx.states,
    };
    'tries: for _ in 0..JOINT_TRIES {
        for &v in &vars {
            let x = match &plans[v] {
                Some(plan) => plan.sample(&own[v], seed)?,
                None => own[v].sample(seed),
            };
            probe.vars[v] = Some(x);
        }
        for c in joint {
            match probe.truth(c) {
                Ok(true) => {}
                // A draw whose expression is NaN or divides by zero is not a
                // solution; try another.
                Ok(false) | Err(RealError::NonFinite(_)) => continue 'tries,
                Err(e) => return Err(e),
            }
        }
        return Ok(vars
            .iter()
            .map(|&v| (v, probe.vars[v].unwrap_or(0.0)))
            .collect());
    }
    Err(RealError::Shape(format!(
        "no draw satisfied the constraints that relate several random reals in {JOINT_TRIES} \
         tries (an equality between random reals, or a region too thin to sample)"
    )))
}

/// Whether `a` is ordered before `b` by the transitive closure of `before`.
fn ordered_before(before: &[(RealVarId, RealVarId)], a: RealVarId, b: RealVarId) -> bool {
    let mut seen = vec![a];
    let mut work = vec![a];
    while let Some(x) = work.pop() {
        for &(p, q) in before {
            if p == x && !seen.contains(&q) {
                if q == b {
                    return true;
                }
                seen.push(q);
                work.push(q);
            }
        }
    }
    false
}

/// Solves the real variables of `p` (module docs; `Ok` has `Some` for every
/// variable that was free, i.e. `known[v]` is `None`).
///
/// `known[v]` is the fixed value of an inactive (`rand_mode(0)`) variable.
/// All draws come from `seed`, and none are made when no variable is free.
pub fn solve(
    p: &RealProblem,
    known: &[Option<f64>],
    states: &[Option<BigUint>],
    seed: &mut u64,
) -> Result<Vec<Option<f64>>, RealError> {
    let n = p.kinds.len();
    if known.len() != n {
        return Err(RealError::Shape("one known value per real variable".into()));
    }
    let ctx = Ctx {
        kinds: &p.kinds,
        vars: known.to_vec(),
        states,
    };
    let mut low = Lowered::default();
    for c in &p.conds {
        ctx.lower(c, &mut low)?;
    }
    let mut plans: Vec<Option<DistPlan>> = (0..n).map(|_| None).collect();
    for d in &p.dists {
        let v = d.var as usize;
        if v >= n || plans[v].is_some() {
            return Err(RealError::Shape(
                "a `dist` names an unknown variable, or a variable has two".into(),
            ));
        }
        let plan = dist_plan(&ctx, d, p.kinds[v])?;
        if let Some(x) = known[v]
            && !plan.support.contains(x)
        {
            low.infeasible = true;
        }
        plans[v] = Some(plan);
    }
    if low.infeasible {
        return Err(RealError::Unsatisfiable);
    }

    let active: Vec<usize> = (0..n).filter(|&v| known[v].is_none()).collect();
    let mut own: Vec<RealSet> = (0..n).map(|v| RealSet::full(p.kinds[v])).collect();
    for (v, s) in &low.unary {
        own[*v as usize] = own[*v as usize].intersect(s);
    }
    for &v in &active {
        if let Some(plan) = &plans[v] {
            own[v] = own[v].intersect(&plan.support);
        }
    }

    // Joint conjuncts (several free random reals) are drawn first; their
    // variables may not also take part in a coupling guard.
    let mut out: Vec<Option<f64>> = vec![None; n];
    if !low.joint.is_empty() {
        let drawn = joint_draw(&ctx, &low.joint, &own, &plans, seed)?;
        for g in &low.coupled {
            let mut touched = vec![g.a as usize];
            touched.extend(g.then.iter().chain(&g.other).map(|(b, _)| *b as usize));
            if touched.iter().any(|t| drawn.iter().any(|(v, _)| v == t)) {
                return Err(RealError::Shape(
                    "a random real that a constraint relates to another random real also takes \
                     part in a coupling guard"
                        .into(),
                ));
            }
        }
        for (v, x) in drawn {
            out[v] = Some(x);
        }
    }

    // The coupling forest: each variable has at most one parent.
    let mut parent: Vec<Option<RealVarId>> = vec![None; n];
    for g in &low.coupled {
        for arm in [&g.then, &g.other] {
            for (b, _) in arm {
                match parent[*b as usize] {
                    Some(q) if q != g.a => {
                        return Err(RealError::Shape(format!(
                            "real variable #{b} is constrained by guards on two other variables"
                        )));
                    }
                    _ => parent[*b as usize] = Some(g.a),
                }
                if !ordered_before(&p.before, g.a, *b) {
                    return Err(RealError::Shape(format!(
                        "a guard on real variable #{} constrains #{b}; add `solve` the first \
                         `before` the second",
                        g.a
                    )));
                }
            }
        }
    }
    let depth = |mut v: usize| -> Result<usize, RealError> {
        let mut d = 0;
        while let Some(q) = parent[v] {
            d += 1;
            v = q as usize;
            if d > n {
                return Err(RealError::Shape("a cycle of coupled real guards".into()));
            }
        }
        Ok(d)
    };
    let mut order: Vec<(usize, usize)> = Vec::with_capacity(active.len());
    for &v in &active {
        order.push((depth(v)?, v));
    }
    order.sort_unstable();
    let order: Vec<usize> = order.into_iter().map(|(_, v)| v).collect();

    // Backward: the values of each variable some completion exists for.
    let mut down: Vec<RealSet> = own.clone();
    for &v in order.iter().rev() {
        let guards: Vec<&Coupled> = low.coupled.iter().filter(|g| g.a as usize == v).collect();
        if guards.is_empty() {
            if down[v].is_empty() {
                return Err(RealError::Unsatisfiable);
            }
            continue;
        }
        let mut cells: Vec<(RealSet, Vec<bool>)> = vec![(own[v].clone(), Vec::new())];
        for g in &guards {
            let mut next = Vec::new();
            for (s, sig) in &cells {
                for (then, pset) in [(true, g.pred.clone()), (false, g.pred.complement())] {
                    let t = s.intersect(&pset);
                    if !t.is_empty() {
                        let mut sg = sig.clone();
                        sg.push(then);
                        next.push((t, sg));
                    }
                }
            }
            cells = next;
        }
        let mut children: Vec<usize> = guards
            .iter()
            .flat_map(|g| g.then.iter().chain(&g.other).map(|(b, _)| *b as usize))
            .collect();
        children.sort_unstable();
        children.dedup();
        let mut feasible = RealSet::empty(p.kinds[v]);
        for (s, sig) in &cells {
            let ok = children.iter().all(|&b| {
                let mut need = down[b].clone();
                for (g, &then) in guards.iter().zip(sig) {
                    for (c, set) in g.arm(then) {
                        if *c as usize == b {
                            need = need.intersect(set);
                        }
                    }
                }
                !need.is_empty()
            });
            if ok {
                feasible = feasible.union(s);
            }
        }
        if feasible.is_empty() {
            return Err(RealError::Unsatisfiable);
        }
        down[v] = feasible;
    }

    // Forward: each variable in turn, parents first.
    let mut extra: Vec<RealSet> = (0..n).map(|v| RealSet::full(p.kinds[v])).collect();
    for &v in &order {
        if out[v].is_some() {
            continue;
        }
        let f = down[v].intersect(&extra[v]);
        if f.is_empty() {
            return Err(RealError::Shape(
                "a real variable has no value left after its parents were drawn".into(),
            ));
        }
        let x = match &plans[v] {
            Some(plan) => plan.sample(&f, seed)?,
            None => f.sample(seed),
        };
        out[v] = Some(x);
        for g in low.coupled.iter().filter(|g| g.a as usize == v) {
            for (b, s) in g.arm(g.pred.contains(x)) {
                extra[*b as usize] = extra[*b as usize].intersect(s);
            }
        }
    }
    Ok(out)
}

impl RealProblem {
    /// Checks the shape rules that do not depend on a call's state values or
    /// `rand_mode`s (module docs), so the builder can reject a class whose
    /// real constraints this solver cannot handle at compile time. [`solve`]
    /// re-checks everything it relies on, with the same messages.
    pub fn check(&self) -> Result<(), String> {
        let n = self.kinds.len() as u32;
        let mut edges: Vec<(RealVarId, RealVarId)> = Vec::new();
        let mut joint: Vec<RealVarId> = Vec::new();
        for c in &self.conds {
            check_cond(c, &mut edges, &mut joint)?;
            check_affine(c)?;
        }
        if let Some(&(a, b)) = edges
            .iter()
            .find(|&&(a, b)| joint.contains(&a) || joint.contains(&b))
        {
            return Err(format!(
                "real variable #{a} or #{b} is both related to another random real and part of \
                 a coupling guard"
            ));
        }
        let mut seen = vec![false; self.kinds.len()];
        for d in &self.dists {
            if d.var >= n || std::mem::replace(&mut seen[d.var as usize], true) {
                return Err("a real variable has two `dist` constraints".into());
            }
        }
        if let Some(&(a, b)) = self.before.iter().find(|&&(a, b)| a >= n || b >= n) {
            return Err(format!(
                "`solve` names real variable #{a} or #{b} that does not exist"
            ));
        }
        let mut parent: Vec<Option<RealVarId>> = vec![None; self.kinds.len()];
        for &(a, b) in &edges {
            match parent[b as usize] {
                Some(q) if q != a => {
                    return Err(format!(
                        "real variable #{b} is constrained by guards on two other variables"
                    ));
                }
                _ => parent[b as usize] = Some(a),
            }
            if !ordered_before(&self.before, a, b) {
                return Err(format!(
                    "a guard on real variable #{a} constrains #{b}; add `solve` the first \
                     `before` the second"
                ));
            }
        }
        Ok(())
    }

    /// Validates every integral expression embedded in the problem (the
    /// builder's last line of defence against building a malformed tree).
    pub fn validate_ir(&self) -> Result<(), super::ir::IrError> {
        fn expr(e: &RealExpr) -> Result<(), super::ir::IrError> {
            match e {
                RealExpr::FromInt(x) => x.validate(),
                RealExpr::Const(..) | RealExpr::Var(_) | RealExpr::State(..) => Ok(()),
                RealExpr::Neg(x) | RealExpr::ToShort(x) => expr(x),
                RealExpr::Add(a, b)
                | RealExpr::Sub(a, b)
                | RealExpr::Mul(a, b)
                | RealExpr::Div(a, b) => {
                    expr(a)?;
                    expr(b)
                }
            }
        }
        fn cond(c: &RealCond) -> Result<(), super::ir::IrError> {
            match c {
                RealCond::Const(_) => Ok(()),
                RealCond::Int(x) => x.validate(),
                RealCond::Cmp(_, l, r) => {
                    expr(l)?;
                    expr(r)
                }
                RealCond::Not(x) => cond(x),
                RealCond::And(xs) | RealCond::Or(xs) => xs.iter().try_for_each(cond),
                RealCond::Guard(p, t, e) => {
                    cond(p)?;
                    cond(t)?;
                    cond(e)
                }
            }
        }
        self.conds.iter().try_for_each(cond)?;
        for d in &self.dists {
            for it in &d.items {
                it.weight.validate()?;
                match &it.span {
                    RealSpan::Point(x) => expr(x)?,
                    RealSpan::Range { lo, hi } => {
                        lo.iter().chain(hi).try_for_each(expr)?;
                    }
                }
            }
            if let Some(w) = &d.default {
                w.validate()?;
            }
        }
        Ok(())
    }
}

fn check_cond(
    c: &RealCond,
    edges: &mut Vec<(RealVarId, RealVarId)>,
    joint: &mut Vec<RealVarId>,
) -> Result<(), String> {
    if c.vars().len() <= 1 {
        return Ok(());
    }
    match c {
        RealCond::And(xs) => xs.iter().try_for_each(|x| check_cond(x, edges, joint)),
        RealCond::Guard(p, t, e) => match p.vars()[..] {
            [] => {
                check_cond(t, edges, joint)?;
                check_cond(e, edges, joint)
            }
            [a] => {
                check_arm(t, a, edges)?;
                check_arm(e, a, edges)
            }
            _ => {
                joint.extend(c.vars());
                Ok(())
            }
        },
        _ => {
            joint.extend(c.vars());
            Ok(())
        }
    }
}

/// Rejects, at build time, an atom that reads one random variable
/// non-affinely (`x * x`, `1.0 / x`): the structural form of the rule
/// [`Ctx::lin`] applies at solve time. The arms of a guard are skipped: the
/// variable they read may already be bound by then, which makes the product
/// affine in the one that is not.
fn check_affine(c: &RealCond) -> Result<(), String> {
    fn reads(e: &RealExpr, v: RealVarId) -> bool {
        let mut any = false;
        e.for_each_var(&mut |x| any |= x == v);
        any
    }
    fn affine(e: &RealExpr, v: RealVarId) -> Result<(), String> {
        match e {
            RealExpr::Const(..) | RealExpr::Var(_) | RealExpr::State(..) | RealExpr::FromInt(_) => {
                Ok(())
            }
            RealExpr::Neg(x) | RealExpr::ToShort(x) => affine(x, v),
            RealExpr::Add(a, b) | RealExpr::Sub(a, b) => {
                affine(a, v)?;
                affine(b, v)
            }
            RealExpr::Mul(a, b) => {
                if reads(a, v) && reads(b, v) {
                    return Err(format!(
                        "a product of two factors that both read real variable #{v} is not \
                         affine in it"
                    ));
                }
                affine(a, v)?;
                affine(b, v)
            }
            RealExpr::Div(a, b) => {
                if reads(b, v) {
                    return Err(format!(
                        "a quotient whose denominator reads real variable #{v} is not affine \
                         in it"
                    ));
                }
                affine(a, v)?;
                affine(b, v)
            }
        }
    }
    match c {
        RealCond::Cmp(_, l, r) => match c.vars()[..] {
            [v] => {
                affine(l, v)?;
                affine(r, v)
            }
            _ => Ok(()),
        },
        RealCond::Not(x) => check_affine(x),
        RealCond::And(xs) | RealCond::Or(xs) => xs.iter().try_for_each(check_affine),
        RealCond::Guard(p, _, _) => check_affine(p),
        RealCond::Const(_) | RealCond::Int(_) => Ok(()),
    }
}

fn check_arm(
    c: &RealCond,
    a: RealVarId,
    edges: &mut Vec<(RealVarId, RealVarId)>,
) -> Result<(), String> {
    let vs = c.vars();
    if vs.len() <= 1 {
        if let [b] = vs[..]
            && b != a
        {
            edges.push((a, b));
        }
        return Ok(());
    }
    match c {
        RealCond::And(xs) => xs.iter().try_for_each(|x| check_arm(x, a, edges)),
        RealCond::Guard(p, t, e) if p.vars().is_empty() => {
            check_arm(t, a, edges)?;
            check_arm(e, a, edges)
        }
        RealCond::Guard(..) => {
            Err("a guard nested in a coupling guard couples two more variables".into())
        }
        _ => Err("a conjunct reads two random real variables".into()),
    }
}

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

    const D: FloatKind = FloatKind::F64;
    const S: FloatKind = FloatKind::F32;

    #[test]
    fn the_grid_is_order_preserving_and_round_trips() {
        for k in [D, S] {
            let xs = [
                -k.max_val(),
                -1.5,
                -1e-30,
                0.0,
                1e-30,
                1.0,
                2.5,
                k.max_val(),
            ]
            .map(|x| k.round(x));
            for w in xs.windows(2) {
                assert!(k.idx(w[0]) < k.idx(w[1]), "{k:?} {w:?}");
            }
            for x in xs {
                assert_eq!(k.val(k.idx(x)), x, "{k:?} {x}");
            }
            assert_eq!(k.idx(-0.0), 0);
            assert_eq!(k.val(k.idx(1.0)), 1.0);
        }
        // Adjacent floats are adjacent indices.
        assert_eq!(D.idx(f64::from_bits(1.0f64.to_bits() + 1)), D.idx(1.0) + 1);
        assert_eq!(
            S.val(S.idx(1.0) + 1),
            f64::from(f32::from_bits(1.0f32.to_bits() + 1))
        );
    }

    #[test]
    fn ceil_and_floor_round_the_right_way_for_shortreal() {
        let up = f64::from(f32::from_bits(1.0f32.to_bits() + 1));
        let mid = (1.0 + up) / 2.0;
        assert_eq!(S.ceil_idx(mid), Some(S.idx(1.0) + 1));
        assert_eq!(S.floor_idx(mid), Some(S.idx(1.0)));
        assert_eq!(S.ceil_idx(1.0), Some(S.idx(1.0)));
        assert_eq!(S.floor_idx(1.0), Some(S.idx(1.0)));
        assert_eq!(S.ceil_idx(1e300), None);
        assert_eq!(S.floor_idx(1e300), Some(S.max_idx()));
        assert_eq!(S.ceil_idx(f64::NEG_INFINITY), Some(-S.max_idx()));
        assert_eq!(D.floor_idx(f64::NAN), None);
    }

    #[test]
    fn set_algebra() {
        let a = RealSet::range(D, 0.0, 2.0);
        let b = RealSet::range(D, 1.0, 3.0);
        assert!(a.intersect(&b).contains(1.5) && !a.intersect(&b).contains(2.5));
        assert!(a.union(&b).contains(2.5) && !a.union(&b).contains(3.5));
        let c = a.complement();
        assert!(!c.contains(1.0) && c.contains(-1.0) && c.contains(2.5));
        assert_eq!(c.complement(), a);
        assert!(RealSet::full(D).complement().is_empty());
        assert_eq!(RealSet::empty(D).complement(), RealSet::full(D));
        // Adjacent floats merge into one interval.
        let p = RealSet::point(D, 1.0);
        let q = RealSet::point(D, f64::from_bits(1.0f64.to_bits() + 1));
        assert_eq!(p.union(&q).intervals().len(), 1);
        // A hole.
        let hole = RealSet::point(D, 1.0).complement();
        assert!(!hole.contains(1.0) && hole.contains(0.999) && hole.contains(1.001));
        // A shortreal never equals a real it cannot represent.
        assert!(RealSet::point(S, 0.1).is_empty());
        assert!(!RealSet::point(S, 0.5).is_empty());
        assert!(RealSet::range(D, 2.0, 1.0).is_empty());
    }

    #[test]
    fn measure_does_not_overflow_on_the_full_range() {
        let m = RealSet::full(D).measure();
        assert!(m.is_finite() && m > 1e308);
        assert_eq!(RealSet::point(D, 3.0).measure(), 0.0);
        assert_eq!(RealSet::range(D, 1.0, 3.0).measure(), 1.0);
    }

    #[test]
    fn sampling_stays_in_the_set_and_is_uniform_by_length() {
        let set = RealSet::range(D, 0.0, 2.0).union(&RealSet::range(D, 10.0, 12.0));
        let mut seed = 7u64;
        let (mut lo, mut hi, mut n_a) = (0u32, 0u32, 0u32);
        for _ in 0..4000 {
            let x = set.sample(&mut seed);
            assert!(set.contains(x), "{x}");
            if x < 5.0 {
                n_a += 1;
            }
            if x % 10.0 < 1.0 {
                lo += 1;
            } else {
                hi += 1;
            }
        }
        // Two equal-length intervals, and each is uniform inside.
        assert!((1850..=2150).contains(&n_a), "{n_a}");
        assert!((1850..=2150).contains(&lo), "{lo} vs {hi}");
        // Full range: finite and both signs.
        let full = RealSet::full(D);
        let xs: Vec<f64> = (0..200).map(|_| full.sample(&mut seed)).collect();
        assert!(xs.iter().all(|x| x.is_finite()));
        assert!(xs.iter().any(|&x| x < 0.0) && xs.iter().any(|&x| x > 0.0));
    }

    #[test]
    fn a_point_set_is_sampled_over_its_points() {
        let set = RealSet::point(D, 3.0).union(&RealSet::point(D, 5.0));
        let mut seed = 1u64;
        let mut n3 = 0;
        for _ in 0..1000 {
            let x = set.sample(&mut seed);
            assert!(x == 3.0 || x == 5.0);
            n3 += usize::from(x == 3.0);
        }
        assert!((400..=600).contains(&n3), "{n3}");
        // A shortreal range is sampled on the f32 grid.
        let sr = RealSet::range(S, 0.0, 1.0);
        for _ in 0..200 {
            let x = sr.sample(&mut seed);
            assert_eq!(f64::from(x as f32), x);
            assert!((0.0..=1.0).contains(&x));
        }
    }

    // ---- the constraint solver ------------------------------------------

    fn v(i: u32) -> RealExpr {
        RealExpr::Var(i)
    }
    fn k(x: f64) -> RealExpr {
        RealExpr::Const(Real::new(x), D)
    }
    fn bx(e: RealExpr) -> Box<RealExpr> {
        Box::new(e)
    }
    fn cmp(op: CmpOp, l: RealExpr, r: RealExpr) -> RealCond {
        RealCond::Cmp(op, l, r)
    }
    fn imp(p: RealCond, t: RealCond) -> RealCond {
        RealCond::Guard(Box::new(p), Box::new(t), Box::new(RealCond::Const(true)))
    }
    fn problem(n: usize, conds: Vec<RealCond>) -> RealProblem {
        RealProblem {
            kinds: vec![D; n],
            conds,
            ..RealProblem::default()
        }
    }
    fn run(p: &RealProblem, seed: &mut u64) -> Result<Vec<f64>, RealError> {
        let known = vec![None; p.kinds.len()];
        solve(p, &known, &[], seed).map(|v| v.into_iter().map(|x| x.unwrap()).collect())
    }
    fn set_of(p: &RealProblem, c: &RealCond) -> RealSet {
        let ctx = Ctx {
            kinds: &p.kinds,
            vars: vec![None; p.kinds.len()],
            states: &[],
        };
        ctx.set_of(c, 0).unwrap()
    }

    #[test]
    fn an_atom_is_exactly_the_floats_for_which_the_rounded_comparison_holds() {
        // x / 3.0 < 1.0: the boundary is where the *rounded* quotient reaches
        // 1.0, not at 3.0.
        let p = problem(1, vec![]);
        let c = cmp(CmpOp::Lt, RealExpr::Div(bx(v(0)), bx(k(3.0))), k(1.0));
        let s = set_of(&p, &c);
        let (_, top) = *s.intervals().last().unwrap();
        assert!(D.val(top) / 3.0 < 1.0);
        assert!(D.val(top + 1) / 3.0 >= 1.0);
        assert_eq!(s.intervals().len(), 1);
        // A decreasing form: 5.0 - x <= 2.0 holds for x >= 3.0.
        let c = cmp(CmpOp::Le, RealExpr::Sub(bx(k(5.0)), bx(v(0))), k(2.0));
        let s = set_of(&p, &c);
        assert!(s.contains(3.0) && !s.contains(f64::from_bits(3.0f64.to_bits() - 1)));
        // Eq is a single point, Ne its complement, Gt/Ge the open/closed tails.
        let eq = set_of(&p, &cmp(CmpOp::Eq, v(0), k(2.5)));
        assert_eq!(eq, RealSet::point(D, 2.5));
        assert_eq!(set_of(&p, &cmp(CmpOp::Ne, v(0), k(2.5))), eq.complement());
        let gt = set_of(&p, &cmp(CmpOp::Gt, v(0), k(2.5)));
        let ge = set_of(&p, &cmp(CmpOp::Ge, v(0), k(2.5)));
        assert!(!gt.contains(2.5) && gt.contains(f64::from_bits(2.5f64.to_bits() + 1)));
        assert!(ge.contains(2.5) && !ge.contains(f64::from_bits(2.5f64.to_bits() - 1)));
        // A negative coefficient flips the direction: -2.0 * x > 4.0 is x < -2.0.
        let c = cmp(CmpOp::Gt, RealExpr::Mul(bx(k(-2.0)), bx(v(0))), k(4.0));
        let s = set_of(&p, &c);
        assert!(s.contains(-2.0000000000000004) && !s.contains(-2.0));
    }

    #[test]
    fn a_shortreal_atom_lives_on_the_binary32_grid() {
        let p = RealProblem {
            kinds: vec![S],
            ..RealProblem::default()
        };
        let ctx = Ctx {
            kinds: &p.kinds,
            vars: vec![None],
            states: &[],
        };
        let s = ctx.set_of(&cmp(CmpOp::Lt, v(0), k(1.0)), 0).unwrap();
        let below = f64::from(f32::from_bits(1.0f32.to_bits() - 1));
        assert!(s.contains(below) && !s.contains(1.0));
        // Values handed out by the sampler are binary32 numbers.
        let mut seed = 3;
        for _ in 0..50 {
            let x = s.sample(&mut seed);
            assert_eq!(f64::from(x as f32), x);
        }
    }

    #[test]
    fn an_open_interval_is_uniform() {
        // 0.0 < v < 2.0
        let p = problem(
            1,
            vec![RealCond::And(vec![
                cmp(CmpOp::Gt, v(0), k(0.0)),
                cmp(CmpOp::Lt, v(0), k(2.0)),
            ])],
        );
        let mut seed = 42u64;
        let mut bins = [0u32; 8];
        for _ in 0..8000 {
            let x = run(&p, &mut seed).unwrap()[0];
            assert!(x > 0.0 && x < 2.0, "{x}");
            bins[(x * 4.0) as usize] += 1;
        }
        for (i, b) in bins.iter().enumerate() {
            assert!((850..=1150).contains(b), "bin {i}: {bins:?}");
        }
    }

    #[test]
    fn an_unsatisfiable_set_is_reported_and_draws_nothing() {
        let p = problem(
            1,
            vec![cmp(CmpOp::Gt, v(0), k(2.0)), cmp(CmpOp::Lt, v(0), k(1.0))],
        );
        let mut seed = 5;
        assert_eq!(run(&p, &mut seed), Err(RealError::Unsatisfiable));
        assert_eq!(seed, 5);
        // No free variable: nothing is drawn and the result is empty.
        let q = problem(1, vec![cmp(CmpOp::Lt, v(0), k(9.0))]);
        let r = solve(&q, &[Some(1.0)], &[], &mut seed).unwrap();
        assert_eq!(r, vec![None]);
        assert_eq!(seed, 5);
        // ... and a fixed value that breaks the constraint is unsatisfiable.
        assert_eq!(
            solve(&q, &[Some(10.0)], &[], &mut seed),
            Err(RealError::Unsatisfiable)
        );
    }

    #[test]
    fn a_state_bound_is_read_per_call_and_an_unreadable_one_is_an_error() {
        let p = problem(1, vec![cmp(CmpOp::Lt, v(0), RealExpr::State(0, D))]);
        let mut seed = 1;
        let st = [Some(BigUint::from((-4.5f64).to_bits()))];
        let r = solve(&p, &[None], &st, &mut seed).unwrap();
        assert!(r[0].unwrap() < -4.5);
        assert_eq!(
            solve(&p, &[None], &[None], &mut seed),
            Err(RealError::StateUnreadable(0))
        );
    }

    #[test]
    fn nan_and_division_by_zero_are_errors_not_silent_values() {
        let p = problem(
            1,
            vec![cmp(CmpOp::Lt, RealExpr::Div(bx(v(0)), bx(k(0.0))), k(1.0))],
        );
        let mut seed = 1;
        assert!(matches!(run(&p, &mut seed), Err(RealError::NonFinite(_))));
    }

    #[test]
    fn a_non_affine_or_two_variable_conjunct_is_rejected() {
        let sq = cmp(CmpOp::Lt, RealExpr::Mul(bx(v(0)), bx(v(0))), k(4.0));
        let mut seed = 1;
        assert!(matches!(
            run(&problem(1, vec![sq.clone()]), &mut seed),
            Err(RealError::Shape(_))
        ));
        // `check` (the build-time gate) rejects it too, by name.
        let why = problem(1, vec![sq]).check().unwrap_err();
        assert!(why.contains("not affine") && why.contains("#0"), "{why}");
        let recip = cmp(CmpOp::Gt, RealExpr::Div(bx(k(1.0)), bx(v(0))), k(2.0));
        assert!(
            problem(1, vec![recip])
                .check()
                .unwrap_err()
                .contains("denominator")
        );
        // A product with one constant factor, and one whose other factor is
        // another variable (bound by the time it is read), stay legal.
        let scaled = cmp(CmpOp::Lt, RealExpr::Mul(bx(k(2.0)), bx(v(0))), k(4.0));
        assert!(problem(1, vec![scaled]).check().is_ok());
    }

    #[test]
    fn conjuncts_relating_two_random_reals_are_solved_jointly() {
        // LRM 18.5: `f < g` and `f + g < 1.0` over unit ranges. Every draw
        // satisfies both; the draws spread over the whole feasible triangle
        // (f gets near 0, g gets near 1).
        let unit_range = |i: u32| [cmp(CmpOp::Ge, v(i), k(0.0)), cmp(CmpOp::Le, v(i), k(1.0))];
        let mut conds: Vec<RealCond> = unit_range(0).into_iter().chain(unit_range(1)).collect();
        conds.push(cmp(CmpOp::Lt, v(0), v(1)));
        conds.push(cmp(CmpOp::Lt, RealExpr::Add(bx(v(0)), bx(v(1))), k(1.0)));
        let p = problem(2, conds);
        assert!(p.check().is_ok());
        let mut seed = 5;
        let (mut lo_f, mut hi_g) = (f64::MAX, f64::MIN);
        for _ in 0..500 {
            let r = run(&p, &mut seed).unwrap();
            assert!(r[0] < r[1] && r[0] + r[1] < 1.0, "{r:?}");
            assert!((0.0..=1.0).contains(&r[0]) && (0.0..=1.0).contains(&r[1]));
            lo_f = lo_f.min(r[0]);
            hi_g = hi_g.max(r[1]);
        }
        assert!(lo_f < 0.05 && hi_g > 0.8, "{lo_f} {hi_g}");
        // An equality between random reals has measure zero: refused by name.
        let eq = problem(2, vec![cmp(CmpOp::Eq, v(0), v(1))]);
        assert!(eq.check().is_ok());
        let mut seed = 1;
        assert!(matches!(run(&eq, &mut seed), Err(RealError::Shape(_))));
        // A contradiction is unsatisfiable only through exhausting the draws,
        // which is also a loud Shape error, never a solution.
        let bad = problem(
            2,
            vec![cmp(CmpOp::Lt, v(0), v(1)), cmp(CmpOp::Gt, v(0), v(1))],
        );
        assert!(matches!(run(&bad, &mut seed), Err(RealError::Shape(_))));
    }

    #[test]
    fn a_guard_with_a_state_predicate_follows_figure_18_3() {
        // (state == 1) -> v > 5.0
        let g = |st: Option<u64>| {
            let p = problem(
                1,
                vec![imp(
                    RealCond::Int(Expr::Binary(
                        super::super::ir::BinOp::Eq,
                        Box::new(Expr::State(0, super::super::ir::Ty::unsigned(8))),
                        Box::new(Expr::constant(super::super::ir::Ty::unsigned(8), 1u64)),
                        super::super::ir::Ty::BIT,
                    )),
                    cmp(CmpOp::Gt, v(0), k(5.0)),
                )],
            );
            let mut seed = 9;
            solve(&p, &[None], &[st.map(BigUint::from)], &mut seed)
        };
        assert!(g(Some(1)).unwrap()[0].unwrap() > 5.0);
        assert!(g(Some(2)).is_ok());
        assert_eq!(g(None), Err(RealError::GuardError));
    }

    /// `a` in [0, 10]; `(a < 5) -> b > 20`; `b < 10`.
    fn coupled() -> RealProblem {
        RealProblem {
            kinds: vec![D, D],
            conds: vec![
                RealCond::And(vec![
                    cmp(CmpOp::Ge, v(0), k(0.0)),
                    cmp(CmpOp::Le, v(0), k(10.0)),
                ]),
                imp(cmp(CmpOp::Lt, v(0), k(5.0)), cmp(CmpOp::Gt, v(1), k(20.0))),
                cmp(CmpOp::Lt, v(1), k(10.0)),
            ],
            dists: vec![],
            before: vec![(0, 1)],
        }
    }

    #[test]
    fn a_coupled_guard_needs_its_order_and_is_solved_staged() {
        let mut p = coupled();
        assert!(p.check().is_ok());
        let mut seed = 11;
        // `a < 5` would force `b > 20` against `b < 10`, so only `a >= 5`
        // completes: `a` is drawn from exactly those values.
        for _ in 0..300 {
            let r = run(&p, &mut seed).unwrap();
            assert!((5.0..=10.0).contains(&r[0]), "{r:?}");
            assert!(r[1] < 10.0);
        }
        p.before.clear();
        let e = p.check().unwrap_err();
        assert!(e.contains("before"), "{e}");
        assert!(matches!(run(&p, &mut seed), Err(RealError::Shape(_))));
    }

    #[test]
    fn a_coupled_guard_draws_the_child_from_the_arm_the_parent_chose() {
        // a in [0,10]; (a < 5) -> b in [20,30]; !(a < 5) -> b in [-30,-20]
        let both = |lo: f64, hi: f64| {
            RealCond::And(vec![
                cmp(CmpOp::Ge, v(1), k(lo)),
                cmp(CmpOp::Le, v(1), k(hi)),
            ])
        };
        let p = RealProblem {
            kinds: vec![D, D],
            conds: vec![
                RealCond::And(vec![
                    cmp(CmpOp::Ge, v(0), k(0.0)),
                    cmp(CmpOp::Le, v(0), k(10.0)),
                ]),
                RealCond::Guard(
                    Box::new(cmp(CmpOp::Lt, v(0), k(5.0))),
                    Box::new(both(20.0, 30.0)),
                    Box::new(both(-30.0, -20.0)),
                ),
            ],
            dists: vec![],
            before: vec![(0, 1)],
        };
        assert!(p.check().is_ok());
        let mut seed = 77;
        let mut low = 0;
        for _ in 0..2000 {
            let r = run(&p, &mut seed).unwrap();
            assert_eq!(r[0] < 5.0, r[1] > 0.0, "{r:?}");
            assert!(r[1].abs() >= 20.0 && r[1].abs() <= 30.0);
            low += usize::from(r[0] < 5.0);
        }
        // Table 18-2's rule: `a` is drawn first, so it is not skewed by `b`.
        assert!((900..=1100).contains(&low), "{low}");
    }

    #[test]
    fn coupling_shapes_the_solver_cannot_do_are_named() {
        // a -> b and c -> b: b has two parents.
        let p = RealProblem {
            kinds: vec![D, D, D],
            conds: vec![
                imp(cmp(CmpOp::Lt, v(0), k(5.0)), cmp(CmpOp::Gt, v(1), k(0.0))),
                imp(cmp(CmpOp::Lt, v(2), k(5.0)), cmp(CmpOp::Gt, v(1), k(1.0))),
            ],
            dists: vec![],
            before: vec![(0, 1), (2, 1)],
        };
        assert!(p.check().unwrap_err().contains("two other variables"));
        // A predicate over two variables is a joint conjunct (rejection
        // sampling), not a coupling guard.
        let q = problem(
            3,
            vec![imp(
                RealCond::And(vec![
                    cmp(CmpOp::Lt, v(0), k(5.0)),
                    cmp(CmpOp::Lt, v(1), k(5.0)),
                ]),
                cmp(CmpOp::Gt, v(2), k(0.0)),
            )],
        );
        assert!(q.check().is_ok());
        let mut seed = 3;
        for _ in 0..100 {
            let s = run(&q, &mut seed).unwrap();
            assert!(!(s[0] < 5.0 && s[1] < 5.0) || s[2] > 0.0, "{s:?}");
        }
        // A nested coupling guard.
        let r = RealProblem {
            before: vec![(0, 1), (1, 2)],
            ..problem(
                3,
                vec![imp(
                    cmp(CmpOp::Lt, v(0), k(5.0)),
                    imp(cmp(CmpOp::Lt, v(1), k(5.0)), cmp(CmpOp::Gt, v(2), k(0.0))),
                )],
            )
        };
        assert!(r.check().unwrap_err().contains("nested"));
    }

    fn weight(n: u64) -> Expr {
        Expr::constant(super::super::ir::Ty::unsigned(8), n)
    }

    #[test]
    fn a_real_dist_weights_points_and_spreads_ranges_over_their_length() {
        // a dist { 1.0 := 1, [2.0:4.0] :/ 3 }, a < 3.0: the range keeps half
        // its length, so P(a == 1.0) = 1 / (1 + 1.5).
        let dist = RealDist {
            var: 0,
            items: vec![
                RealDistItem {
                    span: RealSpan::Point(k(1.0)),
                    weight: weight(1),
                },
                RealDistItem {
                    span: RealSpan::Range {
                        lo: Some(k(2.0)),
                        hi: Some(k(4.0)),
                    },
                    weight: weight(3),
                },
            ],
            default: None,
        };
        let mut p = problem(1, vec![cmp(CmpOp::Lt, v(0), k(3.0))]);
        p.dists.push(dist);
        let mut seed = 21;
        let (mut ones, n) = (0, 6000);
        for _ in 0..n {
            let x = run(&p, &mut seed).unwrap()[0];
            assert!(x == 1.0 || (2.0..3.0).contains(&x), "{x}");
            ones += usize::from(x == 1.0);
        }
        let want = f64::from(n as u32) * 0.4;
        assert!((ones as f64 - want).abs() < 150.0, "{ones} vs {want}");
    }

    #[test]
    fn a_dist_default_covers_the_complement_and_zero_weights_are_excluded() {
        // a dist { [0.0:1.0] :/ 0, default :/ 1 }, a in [0, 2]: never in [0, 1].
        let dist = RealDist {
            var: 0,
            items: vec![RealDistItem {
                span: RealSpan::Range {
                    lo: Some(k(0.0)),
                    hi: Some(k(1.0)),
                },
                weight: weight(0),
            }],
            default: Some(weight(1)),
        };
        let mut p = problem(
            1,
            vec![RealCond::And(vec![
                cmp(CmpOp::Ge, v(0), k(0.0)),
                cmp(CmpOp::Le, v(0), k(2.0)),
            ])],
        );
        p.dists.push(dist);
        let mut seed = 2;
        for _ in 0..500 {
            let x = run(&p, &mut seed).unwrap()[0];
            assert!(x > 1.0 && x <= 2.0, "{x}");
        }
        // A `$` bound is the type's extreme.
        let dollar = RealDist {
            var: 0,
            items: vec![RealDistItem {
                span: RealSpan::Range {
                    lo: Some(k(5.0)),
                    hi: None,
                },
                weight: weight(1),
            }],
            default: None,
        };
        let mut q = problem(1, vec![]);
        q.dists.push(dollar);
        for _ in 0..50 {
            assert!(run(&q, &mut seed).unwrap()[0] >= 5.0);
        }
        // A fixed value outside every positive-weight item breaks the dist.
        assert_eq!(
            solve(&q, &[Some(1.0)], &[], &mut seed),
            Err(RealError::Unsatisfiable)
        );
    }

    #[test]
    fn a_point_in_a_range_item_falls_back_to_the_nominal_weights() {
        // a == 3.0 with a dist over [0:10]: the range keeps zero length but
        // its one point, so the nominal weight decides and the call succeeds.
        let dist = RealDist {
            var: 0,
            items: vec![RealDistItem {
                span: RealSpan::Range {
                    lo: Some(k(0.0)),
                    hi: Some(k(10.0)),
                },
                weight: weight(1),
            }],
            default: None,
        };
        let mut p = problem(1, vec![cmp(CmpOp::Eq, v(0), k(3.0))]);
        p.dists.push(dist);
        let mut seed = 4;
        assert_eq!(run(&p, &mut seed).unwrap(), vec![3.0]);
    }

    #[test]
    fn a_negative_dist_weight_is_an_error() {
        let neg = Expr::constant(super::super::ir::Ty::signed(8), 255u64);
        let mut p = problem(1, vec![]);
        p.dists.push(RealDist {
            var: 0,
            items: vec![RealDistItem {
                span: RealSpan::Point(k(1.0)),
                weight: neg,
            }],
            default: None,
        });
        let mut seed = 1;
        assert_eq!(run(&p, &mut seed), Err(RealError::BadWeight));
    }

    #[test]
    fn an_integral_state_converts_to_real_signed_or_unsigned() {
        use super::super::ir::Ty;
        let zstate = RealExpr::FromInt(Box::new(Expr::State(0, Ty::signed(32))));
        let ctx = Ctx {
            kinds: &[],
            vars: vec![],
            states: &[Some(BigUint::from(0xffff_ff9cu64))],
        };
        assert_eq!(ctx.eval(&zstate).unwrap(), -100.0);
        let unsigned = RealExpr::FromInt(Box::new(Expr::State(0, Ty::unsigned(32))));
        assert_eq!(ctx.eval(&unsigned).unwrap(), f64::from(0xffff_ff9cu32));
    }

    #[test]
    fn the_lrm_real_constraint_example_is_solved_with_its_ordering() {
        // IEEE 1800-2023 18.5.13: with ZSTATE=-100, LOW=0.7, MIN=1.43,
        // NOM=3.3, MAX=3.65.
        let (zs, low, min, nom, max) = (-100.0, 0.7, 1.43, 3.3, 3.65);
        let range = |lo: f64, hi: f64, w: u64| RealDistItem {
            span: RealSpan::Range {
                lo: Some(k(lo)),
                hi: Some(k(hi)),
            },
            weight: weight(w),
        };
        let point = |x: f64, w: u64| RealDistItem {
            span: RealSpan::Point(k(x)),
            weight: weight(w),
        };
        let p = RealProblem {
            kinds: vec![D, D],
            conds: vec![imp(
                RealCond::And(vec![
                    cmp(CmpOp::Ge, v(0), k(low)),
                    cmp(CmpOp::Le, v(0), k(min)),
                ]),
                cmp(CmpOp::Eq, v(1), k(zs)),
            )],
            dists: vec![
                RealDist {
                    var: 0,
                    items: vec![
                        point(zs, 5),
                        range(low, min, 1),
                        range(nom - nom * 0.01, nom + nom * 0.01, 13),
                        range(min, max, 1),
                    ],
                    default: None,
                },
                RealDist {
                    var: 1,
                    items: vec![point(zs, 1), range(min, max, 20)],
                    default: None,
                },
            ],
            before: vec![(0, 1)],
        };
        assert!(p.check().is_ok());
        let mut seed = 99;
        let n = 12_000;
        let (mut a_z, mut b_z_when_low) = (0, (0, 0));
        for _ in 0..n {
            let r = run(&p, &mut seed).unwrap();
            let (a, b) = (r[0], r[1]);
            assert!(a == zs || (low..=max).contains(&a), "{a}");
            assert!(b == zs || (min..=max).contains(&b), "{b}");
            if (low..=min).contains(&a) {
                assert_eq!(b, zs, "{a} {b}");
                b_z_when_low.0 += 1;
            } else if b == zs {
                b_z_when_low.1 += 1;
            }
            a_z += usize::from(a == zs);
        }
        // a == ZSTATE has weight 5 of 5 + (1.43-0.7)/... : lengths enter as
        // shares, and the items overlap, so only the point's weight is exact.
        let total = 5.0 + 1.0 + 13.0 + 1.0;
        let want = f64::from(n as u32) * 5.0 / total;
        assert!((a_z as f64 - want).abs() < 250.0, "{a_z} vs {want}");
        // Outside the guarded range b is ZSTATE with weight 1 against 20.
        assert!(b_z_when_low.0 > 0 && b_z_when_low.1 > 0);
    }

    #[test]
    fn check_rejects_a_second_dist_on_one_variable() {
        let d = |var| RealDist {
            var,
            items: vec![],
            default: None,
        };
        let mut p = problem(1, vec![]);
        p.dists = vec![d(0), d(0)];
        assert!(p.check().is_err());
        p.dists = vec![d(0)];
        assert!(p.check().is_ok());
    }
}