probl_engine/

interp.rs

//! The world-set interpreter.
//!
//! Every statement maps a set of worlds to a set of worlds. All the worlds in
//! a set are at the same point of the program, so control flow is shared and
//! only the data differs. Splitting happens in statements; expressions are
//! evaluated in one world at a time. See docs/semantics.md for the rules this
//! implements.

use crate::analytic::{self, Analytic};
use crate::builtins;
use crate::chain::{Chain, Solution};
use crate::conjugate::{self, Seen};
use crate::continuous::{Family, Rng};
use crate::dist::{Budget, Counts, Dist};
use crate::error::{ErrorKind, Fault, OpError, OpResult, Result, RuntimeError};
use crate::failure::Failures;
use crate::ops::{self, Truth};
use crate::report::Sink;
use crate::value::{Closure, Delayed, Value, fmt_prob};
use crate::weight::Weight;
use crate::world::Returned;
use crate::world::{Flow, World, clear, clear_dead, live_slots, merge, merge_values, state_hash, total_weight};
use probl_sema::builtins::Lifting;
use probl_sema::conjugate::{Conjugacy, Likelihood, Update};
use probl_sema::effects::EvidenceOrder;
use probl_sema::ir::*;
use probl_sema::{Builtin, Liveness};
use probl_syntax::Span;
use probl_syntax::ast::BinOp;
use rustc_hash::{FxHashMap, FxHashSet};
use std::collections::BTreeMap;
use std::sync::Arc;
use std::sync::atomic::{AtomicBool, Ordering};

#[derive(Clone, Debug, Default)]
pub struct Stats {
    /// Largest number of worlds any statement ran on.
    pub peak_worlds: usize,
    /// Statement executions, counting each world separately.
    pub world_steps: u64,
    pub calls: u64,
    pub memo_hits: u64,
    /// Loops solved as Markov chains, and their states (section 10).
    pub solved_loops: u64,
    pub chain_states: u64,
    /// Calls that came back to themselves, solved by iteration, and the
    /// rounds that took (section 6).
    pub solved_calls: u64,
    pub call_rounds: u64,
    /// When sampling: what happened to each variable whose draws may be
    /// delayed (`Conjugacy::variables`).
    pub updates: Vec<Updates>,
    /// When the program was answered with decision diagrams instead of
    /// worlds (`Options::symbolic`): their size.
    pub diagrams: Option<crate::symbolic::Diagrams>,
}

/// How often a variable's draws were delayed and updated exactly
/// (docs/semantics.md, section 14), counted over all runs.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct Updates {
    /// Draws delayed.
    pub delayed: u64,
    /// Observations that updated it exactly.
    pub exact: u64,
    /// Times it was drawn when first needed, from its updated distribution.
    pub drawn: u64,
}

impl Stats {
    /// Add up another batch's statistics.
    pub fn absorb(&mut self, other: &Stats) {
        self.peak_worlds = self.peak_worlds.max(other.peak_worlds);
        self.world_steps += other.world_steps;
        self.calls += other.calls;
        self.memo_hits += other.memo_hits;
        self.solved_loops += other.solved_loops;
        self.chain_states += other.chain_states;
        self.solved_calls += other.solved_calls;
        self.call_rounds += other.call_rounds;
        if self.updates.len() < other.updates.len() {
            self.updates.resize(other.updates.len(), Updates::default());
        }
        for (mine, theirs) in self.updates.iter_mut().zip(&other.updates) {
            mine.delayed += theirs.delayed;
            mine.exact += theirs.exact;
            mine.drawn += theirs.drawn;
        }
    }
}

/// How the engine runs, with the host's limits already applied.
#[derive(Clone, Debug)]
pub struct Config {
    pub today: Option<i32>,
    pub epsilon: f64,
    pub merging: bool,
    pub memoizing: bool,
    pub max_worlds: usize,
    pub max_iterations: u64,
    pub max_call_depth: usize,
    pub max_cached_calls: usize,
    pub max_output: usize,
    pub budget: Budget,
    pub cancel: Option<Arc<AtomicBool>>,
    /// Sample with this seed, instead of enumerating (docs/semantics.md,
    /// section 14).
    pub sample_seed: Option<u64>,
    /// When sampling, delay the draws of conjugate priors and update them
    /// exactly (section 14).
    pub conjugate: bool,
    /// When enumerating, solve loops that cycle as Markov chains (section 10).
    pub solving: bool,
    /// The most states a loop's chain may have to be solved.
    pub max_chain_states: usize,
    /// A fault ends only the world it happens in (docs/semantics.md,
    /// section 11), instead of the run.
    pub partial: bool,
    /// Compute what has no formula by numerical integration, to this
    /// relative tolerance (section 13).
    pub integrate: Option<f64>,
}

/// The most steps eliminating one loop's chain may take before the loop is
/// unrolled instead.
const MAX_ELIMINATION: u64 = 200_000_000;

/// Runs are executed this many at a time, and each batch has a random
/// stream of its own. The size is fixed: changing it changes the output for
/// a given seed.
pub const BATCH: u64 = 1000;

/// What a batch of sampled runs produced, to be combined with the other
/// batches in order.
#[derive(Debug)]
pub struct Batch {
    pub sinks: Vec<Sink>,
    pub totals: SampleTotals,
    pub observed: bool,
    pub densities: bool,
    pub unresolved: Weight,
    pub last_ruling_out: Option<Span>,
    pub stats: Stats,
    pub failures: Failures,
}

/// Lines printed by `print(…)`, with where.
pub type Printed = Vec<(Span, String)>;

/// The final weights of all the runs, when sampling.
#[derive(Clone, Copy, Debug)]
pub struct SampleTotals {
    pub weight: Weight,
    /// Sum of the squared weights, for the effective sample size.
    pub squares: Weight,
}

/// The distribution of a function's return value, unnormalized.
#[derive(Clone, Debug)]
pub struct CallResult {
    pub outcomes: Vec<Returned>,
    /// Weight the call left unresolved.
    pub unresolved: Weight,
    /// Weight its observations ruled out, when enumerating.
    pub lost: Weight,
    /// Whether the call ran `observe` (outside any `simulate`).
    pub observed: bool,
    /// The part of `unresolved` waiting on calls that came back to
    /// themselves and are still being solved (section 6), for each by its
    /// depth in the stack. A result waiting on one isn't final.
    pub pending: Vec<(usize, Weight)>,
    /// In partial mode: its worlds that failed.
    pub failures: Failures,
}

impl CallResult {
    /// The result of a call that came back to itself, before the first
    /// round: nothing resolved, all of it waiting on that call.
    fn waiting(depth: usize) -> CallResult {
        CallResult {
            outcomes: Vec::new(),
            unresolved: Weight::ONE,
            lost: Weight::ZERO,
            observed: false,
            pending: vec![(depth, Weight::ONE)],
            failures: Failures::default(),
        }
    }

    /// Take out the weight waiting on the call at `depth`.
    fn take_pending(&mut self, depth: usize) -> Weight {
        match self.pending.iter().position(|&(d, _)| d == depth) {
            Some(i) => self.pending.remove(i).1,
            None => Weight::ZERO,
        }
    }
}

/// A call's result so far, while the call it waits on is being solved
/// (section 6).
struct Approx {
    result: Arc<CallResult>,
    /// The depth of that call, and its frame and round when this was worked
    /// out. It stands for the call as long as that frame runs, and can be
    /// reused as it is until that frame's next round.
    head: usize,
    frame: u64,
    round: u64,
}

/// Add `weight` times what's waiting in `from` to `into`.
fn add_pending(into: &mut Vec<(usize, Weight)>, from: &[(usize, Weight)], weight: Weight) {
    for &(depth, w) in from {
        match into.iter_mut().find(|(d, _)| *d == depth) {
            Some((_, total)) => *total += weight * w,
            None => into.push((depth, weight * w)),
        }
    }
}

type CallKey = (FnId, Vec<Value>);

pub struct Engine<'p> {
    prog: &'p Program,
    live: &'p Liveness,
    /// Which draws may be delayed, and what must draw them.
    conj: &'p Conjugacy,
    config: Config,
    budget: Budget,
    memo: FxHashMap<CallKey, Arc<CallResult>>,
    /// `simulate` results that are distributions whose outcomes own their
    /// draws: pure values of the block and its captures, which a call's
    /// result, holding analytic outcomes, can't be cached as.
    simulated: FxHashMap<CallKey, Value>,
    /// The calls running, with their depth in the stack.
    active: FxHashMap<CallKey, usize>,
    /// Results so far of calls that wait on a call being solved by
    /// iteration (section 6).
    approx: FxHashMap<CallKey, Approx>,
    /// For each depth in the stack: the frame running there, and its round,
    /// each a number never used before.
    frames_at: Vec<u64>,
    rounds_at: Vec<u64>,
    last_number: u64,
    /// Weight waiting on calls being solved, in the current call.
    pending: Vec<(usize, Weight)>,
    depth: usize,
    /// Weight left unaccounted for (loops cut short, truncated tails).
    pub unresolved: Weight,
    /// When enumerating: weight ruled out by observations, which solving a
    /// loop counts as a way out of it.
    lost: Weight,
    /// For each statement: whether it's an unbounded loop that may be
    /// solved as a Markov chain (section 10).
    solvable: Vec<bool>,
    /// For each loop: the variables it only adds to, which the chain carries
    /// as rewards rather than in its states.
    rewards: FxHashMap<StmtId, Vec<SlotId>>,
    /// Loops whose chains were found too large to solve: they're unrolled
    /// from then on.
    too_large: FxHashSet<StmtId>,
    /// Whether the current inference scope has run `observe`.
    pub observed: bool,
    /// Whether it has observed a value with a density, which makes the
    /// evidence a density too (docs/semantics.md, section 14).
    pub densities: bool,
    /// The last observation that ruled out a world, for impossible evidence.
    pub last_ruling_out: Option<Span>,
    /// When enumerating, the first observation of a density.
    pub density_at: Option<Span>,
    /// How many densities the weights of the worlds that reported or
    /// finished include. It must be as many for each.
    pub reported_densities: Option<u32>,
    /// In the innermost `simulate` running: its first observation of a
    /// density, and how many its returning worlds' weights include, which
    /// must be as many for each, since its result is normalized.
    local_density_at: Option<Span>,
    returned_densities: Option<u32>,
    /// The tolerance, if a result is computed by numerical integration.
    pub integrated: Option<f64>,
    pub sinks: Vec<Sink>,
    dice: FxHashMap<(u32, u32), Value>,
    pools: FxHashMap<(u32, Value), Value>,
    record_names: Vec<Arc<str>>,
    enum_values: Vec<Vec<Value>>,
    print: &'p mut (dyn FnMut(&str) + Send),
    /// The values of the program's inputs (`read`).
    inputs: &'p [Value],
    /// Bytes printed, counting the batches printed before this one.
    printed: usize,
    /// When sampling a batch: what it printed, kept to be printed in batch
    /// order.
    lines: Option<Printed>,
    pub stats: Stats,
    /// When sampling: the random numbers. `None` when enumerating, including
    /// inside `simulate` while sampling.
    sampler: Option<Rng>,
    /// Nesting depth of local `simulate` inference (in either outer mode).
    nested: usize,
    next_latent: u64,
    /// Dynamic effect boundary for collection callbacks. Effects inside an
    /// explicit `simulate` are local to that computation instead.
    callback: Option<&'static str>,
    /// In partial mode: the worlds that failed, in the current call.
    pub failures: Failures,
    /// The `try`s running, innermost last, each with the depth of the call
    /// it's in, and its catches.
    handlers: Vec<(usize, &'p [Catch])>,
    /// Worlds of the statement running that faulted inside a `try` of this
    /// call that catches the fault, on their way to it.
    faulted: Vec<(World, RuntimeError)>,
    /// Where evidence can still be applied, to tell whether a failed
    /// world's weight can be added to the finished worlds'.
    evidence: EvidenceOrder,
}

/// Where a world was when it failed.
#[derive(Clone, Copy)]
struct At {
    f: FnId,
    stmt: StmtId,
    weight: Weight,
    run: u32,
    densities: u32,
}

impl At {
    fn new(f: FnId, stmt: &Stmt, w: &World) -> At {
        At {
            f,
            stmt: stmt.id,
            weight: w.weight,
            run: w.run,
            densities: w.densities,
        }
    }
}

/// The value of `$e` in one world. If it fails with a fault that a `catch`
/// or partial mode takes, the world leaves the statement: `$undo` runs, and
/// the statement goes on with the next world (see `Engine::fail`). A catch
/// gets the world as its statement began: `world w` when `$e` leaves `w` as
/// it was, or `saved s` with a copy `s` taken before `$e` changed it.
macro_rules! each {
    ($self:ident, $at:expr, world $w:expr, $e:expr $(, $undo:block)?) => {
        match $e {
            Ok(v) => v,
            Err(e) => {
                let saved = $self.snapshot(&$w);
                $self.fail(e, $at, saved)?;
                $($undo)?
                continue;
            }
        }
    };
    ($self:ident, $at:expr, saved $saved:ident, $e:expr $(, $undo:block)?) => {
        match $e {
            Ok(v) => v,
            Err(e) => {
                $self.fail(e, $at, $saved.take())?;
                $($undo)?
                continue;
            }
        }
    };
}

impl<'p> Engine<'p> {
    pub fn new(
        prog: &'p Program,
        live: &'p Liveness,
        conj: &'p Conjugacy,
        config: Config,
        inputs: &'p [Value],
        print: &'p mut (dyn FnMut(&str) + Send),
    ) -> Engine<'p> {
        let sampler = config.sample_seed.map(Rng::new);
        Engine {
            prog,
            live,
            conj,
            budget: config.budget.clone(),
            config,
            memo: FxHashMap::default(),
            simulated: FxHashMap::default(),
            active: FxHashMap::default(),
            approx: FxHashMap::default(),
            frames_at: Vec::new(),
            rounds_at: Vec::new(),
            last_number: 0,
            pending: Vec::new(),
            depth: 0,
            unresolved: Weight::ZERO,
            lost: Weight::ZERO,
            solvable: probl_sema::effects::solvable_loops(prog),
            rewards: probl_sema::rewards::loop_rewards(prog),
            too_large: FxHashSet::default(),
            observed: false,
            densities: false,
            last_ruling_out: None,
            density_at: None,
            reported_densities: None,
            local_density_at: None,
            returned_densities: None,
            integrated: None,
            sinks: vec![Sink::default(); prog.reports.len()],
            dice: FxHashMap::default(),
            pools: FxHashMap::default(),
            record_names: prog.records.iter().map(|r| Arc::from(r.name.as_str())).collect(),
            enum_values: prog
                .enums
                .iter()
                .enumerate()
                .map(|(t, e)| {
                    e.variants
                        .iter()
                        .enumerate()
                        .map(|(v, name)| ops::enum_value(t as u32, v as u32, name))
                        .collect()
                })
                .collect(),
            print,
            inputs,
            printed: 0,
            lines: None,
            stats: Stats::default(),
            sampler,
            nested: 0,
            next_latent: 0,
            callback: None,
            failures: Failures::default(),
            evidence: probl_sema::effects::evidence_order(prog),
            handlers: Vec::new(),
            faulted: Vec::new(),
        }
    }

    /// Run the top level. Returns the total weight of the worlds that finish.
    pub fn run_main(&mut self) -> Result<Weight> {
        let prog = self.prog;
        let main = prog.main();
        let world = World {
            slots: vec![Value::Dead; main.n_slots()],
            constraints: Default::default(),
            inherited: Default::default(),
            weight: Weight::ONE,
            run: 0,
            densities: 0,
        };
        let flow = self.exec_block(MAIN, &main.body, vec![world])?;
        uncaught(&flow)?;
        for w in &flow.next {
            self.same_units(w, main.span)?;
        }
        Ok(total_weight(&flow.next))
    }

    /// Weights with different numbers of densities aren't in the same
    /// units, so they can't be added up or compared.
    fn same_units(&mut self, w: &World, span: Span) -> Result<()> {
        match self.reported_densities {
            Some(n) if n != w.densities => Err(OpError::unsupported(
                "worlds that observe different numbers of continuous values aren't supported when enumerating",
            )
            .help("a density is per unit of what it observes, so worlds with more of them can't be weighed against the others; observe the same values in every world")
            .at(self.density_at.unwrap_or(span))),
            _ => {
                self.reported_densities = Some(w.densities);
                Ok(())
            }
        }
    }

    /// Sample the runs `first..first + n` with the random numbers `rng`
    /// (docs/semantics.md, section 14), after earlier batches printed
    /// `printed` bytes. A batch starts afresh: its reports and counters are
    /// its own, so batches can run on any thread; only caches carry over.
    /// What it prints comes back with its result, even when it fails.
    pub fn run_batch(&mut self, rng: Rng, first: u64, n: u64, printed: usize) -> (Result<Batch>, Printed) {
        self.sampler = Some(rng);
        self.sinks = vec![Sink::default(); self.prog.reports.len()];
        self.observed = false;
        self.densities = false;
        self.unresolved = Weight::ZERO;
        self.last_ruling_out = None;
        self.stats = Stats::default();
        self.failures = Failures::default();
        self.printed = printed;
        self.lines = Some(Vec::new());
        self.budget = self.config.budget.clone();
        let result = self.sample_runs(first, n);
        self.budget.give_back();
        (result, self.lines.take().unwrap_or_default())
    }

    fn sample_runs(&mut self, first: u64, n: u64) -> Result<Batch> {
        let main = self.prog.main();
        let worlds = (first..first + n)
            .map(|run| World {
                slots: vec![Value::Dead; main.n_slots()],
                constraints: Default::default(),
                inherited: Default::default(),
                weight: Weight::ONE,
                run: run as u32,
                densities: 0,
            })
            .collect();
        let flow = self.exec_block(MAIN, &main.body, worlds)?;
        uncaught(&flow)?;
        let mut totals = SampleTotals {
            weight: Weight::ZERO,
            squares: Weight::ZERO,
        };
        for w in &flow.next {
            totals.weight += w.weight;
            totals.squares += w.weight * w.weight;
        }
        for sink in &mut self.sinks {
            sink.end_batch();
        }
        Ok(Batch {
            sinks: std::mem::take(&mut self.sinks),
            totals,
            observed: self.observed,
            densities: self.densities,
            unresolved: self.unresolved,
            last_ruling_out: self.last_ruling_out,
            stats: std::mem::take(&mut self.stats),
            failures: std::mem::take(&mut self.failures),
        })
    }

    /// When sampling: one draw from a distribution. A finite distribution
    /// gives one of its outcomes, chosen with its probability (missing mass
    /// is ignored); a continuous one, a number.
    fn sample(&mut self, v: &Value) -> Value {
        let rng = self.sampler.as_mut().expect("only called when sampling");
        match v {
            Value::Dist(d) if !d.outcomes.is_empty() => {
                let i = rng.choose(d.outcomes.iter().map(|(_, p)| *p)).unwrap_or(0);
                match &d.outcomes[i].0 {
                    Value::Continuous(f) => Value::Float(f.sample(rng)),
                    Value::Joint(j) => crate::joint::sample(j, rng),
                    x => x.clone(),
                }
            }
            Value::Continuous(f) => Value::Float(f.sample(rng)),
            Value::Counts(c) => Value::Int(c.sample(rng).into()),
            Value::Joint(j) => crate::joint::sample(j, rng),
            other => other.clone(),
        }
    }

    /// When sampling, a value to report: continuous parts are drawn; a
    /// finite distribution is kept whole, so each outcome counts with its
    /// probability.
    #[inline(never)]
    fn sample_continuous(&mut self, v: Value) -> Value {
        match &v {
            Value::Continuous(_) | Value::Joint(_) => self.sample(&v),
            Value::Dist(d)
                if d.outcomes
                    .iter()
                    .any(|(x, _)| matches!(x, Value::Continuous(_) | Value::Joint(_))) =>
            {
                let pairs = d.outcomes.iter().map(|(x, p)| (self.sample(x), *p)).collect();
                Dist::from_pairs(pairs, d.missing).into_value()
            }
            _ => v,
        }
    }

    /// When sampling, `binomial(…)`, `poisson(…)` or `geometric(…)` with plain
    /// arguments is drawn from, or observed, without listing its outcomes:
    /// the distribution is the same, and it's much faster.
    fn direct_counts(&mut self, f: FnId, e: &'p Expr, w: &World) -> Result<Option<Counts>> {
        let ExprKind::Builtin {
            func: b @ (Builtin::Binomial | Builtin::Poisson | Builtin::Geometric),
            args,
            ..
        } = &e.kind
        else {
            return Ok(None);
        };
        if self.sampler.is_none() {
            return Ok(None);
        }
        let mut values = Vec::with_capacity(args.len());
        for a in args {
            values.push(self.eval(f, a, w)?);
        }
        let counts = builtins::counts(*b, &values, &mut self.budget).map_err(|err| err.at(e.span))?;
        Ok(counts.filter(Counts::direct))
    }

    /// Whether draws of conjugate priors are delayed now: when sampling,
    /// but not inside `simulate` (section 14).
    fn delaying(&self) -> bool {
        self.config.conjugate && self.sampler.is_some()
    }

    fn updates(&mut self, variable: u32) -> &mut Updates {
        let i = variable as usize;
        if self.stats.updates.len() <= i {
            let n = self.conj.variables.len().max(i + 1);
            self.stats.updates.resize(n, Updates::default());
        }
        &mut self.stats.updates[i]
    }

    /// Draw a delayed variable from its distribution, updated by the
    /// observations so far, now that its value is needed.
    fn draw_delayed(&mut self, w: &mut World, slot: SlotId) {
        if let Value::Delayed(d) = &w.slots[slot as usize] {
            let d = **d;
            let x = d
                .family
                .sample(self.sampler.as_mut().expect("delayed only when sampling"));
            w.slots[slot as usize] = Value::Float(x);
            self.updates(d.variable).drawn += 1;
        }
    }

    /// `observe` of a delayed variable through a conjugate form: the log of
    /// the observation's probability (or density), with the variable
    /// updated to its distribution after it. `None` when it isn't an exact
    /// update: the variable isn't delayed, or doesn't have the right family,
    /// or the rest of the observation isn't plain. Then the variable is
    /// drawn, and the observation is made as usual. `from` is where
    /// problems with the observed distribution's parameters are reported,
    /// as when it's evaluated.
    fn observe_exactly(&mut self, f: FnId, u: &Update<'p>, w: &mut World, from: Span) -> Result<Option<f64>> {
        let Value::Delayed(d) = &w.slots[u.slot as usize] else {
            return Ok(None);
        };
        let d = **d;
        let pair = matches!(
            (d.family, u.likelihood),
            (
                Family::Beta { .. },
                Likelihood::Binomial { .. } | Likelihood::Bernoulli { .. }
            ) | (Family::Gamma { .. }, Likelihood::Poisson { .. })
                | (Family::Normal { .. }, Likelihood::Normal { .. })
        );
        let seen = if pair { self.seen(f, u, w, from)? } else { None };
        let Some(seen) = seen else {
            self.draw_delayed(w, u.slot);
            return Ok(None);
        };
        let (ln, posterior) = conjugate::update(&d.family, seen).expect("a conjugate pair");
        w.slots[u.slot as usize] = Value::Delayed(Arc::new(Delayed {
            family: posterior,
            variable: d.variable,
        }));
        self.updates(d.variable).exact += 1;
        self.densities |= matches!(seen, Seen::Normal { .. });
        Ok(Some(ln))
    }

    /// Assigning a discrete outcome, like `floor(x)`, splits its world: one
    /// for each of its values, with the draw restricted to where it's that
    /// value, so that the variable is a plain number. The split worlds go
    /// to `failed`, to be dropped, and theirs to the end. (A place with a
    /// path keeps the outcome.)
    #[cold]
    #[inline(never)]
    fn split_discrete(
        &mut self,
        f: FnId,
        stmt: &'p Stmt,
        place: &Place,
        worlds: &mut Vec<World>,
        failed: &mut Vec<usize>,
    ) -> Result<()> {
        if !place.path.is_empty() {
            return Ok(());
        }
        let mut split = Vec::new();
        let mut gone = Vec::new();
        for (i, w) in worlds.iter().enumerate() {
            let Value::Analytic(a) = &w.slots[place.slot as usize] else {
                continue;
            };
            // `failed` is in order: the worlds that failed to assign.
            if !a.is_discrete() || failed.binary_search(&i).is_ok() {
                continue;
            }
            if a.form.is_some() {
                return Err(
                    analytic::joint_condition("assigning a rounded outcome of several normal draws").at(stmt.span),
                );
            }
            for (v, region, share) in a.atoms() {
                if share <= 0.0 {
                    continue;
                }
                let mut nw = w.clone().scaled(share);
                let latent = analytic::Latent {
                    family: a.family,
                    domain: a.domain.intersect(&region),
                    form: None,
                    weight: a.weight.clone(),
                };
                Arc::make_mut(&mut nw.constraints).insert(a.id, latent);
                if let Err(e) = self.assign(f, place, v, &mut nw, stmt.span) {
                    self.fail(e, At::new(f, stmt, &nw), None)?;
                    continue;
                }
                split.push(nw);
            }
            gone.push(i);
        }
        failed.extend(gone);
        failed.sort_unstable();
        worlds.extend(split);
        self.check_worlds(worlds.len() - failed.len(), stmt.span)
    }

    /// `if p` for a probability `p` drawn from a beta: each branch updates
    /// it, as observing that side from `bernoulli(p)` would. Whether it was
    /// such a `p`. (Out of `exec_stmt_kind`, like the next two, which every
    /// statement runs: inlined, they'd make its frame larger.)
    #[inline(never)]
    fn branch_on_draw(
        &mut self,
        x: &Analytic,
        w: &World,
        yes: &mut Vec<World>,
        no: &mut Vec<World>,
        span: Span,
    ) -> Result<bool> {
        let sides = [Seen::Bernoulli(true), Seen::Bernoulli(false)];
        let mut sides = self.outcomes_of(x, w, &sides, span)?.into_iter();
        let (Some(y), Some(n)) = (sides.next(), sides.next()) else {
            return Ok(false);
        };
        self.check_callback_effect(span)?;
        yes.extend(y);
        no.extend(n);
        Ok(true)
    }

    /// A draw whose distribution depends on an analytic draw, into `out`:
    /// from `normal(mean, sd)` with an outcome of normal draws as its mean
    /// (see `gaussian_draw`), or from `bernoulli(p)` or `binomial(n, p)`
    /// with a probability drawn from a beta (see `dependent_draw`). Whether
    /// it was one.
    #[inline(never)]
    fn draw_dependent(
        &mut self,
        f: FnId,
        place: &Place,
        dist: &'p Expr,
        w: &World,
        out: &mut Vec<World>,
    ) -> Result<bool> {
        if let Some(v) = self.gaussian_draw(f, dist, w)? {
            let mut w = w.clone();
            self.assign(f, place, v, &mut w, dist.span)?;
            out.push(w);
            return Ok(true);
        }
        let Some(outcomes) = self.dependent_draw(f, dist, w)? else {
            return Ok(false);
        };
        let mark = out.len();
        for (v, mut w) in outcomes {
            if let Err(e) = self.assign(f, place, v, &mut w, dist.span) {
                out.truncate(mark);
                return Err(e);
            }
            out.push(w);
        }
        Ok(true)
    }

    /// Update an analytic draw exactly when enumerating: the observation's
    /// probability, as a logarithm, or `None` if it isn't an update this
    /// can make.
    #[inline(never)]
    fn observe_posterior(&mut self, f: FnId, u: &Update<'p>, w: &mut World, from: Span) -> Result<Option<f64>> {
        let Some(x) = self.draw_in(w, u.slot, from)? else {
            return Ok(None);
        };
        let Some(seen) = self.seen(f, u, w, from)? else {
            return Ok(None);
        };
        let Some((ln, latent)) = analytic::update(&x, seen).map_err(|e| e.at(from))? else {
            return Ok(None);
        };
        if matches!(seen, Seen::Normal { .. }) {
            self.check_density(f, from)?;
            self.densities = true;
            w.densities += 1;
        }
        self.know(w, x.id, latent, from)?;
        Ok(Some(ln))
    }

    /// A draw from `normal(mean, sd)` when enumerating, with `mean` an
    /// outcome of normal draws: the mean plus a fresh normal draw. `None`
    /// for another draw.
    fn gaussian_draw(&mut self, f: FnId, dist: &'p Expr, w: &World) -> Result<Option<Value>> {
        let ExprKind::Builtin {
            func: Builtin::Normal,
            args,
            named,
        } = &dist.kind
        else {
            return Ok(None);
        };
        let ([mean, sd], true) = (args.as_slice(), named.is_empty()) else {
            return Ok(None);
        };
        let Value::Analytic(m) = self.eval(f, mean, w)? else {
            return Ok(None);
        };
        let Some(sd) = self.plain_sd(f, sd, w, dist.span)? else {
            return Ok(None);
        };
        self.next_latent += 1;
        match analytic::normal_draw(&m, sd, self.next_latent).map_err(|e| e.at(dist.span))? {
            Some(v) => Ok(Some(v)),
            None => Err(
                analytic::unsupported("drawing from a normal distribution whose mean is this outcome").at(mean.span),
            ),
        }
    }

    /// `observe y from normal(mean, sd)` when enumerating, with `mean` an
    /// outcome of normal draws: the observation's density, as a logarithm,
    /// with the world updated exactly. `None` if the mean isn't an outcome.
    #[inline(never)]
    fn observe_gaussian(&mut self, f: FnId, v: &Value, d: &'p Expr, w: &mut World, span: Span) -> Result<Option<f64>> {
        let ExprKind::Builtin {
            func: Builtin::Normal,
            args,
            named,
        } = &d.kind
        else {
            return Ok(None);
        };
        let ([mean, sd], true) = (args.as_slice(), named.is_empty()) else {
            return Ok(None);
        };
        let Value::Analytic(m) = self.eval(f, mean, w)? else {
            return Ok(None);
        };
        let Some(sd) = self.plain_sd(f, sd, w, d.span)? else {
            return Ok(None);
        };
        let y = match v {
            Value::Int(_) | Value::Float(_) | Value::Prob(_) => v.as_f64(),
            _ => None,
        };
        let Some(y) = y else {
            return Ok(None);
        };
        self.check_density(f, span)?;
        let observed = analytic::observe_normal(&m, y, sd, &mut self.next_latent).map_err(|e| e.at(span))?;
        let Some((ln, known)) = observed else {
            if self.config.integrate.is_some() {
                return Ok(None);
            }
            return Err(analytic::unsupported(
                "observing a value from a normal distribution whose mean is this outcome",
            )
            .at(mean.span));
        };
        self.densities = true;
        w.densities += 1;
        for (id, latent) in known {
            Arc::make_mut(&mut w.constraints).insert(id, latent);
        }
        self.reread(w, span)?;
        Ok(Some(ln))
    }

    /// `observe v from D(…)` with `@integrate`, when parameters of `D` are
    /// outcomes of one draw: what it says about the draw, as a factor of
    /// its likelihood, and its probability or density, as a logarithm.
    /// `None` if it isn't such an observation.
    #[inline(never)]
    fn observe_numerically(
        &mut self,
        f: FnId,
        v: &Value,
        d: &'p Expr,
        w: &mut World,
        span: Span,
    ) -> Result<Option<f64>> {
        use crate::likelihood::{Factor, Law, Likelihood, Param};
        let Some(tolerance) = self.config.integrate else {
            return Ok(None);
        };
        let ExprKind::Builtin { func, args, named } = &d.kind else {
            return Ok(None);
        };
        // Each law's parameters, with where each is valid: from, to, and
        // whether each end is included.
        type Range = (Option<(f64, bool)>, Option<(f64, bool)>);
        let positive: Range = (Some((0.0, false)), None);
        let any: Range = (None, None);
        let probability: Range = (Some((0.0, true)), Some((1.0, true)));
        let (law, ranges): (Law, Vec<Range>) = match func {
            Builtin::Normal => (Law::Normal, vec![any, positive]),
            Builtin::Lognormal => (Law::Lognormal, vec![any, positive]),
            Builtin::Uniform => (Law::Uniform, vec![any, any]),
            Builtin::Beta => (Law::Beta, vec![positive, positive]),
            Builtin::Gamma => (Law::Gamma, vec![positive, positive]),
            Builtin::Exponential => (Law::Exponential, vec![positive]),
            Builtin::Binomial => (Law::Binomial, vec![any, probability]),
            Builtin::Bernoulli => (Law::Bernoulli, vec![probability]),
            Builtin::Poisson => (Law::Poisson, vec![(Some((0.0, true)), None)]),
            Builtin::Geometric => (Law::Geometric, vec![(Some((0.0, false)), Some((1.0, true)))]),
            _ => return Ok(None),
        };
        if !named.is_empty() || args.len() != ranges.len() {
            return Ok(None);
        }
        let mut values = Vec::with_capacity(args.len());
        for a in args {
            values.push(self.eval(f, a, w)?);
        }
        let mut draw: Option<Arc<Analytic>> = None;
        for value in &values {
            match value {
                Value::Analytic(a) => match &draw {
                    None => draw = Some(a.clone()),
                    Some(x) if x.same_axis(a) => {}
                    Some(_) => {
                        return Err(analytic::unsupported(
                            "an observation whose parameters are outcomes of different draws",
                        )
                        .at(d.span));
                    }
                },
                v if v.as_f64().is_some() => {}
                _ => return Ok(None),
            }
        }
        let Some(x) = draw else {
            return Ok(None);
        };
        if x.form.is_some() {
            return Err(
                analytic::unsupported("integrating over a sum of draws, which has several dimensions").at(d.span),
            );
        }
        let value = match (law, v) {
            (Law::Bernoulli, Value::Bool(b)) => f64::from(u8::from(*b)),
            (_, Value::Int(_) | Value::Float(_) | Value::Prob(_)) => match v.as_f64() {
                Some(y) => y,
                None => return Ok(None),
            },
            _ => return Ok(None),
        };
        // A whole number of trials, as `binomial` takes it.
        if law == Law::Binomial {
            let n = &values[0];
            if !matches!(n, Value::Int(_)) {
                return Err(analytic::unsupported("a number of trials that's an outcome of a draw").at(d.span));
            }
            builtins::call_plain(Builtin::Binomial, &[n.clone(), Value::Prob(0.5)], &mut self.budget)
                .map_err(|e| e.at(d.span))?;
        }
        // Parameters outside their range: everywhere, the fault a number
        // makes; on part of the draw's range, not supported yet.
        let mut params = Vec::with_capacity(values.len());
        for (value, &(lo, hi)) in values.iter().zip(&ranges) {
            match value {
                Value::Analytic(a) => {
                    let bad = analytic::outside(&x, &a.pieces, lo, hi).map_err(|e| e.at(d.span))?;
                    let mass = bad.mass();
                    if mass > 0.0 {
                        let name = law.name();
                        return Err(if mass >= x.domain.mass() {
                            RuntimeError::new(d.span, format!("`{name}` isn't defined for these parameters"))
                                .as_fault(Fault::DomainError)
                        } else {
                            analytic::unsupported(&format!(
                                "`{name}` with a parameter outside its range on part of the draw's range"
                            ))
                            .at(d.span)
                        });
                    }
                    params.push(Param::Outcome(a.pieces.clone()));
                }
                v => params.push(Param::Plain(v.as_f64().expect("a number"))),
            }
        }
        if law.is_density() {
            self.check_density(f, span)?;
        }
        let mut factors = x.weight.as_ref().map(|l| l.factors.clone()).unwrap_or_default();
        factors.push(Factor { law, params, value });
        let before = match &x.weight {
            None => libm::log(x.domain.mass()),
            Some(l) => l.shift + libm::log(l.mass(&x.domain)),
        };
        let likelihood = Likelihood::new(x.family, factors, tolerance, &mut self.budget).map_err(|e| e.at(span))?;
        self.integrated = Some(tolerance);
        if law.is_density() {
            self.densities = true;
            w.densities += 1;
        }
        let Some(likelihood) = likelihood else {
            return Ok(Some(f64::NEG_INFINITY));
        };
        let after = likelihood.shift + libm::log(likelihood.mass(&x.domain));
        let latent = analytic::Latent {
            family: x.family,
            domain: x.domain.clone(),
            form: None,
            weight: Some(Arc::new(likelihood)),
        };
        self.know(w, x.id, latent, span)?;
        Ok(Some(after - before))
    }

    /// An observation's distribution that takes outcomes of draws as
    /// parameters, without `@integrate`: its error says that `@integrate`
    /// would compute it.
    #[cold]
    fn integration_hint(&self, d: &Expr, e: RuntimeError) -> RuntimeError {
        let law = matches!(
            &d.kind,
            ExprKind::Builtin {
                func: Builtin::Normal
                    | Builtin::Lognormal
                    | Builtin::Uniform
                    | Builtin::Beta
                    | Builtin::Gamma
                    | Builtin::Exponential
                    | Builtin::Binomial
                    | Builtin::Bernoulli
                    | Builtin::Poisson
                    | Builtin::Geometric,
                ..
            }
        );
        if e.kind != ErrorKind::Unsupported || !law || self.config.integrate.is_some() || self.sampler.is_some() {
            return e;
        }
        e.with_help(
            "`@integrate` computes such an observation of one draw by numerical integration, approximately; or use `@mode sample(runs: 10_000)`",
        )
    }

    /// A normal distribution's standard deviation, checked as `normal`
    /// checks it; `None` if it isn't a plain value.
    fn plain_sd(&mut self, f: FnId, sd: &Expr, w: &World, span: Span) -> Result<Option<f64>> {
        let sd = self.eval(f, sd, w)?;
        if analytic::contains(&sd) || sd.is_uncertain() {
            return Ok(None);
        }
        let checked = builtins::call_plain(Builtin::Normal, &[Value::Float(0.0), sd], &mut self.budget)
            .map_err(|e| e.at(span))?;
        let Value::Continuous(family) = checked else {
            unreachable!("`normal` gives a continuous distribution")
        };
        let Family::Normal { sd, .. } = *family else {
            unreachable!("`normal` gives a normal distribution")
        };
        Ok(Some(sd))
    }

    /// Drawing from `bernoulli(p)` or `binomial(n, p)` when enumerating, with
    /// `p` a variable drawn from a beta: each outcome, in a world where `p`
    /// is updated exactly by it. `None` for another draw.
    fn dependent_draw(&mut self, f: FnId, dist: &'p Expr, w: &World) -> Result<Option<Vec<(Value, World)>>> {
        let ExprKind::Builtin { func, args, named } = &dist.kind else {
            return Ok(None);
        };
        if self.sampler.is_some() || !named.is_empty() {
            return Ok(None);
        }
        let (p, trials) = match (func, args.as_slice()) {
            (Builtin::Bernoulli, [p]) => (p, None),
            (Builtin::Binomial, [n, p]) => (p, Some(n)),
            _ => return Ok(None),
        };
        // The variable, or `prob` of it.
        let slot = match &p.kind {
            ExprKind::Slot(s) => *s,
            ExprKind::Builtin {
                func: Builtin::Prob,
                args,
                ..
            } => match args.as_slice() {
                [
                    Expr {
                        kind: ExprKind::Slot(s),
                        ..
                    },
                ] => *s,
                _ => return Ok(None),
            },
            _ => return Ok(None),
        };
        let Some(x) = self.draw_in(w, slot, p.span)? else {
            return Ok(None);
        };
        let outcomes: Vec<(Value, Seen)> = match trials {
            None => [false, true].map(|b| (Value::Bool(b), Seen::Bernoulli(b))).into(),
            Some(n) => {
                let n = self.eval(f, n, w)?;
                if analytic::contains(&n) {
                    return Ok(None);
                }
                let counts = builtins::counts(Builtin::Binomial, &[n, Value::Prob(0.5)], &mut self.budget)
                    .map_err(|e| e.at(dist.span))?;
                let Some(Counts::Binomial { n, .. }) = counts else {
                    return Ok(None);
                };
                self.budget.outcomes(n as u128 + 1).map_err(|e| e.at(dist.span))?;
                (0..=n)
                    .map(|k| (Value::Int(k.into()), Seen::Binomial { trials: n, k: k as f64 }))
                    .collect()
            }
        };
        let seen: Vec<Seen> = outcomes.iter().map(|(_, s)| *s).collect();
        let worlds = self.outcomes_of(&x, w, &seen, dist.span)?;
        if worlds.is_empty() {
            return Ok(None);
        }
        Ok(Some(
            outcomes
                .into_iter()
                .zip(worlds)
                .filter_map(|((v, _), w)| Some((v, w?)))
                .collect(),
        ))
    }

    /// The worlds where a draw `x`, the parameter of a likelihood, saw each
    /// of `seen`, updated exactly: an empty list if `x` isn't the prior of
    /// that likelihood, and `None` for what `x` can't have seen.
    fn outcomes_of(&mut self, x: &Analytic, w: &World, seen: &[Seen], span: Span) -> Result<Vec<Option<World>>> {
        let mut out = Vec::with_capacity(seen.len());
        for &s in seen {
            let Some((ln, latent)) = analytic::update(x, s).map_err(|e| e.at(span))? else {
                return Ok(Vec::new());
            };
            if ln == f64::NEG_INFINITY {
                out.push(None);
                continue;
            }
            let mut w = w.clone();
            w.weight = w.weight * Weight::from_ln(ln);
            self.know(&mut w, x.id, latent, span)?;
            out.push(Some(w));
        }
        Ok(out)
    }

    /// The analytic outcome in `slot`, as the world reads it.
    fn draw_in(&mut self, w: &World, slot: SlotId, span: Span) -> Result<Option<Analytic>> {
        let v = &w.slots[slot as usize];
        if !matches!(v, Value::Analytic(_)) {
            return Ok(None);
        }
        Ok(
            match analytic::resolve(v, &w.constraints, &mut self.budget).map_err(|e| e.at(span))? {
                Value::Analytic(x) => Some((*x).clone()),
                _ => None,
            },
        )
    }

    /// Record what a world now knows about a latent, and read the values
    /// that hold it again, as they would be read.
    fn know(&mut self, w: &mut World, id: u64, latent: analytic::Latent, span: Span) -> Result<()> {
        Arc::make_mut(&mut w.constraints).insert(id, latent);
        self.reread(w, span)
    }

    /// Read the values that hold latents again, as the world reads them.
    fn reread(&mut self, w: &mut World, span: Span) -> Result<()> {
        for slot in &mut w.slots {
            if analytic::contains(slot) {
                *slot = analytic::resolve(slot, &w.constraints, &mut self.budget).map_err(|e| e.at(span))?;
            }
        }
        Ok(())
    }

    /// What a conjugate observation saw, checked as when its variable is
    /// drawn; `None` if a part of it isn't a plain value.
    fn seen(&mut self, f: FnId, u: &Update<'p>, w: &World, from: Span) -> Result<Option<Seen>> {
        if self.sampler.is_none() {
            for e in u.others() {
                if analytic::contains(&self.eval(f, e, w)?) {
                    return Ok(None);
                }
            }
        }
        // A count, as observing one from a drawn distribution reads it:
        // anything but a whole number of 0 or more is impossible.
        let count = |v: &Value| match v {
            Value::Bool(_) => f64::NAN,
            _ => v.as_f64().unwrap_or(f64::NAN),
        };
        Ok(match u.likelihood {
            Likelihood::Binomial { value, trials } => {
                let v = self.eval(f, value, w)?;
                let n = self.eval(f, trials, w)?;
                match builtins::counts(Builtin::Binomial, &[n, Value::Prob(0.5)], &mut self.budget)
                    .map_err(|e| e.at(from))?
                {
                    Some(Counts::Binomial { n, .. }) => Some(Seen::Binomial {
                        trials: n,
                        k: count(&v),
                    }),
                    _ => None,
                }
            }
            Likelihood::Bernoulli { value } => match self.eval(f, value, w)? {
                Value::Bool(b) => Some(Seen::Bernoulli(b)),
                _ => None,
            },
            Likelihood::Poisson { value } => {
                let v = self.eval(f, value, w)?;
                Some(Seen::Poisson { k: count(&v) })
            }
            Likelihood::Normal { value, sd } => {
                let v = self.eval(f, value, w)?;
                let sd = self.eval(f, sd, w)?;
                if sd.is_uncertain() {
                    return Ok(None);
                }
                let checked = builtins::call_plain(Builtin::Normal, &[Value::Float(0.0), sd], &mut self.budget)
                    .map_err(|e| e.at(from))?;
                let Value::Continuous(family) = checked else {
                    unreachable!("`normal` gives a continuous distribution")
                };
                let Family::Normal { sd, .. } = *family else {
                    unreachable!("`normal` gives a normal distribution")
                };
                match v {
                    Value::Int(_) | Value::Float(_) | Value::Prob(_) if v.as_f64().is_some() => Some(Seen::Normal {
                        y: v.as_f64().expect("a number"),
                        sd,
                    }),
                    _ => None,
                }
            }
        })
    }

    fn merge(&self, worlds: Vec<World>, stmt: StmtId) -> Vec<World> {
        merge(worlds, &self.live.after[stmt as usize], self.merging())
    }

    /// Whether equal worlds merge: never when sampling, since every run
    /// counts separately (section 14), but inside `simulate` they do.
    fn merging(&self) -> bool {
        self.config.merging && self.sampler.is_none()
    }

    /// Whether a fault ends only its world: in partial mode, except inside
    /// `simulate` and collection callbacks, which fail as a whole, in the
    /// world that runs them (docs/semantics.md, section 11).
    fn partial(&self) -> bool {
        self.config.partial && self.nested == 0 && self.callback.is_none()
    }

    /// Whether a `try` running in this call catches `fault`.
    fn catches_here(&self, fault: Option<Fault>) -> bool {
        self.handlers
            .iter()
            .any(|(depth, catches)| *depth == self.depth && catches.iter().any(|c| c.catches(fault)))
    }

    /// Whether a `try` running in a call further up catches `fault`.
    fn catches_above(&self, fault: Option<Fault>) -> bool {
        self.handlers
            .iter()
            .any(|(depth, catches)| *depth < self.depth && catches.iter().any(|c| c.catches(fault)))
    }

    /// Inside a `try` of this call: a copy of `w`, for the catch it may
    /// fault to. Elsewhere nothing is copied.
    #[inline]
    fn snapshot(&self, w: &World) -> Option<World> {
        if self.handlers.is_empty() {
            return None;
        }
        self.handlers
            .iter()
            .any(|(depth, _)| *depth == self.depth)
            .then(|| w.clone())
    }

    /// A world at `at` failed with `e`. A fault that a `try` of this call
    /// catches sends `saved`, the world as its statement began, on to the
    /// catch. One that a `try` further up catches, or any fault in partial
    /// mode, ends the world here: it's recorded, for the call's caller to
    /// catch, or as a failure. Either way the statement goes on with the
    /// other worlds. Anything else stops the run.
    fn fail(&mut self, e: RuntimeError, at: At, saved: Option<World>) -> Result<()> {
        if e.fault.is_none() {
            return Err(e);
        }
        if let Some(w) = saved.filter(|_| self.catches_here(e.fault)) {
            self.faulted.push((w, e));
            return Ok(());
        }
        if !self.partial() && !self.catches_above(e.fault) {
            return Err(e);
        }
        let fun = &self.prog.functions[at.f as usize];
        let e = match fun.kind {
            FnKind::Named => e.with_note(format!("in a call to `{}`", fun.name)),
            FnKind::Lambda => e.with_note("inside a lambda"),
            FnKind::Simulate => e.with_note("inside a `simulate` block"),
            FnKind::Main => e,
        };
        let stmt = at.stmt as usize;
        let before_evidence = self.evidence.at[stmt] || self.evidence.after[stmt];
        let run = self.sampler.is_some().then_some(at.run);
        self.failures.record(e, at.weight, run, before_evidence, at.densities);
        Ok(())
    }

    /// Stop if a statement produced more worlds than allowed. (When sampling,
    /// worlds never multiply: there's one per run in the batch.)
    fn check_worlds(&self, n: usize, span: Span) -> Result<()> {
        if n > self.config.max_worlds && self.sampler.is_none() {
            return Err(RuntimeError::limit(
                span,
                format!("more than {} worlds are too many to follow", self.config.max_worlds),
            )
            .with_help("simplify the model, raise the limit, or sample it with `@mode sample(runs: 10_000)`"));
        }
        Ok(())
    }

    fn spend(&mut self, n: u64, span: Span) -> Result<()> {
        self.budget.work(n).map_err(|e| e.at(span))?;
        if let Some(cancel) = &self.config.cancel {
            if cancel.load(Ordering::Relaxed) {
                return Err(RuntimeError::limit(span, "the run was cancelled"));
            }
        }
        Ok(())
    }

    // ── Statements ───────────────────────────────────────────────────────

    fn exec_block(&mut self, f: FnId, block: &'p Block, worlds: Vec<World>) -> Result<Flow> {
        let mut flow = Flow::next(worlds);
        for stmt in &block.stmts {
            if flow.next.is_empty() {
                break;
            }
            let input = std::mem::take(&mut flow.next);
            let out = self.exec_stmt(f, stmt, input)?;
            flow.next = out.next;
            clear(&mut flow.next, &self.live.dies[stmt.id as usize]);
            flow.broke.extend(out.broke);
            flow.continued.extend(out.continued);
            flow.returned.extend(out.returned);
            flow.faulted.extend(out.faulted);
        }
        Ok(flow)
    }

    /// Weighs an enumerated world by the density `exp(ln)` of what an
    /// observation at `span` saw.
    fn observe_density(&mut self, f: FnId, ln: f64, w: &mut World, span: Span) -> Result<()> {
        self.check_density(f, span)?;
        if ln == f64::INFINITY {
            return Err(RuntimeError::new(span, "the density at the observed value is infinite"));
        }
        if ln == f64::NEG_INFINITY {
            self.lost += w.weight;
            w.weight = Weight::ZERO;
            self.last_ruling_out = Some(span);
        } else {
            w.weight = w.weight * Weight::from_ln(ln);
            w.densities += 1;
        }
        Ok(())
    }

    /// An enumerated world's weight can include a density only where the
    /// world counts it (see `World::densities`).
    /// A `simulate` block's own observations can be densities too: its
    /// result is normalized, as long as its worlds include as many.
    #[inline(never)]
    fn check_density(&mut self, f: FnId, span: Span) -> Result<()> {
        if self.depth == 0 {
            self.density_at.get_or_insert(span);
            return Ok(());
        }
        if self.prog.functions[f as usize].kind == FnKind::Simulate {
            self.local_density_at.get_or_insert(span);
            return Ok(());
        }
        Err(OpError::unsupported(
            "observing a value from a continuous distribution inside a function isn't supported when enumerating yet",
        )
        .help("observe it in the main program or directly inside `simulate`, or sample the model with `@mode sample(runs: 10_000)`")
        .at(span))
    }

    /// A world returning from a `simulate` block: its weight must include as
    /// many densities as the others'.
    #[inline(never)]
    fn simulate_returns(&mut self, w: &World, span: Span) -> Result<()> {
        match self.returned_densities {
            Some(n) if n != w.densities => Err(OpError::unsupported(
                "a `simulate` whose worlds observe different numbers of continuous values isn't supported when enumerating",
            )
            .help("a density is per unit of what it observes, so worlds with more of them can't be weighed against the others; observe the same values in every world")
            .at(self.local_density_at.unwrap_or(span))),
            _ => {
                self.returned_densities = Some(w.densities);
                Ok(())
            }
        }
    }

    /// The parts of a world where an event holds, or doesn't: each with its
    /// probability given what the world knew, and the restrictions there.
    /// An event of several independent draws has a part for each box.
    fn event_parts(
        &mut self,
        event: &analytic::Event,
        yes: bool,
        constraints: &analytic::Constraints,
        span: Span,
    ) -> Result<Vec<(f64, analytic::Constraints)>> {
        if !event.restrictable() {
            return Err(analytic::joint_condition("a condition on several normal draws together").at(span));
        }
        let boxes = event.side(yes).map_err(|e| e.at(span))?;
        let n = boxes.len().max(1);
        self.budget
            .collection(((constraints.len() + event.draws.len()) * n) as u128)
            .map_err(|e| e.at(span))?;
        self.budget
            .work(((constraints.len() + event.size()) * n) as u64)
            .map_err(|e| e.at(span))?;
        let mut parts = Vec::with_capacity(boxes.len());
        for part in &boxes {
            let mut c = constraints.clone();
            let p = event.restrict(part, &mut c);
            parts.push((p, c));
        }
        Ok(parts)
    }

    /// `if` on an event: each part of the world where it holds goes to
    /// `yes`, and each where it doesn't to `no`, restricted to that part.
    #[inline(never)]
    fn branch_on_event(
        &mut self,
        event: &analytic::Event,
        w: World,
        yes: &mut Vec<World>,
        no: &mut Vec<World>,
        span: Span,
    ) -> Result<()> {
        let holds = self.event_parts(event, true, &w.constraints, span)?;
        let fails = self.event_parts(event, false, &w.constraints, span)?;
        let parts: Vec<_> = (holds.into_iter().map(|(p, c)| (true, p, c)))
            .chain(fails.into_iter().map(|(p, c)| (false, p, c)))
            .filter(|(_, p, _)| *p > 0.0)
            .collect();
        let last = parts.len().saturating_sub(1);
        let mut w = Some(w);
        for (i, (side, p, c)) in parts.into_iter().enumerate() {
            // The last part takes the world itself.
            let mut nw = if i == last { w.take() } else { w.clone() }.expect("the world");
            nw.constraints = c;
            let out = if side { &mut *yes } else { &mut *no };
            out.push(nw.scaled(p));
        }
        Ok(())
    }

    /// `observe` of an event, or of one side of it: each part of the world
    /// where it holds goes on, restricted to that part.
    #[inline(never)]
    fn observe_event(
        &mut self,
        event: &analytic::Event,
        yes: bool,
        w: World,
        out: &mut Vec<World>,
        span: Span,
    ) -> Result<()> {
        let parts = self.event_parts(event, yes, &w.constraints, span)?;
        let p: f64 = parts.iter().map(|(p, _)| p).sum();
        if p < 1.0 && self.sampler.is_none() {
            self.lost += w.weight.scale(1.0 - p);
        }
        if p <= 0.0 {
            self.last_ruling_out = Some(span);
        }
        let parts: Vec<_> = parts.into_iter().filter(|(p, _)| *p > 0.0).collect();
        let last = parts.len().saturating_sub(1);
        let mut w = Some(w);
        for (i, (p, c)) in parts.into_iter().enumerate() {
            let mut nw = if i == last { w.take() } else { w.clone() }.expect("the world");
            nw.constraints = c;
            out.push(nw.scaled(p));
        }
        Ok(())
    }

    fn exec_stmt(&mut self, f: FnId, stmt: &'p Stmt, worlds: Vec<World>) -> Result<Flow> {
        // Outside every `try`, no world can fault to a catch.
        if self.handlers.is_empty() {
            return self.exec_stmt_kind(f, stmt, worlds);
        }
        // The worlds that fault in this statement, on their way to a catch,
        // leave with its flow.
        let outer = std::mem::take(&mut self.faulted);
        let flow = self.exec_stmt_kind(f, stmt, worlds);
        let faulted = std::mem::replace(&mut self.faulted, outer);
        let mut flow = flow?;
        flow.faulted.extend(faulted);
        Ok(flow)
    }

    fn exec_stmt_kind(&mut self, f: FnId, stmt: &'p Stmt, worlds: Vec<World>) -> Result<Flow> {
        let span = stmt.span;
        if matches!(
            stmt.kind,
            StmtKind::Draw { .. } | StmtKind::Take { .. } | StmtKind::Observe { .. } | StmtKind::Chance { .. }
        ) {
            self.check_callback_effect(span)?;
        }
        let n = worlds.len();
        self.stats.world_steps += n as u64;
        self.stats.peak_worlds = self.stats.peak_worlds.max(n);
        self.spend(n as u64, span)?;
        self.check_worlds(n, span)?;
        let mut worlds = worlds;
        if self.delaying() {
            // Delayed variables this statement reads are drawn first.
            let conj = self.conj;
            let first = &conj.draws_first[stmt.id as usize];
            if !first.is_empty() {
                for w in &mut worlds {
                    for &slot in first {
                        self.draw_delayed(w, slot);
                    }
                }
            }
        }
        match &stmt.kind {
            StmtKind::Set { place, value } | StmtKind::Update { place, value, .. } => {
                let update = match &stmt.kind {
                    StmtKind::Update { check, .. } => Some(check.as_ref()),
                    _ => None,
                };
                let mut worlds = worlds;
                let mut failed = Vec::new();
                let mut discrete = false;
                for (i, w) in worlds.iter_mut().enumerate() {
                    let mut saved = self.snapshot(w);
                    if let Some(check) = update {
                        let updated = self.update(f, place, value, check, w, span);
                        each!(self, At::new(f, stmt, w), saved saved, updated, { failed.push(i) });
                        continue;
                    }
                    let v = match self.record_place_type(place, w) {
                        Some(ty) => self.eval_expected(f, value, w, &ty),
                        None => self.eval(f, value, w),
                    };
                    let v = each!(self, At::new(f, stmt, w), saved saved, v, { failed.push(i) });
                    if let Value::Analytic(_) = v {
                        discrete |= is_discrete(&v);
                    }
                    each!(self, At::new(f, stmt, w), saved saved, self.assign(f, place, v, w, span), {
                        failed.push(i)
                    });
                }
                if discrete {
                    self.split_discrete(f, stmt, place, &mut worlds, &mut failed)?;
                }
                drop_failed(&mut worlds, &failed);
                // When the live values say less than before, as when a
                // call's result is used up (`f = f - $5`), worlds may now be
                // the same: merging spares the statements after from
                // following each of them (a call next would multiply them).
                if self.live.narrows[stmt.id as usize] && self.merging() {
                    return Ok(Flow::next(self.merge(worlds, stmt.id)));
                }
                Ok(Flow::next(worlds))
            }
            StmtKind::Draw { place, dist } => {
                // A conjugate prior's draw may be delayed (section 14).
                let delay = if self.delaying() {
                    self.conj.delays[stmt.id as usize]
                } else {
                    None
                };
                let mut out = Vec::with_capacity(worlds.len());
                for w in worlds {
                    let at = At::new(f, stmt, &w);
                    let mut saved = self.snapshot(&w);
                    if self.sampler.is_some() {
                        if let Some(counts) = each!(self, at, saved saved, self.direct_counts(f, dist, &w)) {
                            let k = counts.sample(self.sampler.as_mut().expect("sampling"));
                            let mut w = w;
                            each!(self, at, saved saved, self.assign(f, place, Value::Int(k.into()), &mut w, dist.span));
                            out.push(w);
                            continue;
                        }
                    }
                    if self.sampler.is_none()
                        && each!(self, at, saved saved, self.draw_dependent(f, place, dist, &w, &mut out))
                    {
                        self.check_worlds(out.len(), span)?;
                        continue;
                    }
                    let d = each!(self, at, saved saved, self.eval(f, dist, &w));
                    if let (Some(variable), Value::Continuous(family)) = (delay, &d) {
                        if conjugate::is_prior(family) {
                            let mut w = w;
                            w.slots[place.slot as usize] = Value::Delayed(Arc::new(Delayed {
                                family: **family,
                                variable,
                            }));
                            self.updates(variable).delayed += 1;
                            out.push(w);
                            continue;
                        }
                    }
                    let mark = out.len();
                    each!(self, at, saved saved, self.split_by(f, place, d, w, dist.span, &mut out), {
                        out.truncate(mark)
                    });
                    self.check_worlds(out.len(), span)?;
                }
                Ok(Flow::next(self.merge(out, stmt.id)))
            }
            StmtKind::Take { place, bag } => {
                let mut out = Vec::new();
                for w in worlds {
                    let at = At::new(f, stmt, &w);
                    let current = each!(self, at, world w, self.read_place(f, bag, &w, span));
                    let Value::Bag(cards) = &current else {
                        return Err(RuntimeError::new(
                            span,
                            format!("`take` needs a bag, found {}", ops::article(&current.kind())),
                        )
                        .with_help("make one with `bag([card: count, …])`"));
                    };
                    let total: u128 = cards.values().map(|n| *n as u128).sum();
                    if total == 0 {
                        let empty = RuntimeError::new(span, "can't take a card from an empty bag")
                            .as_fault(Fault::EmptyCollection);
                        each!(self, at, world w, Err::<(), _>(empty));
                    }
                    let chosen = match &mut self.sampler {
                        Some(rng) => rng.choose(cards.values().map(|n| *n as f64)),
                        None => None,
                    };
                    for (i, (card, count)) in cards.iter().enumerate() {
                        if self.sampler.is_some() && chosen != Some(i) {
                            continue;
                        }
                        let rest = cards.without_nth(i);
                        let mut nw = if self.sampler.is_some() {
                            w.clone()
                        } else {
                            w.clone().scaled(*count as f64 / total as f64)
                        };
                        self.assign(f, bag, Value::multiset(rest), &mut nw, span)?;
                        self.assign(f, place, card.clone(), &mut nw, span)?;
                        out.push(nw);
                    }
                    self.check_worlds(out.len(), span)?;
                }
                Ok(Flow::next(self.merge(out, stmt.id)))
            }
            StmtKind::Call { dest, callee, args } => {
                let mut out = Vec::with_capacity(worlds.len());
                for w in worlds {
                    let at = At::new(f, stmt, &w);
                    let mut key = Vec::with_capacity(args.len() + 4);
                    let evaluated: Result<()> = args.iter().try_for_each(|a| {
                        key.push(self.eval(f, a, &w)?);
                        Ok(())
                    });
                    each!(self, at, world w, evaluated);
                    let func = match callee {
                        Callee::Fn { func, capture_args } => {
                            self.capture(&mut key, f, *func, capture_args, &w, span, false)?;
                            *func
                        }
                        Callee::Value(e, named) => match each!(self, at, world w, self.eval(f, e, &w)) {
                            Value::Builtin(b) => {
                                let value =
                                    each!(self, at, world w, self.builtin_values(b, &key, named, w.weight, span));
                                let mut saved = self.snapshot(&w);
                                let mut nw = w;
                                each!(self, at, saved saved, self.assign(f, dest, value, &mut nw, span));
                                out.push(nw);
                                continue;
                            }
                            Value::Closure(c) => {
                                if !named.is_empty() {
                                    return Err(RuntimeError::new(
                                        span,
                                        "only builtin minimum/maximum accept named defaults",
                                    ));
                                }
                                self.check_arity(&c, args.len(), span)?;
                                key.extend(c.captured.iter().cloned());
                                c.func
                            }
                            other => {
                                return Err(RuntimeError::new(
                                    e.span,
                                    format!("can't call {}", ops::article(&other.kind())),
                                ));
                            }
                        },
                    };
                    let result = each!(self, at, world w, self.call(func, key, span));
                    self.unresolved += w.weight * result.unresolved;
                    self.lost += w.weight * result.lost;
                    // What faulted in the call faulted in this world: a catch
                    // here takes it with the world as it was, and otherwise
                    // the world fails, for a `try` further up or as a
                    // failure, and may still meet evidence after the call.
                    let run = self.sampler.is_some().then_some(w.run);
                    let after = self.evidence.after[stmt.id as usize];
                    for g in &result.failures.groups {
                        if self.catches_here(g.error.fault) {
                            let mut caught = w.clone();
                            caught.weight = caught.weight * g.weight;
                            self.faulted.push((caught, g.error.clone()));
                        } else if self.partial() || self.catches_above(g.error.fault) {
                            self.failures.absorb_group(g, w.weight, run, after, w.densities);
                        } else {
                            return Err(g.error.clone());
                        }
                    }
                    add_pending(&mut self.pending, &result.pending, w.weight);
                    let last = result.outcomes.len().saturating_sub(1);
                    let mut w = Some(w);
                    for (i, (v, p, restrictions)) in result.outcomes.iter().enumerate() {
                        let mut nw = if i == last {
                            w.take().unwrap()
                        } else {
                            w.clone().unwrap()
                        };
                        nw.weight = nw.weight * *p;
                        // Only copy the constraints (shared with the other
                        // outcomes' worlds) when the call added some.
                        if !restrictions.is_empty() {
                            Arc::make_mut(&mut nw.constraints)
                                .extend(restrictions.iter().map(|(id, d)| (*id, d.clone())));
                        }
                        self.assign(f, dest, v.clone(), &mut nw, span)?;
                        out.push(nw);
                    }
                    self.check_worlds(out.len(), span)?;
                }
                Ok(Flow::next(self.merge(out, stmt.id)))
            }
            StmtKind::If { cond, then, otherwise } => {
                let (mut yes, mut no) = (Vec::new(), Vec::new());
                for w in worlds {
                    let condition = each!(self, At::new(f, stmt, &w), world w, self.eval(f, cond, &w));
                    if let Value::Event(event) = condition {
                        self.check_callback_effect(cond.span)?;
                        self.branch_on_event(&event, w, &mut yes, &mut no, cond.span)?;
                        continue;
                    }
                    if let Value::Analytic(x) = &condition {
                        if self.branch_on_draw(x, &w, &mut yes, &mut no, cond.span)? {
                            continue;
                        }
                    }
                    let c = ops::condition(&condition).map_err(|e| e.at(cond.span))?;
                    if (c.yes > 0.0 && c.no > 0.0) || c.missing > 0.0 {
                        self.check_callback_effect(cond.span)?;
                    }
                    if let Some(rng) = &mut self.sampler {
                        // One branch, chosen with its probability.
                        match rng.choose([c.yes, c.no].into_iter()) {
                            Some(0) => yes.push(w),
                            Some(_) => no.push(w),
                            None => self.unresolved += w.weight,
                        }
                        continue;
                    }
                    let (p, q) = (c.yes, c.no);
                    if c.missing > 0.0 {
                        self.unresolved += w.weight.scale(c.missing);
                    }
                    if q <= 0.0 {
                        yes.push(w.scaled(p));
                    } else if p <= 0.0 {
                        no.push(w.scaled(q));
                    } else {
                        no.push(w.clone().scaled(q));
                        yes.push(w.scaled(p));
                    }
                }
                self.check_worlds(yes.len() + no.len(), span)?;
                let mut flow = Flow::default();
                if !yes.is_empty() {
                    flow.join(self.exec_block(f, then, yes)?);
                }
                if !no.is_empty() {
                    flow.join(self.exec_block(f, otherwise, no)?);
                }
                flow.next = self.merge(flow.next, stmt.id);
                Ok(flow)
            }
            StmtKind::Chance {
                arms,
                otherwise,
                exhaustive,
            } => {
                let mut buckets: Vec<Vec<World>> = vec![Vec::new(); arms.len()];
                let mut rest = Vec::new();
                for w in worlds {
                    let at = At::new(f, stmt, &w);
                    // Chance weights are explicit probabilities.
                    let chances: Result<Vec<f64>> = arms
                        .iter()
                        .map(|(weight, _)| {
                            let value = self.eval(f, weight, &w)?;
                            ops::to_prob(&value).map_err(|e| e.at(weight.span))
                        })
                        .collect();
                    let chances = each!(self, at, world w, chances);
                    let sum: f64 = chances.iter().sum();
                    if sum > 1.0 + 1e-9 {
                        let over = RuntimeError::new(
                            span,
                            format!(
                                "the chances in this `chance` add up to {}, more than 100%",
                                fmt_prob(sum)
                            ),
                        );
                        each!(self, at, world w, Err::<(), _>(over.as_fault(Fault::DomainError)));
                    }
                    let remainder = (1.0 - sum).max(0.0);
                    if remainder > 1e-9 && *exhaustive && otherwise.is_none() {
                        let short = RuntimeError::new(
                            span,
                            format!("the chances add up to {} and there's no `else` arm", fmt_prob(sum)),
                        )
                        .with_help("a `chance` used as a value needs its chances to add up to 100%, or an `else`");
                        each!(self, at, world w, Err::<(), _>(short.as_fault(Fault::DomainError)));
                    }
                    if let Some(rng) = &mut self.sampler {
                        // One arm, chosen with its probability.
                        let options = chances.iter().copied().chain(std::iter::once(remainder));
                        match rng.choose(options) {
                            Some(i) if i < arms.len() => buckets[i].push(w),
                            Some(_) => rest.push(w),
                            None => self.unresolved += w.weight,
                        }
                        continue;
                    }
                    for (bucket, p) in buckets.iter_mut().zip(&chances) {
                        if *p > 0.0 {
                            bucket.push(w.clone().scaled(*p));
                        }
                    }
                    if remainder > 1e-12 {
                        rest.push(w.scaled(remainder));
                    }
                }
                self.check_worlds(buckets.iter().map(Vec::len).sum::<usize>() + rest.len(), span)?;
                let mut flow = Flow::default();
                for ((_, body), bucket) in arms.iter().zip(buckets) {
                    if !bucket.is_empty() {
                        flow.join(self.exec_block(f, body, bucket)?);
                    }
                }
                if !rest.is_empty() {
                    match otherwise {
                        Some(body) => flow.join(self.exec_block(f, body, rest)?),
                        None if *exhaustive => {}
                        None => flow.next.extend(rest),
                    }
                }
                flow.next = self.merge(flow.next, stmt.id);
                Ok(flow)
            }
            StmtKind::Loop { body, bounded } => self.exec_loop(f, stmt, body, *bounded, worlds),
            StmtKind::Break => Ok(Flow {
                broke: worlds,
                ..Flow::default()
            }),
            StmtKind::Continue => Ok(Flow {
                continued: worlds,
                ..Flow::default()
            }),
            StmtKind::Return(value) => self.exec_return(f, stmt, value, worlds),
            StmtKind::Observe { value, from } => {
                self.observed = true;
                // When sampling, an update of a delayed draw; when
                // enumerating, of an analytic one.
                let exact = if self.delaying() || self.sampler.is_none() {
                    probl_sema::conjugate::update(value, from.as_ref())
                } else {
                    None
                };
                let mut out = Vec::with_capacity(worlds.len());
                for mut w in worlds {
                    let here = At::new(f, stmt, &w);
                    let mut saved = self.snapshot(&w);
                    if let Some(u) = &exact {
                        let at = from.as_ref().map_or(span, |d| d.span);
                        let ln = if self.sampler.is_some() {
                            self.observe_exactly(f, u, &mut w, at)
                        } else {
                            self.observe_posterior(f, u, &mut w, at)
                        };
                        if let Some(ln) = each!(self, here, saved saved, ln) {
                            w.weight = w.weight * Weight::from_ln(ln);
                            if w.weight.is_zero() {
                                self.last_ruling_out = Some(span);
                            } else {
                                out.push(w);
                            }
                            continue;
                        }
                    }
                    let (factor, missing, ruled_out) = match from {
                        None => {
                            let v = each!(self, here, saved saved, self.eval(f, value, &w));
                            if let Value::Event(event) = v {
                                self.observe_event(&event, true, w, &mut out, value.span)?;
                                continue;
                            } else {
                                let b = ops::fact(&v, "observe").map_err(|e| e.at(value.span))?;
                                if b { (1.0, 0.0, 0.0) } else { (0.0, 0.0, 1.0) }
                            }
                        }
                        Some(d) => {
                            let v = each!(self, here, saved saved, self.eval(f, value, &w));
                            match each!(self, here, saved saved, self.direct_counts(f, d, &w)) {
                                Some(counts) if self.sampler.is_some() => {
                                    let x = match v {
                                        Value::Bool(_) => None,
                                        _ => v.as_f64(),
                                    };
                                    (x.map_or(0.0, |x| counts.pmf(x)), 0.0, 0.0)
                                }
                                _ => {
                                    if self.sampler.is_none() {
                                        let mut exact = self.observe_gaussian(f, &v, d, &mut w, span);
                                        if matches!(exact, Ok(None)) {
                                            exact = self.observe_numerically(f, &v, d, &mut w, span);
                                        }
                                        if let Some(ln) = each!(self, here, saved saved, exact) {
                                            w.weight = w.weight * Weight::from_ln(ln);
                                            if w.weight.is_zero() {
                                                self.last_ruling_out = Some(span);
                                            } else {
                                                out.push(w);
                                            }
                                            continue;
                                        }
                                    }
                                    let dist = self.eval(f, d, &w).map_err(|e| self.integration_hint(d, e));
                                    let dist = each!(self, here, saved saved, dist);
                                    if let (Value::Bool(b), Value::Event(event)) = (&v, &dist) {
                                        self.observe_event(event, *b, w, &mut out, span)?;
                                        continue;
                                    } else {
                                        if analytic::contains(&v) || analytic::contains(&dist) {
                                            return Err(analytic::unsupported("this likelihood observation").at(span));
                                        }
                                        self.densities |= is_density(&dist);
                                        match each!(
                                            self,
                                            here,
                                            saved saved,
                                            likelihood(&dist, &v, self.sampler.is_some()).map_err(|e| e.at(span))
                                        ) {
                                            Weighs::Probability(p, missing, other) => (p, missing, other),
                                            Weighs::LogDensity(ln) => {
                                                self.observe_density(f, ln, &mut w, span)?;
                                                if !w.weight.is_zero() {
                                                    out.push(w);
                                                }
                                                continue;
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    };
                    if missing > 0.0 && self.sampler.is_none() {
                        self.unresolved += w.weight.scale(missing);
                    }
                    if ruled_out > 0.0 && self.sampler.is_none() {
                        self.lost += w.weight.scale(ruled_out);
                    }
                    if factor > 0.0 {
                        w.weight = w.weight.scale(factor);
                        out.push(w);
                    } else {
                        self.last_ruling_out = Some(span);
                    }
                }
                Ok(Flow::next(out))
            }
            StmtKind::Report { site, value, key } => {
                let mut failed = Vec::new();
                for (i, w) in worlds.iter().enumerate() {
                    let v = each!(self, At::new(f, stmt, w), world w, self.eval(f, value, w), { failed.push(i) });
                    // A report of a count law lists it.
                    let v = match v {
                        Value::Counts(_) => ops::listed(&v, &mut self.budget)
                            .map_err(|e| e.at(value.span))?
                            .into_owned(),
                        v => v,
                    };
                    let k = match key {
                        Some(k) => {
                            let k_value =
                                each!(self, At::new(f, stmt, w), world w, self.eval(f, k, w), { failed.push(i) });
                            if analytic::contains(&k_value) && !matches!(k_value, Value::Event(_)) {
                                return Err(analytic::unsupported("grouping by a continuous outcome").at(k.span));
                            }
                            if k_value.is_uncertain() {
                                return Err(RuntimeError::new(
                                    k.span,
                                    "a `by` key must be a plain value, not a distribution",
                                )
                                .with_help("draw a value first with `~`"));
                            }
                            k_value
                        }
                        None => Value::Unit,
                    };
                    // A distribution of outcomes, from `max(x, d6)` say, reports
                    // as their mixture, and a record or list of them field by
                    // field.
                    let (v, run) = if self.sampler.is_some() {
                        (self.sample_continuous(v), Some(w.run))
                    } else {
                        (v, None)
                    };
                    if let Value::Event(event) = &k {
                        // An event of a draw is a key with two values: each
                        // group has its part of the world.
                        let span = key.as_ref().map_or(span, |k| k.span);
                        for yes in [true, false] {
                            for (p, constraints) in self.event_parts(event, yes, &w.constraints, span)? {
                                if p <= 0.0 {
                                    continue;
                                }
                                let v = analytic::resolve(&v, &constraints, &mut self.budget)
                                    .map_err(|e| e.at(value.span))?;
                                let k = Value::Bool(yes);
                                self.same_units(w, value.span)?;
                                let sink = &self.sinks[*site as usize];
                                sink.validate_analytic(&k, &v).map_err(|e| e.at(value.span))?;
                                sink.validate_structure(&k, &v).map_err(|e| e.at(value.span))?;
                                self.sinks[*site as usize].add(k, &v, w.weight.scale(p), run);
                            }
                        }
                        continue;
                    }
                    if run.is_none() {
                        self.same_units(w, value.span)?;
                        self.sinks[*site as usize]
                            .validate_structure(&k, &v)
                            .map_err(|e| e.at(value.span))?;
                    }
                    self.sinks[*site as usize]
                        .validate_analytic(&k, &v)
                        .map_err(|e| e.at(value.span))?;
                    self.sinks[*site as usize].add(k, &v, w.weight, run);
                }
                let mut worlds = worlds;
                drop_failed(&mut worlds, &failed);
                Ok(Flow::next(worlds))
            }
            StmtKind::Fail { message } => Err(RuntimeError::new(span, message.clone())),
            StmtKind::Try { body, catches } => {
                self.handlers.push((self.depth, catches.as_slice()));
                let flow = self.exec_block(f, body, worlds);
                self.handlers.pop();
                let mut flow = flow?;
                // Each faulted world goes on in the first catch that takes
                // its fault; the others are for a `try` around this one.
                let mut caught: Vec<Vec<World>> = vec![Vec::new(); catches.len()];
                for (w, e) in std::mem::take(&mut flow.faulted) {
                    match catches.iter().position(|c| c.catches(e.fault)) {
                        Some(i) => caught[i].push(w),
                        None => flow.faulted.push((w, e)),
                    }
                }
                for (c, worlds) in catches.iter().zip(caught) {
                    if !worlds.is_empty() {
                        flow.join(self.exec_block(f, &c.body, worlds)?);
                    }
                }
                flow.next = self.merge(flow.next, stmt.id);
                Ok(flow)
            }
            StmtKind::Check { slot, ty } => {
                for w in &mut worlds {
                    // A delayed variable is drawn only if some of its values
                    // could fail the check.
                    if let Value::Delayed(_) = &w.slots[*slot as usize] {
                        let always = *ty == TypeSpec::Float;
                        if always {
                            continue;
                        }
                        self.draw_delayed(w, *slot);
                    }
                    self.check_slot(f, *slot, ty, w, span)?;
                }
                Ok(Flow::next(worlds))
            }
        }
    }

    /// `return`: each world's value, with the restrictions its caller can
    /// see. (Out of `exec_stmt_kind`, whose frame every statement pays for.)
    #[inline(never)]
    fn exec_return(&mut self, f: FnId, stmt: &'p Stmt, value: &'p Expr, worlds: Vec<World>) -> Result<Flow> {
        let mut flow = Flow::default();
        let simulate = self.prog.functions[f as usize].kind == FnKind::Simulate;
        for w in worlds {
            let v = each!(self, At::new(f, stmt, &w), world w, self.eval(f, value, &w));
            if simulate {
                self.simulate_returns(&w, stmt.span)?;
            }
            let mut constraints = w.constraints;
            if constraints.keys().any(|id| !w.inherited.contains(id)) {
                Arc::make_mut(&mut constraints).retain(|id, _| w.inherited.contains(id));
            }
            flow.returned.push((v, w.weight, constraints));
        }
        Ok(flow)
    }

    fn exec_loop(
        &mut self,
        f: FnId,
        stmt: &'p Stmt,
        body: &'p Block,
        bounded: bool,
        worlds: Vec<World>,
    ) -> Result<Flow> {
        let live = self.live;
        let span = stmt.span;
        let entered = total_weight(&worlds);
        let cutoff = entered.scale(self.config.epsilon);
        let mut inside = worlds;
        let mut out = Flow::default();
        let mut iterations: u64 = 0;
        // The states met at the loop's start in earlier iterations, while it
        // may be solved. Once one comes back, the loop cycles, and it's
        // solved as a Markov chain from there (section 10).
        let mut met = self.may_solve(stmt, bounded).then(FxHashSet::<u64>::default);
        while !inside.is_empty() {
            let cycles = match &mut met {
                Some(met) => {
                    // A loop's rewards don't make its states new.
                    let rewards = self.rewards.get(&stmt.id).map_or(&[][..], Vec::as_slice);
                    let mut head = live_slots(&live.loop_head[stmt.id as usize], inside[0].slots.len());
                    head.retain(|&i| !rewards.contains(&(i as SlotId)));
                    let hashes: Vec<u64> = inside.iter().map(|w| state_hash(w, &head)).collect();
                    let again = hashes.iter().any(|h| met.contains(h));
                    met.extend(hashes);
                    again
                }
                None => false,
            };
            if cycles {
                met = None;
                if let Some(solved) = self.solve_loop(f, stmt, body, &inside, cutoff, &mut iterations)? {
                    out.join(solved);
                    break;
                }
            }
            let mass = total_weight(&inside);
            if !bounded && mass < cutoff && self.sampler.is_none() {
                self.unresolved += mass;
                break;
            }
            iterations += 1;
            if iterations > self.config.max_iterations {
                return Err(RuntimeError::limit(
                    span,
                    format!("this loop ran {} times without finishing", self.config.max_iterations),
                )
                .with_note(format!(
                    "worlds still inside the loop weigh {} of what entered it",
                    fmt_prob(mass.ratio(entered))
                ))
                .with_help("check that every world can leave the loop"));
            }
            let mut flow = self.exec_block(f, body, inside)?;
            clear_dead(&mut flow.broke, &live.after[stmt.id as usize]);
            clear_dead(&mut flow.continued, &live.loop_head[stmt.id as usize]);
            out.next.extend(flow.broke);
            out.returned.extend(flow.returned);
            out.faulted.extend(flow.faulted);
            let mut again = flow.next;
            again.extend(flow.continued);
            inside = merge(again, &live.loop_head[stmt.id as usize], self.merging());
        }
        out.next = self.merge(out.next, stmt.id);
        Ok(out)
    }

    /// Whether a loop may be solved as a Markov chain: an unbounded loop
    /// whose body doesn't report or print, when enumerating, unless its
    /// chain was already found too large.
    fn may_solve(&self, stmt: &Stmt, bounded: bool) -> bool {
        !bounded
            && self.config.solving
            && self.sampler.is_none()
            && self.solvable[stmt.id as usize]
            && !self.too_large.contains(&stmt.id)
    }

    /// Solve a loop as an absorbing Markov chain, from the worlds `inside`
    /// at its start (section 10). Its states are the worlds at its start,
    /// told apart by the variables read later. Running the body once from
    /// each state reachable from `inside` gives the chance of going on to
    /// each state, and the worlds that leave. The expected number of visits
    /// to each state then says how much of each leaves. `None` if the chain
    /// is too large, and the loop should be unrolled instead.
    fn solve_loop(
        &mut self,
        f: FnId,
        stmt: &'p Stmt,
        body: &'p Block,
        inside: &[World],
        cutoff: Weight,
        iterations: &mut u64,
    ) -> Result<Option<Flow>> {
        // The chain's transitions are probabilities: a density isn't one.
        if inside
            .iter()
            .any(|w| w.densities > 0 || !w.constraints.is_empty() || w.slots.iter().any(analytic::contains))
        {
            self.too_large.insert(stmt.id);
            return Ok(None);
        }
        let head_set = &self.live.loop_head[stmt.id as usize];
        let after = &self.live.after[stmt.id as usize];
        let n = inside[0].slots.len();
        let rewards: Vec<usize> = self
            .rewards
            .get(&stmt.id)
            .map_or(&[][..], Vec::as_slice)
            .iter()
            .map(|&r| r as usize)
            .collect();
        let mut head = live_slots(head_set, n);
        head.retain(|i| !rewards.contains(i));
        let dead: Vec<SlotId> = head_set.iter_missing(n).collect();
        if rewards.len() > MAX_REWARDS {
            self.too_large.insert(stmt.id);
            return Ok(None);
        }
        // A world's rewards, where they're plain numbers: whole numbers or
        // floats, by their bits, and which are floats.
        let reward_of = |slots: &[Value]| -> Option<Raw> {
            let (mut bits, mut floats) = ([0u64; MAX_REWARDS], 0u8);
            for (i, &r) in rewards.iter().enumerate() {
                bits[i] = match &slots[r] {
                    Value::Int(n) => n.to_i64()? as u64,
                    Value::Float(x) => {
                        floats |= 1 << i;
                        x.to_bits()
                    }
                    _ => return None,
                };
            }
            Some((bits, floats))
        };
        // Each state: its world, with dead slots cleared and rewards at 0,
        // found by its other live slots.
        let mut states: Vec<Vec<Value>> = Vec::new();
        let mut index: FxHashMap<Vec<Value>, usize> = FxHashMap::default();
        let mut intern = |mut slots: Vec<Value>, states: &mut Vec<Vec<Value>>| -> usize {
            for &d in &dead {
                slots[d as usize] = Value::Dead;
            }
            for &r in &rewards {
                slots[r] = Value::Int(0.into());
            }
            let key: Vec<Value> = head.iter().map(|&i| slots[i].clone()).collect();
            *index.entry(key).or_insert_with(|| {
                states.push(slots);
                states.len() - 1
            })
        };
        let total = total_weight(inside);
        let mut start: Vec<f64> = Vec::new();
        // With rewards: each state and reward the worlds start with.
        let mut starts: Vec<((usize, Raw), f64)> = Vec::new();
        for w in inside {
            let reward = if rewards.is_empty() {
                ([0; MAX_REWARDS], 0)
            } else {
                match reward_of(&w.slots) {
                    Some(r) => r,
                    None => return Ok(None),
                }
            };
            let i = intern(w.slots.clone(), &mut states);
            if start.len() <= i {
                start.resize(i + 1, 0.0);
            }
            start[i] += w.weight.ratio(total);
            starts.push(((i, reward), w.weight.ratio(total)));
        }
        // With rewards: each state's steps, to each state with how much the
        // rewards grow.
        let mut moves: Vec<Vec<(u32, Raw, f64)>> = Vec::new();
        let mut chain = Chain::default();
        let mut exits: Vec<Vec<World>> = Vec::new();
        let mut returns: Vec<Vec<Returned>> = Vec::new();
        let mut unresolved: Vec<Weight> = Vec::new();
        let mut failed: Vec<Failures> = Vec::new();
        let mut faulted: Vec<Vec<(World, RuntimeError)>> = Vec::new();
        // A body that uses a call's result so far can't be solved as a chain:
        // that result changes from round to round.
        let saved_pending = std::mem::take(&mut self.pending);
        let mut k = 0;
        while k < states.len() {
            if !self.pending.is_empty() {
                self.pending = saved_pending;
                return Ok(None);
            }
            if states.len() > self.config.max_chain_states {
                self.pending = saved_pending;
                self.too_large.insert(stmt.id);
                return Ok(None);
            }
            let world = World {
                slots: states[k].clone(),
                constraints: Default::default(),
                inherited: Default::default(),
                weight: Weight::ONE,
                run: 0,
                densities: 0,
            };
            let saved_unresolved = std::mem::replace(&mut self.unresolved, Weight::ZERO);
            let saved_lost = std::mem::replace(&mut self.lost, Weight::ZERO);
            let saved_failures = std::mem::take(&mut self.failures);
            let flow = self.exec_block(f, body, vec![world]);
            let left = std::mem::replace(&mut self.unresolved, saved_unresolved);
            let lost = std::mem::replace(&mut self.lost, saved_lost);
            let failures = std::mem::replace(&mut self.failures, saved_failures);
            let mut flow = flow?;
            if flow
                .next
                .iter()
                .chain(&flow.continued)
                .chain(&flow.broke)
                .chain(flow.faulted.iter().map(|(w, _)| w))
                .any(|w| w.densities > 0 || !w.constraints.is_empty() || w.slots.iter().any(analytic::contains))
                || flow
                    .returned
                    .iter()
                    .any(|(v, _, c)| !c.is_empty() || analytic::contains(v))
            {
                self.pending = saved_pending;
                self.too_large.insert(stmt.id);
                return Ok(None);
            }
            clear_dead(&mut flow.broke, after);
            // Leaving: by `break` or `return`, unresolved, ruled out,
            // failed, or faulted to a catch outside the loop.
            let mut leave = left + lost + failures.weight;
            let faults = std::mem::take(&mut flow.faulted);
            for (w, _) in &faults {
                leave += w.weight;
            }
            for w in &flow.broke {
                leave += w.weight;
            }
            for (_, w, _) in &flow.returned {
                leave += *w;
            }
            let mut next: FxHashMap<usize, f64> = FxHashMap::default();
            let mut steps: Vec<(u32, Raw, f64)> = Vec::new();
            for w in flow.next.into_iter().chain(flow.continued) {
                let p = w.weight.to_f64();
                if p == 0.0 {
                    // Too small for the chain's arithmetic.
                    self.too_large.insert(stmt.id);
                    self.pending = saved_pending;
                    return Ok(None);
                }
                let grows = if rewards.is_empty() {
                    ([0; MAX_REWARDS], 0)
                } else {
                    match reward_of(&w.slots) {
                        Some(r) => r,
                        None => {
                            self.too_large.insert(stmt.id);
                            self.pending = saved_pending;
                            return Ok(None);
                        }
                    }
                };
                let j = intern(w.slots, &mut states);
                *next.entry(j).or_insert(0.0) += p;
                if !rewards.is_empty() {
                    steps.push((j as u32, grows, p));
                }
            }
            moves.push(steps);
            let mut next: Vec<(usize, f64)> = next.into_iter().collect();
            next.sort_by_key(|&(j, _)| j);
            chain.next.push(next);
            chain.leave.push(leave.to_f64());
            exits.push(flow.broke);
            returns.push(flow.returned);
            unresolved.push(left);
            failed.push(failures);
            faulted.push(faults);
            k += 1;
        }
        let waiting = !self.pending.is_empty();
        self.pending = saved_pending;
        if waiting {
            return Ok(None);
        }
        start.resize(states.len(), 0.0);
        // Rewards are iterated, which needs no elimination: only states that
        // can't be left are found first.
        if !rewards.is_empty() {
            if let Some(s) = chain.stuck() {
                return Err(self.never_leaves(f, stmt, &states[s], &head));
            }
            drop(index);
            drop(states);
            let chain = Rewarded {
                rewards: &rewards,
                moves: &moves,
                exits,
                returns,
                unresolved,
                failed,
                faulted,
            };
            let solved = self.iterate_rewards(stmt, chain, starts, total, cutoff, iterations)?;
            if solved.is_none() {
                self.too_large.insert(stmt.id);
            }
            return Ok(solved);
        }
        let mut steps = MAX_ELIMINATION;
        let solution = chain.visits(&start, &mut steps);
        self.spend(MAX_ELIMINATION - steps, stmt.span)?;
        let visits = match solution {
            Solution::Visits(v) => v,
            Solution::TooBig => {
                self.too_large.insert(stmt.id);
                return Ok(None);
            }
            Solution::Stuck(s) => return Err(self.never_leaves(f, stmt, &states[s], &head)),
        };
        self.stats.solved_loops += 1;
        self.stats.chain_states += states.len() as u64;
        let mut flow = Flow::default();
        let mut left = Weight::ZERO;
        for (k, v) in visits.into_iter().enumerate() {
            if v == 0.0 {
                continue;
            }
            let times = total.scale(v);
            for mut w in std::mem::take(&mut exits[k]) {
                w.weight = w.weight * times;
                flow.next.push(w);
            }
            for (value, w, c) in std::mem::take(&mut returns[k]) {
                flow.returned.push((value, w * times, c));
            }
            left += unresolved[k] * times;
            self.failures.absorb(&failed[k], times, None, false);
            for (mut w, e) in std::mem::take(&mut faulted[k]) {
                w.weight = w.weight * times;
                flow.faulted.push((w, e));
            }
        }
        self.unresolved += left;
        Ok(Some(flow))
    }

    /// A loop solved as a chain whose steps carry rewards: the law of each
    /// state and reward, step after step, from the worlds that start it, as
    /// unrolling would follow them, but without running the body again.
    /// Like unrolling, it stops once what's still inside weighs less than
    /// `cutoff`, which is then unresolved, or after the most iterations.
    /// `None` if a reward doesn't fit the numbers it's followed in.
    fn iterate_rewards(
        &mut self,
        stmt: &Stmt,
        chain: Rewarded,
        starts: Vec<((usize, Raw), f64)>,
        total: Weight,
        cutoff: Weight,
        iterations: &mut u64,
    ) -> Result<Option<Flow>> {
        let span = stmt.span;
        let count = chain.rewards.len();
        // Each reward is followed as a whole number while every value it
        // takes is an int, and as a float once one is a float.
        let floats = starts
            .iter()
            .map(|((_, (_, f)), _)| *f)
            .chain(chain.moves.iter().flatten().map(|(_, (_, f), _)| *f))
            .fold(0u8, |a, f| a | f);
        let float: [bool; MAX_REWARDS] = std::array::from_fn(|i| floats & (1 << i) != 0);
        let compact = |(bits, floats): &Raw| -> Option<Reward> {
            let mut out = *bits;
            for i in 0..count {
                if float[i] && floats & (1 << i) == 0 {
                    let n = bits[i] as i64;
                    if n.unsigned_abs() > 1 << 53 {
                        return None;
                    }
                    out[i] = (n as f64).to_bits();
                }
            }
            Some(out)
        };
        let plus = |a: &Reward, b: &Reward| -> Option<Reward> {
            let mut out = [0u64; MAX_REWARDS];
            for i in 0..count {
                out[i] = if float[i] {
                    (f64::from_bits(a[i]) + f64::from_bits(b[i])).to_bits()
                } else {
                    (a[i] as i64).checked_add(b[i] as i64)? as u64
                };
            }
            Some(out)
        };
        let value = |r: &Reward, i: usize| -> Value {
            if float[i] {
                Value::Float(f64::from_bits(r[i]))
            } else {
                Value::Int((r[i] as i64).into())
            }
        };
        let mut moves: Vec<Vec<(u32, Reward, f64)>> = Vec::with_capacity(chain.moves.len());
        for steps in chain.moves {
            let mut out = Vec::with_capacity(steps.len());
            for (j, d, p) in steps {
                let Some(d) = compact(d) else { return Ok(None) };
                out.push((*j, d, *p));
            }
            moves.push(out);
        }
        let mut current: FxHashMap<(u32, Reward), f64> = FxHashMap::default();
        for ((s, r), p) in &starts {
            let Some(r) = compact(r) else { return Ok(None) };
            *current.entry((*s as u32, r)).or_insert(0.0) += p;
        }
        // How much weight ran the body from each state; leaves by each of
        // the distinct worlds that leave, with each reward, as unrolling's
        // merging would add it up; and faulted, by state and reward.
        let mut visits = vec![0.0; moves.len()];
        let mut templates: Vec<&World> = Vec::new();
        let mut template_index: FxHashMap<&[Value], u32> = FxHashMap::default();
        let exits: Vec<Vec<(u32, f64)>> = chain
            .exits
            .iter()
            .map(|worlds| {
                worlds
                    .iter()
                    .map(|w| {
                        let t = *template_index.entry(&w.slots[..]).or_insert_with(|| {
                            templates.push(w);
                            (templates.len() - 1) as u32
                        });
                        (t, w.weight.to_f64())
                    })
                    .collect()
            })
            .collect();
        drop(template_index);
        let mut leaving: FxHashMap<(u32, Reward), f64> = FxHashMap::default();
        let mut ran: FxHashMap<(u32, Reward), f64> = FxHashMap::default();
        let mut next: FxHashMap<(u32, Reward), f64> = FxHashMap::default();
        // When every step adds the same, and every world starts with the
        // same reward, as a counter of rounds does, all the worlds have the
        // same reward at each step: the weights are a vector by state.
        let mut increments = moves.iter().flatten().map(|(_, d, _)| d);
        let same = increments.next().copied().filter(|d| increments.all(|e| e == d));
        let mut rewards = current.keys().map(|(_, r)| r);
        let start = rewards.next().copied().filter(|r| rewards.all(|q| q == r));
        if let (Some(d), Some(mut reward)) = (same, start) {
            let n = moves.len();
            let mut weights = vec![0.0; n];
            for (&(s, _), &m) in &current {
                weights[s as usize] += m;
            }
            let mut after = vec![0.0; n];
            loop {
                let mass: f64 = weights.iter().sum();
                if mass == 0.0 {
                    break;
                }
                if total.scale(mass) < cutoff {
                    self.unresolved += total.scale(mass);
                    break;
                }
                *iterations += 1;
                if *iterations > self.config.max_iterations {
                    return Err(RuntimeError::limit(
                        span,
                        format!("this loop ran {} times without finishing", self.config.max_iterations),
                    )
                    .with_note(format!(
                        "worlds still inside the loop weigh {} of what entered it",
                        fmt_prob(mass)
                    ))
                    .with_help("check that every world can leave the loop"));
                }
                self.spend(n as u64, span)?;
                for (s, &m) in weights.iter().enumerate() {
                    if m == 0.0 {
                        continue;
                    }
                    visits[s] += m;
                    for &(j, _, p) in &moves[s] {
                        after[j as usize] += m * p;
                    }
                    for &(t, p) in &exits[s] {
                        *leaving.entry((t, reward)).or_insert(0.0) += m * p;
                    }
                    if !chain.faulted[s].is_empty() {
                        *ran.entry((s as u32, reward)).or_insert(0.0) += m;
                    }
                }
                std::mem::swap(&mut weights, &mut after);
                after.iter_mut().for_each(|w| *w = 0.0);
                let Some(r) = plus(&reward, &d) else { return Ok(None) };
                reward = r;
            }
            current.clear();
        }
        while !current.is_empty() {
            let mass: f64 = current.values().sum();
            if total.scale(mass) < cutoff {
                self.unresolved += total.scale(mass);
                break;
            }
            *iterations += 1;
            if *iterations > self.config.max_iterations {
                return Err(RuntimeError::limit(
                    span,
                    format!("this loop ran {} times without finishing", self.config.max_iterations),
                )
                .with_note(format!(
                    "worlds still inside the loop weigh {} of what entered it",
                    fmt_prob(mass)
                ))
                .with_help("check that every world can leave the loop"));
            }
            self.spend(current.len() as u64, span)?;
            for (&(s, r), &m) in &current {
                visits[s as usize] += m;
                for (j, d, p) in &moves[s as usize] {
                    let Some(r) = plus(&r, d) else { return Ok(None) };
                    *next.entry((*j, r)).or_insert(0.0) += m * p;
                }
                for &(t, p) in &exits[s as usize] {
                    *leaving.entry((t, r)).or_insert(0.0) += m * p;
                }
                if !chain.faulted[s as usize].is_empty() {
                    *ran.entry((s, r)).or_insert(0.0) += m;
                }
            }
            self.check_worlds(next.len(), span)?;
            std::mem::swap(&mut current, &mut next);
            next.clear();
        }
        self.stats.solved_loops += 1;
        self.stats.chain_states += moves.len() as u64;
        // A world leaving from a state with a reward: what it gained in
        // its last step, added to the reward.
        let add = |a: Value, b: &Value, budget: &mut crate::dist::Budget| -> Result<Value> {
            ops::binary(probl_syntax::ast::BinOp::Add, &a, b, budget).map_err(|e| e.at(span))
        };
        let mut flow = Flow::default();
        let rewarded = |w: &World, r: &Reward, budget: &mut crate::dist::Budget| -> Result<World> {
            let mut w = w.clone();
            for (i, &slot) in chain.rewards.iter().enumerate() {
                if !matches!(w.slots[slot], Value::Dead) {
                    w.slots[slot] = add(value(r, i), &w.slots[slot], budget)?;
                }
            }
            Ok(w)
        };
        let mut leaving: Vec<((u32, Reward), f64)> = leaving.into_iter().collect();
        leaving.sort_by_key(|a| a.0);
        for ((t, r), m) in leaving {
            let mut w = rewarded(templates[t as usize], &r, &mut self.budget)?;
            w.weight = total.scale(m);
            flow.next.push(w);
        }
        let mut ran: Vec<((u32, Reward), f64)> = ran.into_iter().collect();
        ran.sort_by_key(|a| a.0);
        for ((s, r), m) in ran {
            let times = total.scale(m);
            for (w, e) in &chain.faulted[s as usize] {
                let mut w = rewarded(w, &r, &mut self.budget)?;
                w.weight = w.weight * times;
                flow.faulted.push((w, e.clone()));
            }
        }
        let mut left = Weight::ZERO;
        for (s, v) in visits.into_iter().enumerate() {
            if v == 0.0 {
                continue;
            }
            let times = total.scale(v);
            for (value, w, c) in &chain.returns[s] {
                flow.returned.push((value.clone(), *w * times, c.clone()));
            }
            left += chain.unresolved[s] * times;
            self.failures.absorb(&chain.failed[s], times, None, false);
        }
        self.unresolved += left;
        Ok(Some(flow))
    }

    /// The error for a loop that some worlds can never leave, with one of
    /// them.
    fn never_leaves(&self, f: FnId, stmt: &Stmt, slots: &[Value], head: &[usize]) -> RuntimeError {
        let names = &self.prog.functions[f as usize].slots;
        let parts: Vec<String> = head
            .iter()
            .filter(|&&i| names[i].name != TEMP)
            .map(|&i| format!("`{}` is {:?}", names[i].name, slots[i]))
            .collect();
        let err = RuntimeError::new(stmt.span, "some worlds can never leave this loop")
            .with_help("check that every world can leave the loop");
        if parts.is_empty() {
            err
        } else {
            err.with_note(format!("for example, the worlds where {}", parts.join(" and ")))
        }
    }

    /// Store each possible value of `d` into `place`, one world per outcome.
    fn split_by(&mut self, f: FnId, place: &Place, d: Value, w: World, span: Span, out: &mut Vec<World>) -> Result<()> {
        match d {
            Value::Dist(_) | Value::Continuous(_) | Value::Counts(_) | Value::Joint(_) if self.sampler.is_some() => {
                let v = self.sample(&d);
                let mut w = w;
                self.assign(f, place, v, &mut w, span)?;
                out.push(w);
            }
            Value::Continuous(family) => {
                self.next_latent += 1;
                let v = Analytic::new(self.next_latent, *family)
                    .value()
                    .map_err(|e| e.at(span))?;
                let mut w = w;
                self.assign(f, place, v, &mut w, span)?;
                out.push(w);
            }
            Value::Dist(dist) => {
                self.unresolved += w.weight.scale(dist.missing);
                let last = dist.outcomes.len().saturating_sub(1);
                let mut w = Some(w);
                for (i, (v, p)) in dist.outcomes.iter().enumerate() {
                    let base = if i == last {
                        w.take().unwrap()
                    } else {
                        w.clone().unwrap()
                    };
                    let mut nw = base.scaled(*p);
                    if matches!(v, Value::Continuous(_) | Value::Joint(_)) {
                        self.split_by(f, place, v.clone(), nw, span, out)?;
                    } else {
                        self.assign(f, place, v.clone(), &mut nw, span)?;
                        out.push(nw);
                    }
                }
            }
            Value::Prob(p) => {
                return self.split_by(f, place, Dist::bernoulli(p).into_value(), w, span, out);
            }
            Value::Joint(dist) => self.draw_joint(f, place, &dist, w, span, out)?,
            // Enumerating needs its outcomes.
            Value::Counts(_) => {
                let d = ops::listed(&d, &mut self.budget).map_err(|e| e.at(span))?.into_owned();
                return self.split_by(f, place, d, w, span, out);
            }
            other => {
                let mut w = w;
                self.assign(f, place, other, &mut w, span)?;
                out.push(w);
            }
        }
        Ok(())
    }

    /// A draw from a distribution whose outcomes own their draws: a world
    /// for each outcome, with its draws fresh there.
    #[inline(never)]
    fn draw_joint(
        &mut self,
        f: FnId,
        place: &Place,
        dist: &Dist,
        w: World,
        span: Span,
        out: &mut Vec<World>,
    ) -> Result<()> {
        self.unresolved += w.weight.scale(dist.missing);
        for (v, p) in &dist.outcomes {
            let v = crate::joint::fresh(v, &mut self.next_latent);
            let mut nw = w.clone().scaled(*p);
            self.assign(f, place, v, &mut nw, span)?;
            out.push(nw);
        }
        Ok(())
    }

    fn check_arity(&self, c: &Closure, given: usize, span: Span) -> Result<()> {
        let expected = self.prog.functions[c.func as usize].n_params as usize;
        if expected != given {
            return Err(RuntimeError::new(
                span,
                format!("this function takes {expected} argument(s), but {given} were given"),
            ));
        }
        Ok(())
    }

    /// Convert values at declared type boundaries, recursively through
    /// containers. Numeric probability conversions are checked, never clamped.
    fn coerce(&mut self, v: Value, ty: &TypeSpec, span: Span) -> Result<Value> {
        self.coerce_at(v, ty, span, 0)
    }

    fn coerce_at(&mut self, v: Value, ty: &TypeSpec, span: Span, depth: usize) -> Result<Value> {
        self.budget.work(1).map_err(|e| e.at(span))?;
        if depth > 64 {
            return Err(RuntimeError::new(
                span,
                "type conversion nesting exceeds the limit of 64",
            ));
        }
        Ok(match (v, ty) {
            (Value::Analytic(_) | Value::Event(_), TypeSpec::Prob) => {
                return Err(analytic::unsupported("converting this outcome to prob").at(span));
            }
            (v @ Value::Float(_), TypeSpec::Int) => Value::Int(
                ops::integer(&v, "int conversion", &mut self.budget)
                    .map_err(|e| e.at(span))?
                    .into_owned(),
            ),
            // An outcome like `floor(x)` is an int; another is never rounded
            // into one.
            (Value::Analytic(a), TypeSpec::Int) if a.int => Value::Analytic(a),
            (Value::Analytic(_), TypeSpec::Int) => {
                return Err(analytic::unsupported("converting this outcome to int").at(span));
            }
            (Value::Counts(c), TypeSpec::Dist(t)) if **t == TypeSpec::Int => Value::Counts(c),
            (Value::Prob(p), TypeSpec::Float) => Value::Float(p),
            // A declared type is a contract: failing it isn't a fault of the
            // world's values, even in partial mode.
            (v @ (Value::Int(_) | Value::Float(_)), TypeSpec::Prob) => {
                ops::make_prob(&v).map_err(|e| OpError { fault: None, ..e }.at(span))?
            }
            (Value::List(xs), TypeSpec::List(t)) => {
                self.budget.collection(xs.len() as u128).map_err(|e| e.at(span))?;
                let mut out = Vec::with_capacity(xs.len());
                for x in xs.iter() {
                    out.push(self.coerce_at(x.clone(), t, span, depth + 1)?);
                }
                Value::list(out)
            }
            (Value::Map(xs), TypeSpec::Map(kt, vt)) => {
                let mut out = BTreeMap::new();
                for (k, v) in xs.iter() {
                    let key = self.coerce_at(k.clone(), kt, span, depth + 1)?;
                    let value = self.coerce_at(v.clone(), vt, span, depth + 1)?;
                    if out.insert(key, value).is_some() {
                        return Err(RuntimeError::new(span, "map key collision during type conversion"));
                    }
                }
                Value::map(out)
            }
            (Value::Bag(xs), TypeSpec::Bag(t)) => {
                let mut out = BTreeMap::<Value, u64>::new();
                for (x, n) in xs.iter() {
                    let key = self.coerce_at(x.clone(), t, span, depth + 1)?;
                    let count = out.entry(key).or_default();
                    *count = count
                        .checked_add(*n)
                        .ok_or_else(|| RuntimeError::new(span, "bag count overflow during type conversion"))?;
                }
                Value::multiset(crate::value::Multiset::new(out))
            }
            (Value::Dist(d), TypeSpec::Dist(t)) => {
                let mut out = Vec::with_capacity(d.outcomes.len());
                for (x, p) in &d.outcomes {
                    out.push((self.coerce_at(x.clone(), t, span, depth + 1)?, *p));
                }
                ops::combine(out, d.missing, &mut self.budget).map_err(|e| e.at(span))?
            }
            (Value::Joint(d), TypeSpec::Dist(t)) => {
                let mut out = Vec::with_capacity(d.outcomes.len());
                for (x, p) in &d.outcomes {
                    out.push((self.coerce_at(x.clone(), t, span, depth + 1)?, *p));
                }
                crate::joint::of(out, d.missing, &mut self.budget).map_err(|e| e.at(span))?
            }
            (Value::Record(r), TypeSpec::Record(t)) => {
                let types: Vec<_> = self.prog.records[*t as usize]
                    .fields
                    .iter()
                    .map(|f| (f.name.clone(), f.ty.clone()))
                    .collect();
                let mut fields = Vec::with_capacity(r.fields.len());
                for (name, value) in &r.fields {
                    let value = match types.iter().find(|(n, _)| n == &**name) {
                        Some((_, t)) => self.coerce_at(value.clone(), t, span, depth + 1)?,
                        None => value.clone(),
                    };
                    fields.push((name.clone(), value));
                }
                ops::make_record(r.ty.clone(), fields)
            }
            (Value::Record(r), TypeSpec::AnonRecord(types)) => {
                let mut fields = Vec::with_capacity(r.fields.len());
                for (name, value) in &r.fields {
                    let value = match types.iter().find(|(n, _)| n == &**name) {
                        Some((_, t)) => self.coerce_at(value.clone(), t, span, depth + 1)?,
                        None => value.clone(),
                    };
                    fields.push((name.clone(), value));
                }
                ops::make_record(r.ty.clone(), fields)
            }
            (v, _) => v,
        })
    }

    /// Does `v` have the declared type?
    fn conforms(&self, v: &Value, ty: &TypeSpec) -> bool {
        match (ty, v) {
            (TypeSpec::Int, Value::Int(_)) => true,
            (TypeSpec::Float, Value::Float(_) | Value::Int(_) | Value::Analytic(_)) => true,
            (TypeSpec::Complex, Value::Complex(_)) => true,
            (TypeSpec::Prob, Value::Prob(_)) => true,
            (TypeSpec::Bool, Value::Bool(_) | Value::Event(_)) => true,
            (TypeSpec::Str, Value::Str(_)) => true,
            (TypeSpec::Date, Value::Date(_)) => true,
            (TypeSpec::Unit, Value::Unit) => true,
            (TypeSpec::Function, Value::Closure(_) | Value::Builtin(_)) => true,
            (TypeSpec::List(t), Value::List(items)) => items.iter().all(|x| self.conforms(x, t)),
            (TypeSpec::List(t), Value::Range(..)) => **t == TypeSpec::Int,
            (TypeSpec::Map(k, t), Value::Map(m)) => m.iter().all(|(a, b)| self.conforms(a, k) && self.conforms(b, t)),
            (TypeSpec::Bag(t), Value::Bag(b)) => b.keys().all(|x| self.conforms(x, t)),
            (TypeSpec::Dist(t), Value::Dist(d) | Value::Joint(d)) => d.outcomes.iter().all(|(x, _)| match x {
                Value::Continuous(_) => **t == TypeSpec::Float,
                _ => self.conforms(x, t),
            }),
            (TypeSpec::Dist(t), Value::Continuous(_)) => **t == TypeSpec::Float,
            (TypeSpec::Dist(t), Value::Counts(_)) => **t == TypeSpec::Int,
            (TypeSpec::Record(r), Value::Record(rec)) => {
                let declared = &self.prog.records[*r as usize];
                rec.ty.as_deref() == Some(declared.name.as_str())
                    && rec.fields.len() == declared.fields.len()
                    && declared
                        .fields
                        .iter()
                        .all(|f| rec.get(&f.name).is_some_and(|x| self.conforms(x, &f.ty)))
            }
            (TypeSpec::Enum(e), Value::Enum(x)) => x.ty == *e,
            (TypeSpec::AnonRecord(fields), Value::Record(rec)) => {
                rec.fields.len() == fields.len()
                    && rec
                        .fields
                        .iter()
                        .zip(fields)
                        .all(|((n, x), (m, t))| &**n == m.as_str() && self.conforms(x, t))
            }
            _ => false,
        }
    }

    // ── Calls ────────────────────────────────────────────────────────────

    /// Run a function as a sub-simulation and return the unnormalized
    /// distribution of its result. Results are memoized unless the function
    /// prints: its behaviour depends only on its arguments and the outside
    /// values it reads, which are all in `key`.
    ///
    /// When enumerating, a call can come back to itself: the same function,
    /// with the same arguments, while it's still running. Its result is then
    /// solved by iteration (section 6). The call that comes back gets the
    /// result so far, which starts with nothing resolved, and the call it
    /// came back to runs again with each new result, until the weight still
    /// waiting on it is below ε (`run_call`).
    fn call(&mut self, func: FnId, key: Vec<Value>, span: Span) -> Result<Arc<CallResult>> {
        self.stats.calls += 1;
        let prog = self.prog;
        let fun = &prog.functions[func as usize];
        // When sampling, every call makes its own choices (section 14).
        let sampling = self.sampler.is_some();
        // A cached result does not prove a callback's effects were permitted:
        // a draw or observation can collapse to one outcome of weight one.
        let memoizable = self.config.memoizing && !sampling && !fun.effects.prints && self.callback.is_none();
        let call_key = (func, key);
        if memoizable {
            if let Some(r) = self.memo.get(&call_key) {
                self.stats.memo_hits += 1;
                self.observed |= r.observed;
                return Ok(r.clone());
            }
        }
        if !sampling {
            if let Some(&depth) = self.active.get(&call_key) {
                if fun.effects.prints {
                    return Err(RuntimeError::new(
                        span,
                        format!("`{}` comes back to itself, and prints", describe_call(fun, &call_key.1)),
                    )
                    .with_note("a call that comes back to itself runs several times, until its result settles, and would print each time")
                    .with_help("remove the `print`, or print the result where the call is made"));
                }
                let so_far = match self.approx.get(&call_key) {
                    Some(a) if self.frames_at.get(a.head) == Some(&a.frame) => a.result.clone(),
                    _ => {
                        let result = Arc::new(CallResult::waiting(depth));
                        let approx = Approx {
                            result: result.clone(),
                            head: depth,
                            frame: self.frames_at[depth],
                            round: self.rounds_at[depth],
                        };
                        self.approx.insert(call_key, approx);
                        result
                    }
                };
                return Ok(so_far);
            }
            // Worked out already, in this round of the call it waits on.
            if let Some(a) = self.approx.get(&call_key) {
                if self.frames_at.get(a.head) == Some(&a.frame) && self.rounds_at.get(a.head) == Some(&a.round) {
                    self.stats.memo_hits += 1;
                    return Ok(a.result.clone());
                }
            }
        }
        if self.depth >= self.config.max_call_depth {
            return Err(RuntimeError::limit(
                span,
                format!("calls are nested more than {} deep", self.config.max_call_depth),
            )
            .with_help("check for recursion that doesn't stop, or write it as a loop"));
        }
        let mut slots = vec![Value::Dead; fun.n_slots()];
        let n = fun.n_params as usize;
        for (i, v) in call_key.1[..n].iter().enumerate() {
            slots[i] = v.clone();
        }
        for (cap, v) in fun.captures.iter().zip(&call_key.1[n..]) {
            slots[cap.slot as usize] = v.clone();
        }
        let depth = self.depth + 1;
        if !sampling {
            self.active.insert(call_key.clone(), depth);
            if self.frames_at.len() <= depth {
                self.frames_at.resize(depth + 1, 0);
                self.rounds_at.resize(depth + 1, 0);
            }
            self.last_number += 1;
            self.frames_at[depth] = self.last_number;
        }
        let result = self.run_call(func, &call_key, slots, depth, span);
        if !sampling {
            self.active.remove(&call_key);
            // Results that waited on this frame no longer stand for anything.
            self.frames_at[depth] = 0;
        }
        let result = result?;
        // A result still waiting on a call being solved isn't final.
        if memoizable
            && result.pending.is_empty()
            && result
                .outcomes
                .iter()
                .all(|(v, _, c)| c.is_empty() && !analytic::contains(v))
        {
            if self.memo.len() >= self.config.max_cached_calls {
                self.memo.clear();
            }
            self.memo.insert(call_key, result.clone());
        }
        Ok(result)
    }

    /// Run a call's body from `slots`, at `depth` in the stack, and return
    /// its result (section 6).
    ///
    /// If the call came back to itself, and waits on no call further down
    /// the stack, it's solved here: round after round, each using the result
    /// of the round before, until the weight still waiting on it is below ε.
    /// That weight stays unresolved. If it waits on a call further down,
    /// it's part of that call's solving instead: what waits on it waits on
    /// that call, and its result stands for it until that call's next round.
    fn run_call(
        &mut self,
        func: FnId,
        call_key: &CallKey,
        slots: Vec<Value>,
        depth: usize,
        span: Span,
    ) -> Result<Arc<CallResult>> {
        let fun = &self.prog.functions[func as usize];
        let sampling = self.sampler.is_some();
        let mut inherited = std::collections::BTreeSet::new();
        for v in &slots {
            analytic::collect_ids(v, &mut inherited);
        }
        let inherited = Arc::new(inherited);
        let mut rounds: u64 = 0;
        let mut before = f64::INFINITY;
        loop {
            if !sampling {
                self.last_number += 1;
                self.rounds_at[depth] = self.last_number;
            }
            let saved_unresolved = std::mem::replace(&mut self.unresolved, Weight::ZERO);
            let saved_lost = std::mem::replace(&mut self.lost, Weight::ZERO);
            let saved_observed = std::mem::replace(&mut self.observed, false);
            let saved_pending = std::mem::take(&mut self.pending);
            let saved_failures = std::mem::take(&mut self.failures);
            self.depth += 1;
            let flow = self.exec_block(
                func,
                &fun.body,
                vec![World {
                    slots: slots.clone(),
                    constraints: Default::default(),
                    inherited: inherited.clone(),
                    weight: Weight::ONE,
                    run: 0,
                    densities: 0,
                }],
            );
            self.depth -= 1;
            let unresolved = std::mem::replace(&mut self.unresolved, saved_unresolved);
            let lost = std::mem::replace(&mut self.lost, saved_lost);
            let observed = std::mem::replace(&mut self.observed, saved_observed);
            let pending = std::mem::replace(&mut self.pending, saved_pending);
            let failures = std::mem::replace(&mut self.failures, saved_failures);
            self.observed |= observed;
            let flow = flow.map_err(|e| match fun.kind {
                FnKind::Named => e.with_note(format!("in a call to `{}`", fun.name)),
                FnKind::Simulate => e.with_note("inside a `simulate` block"),
                FnKind::Lambda => e.with_note("inside a lambda"),
                FnKind::Main => e,
            })?;
            uncaught(&flow)?;
            let mut result = CallResult {
                outcomes: merge_values(flow.returned),
                unresolved,
                lost,
                observed,
                pending,
                failures,
            };
            if let Some(head) = result.pending.iter().map(|&(d, _)| d).filter(|&d| d < depth).min() {
                // Part of the solving of the call at `head`.
                let waiting = Weight::sum(result.pending.iter().map(|&(_, w)| w));
                result.pending = vec![(head, waiting)];
                let result = Arc::new(result);
                let approx = Approx {
                    result: result.clone(),
                    head,
                    frame: self.frames_at[head],
                    round: self.rounds_at[head],
                };
                self.approx.insert(call_key.clone(), approx);
                return Ok(result);
            }
            let waiting = result.take_pending(depth);
            if waiting.is_zero() {
                return Ok(Arc::new(result));
            }
            // It came back to itself: solved here.
            if rounds == 0 {
                self.stats.solved_calls += 1;
            }
            let w = waiting.to_f64();
            if w <= self.config.epsilon {
                self.settle(depth);
                return Ok(Arc::new(result));
            }
            rounds += 1;
            self.stats.call_rounds += 1;
            if w >= before {
                return Err(RuntimeError::new(
                    span,
                    format!(
                        "`{}` never returns for some of its worlds",
                        describe_call(fun, &call_key.1)
                    ),
                )
                .with_note(format!("{} of its weight keeps coming back to it", fmt_prob(w)))
                .with_help("check that the recursion can end"));
            }
            if rounds >= self.config.max_iterations {
                return Err(RuntimeError::limit(
                    span,
                    format!(
                        "`{}` came back to itself {} times without settling",
                        describe_call(fun, &call_key.1),
                        self.config.max_iterations
                    ),
                )
                .with_note(format!("{} of its weight is still waiting on it", fmt_prob(w)))
                .with_help("check that the recursion can end, or raise `@epsilon`"));
            }
            before = w;
            result.pending.push((depth, waiting));
            let approx = Approx {
                result: Arc::new(result),
                head: depth,
                frame: self.frames_at[depth],
                round: self.rounds_at[depth],
            };
            self.approx.insert(call_key.clone(), approx);
        }
    }

    /// The call at `depth` is solved: the results that waited on it, from
    /// its last round, are final too if the weight still waiting in them is
    /// below ε. The others are dropped.
    fn settle(&mut self, depth: usize) {
        let frame = self.frames_at[depth];
        let epsilon = self.config.epsilon;
        let settled: Vec<(CallKey, CallResult)> = self
            .approx
            .iter()
            .filter(|(_, a)| a.head == depth && a.frame == frame)
            .filter(|(_, a)| a.result.pending.iter().all(|&(_, w)| w.to_f64() <= epsilon))
            .map(|(key, a)| {
                let mut result = (*a.result).clone();
                // What's still waiting stays unresolved.
                result.pending.clear();
                (key.clone(), result)
            })
            .collect();
        self.approx.retain(|_, a| a.head != depth);
        if !self.config.memoizing {
            return;
        }
        for (key, result) in settled {
            if !self.active.contains_key(&key)
                && !self.prog.functions[key.0 as usize].effects.prints
                && result
                    .outcomes
                    .iter()
                    .all(|(v, _, c)| c.is_empty() && !analytic::contains(v))
            {
                self.memo.insert(key, Arc::new(result));
            }
        }
    }

    /// `simulate { … }`: run a block as a separate model and return its
    /// normalized distribution (docs/semantics.md, section 8).
    fn simulate(&mut self, func: FnId, key: Vec<Value>, span: Span) -> Result<Value> {
        if key.iter().any(analytic::contains) {
            return Err(analytic::unsupported("capturing an analytic draw inside `simulate`").at(span));
        }
        if !self.simulated.is_empty() {
            if let Some(v) = self.simulated.get(&(func, key.clone())) {
                self.stats.memo_hits += 1;
                return Ok(v.clone());
            }
        }
        let call_key = self.config.memoizing.then(|| (func, key.clone()));
        let saved_observed = self.observed;
        let saved_densities = self.densities;
        let saved_local = (self.local_density_at.take(), self.returned_densities.take());
        // Enumerated, even when sampling (section 14).
        let sampler = self.sampler.take();
        let callback = self.callback.take();
        self.nested += 1;
        let result = self.call(func, key, span);
        self.nested -= 1;
        self.sampler = sampler;
        self.callback = callback;
        // Observations inside `simulate` condition its result only.
        self.observed = saved_observed;
        self.densities = saved_densities;
        let local_density_at = std::mem::replace(&mut self.local_density_at, saved_local.0);
        self.returned_densities = saved_local.1;
        let result = result?;
        if let Some(span) = local_density_at.filter(|_| !result.unresolved.is_zero()) {
            return Err(OpError::unsupported(
                "observing a value from a continuous distribution while some probability is unresolved isn't supported when enumerating",
            )
            .help("the unresolved worlds could have any density at the value, so nothing would bound the evidence; sample the model with `@mode sample(runs: 10_000)`")
            .at(span));
        }
        if !result.pending.is_empty() {
            return Err(OpError::unsupported(
                "a `simulate` block that comes back to a call still running isn't supported yet",
            )
            .at(span));
        }
        // It's built whole or not at all: a fault on one of its paths, not
        // caught inside it, fails the world that runs it.
        if let Some(g) = result.failures.groups.first() {
            return Err(g.error.clone());
        }
        let resolved = Weight::sum(result.outcomes.iter().map(|(_, w, _)| *w));
        if resolved.is_zero() {
            let message = if result.unresolved.is_zero() {
                "every world in this `simulate` was ruled out by `observe`"
            } else {
                "this `simulate` left all of its weight unresolved"
            };
            return Err(RuntimeError::new(span, message));
        }
        let denom = resolved + result.unresolved;
        let pairs: Vec<(Value, f64)> = result
            .outcomes
            .iter()
            .map(|(v, w, _)| (v.clone(), w.ratio(denom)))
            .collect();
        let missing = result.unresolved.ratio(denom);
        // Outcomes of its own draws make a joint distribution.
        if pairs.iter().any(|(v, _)| analytic::contains(v)) {
            let v = crate::joint::of(pairs, missing, &mut self.budget).map_err(|e| e.at(span))?;
            if let Some(call_key) = call_key {
                if self.simulated.len() >= self.config.max_cached_calls {
                    self.simulated.clear();
                }
                self.simulated.insert(call_key, v.clone());
            }
            return Ok(v);
        }
        ops::combine(pairs, missing, &mut self.budget).map_err(|e| e.at(span))
    }

    /// Call a closure that must not split worlds (used by `map`, `filter`, …).
    fn call_pure(&mut self, closure: &Value, args: Vec<Value>, what: &'static str, span: Span) -> Result<Value> {
        if let Value::Builtin(b) = closure {
            return self.builtin_values(*b, &args, &[], Weight::ONE, span);
        }
        let Value::Closure(c) = closure else {
            return Err(RuntimeError::new(
                span,
                format!(
                    "`{what}` needs a function, like `x -> x * 2`, found {}",
                    ops::article(&closure.kind())
                ),
            ));
        };
        self.check_arity(c, args.len(), span)?;
        let mut key = args;
        key.extend(c.captured.iter().cloned());
        let previous = self.callback.replace(what);
        let result = self.call(c.func, key, span);
        self.callback = previous;
        let result = result?;
        if !result.pending.is_empty() {
            return Err(OpError::unsupported(format!(
                "a function given to `{what}` that comes back to a call still running isn't supported yet"
            ))
            .at(span));
        }
        // A callback that faults fails the whole operation, in the world
        // that runs it.
        if let Some(g) = result.failures.groups.first() {
            return Err(g.error.clone());
        }
        match result.outcomes.as_slice() {
            [(v, p, c)] if c.is_empty() && (p.to_f64() - 1.0).abs() < 1e-12 && result.unresolved.is_zero() => {
                Ok(v.clone())
            }
            _ => Err(RuntimeError::new(
                span,
                format!("the function given to `{what}` can't branch on chances, draw values or observe"),
            )),
        }
    }

    // ── Places ───────────────────────────────────────────────────────────

    fn check_callback_effect(&self, span: Span) -> Result<()> {
        match self.callback {
            Some(what) => Err(RuntimeError::new(
                span,
                format!("the function given to `{what}` can't branch on chances, draw values or observe"),
            )
            .with_help("use a loop for probabilistic traversal, or `simulate` for a local distribution")),
            None => Ok(()),
        }
    }

    fn record_place_type(&self, place: &Place, w: &World) -> Option<TypeSpec> {
        if place.path.is_empty() {
            return None;
        }
        let Value::Record(r) = &w.slots[place.slot as usize] else {
            return None;
        };
        let id = self
            .prog
            .records
            .iter()
            .position(|t| Some(t.name.as_str()) == r.ty.as_deref())?;
        let mut ty = TypeSpec::Record(id as u32);
        for element in &place.path {
            ty = match (element, ty) {
                (PathElem::Field(name), TypeSpec::Record(id)) => self.prog.records[id as usize]
                    .fields
                    .iter()
                    .find(|f| f.name == *name)?
                    .ty
                    .clone(),
                (PathElem::Field(name), TypeSpec::AnonRecord(fields)) => fields.into_iter().find(|(n, _)| n == name)?.1,
                (PathElem::Index(_), TypeSpec::List(t) | TypeSpec::Map(_, t)) => *t,
                _ => return None,
            };
        }
        Some(ty)
    }

    fn assign(&mut self, f: FnId, place: &Place, v: Value, w: &mut World, span: Span) -> Result<()> {
        if place.path.is_empty() {
            w.slots[place.slot as usize] = v;
            return Ok(());
        }
        let v = match self.record_place_type(place, w) {
            Some(ty) => {
                let v = self.coerce(v, &ty, span)?;
                if !self.conforms(&v, &ty) {
                    return Err(RuntimeError::new(
                        span,
                        format!(
                            "assigned field should be {}, found {}",
                            ty.describe(self.prog),
                            v.kind()
                        ),
                    ));
                }
                v
            }
            None => v,
        };
        let mut keys = Vec::with_capacity(place.path.len());
        for elem in &place.path {
            keys.push(match elem {
                PathElem::Field(name) => PathKey::Field(name.as_str()),
                PathElem::Index(e) => {
                    let key = self.eval(f, e, w)?;
                    if analytic::contains(&key) {
                        return Err(analytic::unsupported("an analytic assignment index").at(e.span));
                    }
                    PathKey::Index(key)
                }
            });
        }
        let target = &mut w.slots[place.slot as usize];
        if matches!(target, Value::Dead) {
            return Err(self.no_value(f, place.slot, span));
        }
        update(target, &keys, v, &mut self.budget).map_err(|e| e.at(span))
    }

    /// Convert a variable to its declared type, or fail.
    fn check_slot(&mut self, f: FnId, slot: SlotId, ty: &TypeSpec, w: &mut World, span: Span) -> Result<()> {
        w.slots[slot as usize] = self.coerce(w.slots[slot as usize].clone(), ty, span)?;
        if !self.conforms(&w.slots[slot as usize], ty) {
            return Err(self.declared_type_error(f, slot, ty, w, span));
        }
        Ok(())
    }

    fn declared_type_error(&self, f: FnId, slot: SlotId, ty: &TypeSpec, w: &World, span: Span) -> RuntimeError {
        let name = &self.prog.functions[f as usize].slots[slot as usize].name;
        let what = if name == TEMP {
            "the result".to_string()
        } else {
            format!("`{name}`")
        };
        RuntimeError::new(
            span,
            format!(
                "{what} should be {}, but it's {}",
                ops::article(&ty.describe(self.prog)),
                ops::article(&w.slots[slot as usize].kind())
            ),
        )
    }

    /// `place = value`, where `value` applies a method like `push` to
    /// `place`, then the check of a variable with a declared type. Out of
    /// line: every statement runs the statement loop.
    #[inline(never)]
    fn update(
        &mut self,
        f: FnId,
        place: &Place,
        value: &Expr,
        check: Option<&UpdateCheck>,
        w: &mut World,
        span: Span,
    ) -> Result<()> {
        let ExprKind::Builtin { func, args, named } = &value.kind else {
            return Err(RuntimeError::new(span, "internal error: an update that isn't a method"));
        };
        let mut values = Vec::with_capacity(args.len() + named.len());
        for a in args.iter().chain(named.iter().map(|(_, a)| a)) {
            values.push(self.eval(f, a, w)?);
        }
        let field = self.record_place_type(place, w);
        let types = match check {
            Some(c) => c.args.as_deref().map(Some),
            None => Some(None),
        };
        if let (None, Some(types), true) = (&field, types, named.is_empty()) {
            if self.in_place(f, place, *func, &mut values, types, check, w, span)? {
                return Ok(());
            }
        }
        let names: Vec<String> = named.iter().map(|(n, _)| n.clone()).collect();
        let mut v = self.builtin_values(*func, &values, &names, w.weight, span)?;
        if let Some(ty) = field {
            v = self.coerce(v, &ty, value.span)?;
        }
        self.assign(f, place, v, w, span)?;
        if let Some(c) = check {
            self.check_slot(f, place.slot, &c.ty, w, span)?;
        }
        Ok(())
    }

    /// Change the list or map at `place` as `func` would, given the values
    /// of its arguments, the collection first, without copying it unless
    /// something else holds it. What it adds is checked against `types`, as
    /// a check of the whole variable would. False, with nothing changed,
    /// where it can't, as for a distribution, which the method lifts over.
    #[allow(clippy::too_many_arguments)]
    fn in_place(
        &mut self,
        f: FnId,
        place: &Place,
        func: Builtin,
        values: &mut [Value],
        types: Option<&[Option<TypeSpec>]>,
        check: Option<&UpdateCheck>,
        w: &mut World,
        span: Span,
    ) -> Result<bool> {
        let plain = |v: &Value| {
            !v.is_uncertain()
                && !matches!(
                    v,
                    Value::Analytic(_) | Value::Event(_) | Value::Delayed(_) | Value::Dead
                )
        };
        let len = match (func, &values[0]) {
            (Builtin::Push | Builtin::Insert | Builtin::Remove | Builtin::DropLast, Value::List(items)) => items.len(),
            (Builtin::Insert | Builtin::Remove, Value::Map(m)) => m.len(),
            _ => return Ok(false),
        };
        if !values[1..].iter().all(plain) {
            return Ok(false);
        }
        let at = |e: OpError| e.at(span);
        // Whatever can fail does before anything changes.
        let list = matches!(values[0], Value::List(_));
        let index = match func {
            Builtin::Push | Builtin::Insert => {
                self.budget.collection(len as u128 + 1).map_err(at)?;
                match (func, list) {
                    (Builtin::Insert, true) => {
                        Some(ops::as_index(&values[1], len as u128 + 1, &mut self.budget).map_err(at)?)
                    }
                    _ => None,
                }
            }
            Builtin::Remove if list => Some(ops::as_index(&values[1], len as u128, &mut self.budget).map_err(at)?),
            _ => None,
        };
        if let (Some(types), Some(check)) = (types, check) {
            for (i, ty) in types.iter().enumerate() {
                let Some(ty) = ty else { continue };
                let given = &values[i + 1];
                let v = self.coerce(given.clone(), ty, span)?;
                if !self.conforms(&v, ty) {
                    return Err(self.declared_type_error(f, place.slot, &check.ty, w, span));
                }
                // A key that converts could meet one already there, which the
                // check of the whole map reports.
                if !list && i == 0 && v != *given {
                    return Ok(false);
                }
                values[i + 1] = v;
            }
        }
        let mut keys = Vec::with_capacity(place.path.len());
        for elem in &place.path {
            keys.push(match elem {
                PathElem::Field(name) => PathKey::Field(name.as_str()),
                PathElem::Index(e) => {
                    let key = self.eval(f, e, w)?;
                    if analytic::contains(&key) {
                        return Ok(false);
                    }
                    PathKey::Index(key)
                }
            });
        }
        // The value read for the method is another reference to the
        // collection, which would make changing it copy it.
        values[0] = Value::Unit;
        let target = place_mut(&mut w.slots[place.slot as usize], &keys, &mut self.budget).map_err(at)?;
        match (func, target) {
            (Builtin::Push, Value::List(items)) => {
                Arc::make_mut(items).push(std::mem::replace(&mut values[1], Value::Unit))
            }
            (Builtin::Insert, Value::List(items)) => Arc::make_mut(items).insert(
                index.unwrap_or_default() as usize,
                std::mem::replace(&mut values[2], Value::Unit),
            ),
            (Builtin::Remove, Value::List(items)) => {
                Arc::make_mut(items).remove(index.unwrap_or_default() as usize);
            }
            (Builtin::DropLast, Value::List(items)) => {
                Arc::make_mut(items).pop();
            }
            (Builtin::Insert, Value::Map(m)) => {
                let value = std::mem::replace(&mut values[2], Value::Unit);
                Arc::make_mut(m).insert(std::mem::replace(&mut values[1], Value::Unit), value);
            }
            (Builtin::Remove, Value::Map(m)) => {
                Arc::make_mut(m).remove(&values[1]);
            }
            _ => {
                return Err(RuntimeError::new(
                    span,
                    "internal error: an update found another collection",
                ));
            }
        }
        Ok(true)
    }

    fn read_place(&mut self, f: FnId, place: &Place, w: &World, span: Span) -> Result<Value> {
        let mut v = self.slot(f, place.slot, w, span)?;
        for elem in &place.path {
            v = match elem {
                PathElem::Field(name) => ops::field(&v, name, &mut self.budget),
                PathElem::Index(e) => {
                    let i = self.eval(f, e, w)?;
                    ops::index(&v, &i, &mut self.budget)
                }
            }
            .map_err(|e| e.at(span))?;
        }
        Ok(v)
    }

    /// Add the values `callee` captures to `out`, from the slots of `f` that
    /// hold them. A global it holds only for the function values it may call
    /// can have no value yet, as can any, with `lenient`: that's an error
    /// only where the global is read.
    #[allow(clippy::too_many_arguments)]
    fn capture(
        &mut self,
        out: &mut Vec<Value>,
        f: FnId,
        callee: FnId,
        args: &[SlotId],
        w: &World,
        span: Span,
        lenient: bool,
    ) -> Result<()> {
        let prog = self.prog;
        let captures = &prog.functions[callee as usize].captures;
        for (&s, capture) in args.iter().zip(captures) {
            out.push(match &w.slots[s as usize] {
                Value::Dead if lenient || capture.for_values => Value::Dead,
                _ => self.slot(f, s, w, span)?,
            });
        }
        Ok(())
    }

    fn slot(&mut self, f: FnId, s: SlotId, w: &World, span: Span) -> Result<Value> {
        match &w.slots[s as usize] {
            Value::Dead => Err(self.no_value(f, s, span)),
            Value::Delayed(_) => {
                let name = &self.prog.functions[f as usize].slots[s as usize].name;
                Err(
                    OpError::internal("internal error: a variable was read before it was drawn")
                        .at(span)
                        .with_note(format!(
                            "`{name}`'s draw was delayed for an exact update (docs/semantics.md, section 14)"
                        )),
                )
            }
            v => analytic::resolve(v, &w.constraints, &mut self.budget).map_err(|e| e.at(span)),
        }
    }

    fn no_value(&self, f: FnId, s: SlotId, span: Span) -> RuntimeError {
        let name = &self.prog.functions[f as usize].slots[s as usize].name;
        RuntimeError::new(span, format!("`{name}` has no value here")).with_help(
            "it's used before it's given a value; if a function reads it, call the function after the variable is set",
        )
    }

    // ── Expressions ──────────────────────────────────────────────────────

    fn eval_expected(&mut self, f: FnId, e: &Expr, w: &World, ty: &TypeSpec) -> Result<Value> {
        match (&e.kind, ty) {
            (ExprKind::List(xs), TypeSpec::List(t)) => {
                self.budget.collection(xs.len() as u128).map_err(|err| err.at(e.span))?;
                let mut out = Vec::with_capacity(xs.len());
                for x in xs {
                    out.push(self.eval_expected(f, x, w, t)?);
                }
                Ok(Value::list(out))
            }
            (ExprKind::Map(xs), TypeSpec::Map(kt, vt)) => {
                let mut out = BTreeMap::new();
                for (k, v) in xs {
                    let key = self.eval_expected(f, k, w, kt)?;
                    if analytic::contains(&key) {
                        return Err(analytic::unsupported("analytic map keys").at(k.span));
                    }
                    if key.is_uncertain() {
                        return Err(RuntimeError::new(k.span, "map keys can't be distributions"));
                    }
                    out.insert(key, self.eval_expected(f, v, w, vt)?);
                }
                Ok(Value::map(out))
            }
            (ExprKind::Record { ty: None, fields }, TypeSpec::AnonRecord(types)) => {
                let mut out = Vec::with_capacity(fields.len());
                for (n, e) in fields {
                    let v = match types.iter().find(|(name, _)| name == n) {
                        Some((_, ty)) => self.eval_expected(f, e, w, ty)?,
                        None => self.eval(f, e, w)?,
                    };
                    out.push((Arc::from(n.as_str()), v));
                }
                Ok(ops::make_record(None, out))
            }
            _ => {
                let v = self.eval(f, e, w)?;
                self.coerce(v, ty, e.span)
            }
        }
    }

    /// A shared field type supplies conversion context even when the record is
    /// represented by a finite distribution. Heterogeneous records have no
    /// shared context, and their individual declarations are checked below.
    fn record_field_type(&self, value: &Value, name: &str) -> Option<TypeSpec> {
        match value {
            Value::Record(r) => self
                .prog
                .records
                .iter()
                .find(|t| Some(t.name.as_str()) == r.ty.as_deref())?
                .fields
                .iter()
                .find(|f| f.name == name)
                .map(|f| f.ty.clone()),
            Value::Dist(d) => {
                let first = self.record_field_type(&d.outcomes.first()?.0, name)?;
                d.outcomes
                    .iter()
                    .all(|(v, _)| self.record_field_type(v, name).as_ref() == Some(&first))
                    .then_some(first)
            }
            _ => None,
        }
    }

    fn check_updated_records(&mut self, value: Value, span: Span) -> Result<Value> {
        match value {
            Value::Record(ref r) => {
                let Some(id) = self
                    .prog
                    .records
                    .iter()
                    .position(|t| Some(t.name.as_str()) == r.ty.as_deref())
                else {
                    return Ok(value);
                };
                let ty = TypeSpec::Record(id as u32);
                let value = self.coerce(value, &ty, span)?;
                if !self.conforms(&value, &ty) {
                    return Err(RuntimeError::new(
                        span,
                        format!("updated record doesn't conform to {}", ty.describe(self.prog)),
                    ));
                }
                Ok(value)
            }
            Value::Dist(d) => {
                let mut outcomes = Vec::with_capacity(d.outcomes.len());
                for (v, p) in &d.outcomes {
                    outcomes.push((self.check_updated_records(v.clone(), span)?, *p));
                }
                ops::combine(outcomes, d.missing, &mut self.budget).map_err(|e| e.at(span))
            }
            v => Ok(v),
        }
    }

    fn eval(&mut self, f: FnId, e: &Expr, w: &World) -> Result<Value> {
        let span = e.span;
        let at = |err: OpError| err.at(span);
        match &e.kind {
            ExprKind::Lit(l) => self.literal(l).map_err(at),
            ExprKind::Slot(s) => self.slot(f, *s, w, span),
            ExprKind::Unary(op, x) => {
                let v = self.eval(f, x, w)?;
                ops::unary(*op, &v, &mut self.budget).map_err(at)
            }
            ExprKind::Binary(op @ (BinOp::And | BinOp::Or), a, b) => {
                let and = *op == BinOp::And;
                let word = if and { "and" } else { "or" };
                let va = self.eval(f, a, w)?;
                let ta = ops::truth(&va, word).map_err(|err| err.at(a.span))?;
                if let Truth::Fact(x) = ta {
                    if x != and {
                        // `false and …` is false; `true or …` is true.
                        return Ok(Value::Bool(x));
                    }
                }
                let vb = self.eval(f, b, w)?;
                let tb = ops::truth(&vb, word).map_err(|err| err.at(b.span))?;
                ops::logic(and, ta, tb, &mut self.budget).map_err(at)
            }
            ExprKind::Binary(op, a, b) => {
                let va = self.eval(f, a, w)?;
                let vb = self.eval(f, b, w)?;
                ops::binary(*op, &va, &vb, &mut self.budget).map_err(at)
            }
            ExprKind::List(items) => {
                self.budget.collection(items.len() as u128).map_err(at)?;
                let mut values = Vec::with_capacity(items.len());
                for item in items {
                    values.push(self.eval(f, item, w)?);
                }
                Ok(Value::list(values))
            }
            ExprKind::Map(entries) => {
                let mut map = BTreeMap::new();
                for (k, v) in entries {
                    let key = self.eval(f, k, w)?;
                    if analytic::contains(&key) {
                        return Err(analytic::unsupported("analytic map keys").at(k.span));
                    }
                    if key.is_uncertain() {
                        return Err(RuntimeError::new(k.span, "map keys can't be distributions"));
                    }
                    let value = self.eval(f, v, w)?;
                    map.insert(key, value);
                }
                Ok(Value::map(map))
            }
            ExprKind::Record { ty, fields } => {
                let mut values = Vec::with_capacity(fields.len());
                for (name, v) in fields {
                    let field_ty = ty
                        .and_then(|t| self.prog.records[t as usize].fields.iter().find(|d| d.name == *name))
                        .map(|d| d.ty.clone());
                    let value = match field_ty {
                        Some(t) => {
                            let value = self.eval(f, v, w)?;
                            let value = self.coerce(value, &t, v.span)?;
                            if !self.conforms(&value, &t) {
                                return Err(RuntimeError::new(
                                    v.span,
                                    format!(
                                        "field `{name}` should be {}, found {}",
                                        t.describe(self.prog),
                                        value.kind()
                                    ),
                                ));
                            }
                            value
                        }
                        None => self.eval(f, v, w)?,
                    };
                    values.push((Arc::from(name.as_str()), value));
                }
                Ok(ops::make_record(
                    ty.map(|t| self.record_names[t as usize].clone()),
                    values,
                ))
            }
            ExprKind::Field(x, name) => {
                let v = self.eval(f, x, w)?;
                ops::field(&v, name, &mut self.budget).map_err(at)
            }
            ExprKind::Index(x, i) => {
                let v = self.eval(f, x, w)?;
                let i = self.eval(f, i, w)?;
                ops::index(&v, &i, &mut self.budget).map_err(at)
            }
            ExprKind::With(x, fields) => {
                let base = self.eval(f, x, w)?;
                let mut updates = Vec::with_capacity(fields.len());
                for (name, v) in fields {
                    let field_ty = self.record_field_type(&base, name);
                    let value = match field_ty {
                        Some(ty) => {
                            let value = self.eval_expected(f, v, w, &ty)?;
                            if !self.conforms(&value, &ty) {
                                return Err(RuntimeError::new(
                                    v.span,
                                    format!(
                                        "field `{name}` should be {}, found {}",
                                        ty.describe(self.prog),
                                        value.kind()
                                    ),
                                ));
                            }
                            value
                        }
                        None => self.eval(f, v, w)?,
                    };
                    updates.push((Arc::from(name.as_str()), value));
                }
                let value = ops::lift1(&base, &mut self.budget, |b, _| ops::with_fields(b, &updates)).map_err(at)?;
                self.check_updated_records(value, span)
            }
            ExprKind::Builtin { func, args, named } => self.builtin(f, *func, args, named, w, span),
            ExprKind::Closure { func, capture_args } => {
                // A named function reads its globals when it's called, so a
                // reference to it doesn't need them yet.
                let named = self.prog.functions[*func as usize].kind == FnKind::Named;
                let mut captured = Vec::with_capacity(capture_args.len());
                self.capture(&mut captured, f, *func, capture_args, w, span, named)?;
                Ok(Value::Closure(Arc::new(Closure {
                    func: *func,
                    captured,
                    named,
                })))
            }
            ExprKind::Bind { value, globals } => {
                let v = self.eval(f, value, w)?;
                let Value::Closure(c) = &v else {
                    return Ok(v);
                };
                if !c.named {
                    return Ok(v);
                }
                // A named function reads its globals now, as a call by name
                // does.
                let prog = self.prog;
                let captures = &prog.functions[c.func as usize].captures;
                let mut captured = Vec::with_capacity(captures.len());
                for capture in captures {
                    let CaptureSource::Global(g) = capture.source else {
                        return Err(RuntimeError::new(
                            span,
                            "internal error: a named function captures a local",
                        ));
                    };
                    let Ok(i) = globals.binary_search_by_key(&g, |&(g, _)| g) else {
                        return Err(RuntimeError::new(
                            span,
                            "internal error: a function value's globals aren't available here",
                        ));
                    };
                    captured.push(self.slot(f, globals[i].1, w, span)?);
                }
                Ok(Value::Closure(Arc::new(Closure {
                    func: c.func,
                    captured,
                    named: true,
                })))
            }
            ExprKind::Simulate { func, capture_args } => {
                let mut key = Vec::with_capacity(capture_args.len());
                self.capture(&mut key, f, *func, capture_args, w, span, false)?;
                self.simulate(*func, key, span)
            }
            ExprKind::Interp(parts) => {
                let mut text = String::new();
                for part in parts {
                    match part {
                        InterpPart::Lit(s) => crate::text::push(&mut text, s, &mut self.budget).map_err(at)?,
                        InterpPart::Expr(e) => {
                            let value = self.eval(f, e, w)?;
                            crate::text::push_value(&mut text, &value, &mut self.budget).map_err(at)?;
                        }
                    }
                }
                Ok(Value::str(&text))
            }
            ExprKind::Input(i) => Ok(self.inputs[*i as usize].clone()),
        }
    }

    fn literal(&mut self, l: &Lit) -> OpResult<Value> {
        Ok(match l {
            Lit::Builtin(b) => Value::Builtin(*b),
            Lit::Unit => Value::Unit,
            Lit::Bool(b) => Value::Bool(*b),
            Lit::Int(i) => {
                self.budget.integer_bits(i.bits())?;
                self.budget.work(i.bits().div_ceil(64).max(1))?;
                Value::Int(i.clone())
            }
            Lit::Float(x) | Lit::FloatConstant(x) => Value::Float(*x),
            Lit::Prob(p) => Value::Prob(*p),
            Lit::Str(s) => crate::text::value(s, &mut self.budget)?,
            Lit::Dice { count, sides } => {
                if let Some(d) = self.dice.get(&(*count, *sides)) {
                    return Ok(d.clone());
                }
                let d = Dist::dice(*count, *sides, &mut self.budget)?.into_value();
                self.dice.insert((*count, *sides), d.clone());
                d
            }
            Lit::Enum { ty, variant } => self.enum_values[*ty as usize][*variant as usize].clone(),
        })
    }

    fn builtin(
        &mut self,
        f: FnId,
        b: Builtin,
        args: &[Expr],
        named: &[(String, Expr)],
        w: &World,
        span: Span,
    ) -> Result<Value> {
        if b == Builtin::Typeof {
            // A direct inspection knows a delayed variable's outcome type
            // without drawing it and losing subsequent conjugate updates.
            let v = match args[0].kind {
                ExprKind::Slot(s) if matches!(w.slots[s as usize], Value::Delayed(_)) => w.slots[s as usize].clone(),
                _ => self.eval(f, &args[0], w)?,
            };
            return crate::type_name::of(&v, &self.prog.enums, &mut self.budget).map_err(|e| e.at(span));
        }
        let mut values = Vec::with_capacity(args.len());
        for a in args {
            values.push(self.eval(f, a, w)?);
        }
        let mut names = Vec::with_capacity(named.len());
        for (name, value) in named {
            names.push(name.clone());
            values.push(self.eval(f, value, w)?);
        }
        self.builtin_values(b, &values, &names, w.weight, span)
    }

    fn check_function_arity(&self, value: &Value, given: usize, span: Span) -> Result<()> {
        match value {
            Value::Closure(c) => self.check_arity(c, given, span),
            Value::Builtin(b) => {
                let (min, max) = b.arity();
                if given < min || given > max || (*b == Builtin::Date && given == 2) {
                    return Err(RuntimeError::new(
                        span,
                        format!(
                            "`{}` takes {}, but {given} were given",
                            b.name(),
                            if *b == Builtin::Date {
                                "one ISO string or three integers".into()
                            } else if min == max {
                                format!("{min} argument(s)")
                            } else if max == usize::MAX {
                                format!("at least {min} arguments")
                            } else {
                                format!("{min} to {max} arguments")
                            }
                        ),
                    ));
                }
                Ok(())
            }
            _ => Err(RuntimeError::new(
                span,
                "expected a comparator function, like `(a, b) -> a - b`",
            )),
        }
    }

    /// Shared dispatch for direct calls and first-class builtin values.
    fn builtin_values(
        &mut self,
        b: Builtin,
        values: &[Value],
        named: &[String],
        weight: Weight,
        span: Span,
    ) -> Result<Value> {
        let extrema = matches!(b, Builtin::Minimum | Builtin::Maximum);
        if !named.is_empty() && (!extrema || named != ["default"]) {
            return Err(RuntimeError::new(
                span,
                format!(
                    "`{}` does not accept these named arguments; only minimum/maximum support `default`",
                    b.name()
                ),
            ));
        }
        let positional = values.len() - named.len();
        if !named.is_empty() && positional == 0 {
            return Err(RuntimeError::new(
                span,
                "a collection argument is required before `default`",
            ));
        }
        if !named.is_empty() && positional >= 3 {
            return Err(RuntimeError::new(
                span,
                "the default was supplied both positionally and by name",
            ));
        }
        self.check_function_arity(&Value::Builtin(b), values.len(), span)?;
        let (values, default) = if extrema && (!named.is_empty() || values.len() == 3) {
            (&values[..values.len() - 1], values.last())
        } else {
            (values, None)
        };
        let at = |err: OpError| err.at(span);
        builtins::check_query_input(b, values).map_err(at)?;
        // Asking the built-in first spares a search of every collection
        // argument, on every call, for those that move values unseen.
        if !handles_analytic(b, values) && values.iter().any(analytic::contains) {
            // Internal built-ins stand for the syntax that uses them.
            let name = match b {
                Builtin::ScoreLaw => "score",
                Builtin::RepeatCount => "repeat",
                Builtin::Typeof => "typeof",
                b => b.name(),
            };
            return Err(analytic::unsupported(&format!("`{name}` on this outcome")).at(span));
        }
        match b {
            Builtin::Typeof => crate::type_name::of(&values[0], &self.prog.enums, &mut self.budget).map_err(at),
            Builtin::RunDate => self.config.today.map(Value::Date).ok_or_else(|| {
                RuntimeError::new(span, "the host did not supply an execution date for `today`")
                    .with_help("set Options.today, or supply today in the WASM request")
            }),
            Builtin::Print => {
                let mut text = String::new();
                for (i, value) in values.iter().enumerate() {
                    if i != 0 {
                        crate::text::push(&mut text, " ", &mut self.budget).map_err(at)?;
                    }
                    crate::text::push_value(&mut text, value, &mut self.budget).map_err(at)?;
                }
                let line = if weight == Weight::ONE {
                    text
                } else {
                    format!("[{}] {text}", fmt_weight(weight))
                };
                self.printed += line.len() + 1;
                if self.printed > self.config.max_output {
                    return Err(too_much_output(span));
                }
                match &mut self.lines {
                    Some(lines) => lines.push((span, line)),
                    None => (self.print)(&line),
                }
                Ok(Value::Unit)
            }
            Builtin::Map | Builtin::Filter | Builtin::Reduce => self.higher_order(b, values, span),
            Builtin::Count if values.len() == 2 => self.higher_order(b, values, span),
            Builtin::Minimum | Builtin::Maximum => {
                if values.len() == 2 {
                    self.collection_extreme(b, values, default, span)
                } else {
                    builtins::population_extreme(&values[0], b == Builtin::Maximum, default, &mut self.budget)
                        .map_err(at)
                }
            }
            Builtin::Sort | Builtin::SortDesc if values.len() == 2 => self.higher_order(b, values, span),
            Builtin::Highest | Builtin::Lowest if values.len() == 3 => self.higher_order(b, values, span),
            Builtin::Min | Builtin::Max => builtins::call_plain(b, values, &mut self.budget).map_err(at),
            Builtin::Roll => self.roll(values).map_err(at),
            Builtin::Take => Err(RuntimeError::new(
                span,
                "use `deck.take()` to draw and remove an item from a mutable bag",
            )),
            _ if b.lifting() == Lifting::Raw => builtins::call_raw(b, values, &mut self.budget).map_err(at),
            _ if values.iter().any(|v| matches!(v, Value::Continuous(_))) => {
                let kind = values
                    .iter()
                    .find(|v| matches!(v, Value::Continuous(_)))
                    .unwrap()
                    .kind();
                Err(
                    RuntimeError::new(span, format!("`{}` needs a value, not a {kind}", b.name()))
                        .with_help("draw a value first, like `let x ~ normal(0, 1)`"),
                )
            }
            _ => ops::lift_n(values, &mut self.budget, &|a, budget| {
                builtins::call_plain(b, a, budget)
            })
            .map_err(at),
        }
    }

    fn roll(&mut self, values: &[Value]) -> OpResult<Value> {
        if values[0].is_uncertain() {
            return Err(OpError::new("the number of dice to roll must be a plain number")
                .help("draw it first, like `let n ~ d4`"));
        }
        let count = ops::integer(&values[0], "roll count", &mut self.budget)?;
        let count = count
            .to_u64()
            .filter(|n| *n <= 1000)
            .ok_or_else(|| OpError::new("roll needs between 0 and 1000 dice"))? as u32;
        let die = match &values[1] {
            v @ (Value::Int(_) | Value::Float(_)) => {
                let sides = ops::integer(v, "roll sides", &mut self.budget)?;
                let sides = sides
                    .to_u64()
                    .filter(|n| *n >= 1 && *n <= u32::MAX as u64)
                    .ok_or_else(|| OpError::new("roll needs between 1 and 4294967295 sides"))?;
                Dist::dice(1, sides as u32, &mut self.budget)?
            }
            Value::Dist(d) => (**d).clone(),
            v => {
                return Err(OpError::new(format!(
                    "roll needs a die, like `d6`, found {}",
                    ops::article(&v.kind())
                )));
            }
        };
        let key = (count, die.clone().into_value());
        if let Some(pool) = self.pools.get(&key) {
            return Ok(pool.clone());
        }
        let pool = Dist::pool(count, &die, &mut self.budget)?.into_value();
        self.pools.insert(key, pool.clone());
        Ok(pool)
    }

    fn compare_callback(
        &mut self,
        comparator: &Value,
        a: &Value,
        other: &Value,
        b: Builtin,
        span: Span,
    ) -> Result<std::cmp::Ordering> {
        self.budget.work(1).map_err(|e| e.at(span))?;
        match self.call_pure(comparator, vec![a.clone(), other.clone()], b.name(), span)? {
            Value::Int(n) => Ok(n.cmp(&probl_number::Integer::ZERO)),
            Value::Float(x) if x.is_finite() => Ok(x.partial_cmp(&0.0).expect("finite comparator result")),
            other => Err(RuntimeError::new(
                span,
                format!(
                    "the comparator given to `{}` must return a finite int or float (negative, zero, or positive), found {}",
                    b.name(),
                    ops::article(&other.kind())
                ),
            )),
        }
    }

    fn collection_extreme(
        &mut self,
        b: Builtin,
        values: &[Value],
        default: Option<&Value>,
        span: Span,
    ) -> Result<Value> {
        if matches!(values[0], Value::Dist(_) | Value::Continuous(_)) {
            return Err(RuntimeError::new(
                span,
                "a comparator is supported only for collection extrema, not distribution support bounds",
            ));
        }
        let comparator = &values[1];
        self.check_function_arity(comparator, 2, span)?;
        let items = builtins::items(&values[0], b.name(), &mut self.budget).map_err(|e| e.at(span))?;
        let mut iter = items.into_iter();
        let Some(mut best) = iter.next() else {
            return builtins::empty_extreme(b.name(), default).map_err(|e| e.at(span));
        };
        for x in iter {
            let order = self.compare_callback(comparator, &x, &best, b, span)?;
            if (b == Builtin::Maximum && order.is_gt()) || (b == Builtin::Minimum && order.is_lt()) {
                best = x;
            }
        }
        Ok(best)
    }

    /// Collection operations which call a function, including custom sorting.
    fn higher_order(&mut self, b: Builtin, values: &[Value], span: Span) -> Result<Value> {
        if let Value::Dist(d) = &values[0] {
            let mut results = Vec::with_capacity(d.outcomes.len());
            for (coll, p) in &d.outcomes {
                let mut args = values.to_vec();
                args[0] = coll.clone();
                results.push((self.higher_order(b, &args, span)?, *p));
            }
            return ops::combine(results, d.missing, &mut self.budget).map_err(|e| e.at(span));
        }
        if matches!(b, Builtin::Highest | Builtin::Lowest) {
            if let Value::Dist(d) = &values[1] {
                let mut results = Vec::with_capacity(d.outcomes.len());
                for (count, p) in &d.outcomes {
                    let mut args = values.to_vec();
                    args[1] = count.clone();
                    results.push((self.higher_order(b, &args, span)?, *p));
                }
                return ops::combine(results, d.missing, &mut self.budget).map_err(|e| e.at(span));
            }
        }
        let items = builtins::items(&values[0], b.name(), &mut self.budget).map_err(|e| e.at(span))?;
        match b {
            Builtin::Sort | Builtin::SortDesc | Builtin::Highest | Builtin::Lowest => {
                let comparator = values.last().unwrap();
                // Validate even empty/singleton collections without invoking the callback.
                self.check_function_arity(comparator, 2, span)?;
                let count = if matches!(b, Builtin::Highest | Builtin::Lowest) {
                    Some(builtins::extreme_count(&values[1], &mut self.budget).map_err(|e| e.at(span))?)
                } else {
                    None
                };
                crate::ordering::reserve_sort(items.len(), &mut self.budget).map_err(|e| e.at(span))?;
                let mut items = items;
                crate::ordering::try_sort_by(&mut items, |a, other| {
                    let order = self.compare_callback(comparator, a, other, b, span)?;
                    Ok(if matches!(b, Builtin::SortDesc | Builtin::Highest) {
                        order.reverse()
                    } else {
                        order
                    })
                })?;
                if let Some(n) = count {
                    items.truncate(n);
                }
                Ok(Value::list(items))
            }
            Builtin::Map => {
                let mut out = Vec::with_capacity(items.len());
                for x in items {
                    out.push(self.call_pure(&values[1], vec![x], "map", span)?);
                }
                Ok(Value::list(out))
            }
            Builtin::Filter | Builtin::Count => {
                let mut kept = Vec::new();
                for x in items {
                    let test = self.call_pure(&values[1], vec![x.clone()], b.name(), span)?;
                    match test {
                        Value::Bool(true) => kept.push(x),
                        Value::Bool(false) => {}
                        other => {
                            return Err(RuntimeError::new(
                                span,
                                format!("the test given to `{}` must give a fact (true or false)", b.name()),
                            )
                            .with_help(format!("it gave {}", ops::article(&other.kind()))));
                        }
                    }
                }
                Ok(if b == Builtin::Count {
                    Value::Int((kept.len() as i64).into())
                } else {
                    Value::list(kept)
                })
            }
            Builtin::Reduce => {
                let callback = &values[1];
                if !matches!(callback, Value::Closure(_) | Value::Builtin(_)) {
                    return Err(
                        RuntimeError::new(span, "`reduce` needs a function as its second argument")
                            .with_help("write `reduce(xs, f)` or `reduce(xs, f, initial)`"),
                    );
                }
                self.check_function_arity(callback, 2, span)?;
                let mut items = items.into_iter();
                let mut acc = match values.get(2) {
                    Some(initial) => initial.clone(),
                    None => items.next().ok_or_else(|| {
                        RuntimeError::new(span, "`reduce` of an empty collection needs an initial value")
                            .with_help("supply an initial value: `reduce(xs, f, initial)`")
                    })?,
                };
                for x in items {
                    acc = self.call_pure(callback, vec![acc, x], "reduce", span)?;
                }
                Ok(acc)
            }
            _ => unreachable!(),
        }
    }
}

/// A call as written: `f(3, true)`.
/// The most rewards a loop's chain follows.
const MAX_REWARDS: usize = 4;

/// A loop's rewards, as whole numbers or floats by their bits.
type Reward = [u64; MAX_REWARDS];

/// Rewards as a world has them: their bits, and which are floats.
type Raw = (Reward, u8);

/// A loop's chain with rewards, as `solve_loop` found it: for each state,
/// its steps with how much the rewards grow, and what leaves the loop from
/// it, in one run of the body from the state with its rewards at 0.
struct Rewarded<'a> {
    rewards: &'a [usize],
    moves: &'a [Vec<(u32, Raw, f64)>],
    exits: Vec<Vec<World>>,
    returns: Vec<Vec<Returned>>,
    unresolved: Vec<Weight>,
    failed: Vec<Failures>,
    faulted: Vec<Vec<(World, RuntimeError)>>,
}

fn describe_call(fun: &Function, key: &[Value]) -> String {
    let args: Vec<String> = key[..fun.n_params as usize].iter().map(|v| format!("{v:?}")).collect();
    format!("{}({})", fun.name, args.join(", "))
}

/// The error when `print` at `span` goes over the output limit.
/// A world faults to a catch only inside a `try` of its own call, which takes
/// it: none can leave the call that way.
fn uncaught(flow: &Flow) -> Result<()> {
    match flow.faulted.first() {
        Some((_, e)) => Err(OpError::internal("a fault went past the `try` meant to catch it").at(e.span)),
        None => Ok(()),
    }
}

/// Remove the worlds at the indices `failed`, which are in increasing order.
fn drop_failed(worlds: &mut Vec<World>, failed: &[usize]) {
    if failed.is_empty() {
        return;
    }
    let mut failed = failed.iter().copied().peekable();
    let mut i = 0;
    worlds.retain(|_| {
        let drop = failed.next_if_eq(&i).is_some();
        i += 1;
        !drop
    });
}

pub fn too_much_output(span: Span) -> RuntimeError {
    RuntimeError::limit(span, "the program printed more than the output limit")
}

/// A weight as a percentage, even when it's too small for an `f64`.
pub fn fmt_weight(w: Weight) -> String {
    let x = w.to_f64();
    if x == 0.0 && !w.is_zero() {
        let l = w.log10() + 2.0;
        return format!("{:.1}e{}%", libm::pow(10.0, l - l.floor()), l.floor());
    }
    fmt_prob(x)
}

enum PathKey<'a> {
    Field(&'a str),
    Index(Value),
}

fn update(target: &mut Value, keys: &[PathKey], v: Value, budget: &mut Budget) -> OpResult<()> {
    let Some((first, rest)) = keys.split_first() else {
        *target = v;
        return Ok(());
    };
    match (first, target) {
        (PathKey::Field(name), Value::Record(r)) => {
            let r = Arc::make_mut(r);
            let slot = r
                .get_mut(name)
                .ok_or_else(|| OpError::new(format!("this record has no field `{name}`")))?;
            update(slot, rest, v, budget)
        }
        (PathKey::Index(i), Value::List(items)) => {
            let items = Arc::make_mut(items);
            let k = ops::as_index(i, items.len() as u128, budget)? as usize;
            update(&mut items[k], rest, v, budget)
        }
        (PathKey::Index(k), Value::Map(m)) => {
            let m = Arc::make_mut(m);
            if rest.is_empty() {
                m.insert(k.clone(), v);
                return Ok(());
            }
            let slot = m
                .get_mut(k)
                .ok_or_else(|| OpError::new(format!("the key {k:?} isn't in the map")))?;
            update(slot, rest, v, budget)
        }
        (PathKey::Field(name), other) => Err(OpError::new(format!(
            "can't set the field `{name}` of {}",
            ops::article(&other.kind())
        ))),
        (PathKey::Index(_), other) => Err(OpError::new(format!(
            "can't set an element of {}",
            ops::article(&other.kind())
        ))),
    }
}

/// The value at `keys` in `target`, to change in place: each level is
/// copied only if something else holds it.
fn place_mut<'v>(target: &'v mut Value, keys: &[PathKey], budget: &mut Budget) -> OpResult<&'v mut Value> {
    let Some((first, rest)) = keys.split_first() else {
        return Ok(target);
    };
    match (first, target) {
        (PathKey::Field(name), Value::Record(r)) => {
            let slot = Arc::make_mut(r)
                .get_mut(name)
                .ok_or_else(|| OpError::new(format!("this record has no field `{name}`")))?;
            place_mut(slot, rest, budget)
        }
        (PathKey::Index(i), Value::List(items)) => {
            let items = Arc::make_mut(items);
            let k = ops::as_index(i, items.len() as u128, budget)? as usize;
            place_mut(&mut items[k], rest, budget)
        }
        (PathKey::Index(k), Value::Map(m)) => {
            let slot = Arc::make_mut(m)
                .get_mut(k)
                .ok_or_else(|| OpError::new(format!("the key {k:?} isn't in the map")))?;
            place_mut(slot, rest, budget)
        }
        (PathKey::Field(name), other) => Err(OpError::new(format!(
            "can't set the field `{name}` of {}",
            ops::article(&other.kind())
        ))),
        (PathKey::Index(_), other) => Err(OpError::new(format!(
            "can't set an element of {}",
            ops::article(&other.kind())
        ))),
    }
}

/// Whether observing a value from `d` uses a density.
fn is_density(d: &Value) -> bool {
    match d {
        Value::Continuous(_) => true,
        Value::Dist(dist) => dist.outcomes.iter().any(|(x, _)| matches!(x, Value::Continuous(_))),
        _ => false,
    }
}

/// Whether a value is a discrete outcome of a draw, like `floor(x)`. Out of
/// line: assignments run it on every analytic value they store.
#[cold]
#[inline(never)]
fn is_discrete(v: &Value) -> bool {
    matches!(v, Value::Analytic(a) if a.is_discrete())
}

/// Whether a built-in deals with these analytic arguments deliberately.
/// The others would look inside them as if they were plain values.
fn handles_analytic(b: Builtin, values: &[Value]) -> bool {
    use Builtin as B;
    let plain = |v: &Value| !analytic::contains(v);
    let number = |v: &Value| matches!(v, Value::Analytic(_) | Value::Int(_) | Value::Float(_) | Value::Prob(_));
    let numbers = |v: &Value| matches!(v, Value::List(xs) if xs.iter().all(number));
    match b {
        B::Map | B::Filter | B::Reduce | B::Len | B::IterItems | B::Settled | B::BooleanLaw | B::IsListOfLen => true,
        // They move values without looking at them.
        B::Slice
        | B::Reverse
        | B::Keys
        | B::Values
        | B::Enumerate
        | B::Zip
        | B::Push
        | B::Pop
        | B::Last
        | B::DropLast => true,
        B::Count => values.len() == 1,
        // A key is compared with others, so it can't be an outcome; an index
        // is a number, which an outcome isn't.
        B::Get | B::Insert | B::Remove => matches!(values[0], Value::List(_)) || plain(&values[1]),
        // Affine in one draw, piecewise affine.
        B::Sum | B::Mean => numbers(&values[0]),
        B::Floor | B::Ceil | B::Trunc => matches!(values[0], Value::Analytic(_)),
        B::Sqrt | B::Exp | B::Ln | B::Exp2 | B::Log10 | B::Log2 | B::Log1p | B::Expm1 => {
            matches!(values[0], Value::Analytic(_))
        }
        B::Round => matches!(values[0], Value::Analytic(_)) && values[1..].iter().all(plain),
        B::Abs | B::Min | B::Max | B::Clamp => values
            .iter()
            .all(|v| number(v) || matches!(v, Value::Dist(_)) && plain(v)),
        B::Minimum | B::Maximum => values.len() == 1 && numbers(&values[0]),
        _ => false,
    }
}

/// What observing `v` from `d` weighs.
enum Weighs {
    /// The probability of `v`, the probability that is missing from `d` (so
    /// the true value may be up to that much higher), and the probability
    /// of the other outcomes, which the observation rules out. When
    /// sampling, a continuous `d` gives its density here (section 13).
    Probability(f64, f64, f64),
    /// When enumerating, the logarithm of a continuous `d`'s density at `v`.
    LogDensity(f64),
}

fn likelihood(d: &Value, v: &Value, sampling: bool) -> OpResult<Weighs> {
    let continuous = match d {
        Value::Continuous(f) => Some(vec![(**f, 1.0)]),
        Value::Dist(dist) if dist.outcomes.iter().any(|(x, _)| matches!(x, Value::Continuous(_))) => {
            let parts: Option<Vec<_>> = dist
                .outcomes
                .iter()
                .map(|(x, p)| match x {
                    Value::Continuous(f) => Some((**f, *p)),
                    _ => None,
                })
                .collect();
            let parts = parts.ok_or_else(|| {
                OpError::new("`observe … from` can't mix a density with the probabilities of single values")
            })?;
            Some(parts)
        }
        _ => None,
    };
    if let Some(parts) = continuous {
        if let Value::Dist(dist) = d {
            if dist.missing > 0.0 && !sampling {
                return Err(OpError::unsupported(
                    "observing a value from a mixture with unresolved probability isn't supported when enumerating",
                )
                .help("the unresolved part could have any density at the value; sample the model with `@mode sample(runs: 10_000)`"));
            }
        }
        let x = match v {
            Value::Bool(_) => None,
            _ => v.as_f64(),
        }
        .ok_or_else(|| {
            OpError::new(format!(
                "a continuous distribution can't produce {}",
                ops::article(&v.kind())
            ))
        })?;
        if !sampling {
            return Ok(Weighs::LogDensity(crate::stats::log_sum_exp(
                parts.iter().map(|(f, p)| libm::log(*p) + f.ln_pdf(x)),
            )));
        }
        return Ok(Weighs::Probability(
            parts.iter().map(|(f, p)| p * f.pdf(x)).sum(),
            0.0,
            0.0,
        ));
    }
    match d {
        // As observing compares values: `2.0` is the outcome 2.
        Value::Counts(law) => {
            let x = match v.as_complex() {
                Some(z) if z.im() == 0.0 => z.re(),
                _ => f64::NAN,
            };
            let p = if x.is_nan() { 0.0 } else { law.pmf(x) };
            Ok(Weighs::Probability(p, 0.0, 1.0 - p))
        }
        Value::Dist(dist) => {
            let (mut seen, mut other) = (0.0, 0.0);
            for (x, p) in &dist.outcomes {
                if ops::equals(x, v) {
                    seen += p;
                } else {
                    other += p;
                }
            }
            Ok(Weighs::Probability(seen, dist.missing, other))
        }
        Value::Prob(_) => Err(OpError::new("`observe … from` needs a distribution, not a probability")
            .help("to observe that a fact with probability p is true, write `observe true from bernoulli(p)`")),
        other => Ok(if ops::equals(other, v) {
            Weighs::Probability(1.0, 0.0, 0.0)
        } else {
            Weighs::Probability(0.0, 0.0, 1.0)
        }),
    }
}