vj_sim/constraints/

solve.rs

//! The solver orchestrator: constraint IR to BDD to exact model count to an
//! exactly uniform sample (plan C3; LRM 18.4.1).
//!
//! # What "uniform" means here
//!
//! [`Solver::sample`] draws one assignment uniformly at random from *all*
//! legal combinations of the random variables: every model of the constraint
//! set has probability `1 / count`. It walks the variables in a caller-supplied
//! order, most significant bit first, and fixes each bit with probability
//! `models-with-that-bit / models-remaining`, so the product along any path
//! is exactly `1 / count` whatever the order is. (This is the probability
//! LRM 18.5.9's Table 18-1 states for the *unordered* case.) `solve ...
//! before` (C4) changes which *marginal* is uniform:
//! [`Solver::sample_groups`] samples each ordered set uniformly over its own
//! feasible values, built on the same per-bit machinery.
//!
//! # Counts
//!
//! The plan sketched a `u128` count. Two 64-bit random variables already have
//! up to 2^128 combinations, so [`SatCount`] is an arbitrary-precision
//! integer; nothing in the sampler can silently wrap.
//!
//! # Randomness
//!
//! Draws come from the `&mut u64` the caller passes (for `randomize()`, the
//! object's own PRNG seed; see [`super::rng`]). An unsatisfiable problem
//! consumes no draws, so a failed `randomize()` does not perturb the object's
//! random stream (LRM 18.6.3).

use super::bdd::{BddBuildError, Compiled, NODE_BUDGET, compile};
use super::ir::{Expr, Ty, VarId};
use super::randc::{self, RandcState, fingerprint_of};
use super::rng;
use super::sat::{SatLimits, SatSolver};
use biodivine_lib_bdd::{Bdd, BddValuation};
use core::fmt;
use num_bigint::BigUint;
use num_traits::{One, Zero};
use std::collections::{BTreeMap, BTreeSet, HashMap};
use std::sync::{Arc, Mutex};

/// A number of models (assignments of every random-variable bit).
pub type SatCount = BigUint;

/// Why a solve failed.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum SolveError {
    /// The constraints could not be turned into a BDD or, for an [`Engine`]
    /// that fell back, into a CNF (including [`BddBuildError::TooLarge`] from
    /// [`Solver`], the SAT-fallback trigger).
    Build(BddBuildError),
    /// No assignment satisfies the constraints (LRM 18.6.3: `randomize()`
    /// returns 0 and changes nothing).
    Unsatisfiable,
    /// The caller's variable order names an unknown or repeated variable.
    BadOrder {
        /// The offending entry.
        var: VarId,
    },
    /// The [`Plan`] is malformed in a way other than naming a bad variable.
    BadPlan {
        /// What is wrong.
        what: &'static str,
    },
    /// The SAT fallback ([`super::sat`]) did not finish within its time limit.
    SatTimeout,
    /// The plan asks the SAT fallback for something it cannot sample
    /// faithfully (a `dist` weighted by a count of models or of distinct
    /// values, a `randc` variable too wide to enumerate); refused instead of
    /// silently mis-weighted.
    SatUnsupported {
        /// What is not supported.
        what: &'static str,
    },
    /// A model the SAT solver returned failed the IR's own evaluator: a bug
    /// in the CNF encoding, never a legal outcome.
    SatModelRejected,
    /// The SAT solver itself failed (out of memory, an internal error) or
    /// contradicted itself; never a legal outcome.
    SatSolverFailure {
        /// What went wrong.
        what: &'static str,
    },
}

/// An exact non-negative rational weight `num / den` (`den` is never 0).
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Weight {
    /// The numerator.
    pub num: BigUint,
    /// The denominator.
    pub den: BigUint,
}

impl Weight {
    /// The integer weight `n`.
    pub fn int(n: impl Into<BigUint>) -> Weight {
        Weight {
            num: n.into(),
            den: BigUint::one(),
        }
    }
}

/// How much of the *rest of the problem* a code of a [`Weighted`] variable
/// stands for, which scales the code's weight.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Measure {
    /// Every feasible code counts once.
    Flat,
    /// A code counts once per feasible *tuple* (model) it leaves.
    Models,
    /// A code counts once per distinct feasible value of this variable (a
    /// member of the same [`Group`]) it leaves, and once the code is drawn that
    /// variable is itself drawn, uniformly over those values.
    Distinct(VarId),
}

/// A variable drawn by weight instead of uniformly (a `dist` selector, C6):
/// code `c` is chosen with probability proportional to
/// `weights[c] * measure(c)` among the codes that are still feasible; a code
/// with weight 0, or beyond `weights.len()`, is never chosen.
///
/// With [`Measure::Distinct`] this is exact for LRM 18.5.3: the selector's
/// code is an item of the `dist`, its weight the item's per-value weight, and
/// the measure the number of values of the subject that are both in the item
/// and still feasible, so the probability of a feasible value is proportional
/// to the sum of the per-value weights of the items that contain it, however
/// the items overlap or the other constraints exclude values.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Weighted {
    /// The variable.
    pub var: VarId,
    /// The weight of each value of the variable, indexed by value.
    pub weights: Vec<Weight>,
    /// What a code stands for.
    pub measure: Measure,
}

/// One set of variables solved together ([`Plan`]).
#[derive(Clone, Debug, Default, PartialEq, Eq)]
pub struct Group {
    /// Variables drawn uniformly over the group's feasible tuples.
    pub vars: Vec<VarId>,
    /// Variables drawn by weight, before `vars`.
    pub weighted: Vec<Weighted>,
}

/// How [`Solver::sample_plan`] solves: the `randc` variables first (LRM
/// 18.4.2), then each ordered [`Group`] in turn (18.5.9).
#[derive(Clone, Debug, Default, PartialEq, Eq)]
pub struct Plan {
    /// The ordered groups.
    pub groups: Vec<Group>,
    /// The `randc` variables, in the order they are drawn.
    pub randc: Vec<VarId>,
}

/// A [`Weighted`] resolved to bit-vector indices.
pub(crate) struct ResolvedWeighted<'a> {
    pub(crate) var: usize,
    pub(crate) weights: &'a [Weight],
    pub(crate) measure: Measure,
}

/// A [`Group`] resolved to bit-vector indices: its uniform variables and its
/// weighted ones.
pub(crate) type Resolved<'a> = (Vec<usize>, Vec<ResolvedWeighted<'a>>);

impl fmt::Display for SolveError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            SolveError::Build(e) => fmt::Display::fmt(e, f),
            SolveError::Unsatisfiable => f.write_str("the constraints have no solution"),
            SolveError::BadOrder { var } => write!(
                f,
                "solve order names random variable #{var} which is unknown or listed twice"
            ),
            SolveError::BadPlan { what } => write!(f, "malformed solve plan: {what}"),
            SolveError::SatTimeout => {
                f.write_str("the SAT fallback did not finish within its time limit")
            }
            SolveError::SatUnsupported { what } => {
                write!(f, "not supported by the SAT fallback: {what}")
            }
            SolveError::SatModelRejected => f.write_str(
                "the SAT fallback returned an assignment the constraint evaluator rejects \
                 (an encoding bug)",
            ),
            SolveError::SatSolverFailure { what } => write!(f, "the SAT fallback failed: {what}"),
        }
    }
}

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

impl From<BddBuildError> for SolveError {
    fn from(e: BddBuildError) -> Self {
        SolveError::Build(e)
    }
}

/// Compiled problems by [`problem_digest`]. Compiling is a pure function of
/// its arguments, and a `randomize()` in a loop (or a soft-constraint ladder,
/// LRM 18.5.13.1, which retries the same hard problem) re-asks the same one
/// over and over; the build was most of a call. Bounded in entries and in the
/// size of what is kept; cleared when full.
static SOLVER_CACHE: Mutex<BTreeMap<(u64, u64), Solver>> = Mutex::new(BTreeMap::new());
const SOLVER_CACHE_CAP: usize = 512;
const SOLVER_CACHE_MAX_NODES: usize = 4096;

/// A compiled, counted constraint problem, ready to sample.
#[derive(Clone, Debug)]
pub struct Solver {
    compiled: Compiled,
    count: SatCount,
    /// Per-bit draw steps already computed, shared by every clone (see [`StepMemo`]).
    steps: Arc<Mutex<StepMemo>>,
}

/// One bit decision of the uniform tuple draw: the models with the bit set
/// and their count, and (once a draw needed it) the models with it clear.
#[derive(Debug)]
struct Step {
    ones: Bdd,
    ones_count: BigUint,
    zeros: Option<Bdd>,
}

/// Memo of [`Solver::sample_plan`]'s per-bit BDD conjunctions. The BDD a draw
/// holds before a bit is a pure function of the plan and of every decision
/// taken earlier in the draw (its `trace`), so a repeated prefix reuses the
/// conjunction and the cardinality that an apply and a model count would
/// recompute. Only the random numbers differ between draws, so the output is
/// unchanged. Bounded: cleared when it holds [`STEP_MEMO_CAP`] steps.
#[derive(Debug, Default)]
struct StepMemo {
    plans: Vec<(Plan, HashMap<Vec<u32>, Step>)>,
}
const STEP_MEMO_CAP: usize = 100_000;
const STEP_MEMO_PLANS: usize = 4;

impl Solver {
    /// Compiles the conjunction of `conjuncts` with the default node budget.
    pub fn new(
        var_tys: &[Ty],
        states: &[Option<BigUint>],
        conjuncts: &[Expr],
    ) -> Result<Solver, SolveError> {
        Solver::with_budget(var_tys, states, conjuncts, NODE_BUDGET)
    }

    /// As [`Solver::new`] with an explicit per-BDD node budget.
    pub fn with_budget(
        var_tys: &[Ty],
        states: &[Option<BigUint>],
        conjuncts: &[Expr],
        budget: usize,
    ) -> Result<Solver, SolveError> {
        let digest = problem_digest(var_tys, states, conjuncts, budget);
        Solver::with_digest(digest, var_tys, states, conjuncts, budget)
    }

    /// [`Solver::with_budget`] for a caller that already hashed the problem
    /// (hashing formats the whole expression tree, so it is done once per call).
    fn with_digest(
        digest: (u64, u64),
        var_tys: &[Ty],
        states: &[Option<BigUint>],
        conjuncts: &[Expr],
        budget: usize,
    ) -> Result<Solver, SolveError> {
        if let Some(s) = SOLVER_CACHE
            .lock()
            .unwrap_or_else(|p| p.into_inner())
            .get(&digest)
        {
            return Ok(s.clone());
        }
        let compiled = compile(var_tys, states, conjuncts, budget)?;
        let count = compiled.bdd.exact_cardinality();
        let solver = Solver {
            compiled,
            count,
            steps: Arc::default(),
        };
        if solver.compiled.bdd.size() <= SOLVER_CACHE_MAX_NODES {
            let mut cache = SOLVER_CACHE.lock().unwrap_or_else(|p| p.into_inner());
            if cache.len() >= SOLVER_CACHE_CAP {
                cache.clear();
            }
            cache.insert(digest, solver.clone());
        }
        Ok(solver)
    }

    /// The exact number of legal assignments.
    pub fn count(&self) -> &SatCount {
        &self.count
    }

    /// The underlying BDD and variable layout.
    pub fn compiled(&self) -> &Compiled {
        &self.compiled
    }

    /// Draws one assignment uniformly from all models (module docs), fixing
    /// variables in `order`, then any variable `order` omits in id order.
    /// Returns one value per variable, indexed by [`VarId`].
    pub fn sample(&self, order: &[VarId], seed: &mut u64) -> Result<Vec<BigUint>, SolveError> {
        // One group: no later group to project away, so this is exactly the
        // per-bit walk over all models.
        self.sample_groups(&[order.to_vec()], seed)
    }

    /// Draws one assignment solving the variables group by group, in the
    /// order of `groups` (`solve ... before`, LRM 18.5.9; [`super::order`]
    /// builds them).
    ///
    /// Each group is sampled *uniformly over its own distinct feasible value
    /// tuples*: the variables of every later group are existentially
    /// projected away, so a group value is feasible when some assignment of
    /// the later groups completes it, and its probability does not depend on
    /// how many completions it has. That is the LRM's `solve s before d`
    /// (Table 18-2: `s` is 50/50 although `s == 1` has a single completion
    /// and `s == 0` has 2^32). The chosen values are then fixed in the full
    /// constraint set and the next group is drawn from what remains. Inside a
    /// group the per-bit walk of [`Solver::sample`] applies, so a group is
    /// uniform over its tuples whatever its internal order.
    ///
    /// Any variable `groups` omits is solved with the *last* group (LRM
    /// 18.5.9: unordered variables are solved with the last set of ordered
    /// variables), in id order. A variable listed twice is
    /// [`SolveError::BadOrder`].
    pub fn sample_groups(
        &self,
        groups: &[Vec<VarId>],
        seed: &mut u64,
    ) -> Result<Vec<BigUint>, SolveError> {
        let plan = Plan {
            groups: groups
                .iter()
                .map(|g| Group {
                    vars: g.clone(),
                    weighted: Vec::new(),
                })
                .collect(),
            randc: Vec::new(),
        };
        self.sample_plan(&plan, &mut [], seed)
    }

    /// Draws one assignment following `plan` (see [`Plan`]): first the `randc`
    /// variables, each from its own cycle state (`randc_states[i]` belongs to
    /// `plan.randc[i]`; [`super::randc`]), then the ordered groups, each with
    /// its weighted variables before its uniform ones. A state is only
    /// advanced for a variable that was drawn, so a caller that must not
    /// commit a failed call clones the states first.
    pub fn sample_plan(
        &self,
        plan: &Plan,
        randc_states: &mut [RandcState],
        seed: &mut u64,
    ) -> Result<Vec<BigUint>, SolveError> {
        let layout = &self.compiled.layout;
        let n = layout.tys.len();
        let widths: Vec<usize> = layout.bits.iter().map(Vec::len).collect();
        let full = resolve_plan(&widths, plan, randc_states.len())?;
        let last = full.len() - 1;

        if self.count.is_zero() {
            return Err(SolveError::Unsatisfiable);
        }
        let mut cur = self.compiled.bdd.clone();
        let mut out = vec![BigUint::zero(); n];
        let mut memo_guard = self.steps.lock().unwrap_or_else(|p| p.into_inner());
        let memo_all = &mut *memo_guard;
        if memo_all.plans.iter().map(|(_, m)| m.len()).sum::<usize>() >= STEP_MEMO_CAP {
            *memo_all = StepMemo::default();
        }
        let slot = match memo_all.plans.iter().position(|(p, _)| p == plan) {
            Some(i) => i,
            None => {
                if memo_all.plans.len() >= STEP_MEMO_PLANS {
                    *memo_all = StepMemo::default();
                }
                memo_all.plans.push((plan.clone(), HashMap::new()));
                memo_all.plans.len() - 1
            }
        };
        let memo = &mut memo_all.plans[slot].1;
        // Every decision of this draw so far: the memo key. A tagged word is a
        // randc or weighted code (tag 1) or a bit (tag 0, `var << 1 | bit`).
        let mut trace: Vec<u32> = Vec::new();

        // The `randc` variables come first, before every `rand` one (LRM
        // 18.4.2), in id order.
        for (k, &v) in plan.randc.iter().enumerate() {
            let v = v as usize;
            let others: Vec<_> = (0..n)
                .filter(|&o| o != v)
                .flat_map(|o| layout.bits[o].iter().copied())
                .collect();
            let base = self.compiled.bdd.exists(&others);
            let now = cur.exists(&others);
            let w = layout.bits[v].len();
            let holds = |f: &Bdd, x: u32| -> bool {
                let mut val = BddValuation::all_false(layout.num_bits() as u16);
                for (i, bv) in layout.bits[v].iter().enumerate() {
                    val.set_value(*bv, (x >> i) & 1 == 1);
                }
                f.eval_in(&val)
            };
            let fp = fingerprint_of(&base.to_bytes());
            let picked = randc::pick(
                &mut randc_states[k],
                fp,
                &mut || (0..1u32 << w).filter(|&x| holds(&base, x)).collect(),
                &mut |x| holds(&now, x),
                seed,
            )
            .ok_or(SolveError::Unsatisfiable)?;
            trace.push((1 << 31) | picked);
            for (i, bv) in layout.bits[v].iter().enumerate() {
                let bit = (picked >> i) & 1 == 1;
                out[v].set_bit(i as u64, bit);
                cur = cur.var_select(*bv, bit);
            }
        }

        for (gi, (group, weighted)) in full.iter().enumerate() {
            trace.push(u32::MAX - gi as u32);
            // Existentially project the later groups away. Variables that
            // become free multiply every count by the same power of two, which
            // cancels in the ratios below.
            let later: Vec<_> = full[gi + 1..]
                .iter()
                .flat_map(|(vars, wtd)| vars.iter().chain(wtd.iter().map(|w| &w.var)))
                .flat_map(|&v| layout.bits[v].iter().copied())
                .collect();
            let mut proj = if later.is_empty() {
                cur.clone()
            } else {
                cur.exists(&later)
            };
            // Weighted variables: a code is chosen with probability
            // proportional to its weight, scaled by its measure, among the
            // codes still feasible.
            let mut drawn: Vec<usize> = Vec::new();
            for wv in weighted {
                let v = wv.var;
                let bits = &layout.bits[v];
                let fixing = |code: usize| -> Vec<_> {
                    bits.iter()
                        .enumerate()
                        .map(|(i, bv)| (*bv, (code >> i) & 1 == 1))
                        .collect()
                };
                // Every bit but the subject's, to project the rest away.
                let all_but = |keep: usize| -> Vec<_> {
                    (0..n)
                        .filter(|&o| o != keep)
                        .flat_map(|o| layout.bits[o].iter().copied())
                        .collect()
                };
                let subject_rest = match wv.measure {
                    Measure::Distinct(s) => Some(all_but(s as usize)),
                    Measure::Flat | Measure::Models => None,
                };
                // (code, numerator * measure, denominator)
                let mut feasible: Vec<(usize, BigUint, &BigUint)> = Vec::new();
                for (c, wt) in wv.weights.iter().enumerate() {
                    if wt.num.is_zero() {
                        continue;
                    }
                    let pc = proj.restrict(&fixing(c));
                    if pc.is_false() {
                        continue;
                    }
                    let m = match (&wv.measure, &subject_rest) {
                        (Measure::Models, _) => pc.exact_cardinality(),
                        (_, Some(rest)) => pc.exists(rest).exact_cardinality(),
                        _ => BigUint::one(),
                    };
                    feasible.push((c, &wt.num * m, &wt.den));
                }
                let code = draw_weighted(&feasible, seed)?;
                trace.push((1 << 31) | code as u32);
                let fix = fixing(code);
                proj = proj.restrict(&fix);
                if gi < last {
                    cur = cur.restrict(&fix);
                }
                for (i, _) in bits.iter().enumerate() {
                    out[v].set_bit(i as u64, (code >> i) & 1 == 1);
                }
                // The subject is drawn now, uniformly over its distinct
                // feasible values (they are all in the chosen item).
                if let (Measure::Distinct(s), Some(rest)) = (wv.measure, &subject_rest) {
                    let s = s as usize;
                    let mut pj = proj.exists(rest);
                    let mut remaining = pj.exact_cardinality();
                    for i in (0..layout.bits[s].len()).rev() {
                        let bv = layout.bits[s][i];
                        let ones = pj.var_select(bv, true);
                        let c1 = ones.exact_cardinality();
                        let bit = rng::below(seed, &remaining) < c1;
                        if bit {
                            pj = ones;
                            remaining = c1;
                        } else {
                            remaining -= &c1;
                            pj = pj.var_select(bv, false);
                        }
                        proj = proj.var_select(bv, bit);
                        trace.push((bv.to_index() as u32) << 1 | u32::from(bit));
                        if gi < last {
                            cur = cur.var_select(bv, bit);
                        }
                        out[s].set_bit(i as u64, bit);
                    }
                    drawn.push(s);
                }
            }
            let mut remaining = proj.exact_cardinality();
            // A bit outside `proj`'s support is unconstrained: exactly half of
            // the `remaining` models have it set, so it is drawn from that
            // count without conjoining a literal to `proj` (one full BDD
            // apply per bit is what dominated a randomize call). `proj` then
            // counts every skipped bit as free, so the models with a later
            // bit fixed are `skipped` doublings too many; both are undone
            // exactly below. The last group only: an earlier one must also
            // keep `cur` fixed for the groups after it.
            let support = (gi == last).then(|| proj.support_set());
            let mut skipped = 0u64;
            for &v in group {
                if drawn.contains(&v) {
                    continue;
                }
                for i in (0..layout.bits[v].len()).rev() {
                    let bv = layout.bits[v][i];
                    if support.as_ref().is_some_and(|s| !s.contains(&bv)) {
                        let half = &remaining >> 1u32;
                        let bit = rng::below(seed, &remaining) < half;
                        remaining = half;
                        skipped += 1;
                        if bit {
                            out[v].set_bit(i as u64, true);
                        }
                        continue;
                    }
                    if !memo.contains_key(trace.as_slice()) {
                        let ones = proj.var_select(bv, true);
                        let ones_count = ones.exact_cardinality();
                        memo.insert(
                            trace.clone(),
                            Step {
                                ones,
                                ones_count,
                                zeros: None,
                            },
                        );
                    }
                    let step = memo.get_mut(trace.as_slice()).expect("just inserted");
                    let c1 = &step.ones_count >> skipped;
                    let zeros_count = &remaining - &c1;
                    let pick = rng::below(seed, &remaining);
                    let bit = pick < c1;
                    if bit {
                        out[v].set_bit(i as u64, true);
                        proj = step.ones.clone();
                        remaining = c1;
                    } else {
                        proj = step
                            .zeros
                            .get_or_insert_with(|| proj.var_select(bv, false))
                            .clone();
                        remaining = zeros_count;
                    }
                    trace.push((bv.to_index() as u32) << 1 | u32::from(bit));
                    // The last group's `proj` is `cur`; nothing follows it.
                    if gi < last {
                        cur = cur.var_select(bv, bit);
                    }
                }
            }
        }
        debug_assert!(self.compiled.accepts(&out), "sampled a non-model");
        Ok(out)
    }
}

/// Validates `plan` against the variables' bit `widths` and resolves it to
/// bit-vector indices: one entry per group, the variables the plan omits
/// appended to the last group (LRM 18.5.9), and a single empty group when
/// the plan has none. `randc_states` is how many cycle states the caller
/// holds. Shared by the exact and the SAT sampler so both reject exactly the
/// same malformed plans.
pub(crate) fn resolve_plan<'a>(
    widths: &[usize],
    plan: &'a Plan,
    randc_states: usize,
) -> Result<Vec<Resolved<'a>>, SolveError> {
    let n = widths.len();
    if randc_states != plan.randc.len() {
        return Err(SolveError::BadPlan {
            what: "one randc state per randc variable",
        });
    }
    for &v in &plan.randc {
        let w = widths.get(v as usize).copied().unwrap_or(0);
        if w > randc::MAX_RANDC_WIDTH as usize {
            return Err(SolveError::BadPlan {
                what: "a randc variable is wider than the permutation engine supports",
            });
        }
    }
    let mut seen = vec![false; n];
    let mut mark = |v: VarId| -> Result<usize, SolveError> {
        let i = v as usize;
        if i >= n || seen[i] {
            return Err(SolveError::BadOrder { var: v });
        }
        seen[i] = true;
        Ok(i)
    };
    for &v in &plan.randc {
        mark(v)?;
    }
    let mut full: Vec<Resolved<'a>> = Vec::with_capacity(plan.groups.len().max(1));
    for g in &plan.groups {
        let mut weighted = Vec::with_capacity(g.weighted.len());
        for w in &g.weighted {
            let i = mark(w.var)?;
            let width = widths[i];
            if width >= 32 || w.weights.len() > 1usize << width {
                return Err(SolveError::BadPlan {
                    what: "a weighted variable has more weights than values",
                });
            }
            if w.weights.iter().any(|x| x.den.is_zero()) {
                return Err(SolveError::BadPlan {
                    what: "a weight has a zero denominator",
                });
            }
            weighted.push(ResolvedWeighted {
                var: i,
                weights: w.weights.as_slice(),
                measure: w.measure,
            });
        }
        let mut members = Vec::with_capacity(g.vars.len());
        for &v in &g.vars {
            members.push(mark(v)?);
        }
        for w in &weighted {
            if let Measure::Distinct(s) = w.measure
                && !members.contains(&(s as usize))
            {
                return Err(SolveError::BadPlan {
                    what: "a weighted variable's subject is not a uniform member of its group",
                });
            }
        }
        full.push((members, weighted));
    }
    if full.is_empty() {
        full.push((Vec::new(), Vec::new()));
    }
    let last = full.len() - 1;
    full[last].0.extend((0..n).filter(|i| !seen[*i]));
    Ok(full)
}

/// Draws one code from `feasible` (`(code, numerator * measure,
/// denominator)`) with probability proportional to its weight, in exact
/// rational arithmetic over the common denominator; no feasible weight is
/// [`SolveError::Unsatisfiable`] and consumes no draw.
pub(crate) fn draw_weighted(
    feasible: &[(usize, BigUint, &BigUint)],
    seed: &mut u64,
) -> Result<usize, SolveError> {
    let common = feasible
        .iter()
        .fold(BigUint::one(), |acc, (_, _, den)| acc * *den);
    let scaled: Vec<(usize, BigUint)> = feasible
        .iter()
        .map(|(c, num, den)| (*c, num * (&common / *den)))
        .collect();
    let total: BigUint = scaled.iter().map(|(_, w)| w.clone()).sum();
    if total.is_zero() {
        return Err(SolveError::Unsatisfiable);
    }
    let mut pick = rng::below(seed, &total);
    let mut code = scaled[scaled.len() - 1].0;
    for (c, w) in &scaled {
        if pick < *w {
            code = *c;
            break;
        }
        pick -= w;
    }
    Ok(code)
}

/// The solver `randomize()` uses: the exact BDD [`Solver`] when the problem
/// fits its node budget, the approximate [`SatSolver`] when it does not (plan
/// section 3.4, ledger XC-063). Both take the same [`Plan`] and obey the same
/// contracts; what differs is exactness of the distribution and what the SAT
/// side refuses ([`SolveError::SatUnsupported`]; [`super::sat`] module docs).
#[derive(Debug)]
pub enum Engine {
    /// The problem fit the BDD node budget: exactly uniform.
    Exact(Solver),
    /// It did not: approximately uniform.
    Sat(SatSolver),
}

/// Problems (as a 128-bit digest) whose BDD build already exceeded its budget
/// in this process, so the next call skips straight to SAT instead of paying
/// for another failing build. Bounded; cleared when full.
static BDD_TOO_LARGE: Mutex<BTreeSet<(u64, u64)>> = Mutex::new(BTreeSet::new());
const BDD_TOO_LARGE_CAP: usize = 4096;

/// The digest identifying one build problem: two independent fixed-key SipHash
/// streams over the structural `Hash` of the inputs (no text is formatted; the
/// digest is only a cache key, never persisted).
fn problem_digest(
    var_tys: &[Ty],
    states: &[Option<BigUint>],
    conjuncts: &[Expr],
    budget: usize,
) -> (u64, u64) {
    use std::hash::{Hash as _, Hasher as _};
    let mut lo = std::collections::hash_map::DefaultHasher::new();
    let mut hi = std::collections::hash_map::DefaultHasher::new();
    0x1234_5678_9abc_def1_u64.hash(&mut hi);
    for h in [&mut lo, &mut hi] {
        budget.hash(h);
        var_tys.hash(h);
        states.hash(h);
        conjuncts.hash(h);
    }
    (lo.finish(), hi.finish())
}

impl Engine {
    /// Compiles the conjunction of `conjuncts` with the default BDD budget and
    /// SAT limits.
    pub fn new(
        var_tys: &[Ty],
        states: &[Option<BigUint>],
        conjuncts: &[Expr],
    ) -> Result<Engine, SolveError> {
        Engine::with_limits(
            var_tys,
            states,
            conjuncts,
            NODE_BUDGET,
            SatLimits::default(),
        )
    }

    /// As [`Engine::new`] with an explicit BDD node budget and SAT limits.
    pub fn with_limits(
        var_tys: &[Ty],
        states: &[Option<BigUint>],
        conjuncts: &[Expr],
        budget: usize,
        limits: SatLimits,
    ) -> Result<Engine, SolveError> {
        let digest = problem_digest(var_tys, states, conjuncts, budget);
        let known_large = BDD_TOO_LARGE
            .lock()
            .unwrap_or_else(|p| p.into_inner())
            .contains(&digest);
        if !known_large {
            match Solver::with_digest(digest, var_tys, states, conjuncts, budget) {
                Ok(s) => return Ok(Engine::Exact(s)),
                Err(SolveError::Build(BddBuildError::TooLarge { .. })) => {
                    let mut seen = BDD_TOO_LARGE.lock().unwrap_or_else(|p| p.into_inner());
                    if seen.len() >= BDD_TOO_LARGE_CAP {
                        seen.clear();
                    }
                    seen.insert(digest);
                }
                Err(e) => return Err(e),
            }
        }
        SatSolver::with_limits(var_tys, states, conjuncts, limits).map(Engine::Sat)
    }

    /// Whether this problem fell back to SAT.
    pub fn is_sat(&self) -> bool {
        matches!(self, Engine::Sat(_))
    }

    /// As [`Solver::sample_plan`] (exact) or [`SatSolver::sample_plan`].
    pub fn sample_plan(
        &self,
        plan: &Plan,
        randc_states: &mut [RandcState],
        seed: &mut u64,
    ) -> Result<Vec<BigUint>, SolveError> {
        match self {
            Engine::Exact(s) => s.sample_plan(plan, randc_states, seed),
            Engine::Sat(s) => s.sample_plan(plan, randc_states, seed),
        }
    }
}

/// One-shot: compile, count and draw a single sample.
pub fn solve(
    var_tys: &[Ty],
    states: &[Option<BigUint>],
    conjuncts: &[Expr],
    order: &[VarId],
    seed: &mut u64,
) -> Result<Vec<BigUint>, SolveError> {
    Solver::new(var_tys, states, conjuncts)?.sample(order, seed)
}

/// One-shot: compile, count and draw a single sample solving `groups` in
/// order ([`Solver::sample_groups`]).
pub fn solve_grouped(
    var_tys: &[Ty],
    states: &[Option<BigUint>],
    conjuncts: &[Expr],
    groups: &[Vec<VarId>],
    seed: &mut u64,
) -> Result<Vec<BigUint>, SolveError> {
    Solver::new(var_tys, states, conjuncts)?.sample_groups(groups, seed)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::constraints::ir::BinOp;

    fn bin(op: BinOp, a: Expr, b: Expr, t: Ty) -> Expr {
        Expr::Binary(op, Box::new(a), Box::new(b), t)
    }

    /// LRM 18.5.9 Table 18-1/18-2 with an 8-bit `d`: `s -> d == 0`.
    fn table_18_2() -> Solver {
        let s = Expr::Var(0, Ty::BIT);
        let d = Expr::Var(1, Ty::unsigned(8));
        let c = bin(
            BinOp::LogImpl,
            s,
            bin(BinOp::Eq, d, Expr::constant(Ty::unsigned(8), 0u8), Ty::BIT),
            Ty::BIT,
        );
        Solver::new(&[Ty::BIT, Ty::unsigned(8)], &[], &[c]).expect("compiles")
    }

    /// `(count of s == 1, count of d == 0)` over `n` draws.
    fn tally(sv: &Solver, groups: &[Vec<VarId>], n: u32, seed0: u64) -> (u32, u32) {
        let mut seed = seed0;
        let (mut s1, mut d0) = (0, 0);
        for _ in 0..n {
            let v = sv.sample_groups(groups, &mut seed).expect("sat");
            s1 += u32::from(!v[0].is_zero());
            d0 += u32::from(v[1].is_zero());
        }
        (s1, d0)
    }

    #[test]
    fn solve_s_before_d_is_table_18_2() {
        let sv = table_18_2();
        assert_eq!(sv.count(), &BigUint::from(257u32));
        let n = 40_000;
        let (s1, d0) = tally(&sv, &[vec![0], vec![1]], n, 11);
        // s is 50/50 (sigma 100): the single completion of s == 1 does not
        // make it unlikely.
        assert!((19_400..20_600).contains(&s1), "s==1: {s1}");
        // P(d == 0) = 1/2 + 1/2 * 2^-8 = 0.50195 (sigma 100).
        assert!((19_500..20_700).contains(&d0), "d==0: {d0}");
    }

    #[test]
    fn unordered_is_table_18_1() {
        let sv = table_18_2();
        let n = 40_000;
        // One group, and no group at all, both mean "unordered".
        for groups in [vec![vec![0, 1]], vec![]] {
            let (s1, d0) = tally(&sv, &groups, n, 12);
            // P(s == 1) = 1/257 = 155.6 (sigma 12.5); P(d == 0) = 2/257.
            assert!((80..235).contains(&s1), "s==1: {s1}");
            assert!((235..390).contains(&d0), "d==0: {d0}");
        }
    }

    #[test]
    fn solving_d_first_makes_d_uniform_over_its_feasible_values() {
        let sv = table_18_2();
        let n = 51_200;
        let (s1, d0) = tally(&sv, &[vec![1], vec![0]], n, 13);
        // Every d is feasible (s == 0 always completes it): P(d == 0) = 1/256
        // = 200 (sigma 14). Given d == 0, s is 50/50, so P(s == 1) = 1/512.
        assert!((130..270).contains(&d0), "d==0: {d0}");
        assert!((60..140).contains(&s1), "s==1: {s1}");
    }

    #[test]
    fn a_variable_the_groups_omit_is_solved_with_the_last_group() {
        let sv = table_18_2();
        // Only `s` is listed: `d` joins the last (only) group, so this is the
        // unordered distribution.
        let (s1, _) = tally(&sv, &[vec![0]], 40_000, 14);
        assert!((80..235).contains(&s1), "s==1: {s1}");
    }

    #[test]
    fn one_group_is_exactly_the_per_bit_sample() {
        let sv = table_18_2();
        let (mut a, mut b) = (5u64, 5u64);
        for _ in 0..50 {
            assert_eq!(
                sv.sample(&[1, 0], &mut a).unwrap(),
                sv.sample_groups(&[vec![1, 0]], &mut b).unwrap()
            );
        }
        assert_eq!(a, b, "the same draws were consumed");
    }

    #[test]
    fn a_repeated_or_unknown_variable_is_a_bad_order() {
        let sv = table_18_2();
        let mut seed = 1u64;
        assert_eq!(
            sv.sample_groups(&[vec![0], vec![0]], &mut seed),
            Err(SolveError::BadOrder { var: 0 })
        );
        assert_eq!(
            sv.sample_groups(&[vec![7]], &mut seed),
            Err(SolveError::BadOrder { var: 7 })
        );
    }

    #[test]
    fn an_unsatisfiable_problem_consumes_no_draws_when_grouped() {
        let x = Expr::Var(0, Ty::unsigned(4));
        let c = bin(BinOp::Lt, x, Expr::constant(Ty::unsigned(4), 0u8), Ty::BIT);
        let sv = Solver::new(&[Ty::unsigned(4)], &[], &[c]).unwrap();
        let mut seed = 9u64;
        assert_eq!(
            sv.sample_groups(&[vec![0]], &mut seed),
            Err(SolveError::Unsatisfiable)
        );
        assert_eq!(seed, 9);
    }

    #[test]
    fn three_levels_are_each_uniform_over_their_feasible_values() {
        // a < b < c over 3-bit values: a in 0..=5, c in 2..=7, b in 1..=6.
        // Solving a, then b, then c: a is uniform over its 6 feasible values
        // (unordered it would favour small a), and every draw is a model.
        let t = Ty::unsigned(3);
        let (a, b, c) = (Expr::Var(0, t), Expr::Var(1, t), Expr::Var(2, t));
        let conj = vec![
            bin(BinOp::Lt, a, b.clone(), Ty::BIT),
            bin(BinOp::Lt, b, c, Ty::BIT),
        ];
        let sv = Solver::new(&[t, t, t], &[], &conj).unwrap();
        let n = 36_000u32;
        let mut seed = 21u64;
        let mut hist = [0u32; 8];
        for _ in 0..n {
            let v = sv
                .sample_groups(&[vec![0], vec![1], vec![2]], &mut seed)
                .unwrap();
            assert!(v[0] < v[1] && v[1] < v[2], "{v:?}");
            hist[usize::try_from(v[0].clone()).unwrap()] += 1;
        }
        // 6 values, 6000 each (sigma ~71).
        for (i, h) in hist.iter().enumerate() {
            if i < 6 {
                assert!((5_500..6_500).contains(h), "a=={i}: {hist:?}");
            } else {
                assert_eq!(*h, 0, "a=={i} is infeasible: {hist:?}");
            }
        }
    }

    fn plan(randc: Vec<VarId>, groups: Vec<Group>) -> Plan {
        Plan { groups, randc }
    }

    fn uniform(vars: Vec<VarId>) -> Group {
        Group {
            vars,
            weighted: Vec::new(),
        }
    }

    #[test]
    fn a_randc_variable_cycles_through_every_value_across_calls() {
        // LRM 18.4.2 / Figure 18-1: v0 is a 2-bit randc, v1 a free rand bit.
        let t = Ty::unsigned(2);
        let sv = Solver::new(&[t, Ty::BIT], &[], &[]).unwrap();
        let p = plan(vec![0], vec![uniform(vec![1])]);
        let mut states = vec![RandcState::default()];
        let mut seed = 31u64;
        for cycle in 0..40 {
            let mut got: Vec<u32> = (0..4)
                .map(|_| {
                    let v = sv.sample_plan(&p, &mut states, &mut seed).unwrap();
                    u32::try_from(v[0].clone()).unwrap()
                })
                .collect();
            got.sort_unstable();
            assert_eq!(got, vec![0, 1, 2, 3], "cycle {cycle}");
        }
    }

    #[test]
    fn randc_variables_are_solved_before_every_rand_group() {
        // v0 randc (2 bits), v1 rand: `v1 -> v0 == 0`. v0 is drawn first
        // (uniformly, every value has a completion), so v1 can only be 1 on
        // the quarter of calls where v0 == 0, and there it is 50/50. Had v1
        // been drawn first, P(v1 == 1) would be 1/2.
        let t = Ty::unsigned(2);
        let c = bin(
            BinOp::LogImpl,
            Expr::Var(1, Ty::BIT),
            bin(BinOp::Eq, Expr::Var(0, t), Expr::constant(t, 0u8), Ty::BIT),
            Ty::BIT,
        );
        let sv = Solver::new(&[t, Ty::BIT], &[], &[c]).unwrap();
        // Even when `solve v1 before ...` style grouping lists v1 first.
        let p = plan(vec![0], vec![uniform(vec![1])]);
        let mut states = vec![RandcState::default()];
        let mut seed = 32u64;
        let (mut v1_ones, mut v0_zero) = (0u32, 0u32);
        for _ in 0..4000 {
            let v = sv.sample_plan(&p, &mut states, &mut seed).unwrap();
            v0_zero += u32::from(v[0].is_zero());
            if !v[1].is_zero() {
                v1_ones += 1;
                assert!(v[0].is_zero(), "{v:?}");
            }
        }
        assert_eq!(v0_zero, 1000, "a randc value occurs once per cycle");
        assert!((400..600).contains(&v1_ones), "v1==1: {v1_ones}");
    }

    #[test]
    fn a_randc_variable_whose_legal_set_changes_restarts_its_cycle() {
        // v0 randc (3 bits), `v0 < s` with s a state: the legal set is the
        // first s values; changing s starts a new permutation (18.4.2).
        let t = Ty::unsigned(3);
        let c = bin(BinOp::Lt, Expr::Var(0, t), Expr::State(0, t), Ty::BIT);
        let p = plan(vec![0], vec![]);
        let mut states = vec![RandcState::default()];
        let mut seed = 33u64;
        let wide =
            Solver::new(&[t], &[Some(BigUint::from(4u8))], std::slice::from_ref(&c)).unwrap();
        let first = wide.sample_plan(&p, &mut states, &mut seed).unwrap();
        assert_eq!(states[0].remaining.len(), 3);
        let narrow = Solver::new(&[t], &[Some(BigUint::from(2u8))], &[c]).unwrap();
        let mut seen = std::collections::BTreeSet::new();
        for _ in 0..2 {
            let v = narrow.sample_plan(&p, &mut states, &mut seed).unwrap();
            seen.insert(u32::try_from(v[0].clone()).unwrap());
        }
        assert_eq!(seen, [0, 1].into_iter().collect());
        assert!(first[0] < BigUint::from(4u8));
    }

    #[test]
    fn equal_functions_serialise_to_equal_bytes() {
        // The randc fingerprint hashes the BDD: two spellings of one
        // constraint must agree, a different one must not.
        let t = Ty::unsigned(3);
        let x = || Expr::Var(0, t);
        let k = |n: u8| Expr::constant(t, n);
        let a = bin(BinOp::Lt, x(), k(2), Ty::BIT);
        let b = bin(
            BinOp::LogOr,
            bin(BinOp::Eq, x(), k(0), Ty::BIT),
            bin(BinOp::Eq, x(), k(1), Ty::BIT),
            Ty::BIT,
        );
        let c = bin(BinOp::Lt, x(), k(3), Ty::BIT);
        let bytes = |e: Expr| {
            Solver::new(&[t], &[], &[e])
                .unwrap()
                .compiled()
                .bdd
                .to_bytes()
        };
        assert_eq!(bytes(a.clone()), bytes(b));
        assert_ne!(bytes(a), bytes(c));
    }

    #[test]
    fn a_weighted_variable_is_drawn_in_proportion_to_its_weights() {
        let t = Ty::unsigned(2);
        let sv = Solver::new(&[t], &[], &[]).unwrap();
        // Weights 1 : 3 : 0, and code 3 has none: only 0 and 1 appear.
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![],
                weighted: vec![Weighted {
                    var: 0,
                    weights: ints(&[1, 3, 0]),
                    measure: Measure::Flat,
                }],
            }],
        );
        let mut seed = 34u64;
        let mut hist = [0u32; 4];
        for _ in 0..40_000 {
            let v = sv.sample_plan(&p, &mut [], &mut seed).unwrap();
            hist[usize::try_from(v[0].clone()).unwrap()] += 1;
        }
        assert_eq!(hist[2] + hist[3], 0, "{hist:?}");
        // 10000 / 30000 (sigma 87).
        assert!((9_600..10_400).contains(&hist[0]), "{hist:?}");
        assert!((29_600..30_400).contains(&hist[1]), "{hist:?}");
    }

    #[test]
    fn a_weighted_variable_only_takes_feasible_codes() {
        // Code 1 carries all the weight but is excluded: weights are
        // renormalised over what is left (18.5.3: unlisted values excluded).
        let t = Ty::unsigned(2);
        let c = bin(BinOp::Ne, Expr::Var(0, t), Expr::constant(t, 1u8), Ty::BIT);
        let sv = Solver::new(&[t], &[], &[c]).unwrap();
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![],
                weighted: vec![Weighted {
                    var: 0,
                    weights: ints(&[1, 100, 1, 0]),
                    measure: Measure::Flat,
                }],
            }],
        );
        let mut seed = 35u64;
        let mut hist = [0u32; 4];
        for _ in 0..4_000 {
            let v = sv.sample_plan(&p, &mut [], &mut seed).unwrap();
            hist[usize::try_from(v[0].clone()).unwrap()] += 1;
        }
        assert_eq!(hist[1] + hist[3], 0, "{hist:?}");
        assert!((1_800..2_200).contains(&hist[0]), "{hist:?}");
    }

    #[test]
    fn weights_with_no_feasible_code_are_unsatisfiable() {
        let t = Ty::unsigned(2);
        let sv = Solver::new(&[t], &[], &[]).unwrap();
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![],
                weighted: vec![Weighted {
                    var: 0,
                    weights: ints(&[0; 4]),
                    measure: Measure::Flat,
                }],
            }],
        );
        let mut seed = 36u64;
        assert_eq!(
            sv.sample_plan(&p, &mut [], &mut seed),
            Err(SolveError::Unsatisfiable)
        );
    }

    #[test]
    fn malformed_plans_are_rejected_loudly() {
        let t = Ty::unsigned(2);
        let sv = Solver::new(&[t, Ty::unsigned(17)], &[], &[]).unwrap();
        let mut seed = 37u64;
        // One randc state is missing.
        assert!(matches!(
            sv.sample_plan(&plan(vec![0], vec![]), &mut [], &mut seed),
            Err(SolveError::BadPlan { .. })
        ));
        // A 17-bit randc exceeds the engine's width.
        assert!(matches!(
            sv.sample_plan(
                &plan(vec![1], vec![]),
                &mut [RandcState::default()],
                &mut seed
            ),
            Err(SolveError::BadPlan { .. })
        ));
        // More weights than the 2-bit variable has values.
        let too_many = Group {
            vars: vec![],
            weighted: vec![Weighted {
                var: 0,
                weights: ints(&[1; 5]),
                measure: Measure::Flat,
            }],
        };
        assert!(matches!(
            sv.sample_plan(&plan(vec![], vec![too_many]), &mut [], &mut seed),
            Err(SolveError::BadPlan { .. })
        ));
        // A variable both randc and in a group is a repeat.
        assert_eq!(
            sv.sample_plan(
                &plan(vec![0], vec![uniform(vec![0])]),
                &mut [RandcState::default()],
                &mut seed
            ),
            Err(SolveError::BadOrder { var: 0 })
        );
    }

    fn ints(ws: &[u32]) -> Vec<Weight> {
        ws.iter().map(|&w| Weight::int(w)).collect()
    }

    /// `K == i` over a selector variable `k` of type `kt`.
    fn is_code(k: u32, kt: Ty, i: u64) -> Expr {
        bin(BinOp::Eq, Expr::Var(k, kt), Expr::constant(kt, i), Ty::BIT)
    }

    /// `lo <= x && x <= hi` for the variable `x`.
    fn in_range(x: &Expr, lo: u64, hi: u64) -> Expr {
        let t = x.ty();
        bin(
            BinOp::LogAnd,
            bin(BinOp::Ge, x.clone(), Expr::constant(t, lo), Ty::BIT),
            bin(BinOp::Le, x.clone(), Expr::constant(t, hi), Ty::BIT),
            Ty::BIT,
        )
    }

    /// The conjunct of a `dist` whose item `i` is `items[i]` (inclusive), over
    /// the subject `x` and the selector `k` (id 1, 3 bits).
    fn dist_conjunct(x: &Expr, items: &[(u64, u64)]) -> (Expr, Ty) {
        let kt = Ty::unsigned(3);
        let mut acc: Option<Expr> = None;
        for (i, &(lo, hi)) in items.iter().enumerate() {
            let term = bin(
                BinOp::LogAnd,
                is_code(1, kt, i as u64),
                in_range(x, lo, hi),
                Ty::BIT,
            );
            acc = Some(match acc {
                None => term,
                Some(p) => bin(BinOp::LogOr, p, term, Ty::BIT),
            });
        }
        (acc.expect("an item"), kt)
    }

    /// Draws `n` samples of `x` (var 0) from the plan "selector 1 weighted
    /// by `weights`, subject 0 distinct" and tallies each value.
    fn dist_hist(
        xt: Ty,
        items: &[(u64, u64)],
        weights: Vec<Weight>,
        extra: &[Expr],
        n: u32,
        seed0: u64,
    ) -> std::collections::BTreeMap<u32, u32> {
        let x = Expr::Var(0, xt);
        let (c, kt) = dist_conjunct(&x, items);
        let mut conj = vec![c];
        conj.extend_from_slice(extra);
        let sv = Solver::new(&[xt, kt], &[], &conj).unwrap();
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![0],
                weighted: vec![Weighted {
                    var: 1,
                    weights,
                    measure: Measure::Distinct(0),
                }],
            }],
        );
        let mut seed = seed0;
        let mut hist = std::collections::BTreeMap::new();
        for _ in 0..n {
            let v = sv.sample_plan(&p, &mut [], &mut seed).unwrap();
            let x = u32::try_from(v[0].clone()).unwrap();
            *hist.entry(x).or_insert(0) += 1;
        }
        hist
    }

    #[test]
    fn a_dist_with_an_excluded_value_renormalises_over_the_rest() {
        // LRM 18.5.3: `x != 200` with `{100 := 1, 200 := 2, 300 := 5}` leaves
        // 100 and 300 in the ratio 1 : 5 (N = 36000, 100 expected 6000,
        // sigma 71).
        let t = Ty::unsigned(9);
        let ne = bin(
            BinOp::Ne,
            Expr::Var(0, t),
            Expr::constant(t, 200u32),
            Ty::BIT,
        );
        let h = dist_hist(
            t,
            &[(100, 100), (200, 200), (300, 300)],
            ints(&[1, 2, 5]),
            &[ne],
            36_000,
            41,
        );
        assert_eq!(
            h.keys().copied().collect::<Vec<_>>(),
            vec![100, 300],
            "{h:?}"
        );
        assert!((5_600..6_400).contains(&h[&100]), "{h:?}");
    }

    #[test]
    fn a_partly_excluded_range_item_keeps_only_its_feasible_values() {
        // `[100:102] := 3, 103 := 1` with `x > 101`: the range keeps only 102
        // (weight 3 per value), against 103 (weight 1): 3 : 1. An item-level
        // weighting (3 * 3 values vs 1) would give 9 : 1.
        let t = Ty::unsigned(8);
        let gt = bin(
            BinOp::Gt,
            Expr::Var(0, t),
            Expr::constant(t, 101u32),
            Ty::BIT,
        );
        let h = dist_hist(
            t,
            &[(100, 102), (103, 103)],
            ints(&[3, 1]),
            &[gt],
            40_000,
            42,
        );
        assert_eq!(
            h.keys().copied().collect::<Vec<_>>(),
            vec![102, 103],
            "{h:?}"
        );
        // 30000 expected, sigma 87.
        assert!((29_600..30_400).contains(&h[&102]), "{h:?}");
    }

    #[test]
    fn overlapping_items_add_their_per_value_weights() {
        // `[0:3] := 1, [2:5] := 1`: 2 and 3 are in both, so they weigh 2
        // against 1 for 0, 1, 4, 5 (LRM 18.5.3: weights of a value in several
        // items add). Total weight 8; N = 48000, expected 6000 and 12000.
        let t = Ty::unsigned(4);
        let h = dist_hist(t, &[(0, 3), (2, 5)], ints(&[1, 1]), &[], 48_000, 43);
        for v in [0, 1, 4, 5] {
            assert!((5_600..6_400).contains(&h[&v]), "{v}: {h:?}");
        }
        for v in [2, 3] {
            assert!((11_500..12_500).contains(&h[&v]), "{v}: {h:?}");
        }
        assert_eq!(h.len(), 6, "{h:?}");
    }

    #[test]
    fn a_range_weight_spread_with_a_denominator_is_per_value() {
        // `[0:3] :/ 4` gives each of its four values 4/4 = 1; with `7 := 1`
        // all five values are equally likely (N = 40000, 8000 each, sigma 80).
        let t = Ty::unsigned(4);
        let w = vec![
            Weight {
                num: 4u32.into(),
                den: 4u32.into(),
            },
            Weight::int(1u32),
        ];
        let h = dist_hist(t, &[(0, 3), (7, 7)], w, &[], 40_000, 44);
        assert_eq!(h.len(), 5, "{h:?}");
        for (v, c) in &h {
            assert!((7_550..8_450).contains(c), "{v}: {h:?}");
        }
    }

    #[test]
    fn the_dist_subject_is_uniform_over_its_values_not_over_completions() {
        // x and y are 2 bits with `y <= x`; x = 3 has four completions and
        // x = 0 only one. The dist `[0:3] := 1` makes x uniform over its four
        // feasible values regardless (N = 40000, 10000 each, sigma 87).
        let t = Ty::unsigned(2);
        let (x, y) = (Expr::Var(0, t), Expr::Var(2, t));
        let (c, kt) = dist_conjunct(&x, &[(0, 3)]);
        let le = bin(BinOp::Le, y, x, Ty::BIT);
        let sv = Solver::new(&[t, kt, t], &[], &[c, le]).unwrap();
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![0, 2],
                weighted: vec![Weighted {
                    var: 1,
                    weights: ints(&[1]),
                    measure: Measure::Distinct(0),
                }],
            }],
        );
        let mut seed = 45u64;
        let (mut hx, mut ymax) = ([0u32; 4], 0u32);
        for _ in 0..40_000 {
            let v = sv.sample_plan(&p, &mut [], &mut seed).unwrap();
            assert!(v[2] <= v[0], "{v:?}");
            hx[usize::try_from(v[0].clone()).unwrap()] += 1;
            ymax += u32::from(v[0] == BigUint::from(3u8) && v[2] == BigUint::from(3u8));
        }
        for h in hx {
            assert!((9_550..10_450).contains(&h), "{hx:?}");
        }
        // Given x == 3, y is uniform over 0..=3: a quarter of 10000.
        assert!((2_250..2_750).contains(&ymax), "{ymax}");
    }

    #[test]
    fn the_models_measure_weights_a_selector_by_its_completions() {
        // The same problem with `Measure::Models`: code 0 stands for each
        // (x, y) tuple, so x is weighted by its completions 1 : 2 : 3 : 4
        // (N = 40000, x == 3 expected 16000, sigma 98).
        let t = Ty::unsigned(2);
        let (x, y) = (Expr::Var(0, t), Expr::Var(2, t));
        let (c, kt) = dist_conjunct(&x, &[(0, 3)]);
        let le = bin(BinOp::Le, y, x, Ty::BIT);
        let sv = Solver::new(&[t, kt, t], &[], &[c, le]).unwrap();
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![0, 2],
                weighted: vec![Weighted {
                    var: 1,
                    weights: ints(&[1]),
                    measure: Measure::Models,
                }],
            }],
        );
        let mut seed = 46u64;
        let mut x3 = 0u32;
        for _ in 0..40_000 {
            let v = sv.sample_plan(&p, &mut [], &mut seed).unwrap();
            x3 += u32::from(v[0] == BigUint::from(3u8));
        }
        assert!((15_500..16_500).contains(&x3), "{x3}");
    }

    #[test]
    fn a_selector_with_no_positive_feasible_weight_is_unsatisfiable() {
        // Every item has weight 0: nothing can be chosen (LRM 18.5.3: a
        // zero-weight value is never selected).
        let t = Ty::unsigned(4);
        let x = Expr::Var(0, t);
        let (c, kt) = dist_conjunct(&x, &[(0, 3), (5, 5)]);
        let sv = Solver::new(&[t, kt], &[], &[c]).unwrap();
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![0],
                weighted: vec![Weighted {
                    var: 1,
                    weights: ints(&[0, 0]),
                    measure: Measure::Distinct(0),
                }],
            }],
        );
        let mut seed = 47u64;
        assert_eq!(
            sv.sample_plan(&p, &mut [], &mut seed),
            Err(SolveError::Unsatisfiable)
        );
        assert_eq!(seed, 47, "no draws consumed");
    }

    #[test]
    fn a_distinct_subject_outside_the_group_or_a_zero_denominator_is_a_bad_plan() {
        let t = Ty::unsigned(4);
        let x = Expr::Var(0, t);
        let (c, kt) = dist_conjunct(&x, &[(0, 3)]);
        let sv = Solver::new(&[t, kt], &[], &[c]).unwrap();
        let mut seed = 48u64;
        let weighted = |w: Vec<Weight>| Weighted {
            var: 1,
            weights: w,
            measure: Measure::Distinct(0),
        };
        // The subject is not listed in the selector's group.
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![],
                weighted: vec![weighted(ints(&[1]))],
            }],
        );
        assert!(matches!(
            sv.sample_plan(&p, &mut [], &mut seed),
            Err(SolveError::BadPlan { .. })
        ));
        // A weight whose denominator is 0.
        let bad = vec![Weight {
            num: 1u32.into(),
            den: 0u32.into(),
        }];
        let p = plan(
            vec![],
            vec![Group {
                vars: vec![0],
                weighted: vec![weighted(bad)],
            }],
        );
        assert!(matches!(
            sv.sample_plan(&p, &mut [], &mut seed),
            Err(SolveError::BadPlan { .. })
        ));
    }

    /// Bits outside the constraints' support are drawn without conjoining a
    /// literal (`sample_plan`); the draw must be exactly the one the per-bit
    /// walk over every bit makes. The expected values were produced by the
    /// walk before that shortcut existed.
    #[test]
    fn unconstrained_bits_draw_as_the_full_walk_does() {
        let x = Expr::Var(0, Ty::unsigned(8));
        let c = bin(
            BinOp::Lt,
            x,
            Expr::constant(Ty::unsigned(8), 100u8),
            Ty::BIT,
        );
        let tys = [Ty::unsigned(8), Ty::unsigned(8), Ty::unsigned(4)];
        let sv = Solver::new(&tys, &[], &[c]).expect("compiles");
        let mut seed = 7;
        let mut got = Vec::new();
        for _ in 0..6 {
            let v = sv.sample(&[1, 0, 2], &mut seed).expect("sat");
            assert!(v[0] < BigUint::from(100u8));
            got.push(
                v.iter()
                    .map(ToString::to_string)
                    .collect::<Vec<_>>()
                    .join(","),
            );
        }
        assert_eq!(
            got.join(" "),
            "21,216,3 52,172,1 61,93,5 42,49,7 18,1,0 31,165,6"
        );
    }

    /// A repeated compile is served from the cache with the same count.
    #[test]
    fn repeated_compile_is_cached_and_equal() {
        let x = Expr::Var(0, Ty::unsigned(6));
        let c = bin(BinOp::Lt, x, Expr::constant(Ty::unsigned(6), 37u8), Ty::BIT);
        let tys = [Ty::unsigned(6)];
        let a = Solver::new(&tys, &[], std::slice::from_ref(&c)).unwrap();
        let b = Solver::new(&tys, &[], std::slice::from_ref(&c)).unwrap();
        assert_eq!(a.count(), &BigUint::from(37u32));
        assert_eq!(a.count(), b.count());
        let digest = problem_digest(&tys, &[], &[c], NODE_BUDGET);
        assert!(SOLVER_CACHE.lock().unwrap().contains_key(&digest));
    }
}