diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1343f66d..8da17c72 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -537,8 +537,20 @@ jobs: # memory safety is hardest to get right. Nightly-only: this shard is # already off the PR critical path, so the extra build cost is free # where it would not be in Tier 1a. + # + # RUST_MIN_STACK for the same reason the `tsan` shard sets it, and it is + # load-bearing here too: the 2026-08-14 nightly — the first ASan run to + # see the parallel code — died with + # `AddressSanitizer: stack-overflow ... T1601`, on a rayon worker rather + # than the main thread. It is not corruption. `simplify::dispatch`'s + # stack governor refills at 512 KiB, ASan's instrumented frames are far + # fatter than the uninstrumented ones that margin was tuned for, and + # rayon workers start from a 2 MiB default rather than the main thread's + # 8 MiB. Raising the worker stack fixes it deterministically and it does + # not reproduce uninstrumented. env: LSAN_OPTIONS: detect_leaks=0 + RUST_MIN_STACK: "33554432" run: | RUSTFLAGS="-Zsanitizer=address" \ cargo +nightly test --workspace --lib --tests \ diff --git a/CHANGELOG.md b/CHANGELOG.md index a107f8c0..f58c947f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -127,6 +127,55 @@ Both are detailed under *Behaviour changes to plan for*. this**: a root that never produces a sign change — a double root like `(x−1)²` on `[0, 2]`, or `(x²−1)²` on `[-2, 2]` — still answers `"undecided"`, because no witness pair exists and none is invented. +- **`verified_sign` could hang on a rational constant with enough digits.** + `sin(x)·D − N·x ≥ 0` on `[0, 3/2]` took 0.05 s at `N/D = 636/1000`, 0.08 s at + nine digits, and **over 300 s at twelve** — three extra digits turned + milliseconds into a hang. The cause was neither exact rational arithmetic nor + repeated conversion but a **non-terminating loop** in the branch-and-bound: + a sub-box bisected down to the width floor was pushed back onto the active + list, immediately re-selected as the smallest key with nothing changed, and + the loop spun *without ever consuming its subdivision budget* — which is why + capping `max_subdivisions` at 64 did not help. It only triggered once `tol` + became unreachable, and `tol` is an **absolute** width: at `D = 10¹²` the + function is of order `10¹¹`, so the default `1e-9` asks for twenty digits and + the floor arrives first. Nine digits happened to converge just above the + floor, twelve just below. Boxes that reach the floor are now retired out of + the active list, their keys still folded into the final bound, so the search + always makes progress. Cost is now flat in the size of the constant: 0.33 s at + three digits through 0.47 s at sixteen. +- **Inequalities that are tight at an endpoint no longer stay `"undecided"`.** + Cusa–Huygens, Mitrinović–Adamović, Wilker and Huygens were `"true"` on + `[0.1, 1.5]` but `"undecided"` on `[0.01, 1.5]` and at `x = 0` — precisely the + point that makes them worth stating. Two independent things were in the way. + First, `tol` was also the wrong *stopping rule* for a sign question: on + `[0.01, 1.5]` the true minimum of the Cusa–Huygens form is `1.7·10⁻¹³`, so the + search met the `1e-9` tolerance and stopped with an enclosure that still + straddled zero. `verified_sign` now re-runs the search with the sign itself as + the goal, refining while the running bound straddles zero instead of to an + absolute width. Second, where the margin genuinely *vanishes* no subdivision + can ever help, so the box is split: a collar `[a, a+δ]` at the endpoint is + handled by a truncated Taylor expansion there, and the rest by the usual + branch-and-bound. The two pieces are closed and share the join point `δ`, so + their union is the original box with no gap. All four are now `"true"` on + `[0, 1.5]`, as is Jordan's inequality stated exactly as + `10¹²·sin x − 636619772368·x ≥ 0`. + **The remainder is proven, not assumed.** Coefficients `c_k = g⁽ᵏ⁾(a)/k!` are + accepted as *zero* only when substitution followed by `simplify` lands on the + literal integer `0` — no numeric enclosure can prove a value is zero, and none + is asked to — cross-checked against ball arithmetic, and the tail is a + Lagrange remainder `|R(t)| ≤ t^m·sup|g⁽ᵐ⁾|/m!` whose sup is a rigorous + enclosure over the whole collar. Analyticity, which Taylor's theorem needs, is + certified by requiring every derivative `g … g⁽ᵐ⁾` to enclose successfully + there. With `c_0 … c_{j−1}` proven zero, `g(a+t) ≥ t^j·[c_j − T(δ)]` and + `t^j ≥ 0` finishes it. **Nothing was traded for the extra reach**: a margin + that vanishes in the *interior* — `(x − 7/10)²(x + 1)` on `[0, 3/2]` — is + still `"undecided"`, because the expansion does not apply there. A leading + coefficient proven *negative* now returns `"false"` rather than `"undecided"`, + which settles cases no sampling could see: `x³ − x²/1000` is negative only on + `(0, 1/1000)`, and each of the four inequalities reversed is refuted at the + same endpoint where the original is certified. A strict `"positive"` query is + `"false"` where the expression is proven to vanish exactly, so `x² > 0` on + `[0, 1]` is `"false"` while `x² ≥ 0` is `"true"`. - **`verified_integral` refused removable singularities.** Taylor-model quadrature raised `E-VALIDATED-003` on any sub-interval where the reciprocal's enclosure contained zero, which put `∫₀¹ ln(1+x)/x dx = π²/12` out of reach @@ -372,6 +421,58 @@ Both are detailed under *Behaviour changes to plan for*. ### Added +- **Validated-bounds coverage is queryable: `bounds_supported(expr)` and a + `taylor_model` bit in `capabilities()["primitives"]`.** The only + per-function coverage flag the agent contract exposed was `numeric_ball`, + and it is not the flag that governs `bound_on_box` / `verified_integral` / + `verified_no_roots` / `verified_sign`. Ball arithmetic is *pointwise*; a + Taylor model needs a rule with a rigorous Lagrange remainder, and ten + primitives have the first without the second — `erf`, `erfc`, `bessel_j0`, + `bessel_j1`, `digamma`, `lambert_w`, `acosh`, `asinh`, `floor`, `ceil`. So + `numeric_ball` said `True` for `bessel_j0` and every bound over a box died + on `E-VALIDATED-001`. The boundary was enforced correctly and could not be + found ahead of time, which is how a planning loop loses a whole designed + workload (Turán-type inequalities for Bessel functions, in the 2026-08-13 + autoresearch run) to a route it could have ruled out for free. + + `taylor_model` reports it per primitive — `True` for the elementary + fragment (`exp`, `log`, `sqrt`, `sin`, `cos`, `tan`, `asin`, `acos`, + `atan`, `sinh`, `cosh`, `tanh`, `abs`) and `False` for every special + function. `ak.bounds_supported(expr)` asks for a whole expression, without + running the bound: it is truthy when nothing in the expression will be + refused as unsupported, and carries `.blocker` (the evaluator's own + description of the first construct it has no rule for) and `.functions` + (every blocking function, so a substitution can be planned in one round + rather than found one at a time). + + **Neither is a maintained list.** Both are derived by running the real + Taylor evaluator on a probe expression and asking whether it refuses with + `E-VALIDATED-001` — a second hand-written table is how `numeric_ball` came + to be read as coverage in the first place, and would have been a worse + outcome than no flag at all. `tests/test_taylor_model_coverage.py` + re-derives the bit the only other way there is, by calling `bound_on_box` + on every registered primitive, and fails if the two ever disagree. + + `numeric_ball` itself is *accurate* and stays as it is: those ten + primitives really do have Arb ball arithmetic. It answers a different + question, and now says so next to a flag that answers this one. A `True` + from either means "not `E-VALIDATED-001`" — a covered function can still be + refused on a particular box for a domain violation (`E-VALIDATED-003`) or a + non-finite enclosure (`-004`), which no box-free predicate can rule out. + + This is deliberately *not* folded into `certifiable`, which asks whether an + operation emits a **Lean** certificate and answers from the certificate + ledger. A rigorous enclosure is not a Lean proof term and the validated + subsystem has no ledger rows; one predicate returning `True` for two kinds + of evidence would be a worse contract than two predicates. + + New in Rust: `alkahest_cas::{taylor_model_refusal, taylor_model_blockers, + taylor_model_supports, taylor_model_supports_call}` and + `Capabilities::TAYLOR_MODEL` (also a `taylor_model` column in + `CoverageReport::to_markdown`). `capabilities()["contract_version"]` stays + `3`: the row gained a key and lost none, which is the same additive rule + the `__all__` freeze check applies. + - **Gröbner results can be read back — `GbPoly.to_expr`, iteration over a `GroebnerBasis`, and `expr_to_gbpoly`.** Everything that returned a basis returned a handle nobody could open. `GbPoly` exposed only `is_zero` and diff --git a/alkahest-core/src/holonomic/zeilberger.rs b/alkahest-core/src/holonomic/zeilberger.rs index fa3ff2ed..c106c48d 100644 --- a/alkahest-core/src/holonomic/zeilberger.rs +++ b/alkahest-core/src/holonomic/zeilberger.rs @@ -806,11 +806,23 @@ mod tests { .expect("Franel must be decided at the default bounds"); let elapsed = start.elapsed(); println!("franel: order {} in {:?}", result.value.order, elapsed); - assert!( - elapsed < std::time::Duration::from_secs(10), - "Franel took {elapsed:?} at the default bounds — the exact Q(n)(k) \ - post-processing has regressed to the coefficient blowup it used to have" - ); + // The wall-clock guard is meaningful only in an uninstrumented release + // build. Under a sanitizer it is not: the nightly `lsan` shard runs a + // debug build with LeakSanitizer, where the whole lib suite takes ~23 + // minutes and this test breached a 10 s bound while the thing it guards + // — the Z[n][k] gcd — was perfectly healthy. A timing assertion that + // fails for the instrumentation rather than the regression is a flaky + // test, so it is skipped there; the correctness assertions below always + // run, in every configuration. + // `debug_assertions` is the discriminator: every sanitizer shard builds + // in debug, and the release run this guard is written for does not. + if !cfg!(debug_assertions) { + assert!( + elapsed < std::time::Duration::from_secs(10), + "Franel took {elapsed:?} at the default bounds — the exact Q(n)(k) \ + post-processing has regressed to the coefficient blowup it used to have" + ); + } let r = &result.value; assert_eq!(r.order, 2, "Σ_k C(n,k)³ satisfies an order-2 recurrence"); diff --git a/alkahest-core/src/lib.rs b/alkahest-core/src/lib.rs index 5155ed2a..a1082d08 100644 --- a/alkahest-core/src/lib.rs +++ b/alkahest-core/src/lib.rs @@ -222,7 +222,10 @@ pub use number_theory::{ discrete_log, factorint, isprime, jacobi_symbol, nextprime, nthroot_mod, totient, NumberTheoryError, QuadraticDirichlet, }; -pub use primitive::{Capabilities, CoverageReport, CoverageRow, Primitive, PrimitiveRegistry}; +pub use primitive::{ + taylor_model_blockers, taylor_model_refusal, taylor_model_supports, taylor_model_supports_call, + Capabilities, CoverageReport, CoverageRow, Primitive, PrimitiveRegistry, +}; #[cfg(feature = "groebner")] pub use solver::{ diophantine, expr_to_gbpoly, extract_regular_chain_from_basis, gbpoly_to_expr, @@ -310,7 +313,10 @@ pub mod stable { ResultantError, RootInterval, SparseGcdError, SparseInterpError, UniPoly, UniPolyFactorModP, UniPolyFactorization, }; - pub use crate::primitive::{Primitive, PrimitiveRegistry}; + pub use crate::primitive::{ + taylor_model_blockers, taylor_model_refusal, taylor_model_supports, + taylor_model_supports_call, Primitive, PrimitiveRegistry, + }; pub use crate::real::{ cad_lift, cad_project, decide, decide_expr, routh_hurwitz, CadError, QeResult, RouthHurwitz, }; diff --git a/alkahest-core/src/primitive/mod.rs b/alkahest-core/src/primitive/mod.rs index a64d82bd..bcc1329a 100644 --- a/alkahest-core/src/primitive/mod.rs +++ b/alkahest-core/src/primitive/mod.rs @@ -50,6 +50,12 @@ use crate::kernel::{ExprId, ExprPool}; use std::collections::HashMap; use std::fmt; +pub mod taylor_support; + +pub use taylor_support::{ + taylor_model_blockers, taylor_model_refusal, taylor_model_supports, taylor_model_supports_call, +}; + // --------------------------------------------------------------------------- // Capability flags // --------------------------------------------------------------------------- @@ -65,6 +71,18 @@ bitflags::bitflags! { const NUMERIC_BALL = 1 << 4; const LOWER_LLVM = 1 << 5; const LEAN_THEOREM = 1 << 6; + /// The validated-bounds subsystem ([`crate::validated`]) has a + /// rigorous Taylor-model rule for this primitive, so + /// `bound_on_box` / `verified_integral` / `verified_no_roots` / + /// `verified_sign` will not refuse it with `E-VALIDATED-001`. + /// + /// **This is not implied by `NUMERIC_BALL`, and does not imply it.** + /// Pointwise ball arithmetic and a Taylor model with a rigorous + /// remainder are different pieces of work: `erf`, `bessel_j0`, + /// `digamma`, `floor`, … have the former and not the latter. The bit + /// is derived by running the evaluator (see + /// [`taylor_support`]), never from a list. + const TAYLOR_MODEL = 1 << 7; } } @@ -78,6 +96,7 @@ impl fmt::Display for Capabilities { (Capabilities::NUMERIC_BALL, "numeric_ball"), (Capabilities::LOWER_LLVM, "lower_llvm"), (Capabilities::LEAN_THEOREM, "lean"), + (Capabilities::TAYLOR_MODEL, "taylor_model"), ]; let present: Vec<&str> = names .iter() @@ -182,8 +201,8 @@ pub struct CoverageReport { impl CoverageReport { /// Render as a Markdown table (suitable for CI PR comments or docs). pub fn to_markdown(&self) -> String { - let header = "| Primitive | simplify | diff_fwd | diff_rev | numeric_f64 | numeric_ball | lower_llvm | lean |\n\ - |---|---|---|---|---|---|---|---|"; + let header = "| Primitive | simplify | diff_fwd | diff_rev | numeric_f64 | numeric_ball | lower_llvm | lean | taylor_model |\n\ + |---|---|---|---|---|---|---|---|---|"; let rows: Vec = self .rows .iter() @@ -196,7 +215,7 @@ impl CoverageReport { } }; format!( - "| {} | {} | {} | {} | {} | {} | {} | {} |", + "| {} | {} | {} | {} | {} | {} | {} | {} | {} |", r.name, tick(Capabilities::SIMPLIFY), tick(Capabilities::DIFF_FORWARD), @@ -205,6 +224,7 @@ impl CoverageReport { tick(Capabilities::NUMERIC_BALL), tick(Capabilities::LOWER_LLVM), tick(Capabilities::LEAN_THEOREM), + tick(Capabilities::TAYLOR_MODEL), ) }) .collect(); @@ -261,7 +281,7 @@ impl PrimitiveRegistry { pub fn capabilities(&self, name: &str) -> Capabilities { self.map .get(name) - .map(|e| e.caps) + .map(|e| with_taylor_model(name, e.caps)) .unwrap_or(Capabilities::empty()) } @@ -273,7 +293,7 @@ impl PrimitiveRegistry { .iter() .map(|(name, e)| CoverageRow { name: name.to_string(), - caps: e.caps, + caps: with_taylor_model(name, e.caps), }) .collect(); rows.sort_by(|a, b| a.name.cmp(&b.name)); @@ -372,7 +392,9 @@ impl PrimitiveRegistry { /// Iterate over all registered (name, capabilities) pairs. pub fn iter(&self) -> impl Iterator { - self.map.iter().map(|(k, e)| (*k, e.caps)) + self.map + .iter() + .map(|(k, e)| (*k, with_taylor_model(k, e.caps))) } } @@ -453,9 +475,35 @@ fn probe_caps(p: &dyn Primitive) -> Capabilities { if p.lean_theorem().is_some() { caps |= Capabilities::LEAN_THEOREM; } + // NB: `TAYLOR_MODEL` is deliberately *not* probed here — see + // `with_taylor_model`. Probing it at registration cost every caller that + // builds a registry, which `default_registry()` does on hot paths. caps } +/// Add [`Capabilities::TAYLOR_MODEL`] to a primitive's stored capabilities. +/// +/// This bit is resolved when the capabilities are *read*, not when the +/// primitive is registered. It is not a slot on the `Primitive` trait: the +/// validated-bounds subsystem keeps its own per-function rules in +/// `validated::taylor`, and a primitive cannot self-report whether one exists +/// without that claim being able to drift — so the answer comes from asking +/// the evaluator (memoised, see `taylor_support`). +/// +/// Asking it at registration time made `PrimitiveRegistry::register` pay a +/// probe per primitive, and `default_registry()` is rebuilt on hot paths such +/// as `diff` and `series` — it cost `series(sin x, 12)` roughly 30% steady +/// state and ~4 ms on the first construction in a process. Reading is rare +/// (`capabilities()`, `bounds_supported`, the coverage report), so the cost +/// belongs here. +fn with_taylor_model(name: &str, caps: Capabilities) -> Capabilities { + if taylor_support::taylor_model_supports(name) { + caps | Capabilities::TAYLOR_MODEL + } else { + caps + } +} + // --------------------------------------------------------------------------- // Built-in primitives // --------------------------------------------------------------------------- diff --git a/alkahest-core/src/primitive/taylor_support.rs b/alkahest-core/src/primitive/taylor_support.rs new file mode 100644 index 00000000..fdfcdff7 --- /dev/null +++ b/alkahest-core/src/primitive/taylor_support.rs @@ -0,0 +1,364 @@ +//! Which constructs the validated-bounds subsystem can actually bound. +//! +//! [`crate::validated`] (`bound_on_box`, `verified_integral`, +//! `verified_no_roots`, `verified_sign`) is *not* driven by the +//! [`Capabilities::NUMERIC_BALL`](super::Capabilities::NUMERIC_BALL) bundle +//! slot. Pointwise ball arithmetic gives an enclosure of `f` at a ball; a +//! Taylor model additionally needs a polynomial expansion with a rigorous +//! Lagrange remainder, which is written per function in +//! [`crate::validated::taylor`]. The two sets differ: `erf`, `bessel_j0`, +//! `digamma`, `floor`, … all have real ball arithmetic (via Arb) and no +//! Taylor-model rule, so `bound_on_box` refuses them with +//! `E-VALIDATED-001`. +//! +//! Before this module the boundary was only discoverable by hitting it. The +//! flag exposed here closes that gap **without introducing a second list to +//! maintain**: every answer is produced by *running* the real evaluator on a +//! probe expression and looking at whether it refuses with +//! [`ValidatedError::Unsupported`]. There is nothing here to keep in sync — +//! adding a rule to `validated::taylor` flips the flag on the next call, and +//! removing one flips it off. +//! +//! # Example +//! +//! ``` +//! use alkahest_cas::kernel::{Domain, ExprPool}; +//! use alkahest_cas::primitive::{taylor_model_refusal, taylor_model_supports}; +//! +//! assert!(taylor_model_supports("sin")); +//! // Real ball arithmetic, no Taylor-model rule: +//! assert!(!taylor_model_supports("bessel_j0")); +//! +//! let pool = ExprPool::new(); +//! let x = pool.symbol("x", Domain::Real); +//! assert!(taylor_model_refusal(pool.func("sin", vec![x]), &pool).is_none()); +//! assert!(taylor_model_refusal(pool.func("erf", vec![x]), &pool).is_some()); +//! ``` + +use crate::kernel::{Domain, ExprData, ExprId, ExprPool}; +use crate::validated::taylor::taylor_range; +use crate::validated::ValidatedError; +use rug::Float; +use std::collections::{HashMap, HashSet}; +use std::sync::{OnceLock, RwLock, RwLockReadGuard, RwLockWriteGuard}; + +/// Probe box. Any non-degenerate box works: whether the evaluator has a +/// *rule* for a construct is a structural question, and the box only decides +/// whether that rule then hits a domain violation — a different error, and one +/// this module deliberately reports as "supported" (the rule exists; the +/// caller picked a bad box). +const PROBE_LO: f64 = 0.25; +const PROBE_HI: f64 = 0.5; +/// Order 1 / 64 bits: dispatch does not depend on either, so probe cheaply. +const PROBE_ORDER: usize = 1; +const PROBE_PREC: u32 = 64; +/// Arities probed by [`taylor_model_supports`]. Today the evaluator only has +/// unary function rules, but probing a range means a future binary rule +/// (`atan2`, `hypot`) is picked up with no edit here. +const MAX_PROBE_ARITY: usize = 3; + +/// The evaluator's own description of the first construct in `expr` that it +/// has no rigorous Taylor-model rule for, or `None` if it has a rule for +/// every one of them. +/// +/// `None` is exactly the condition "`bound_on_box` will not fail with +/// `E-VALIDATED-001`". It is **not** a promise that `bound_on_box` succeeds: +/// a supported function can still hit a pole, a branch cut or an overflow on +/// a particular box (`E-VALIDATED-003` / `E-VALIDATED-004`), which depends on +/// the box and not on the expression. +/// +/// Free symbols are enclosed in an arbitrary probe box, so the answer depends +/// only on the shape of `expr`. +pub fn taylor_model_refusal(expr: ExprId, pool: &ExprPool) -> Option { + let (symbols, _) = walk_expr(expr, pool); + let mut probe_box: Vec<(ExprId, Float, Float)> = symbols + .into_iter() + .map(|s| { + ( + s, + Float::with_val(PROBE_PREC, PROBE_LO), + Float::with_val(PROBE_PREC, PROBE_HI), + ) + }) + .collect(); + if probe_box.is_empty() { + // A constant expression still has to be *evaluated* to learn whether + // its functions are supported, and the evaluator rejects an empty + // box. Add `expr` itself as the box variable: only `Symbol` nodes are + // ever matched against the box, so a non-symbol entry is an unused + // extra dimension — and unlike interning a fresh symbol, it does not + // mutate the caller's pool. + probe_box.push(( + expr, + Float::with_val(PROBE_PREC, PROBE_LO), + Float::with_val(PROBE_PREC, PROBE_HI), + )); + } + match taylor_range(expr, pool, &probe_box, PROBE_ORDER, PROBE_PREC) { + Err(ValidatedError::Unsupported { what }) => Some(what), + _ => None, + } +} + +/// Every function *call* inside `expr` that the Taylor-model evaluator has no +/// rule for, by name, deduplicated and sorted. +/// +/// This is the actionable half of [`taylor_model_refusal`]: the refusal names +/// the first blocking construct, this names all the blocking functions at +/// once, so a caller can decide what to substitute. An expression can be +/// refused with an empty list here — the blocker may be a node kind rather +/// than a function (a symbolic exponent over a non-positive base, a +/// `Piecewise`, …); [`taylor_model_refusal`] stays authoritative. +pub fn taylor_model_blockers(expr: ExprId, pool: &ExprPool) -> Vec { + let (_, calls) = walk_expr(expr, pool); + let mut out: Vec = calls + .into_iter() + .filter(|(name, arity)| !taylor_model_supports_call(name, *arity)) + .map(|(name, _)| name) + .collect(); + out.sort(); + out.dedup(); + out +} + +/// Does the Taylor-model evaluator have a rule for `name` at *any* arity it +/// accepts? This is the per-primitive flag reported as `taylor_model` in +/// `capabilities()["primitives"]`. +pub fn taylor_model_supports(name: &str) -> bool { + (1..=MAX_PROBE_ARITY).any(|arity| taylor_model_supports_call(name, arity)) +} + +/// Does the Taylor-model evaluator have a rule for `name` applied to exactly +/// `arity` arguments? +/// +/// Arity matters: the evaluator's rules are unary today, so `atan2(x, y)` is +/// refused for a reason that has nothing to do with whether `atan2` could be +/// bounded in principle. +pub fn taylor_model_supports_call(name: &str, arity: usize) -> bool { + if arity == 0 || arity > MAX_PROBE_ARITY { + // The evaluator's `Func` arms all destructure at least one argument, + // and every arm above `MAX_PROBE_ARITY` is the catch-all refusal. + return false; + } + let cache = &cache()[arity - 1]; + if let Some(&hit) = read_lock(cache).get(name) { + return hit; + } + // Computed *outside* the lock: probing runs the validated evaluator, which + // must never be able to re-enter this cache while it is held. + let answer = probe_call(name, arity); + write_lock(cache).insert(name.to_string(), answer); + answer +} + +// --------------------------------------------------------------------------- +// Probing +// --------------------------------------------------------------------------- + +/// Build `name(x₁, …, x_arity)` in a scratch pool and ask the real evaluator. +fn probe_call(name: &str, arity: usize) -> bool { + let pool = ExprPool::new(); + let args: Vec = (0..arity) + .map(|i| pool.symbol(format!("__taylor_probe_{i}"), Domain::Real)) + .collect(); + let call = pool.func(name, args); + taylor_model_refusal(call, &pool).is_none() +} + +/// Collect the free symbols and the `(function name, arity)` calls in `expr`. +/// +/// Enumerating the nodes is not the same as knowing which of them the +/// evaluator supports — that question is only ever answered by running it. +fn walk_expr(expr: ExprId, pool: &ExprPool) -> (Vec, Vec<(String, usize)>) { + let mut seen: HashSet = HashSet::new(); + let mut symbols: Vec = Vec::new(); + let mut calls: Vec<(String, usize)> = Vec::new(); + let mut stack = vec![expr]; + while let Some(id) = stack.pop() { + if !seen.insert(id) { + continue; + } + let children: Vec = pool.with(id, |data| match data { + ExprData::Symbol { .. } => { + symbols.push(id); + vec![] + } + ExprData::Integer(_) | ExprData::Rational(_) | ExprData::Float(_) => vec![], + ExprData::Add(args) | ExprData::Mul(args) => args.clone(), + ExprData::Pow { base, exp } => vec![*base, *exp], + ExprData::Func { name, args } => { + calls.push((name.clone(), args.len())); + args.clone() + } + ExprData::Piecewise { branches, default } => { + let mut ids: Vec = branches.iter().flat_map(|(c, v)| [*c, *v]).collect(); + ids.push(*default); + ids + } + ExprData::Predicate { args, .. } => args.clone(), + ExprData::Forall { var, body } | ExprData::Exists { var, body } => vec![*var, *body], + ExprData::BigO(arg) => vec![*arg], + ExprData::RootSum { poly, var, body } => vec![*poly, *var, *body], + }); + stack.extend(children); + } + // Deterministic order: the walk visits through a stack, and a box whose + // dimension order depended on traversal luck would make the probe's + // (irrelevant, but observable) numerics irreproducible. + symbols.sort_unstable(); + symbols.dedup(); + (symbols, calls) +} + +// --------------------------------------------------------------------------- +// Cache +// --------------------------------------------------------------------------- + +/// Probing is pure and deterministic — the same name and arity always get the +/// same answer — but it runs a whole Taylor evaluation, and +/// `PrimitiveRegistry::default_registry()` is rebuilt on hot paths (every +/// `diff` of an unregistered `Func` node). So memoise it process-wide, one +/// map per arity so that a `&str` lookup does not have to allocate a key. +type ProbeCache = [RwLock>; MAX_PROBE_ARITY]; + +fn cache() -> &'static ProbeCache { + static CACHE: OnceLock = OnceLock::new(); + CACHE.get_or_init(|| std::array::from_fn(|_| RwLock::new(HashMap::new()))) +} + +/// A poisoned probe cache is not a correctness problem: entries are +/// deterministic and inserted one at a time, so recovering the map is always +/// better than propagating a panic from an unrelated thread. +fn read_lock(c: &RwLock>) -> RwLockReadGuard<'_, HashMap> { + c.read().unwrap_or_else(|e| e.into_inner()) +} + +fn write_lock(c: &RwLock>) -> RwLockWriteGuard<'_, HashMap> { + c.write().unwrap_or_else(|e| e.into_inner()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::validated::bounds::{bound_on_box, BoundOptions}; + + /// The whole point of deriving the flag: it agrees with the subsystem it + /// describes, for every registered primitive, without a list in between. + #[test] + fn flag_matches_bound_on_box_for_every_primitive() { + let reg = crate::primitive::PrimitiveRegistry::default_registry(); + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let opts = BoundOptions { + order: 2, + prec: 64, + tol: 1e-3, + max_subdivisions: 8, + }; + for (name, caps) in reg.iter() { + let flag = caps.contains(crate::primitive::Capabilities::TAYLOR_MODEL); + let call = pool.func(name, vec![x]); + let refused_as_unsupported = matches!( + bound_on_box(call, &pool, &[(x, 0.25, 0.5)], &opts), + Err(ValidatedError::Unsupported { .. }) + ); + assert_eq!( + flag, !refused_as_unsupported, + "`{name}`: taylor_model flag = {flag} but bound_on_box \ + unsupported = {refused_as_unsupported}" + ); + } + } + + /// `NUMERIC_BALL` is not a stale bit that should have been `false`: every + /// primitive that has it and lacks a Taylor-model rule really does + /// evaluate as a ball. The two flags answer different questions — + /// pointwise enclosure vs. a polynomial model with a rigorous remainder — + /// which is exactly why reading the first as the second is a trap. + #[test] + fn numeric_ball_is_accurate_where_it_differs_from_taylor_model() { + use crate::ball::ArbBall; + let reg = crate::primitive::PrimitiveRegistry::default_registry(); + let mut differing = 0usize; + for (name, caps) in reg.iter() { + if !caps.contains(crate::primitive::Capabilities::NUMERIC_BALL) + || caps.contains(crate::primitive::Capabilities::TAYLOR_MODEL) + { + continue; + } + differing += 1; + let arg = [ArbBall::from_f64(1.0, 128)]; + assert!( + reg.numeric_ball(name, &arg).is_some(), + "`{name}` advertises numeric_ball but has none" + ); + } + assert!( + differing > 0, + "the two flags have become the same question — this test is now vacuous" + ); + } + + #[test] + fn refusal_names_the_blocking_function() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = pool.mul(vec![x, pool.func("bessel_j0", vec![x])]); + let what = taylor_model_refusal(e, &pool).expect("bessel_j0 has no Taylor rule"); + assert!(what.contains("bessel_j0"), "{what}"); + assert_eq!( + taylor_model_blockers(e, &pool), + vec!["bessel_j0".to_string()] + ); + } + + #[test] + fn supported_expression_has_no_refusal_and_no_blockers() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = pool.add(vec![pool.func("sin", vec![x]), pool.func("exp", vec![x])]); + assert!(taylor_model_refusal(e, &pool).is_none()); + assert!(taylor_model_blockers(e, &pool).is_empty()); + } + + /// A domain violation is not a refusal to model: `log` has a rule, the + /// box is just bad. Reporting it as unsupported would send a planner off + /// a perfectly good route. + #[test] + fn domain_violation_is_not_unsupported() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = pool.func("log", vec![x]); + assert!(taylor_model_refusal(e, &pool).is_none()); + let opts = BoundOptions::default(); + assert!(bound_on_box(e, &pool, &[(x, -2.0, -1.0)], &opts).is_err()); + } + + /// Constant expressions have no free symbols; the probe must still reach + /// the function nodes rather than failing on an empty box. + #[test] + fn constant_expression_is_classified() { + let pool = ExprPool::new(); + let two = pool.integer(2_i32); + assert!(taylor_model_refusal(pool.func("sin", vec![two]), &pool).is_none()); + assert!(taylor_model_refusal(pool.func("erf", vec![two]), &pool).is_some()); + } + + #[test] + fn arity_is_part_of_the_question() { + assert!(taylor_model_supports_call("sin", 1)); + assert!(!taylor_model_supports_call("sin", 2)); + assert!(!taylor_model_supports_call("sin", 0)); + assert!(!taylor_model_supports("atan2")); + } + + /// The cached answer is the freshly probed answer. + #[test] + fn cache_is_transparent() { + for name in ["sin", "erf", "sqrt", "digamma"] { + let cached = taylor_model_supports_call(name, 1); + assert_eq!(cached, probe_call(name, 1), "{name}"); + assert_eq!(cached, taylor_model_supports_call(name, 1), "{name}"); + } + } +} diff --git a/alkahest-core/src/validated/bounds.rs b/alkahest-core/src/validated/bounds.rs index 9a67e9a9..8beaafbc 100644 --- a/alkahest-core/src/validated/bounds.rs +++ b/alkahest-core/src/validated/bounds.rs @@ -28,7 +28,9 @@ //! omitting that piece of the domain from the answer. use super::taylor::{taylor_range, TaylorContext, MAX_ORDER}; -use super::{contains_zero, from_bounds, from_float, is_finite, lb, ub, width, ValidatedError}; +use super::{ + contains_zero, from_bounds, from_float, is_finite, lb, mag, ub, width, ValidatedError, +}; use crate::ball::ArbBall; use crate::diff::diff; use crate::kernel::subs::subs; @@ -249,10 +251,22 @@ fn bound_on_fboxes( } let prec = opts.prec; - let (lo_bound, used_lo, exhausted_lo) = - extremum_search(expr, pool, boxes0, opts, Extremum::Min)?; - let (hi_bound, used_hi, exhausted_hi) = - extremum_search(expr, pool, boxes0, opts, Extremum::Max)?; + let (lo_bound, used_lo, exhausted_lo) = extremum_search( + expr, + pool, + boxes0, + opts, + Extremum::Min, + SearchGoal::Tolerance, + )?; + let (hi_bound, used_hi, exhausted_hi) = extremum_search( + expr, + pool, + boxes0, + opts, + Extremum::Max, + SearchGoal::Tolerance, + )?; let enclosure = from_bounds(&lo_bound, &hi_bound, prec); if !is_finite(&enclosure) { @@ -274,6 +288,25 @@ enum Extremum { Max, } +/// What makes an [`extremum_search`] pass stop early. +#[derive(Clone, Copy, PartialEq, Eq)] +enum SearchGoal { + /// Stop once the extremum is pinned down to within `opts.tol`. + Tolerance, + /// Stop once the running rigorous bound is proven strictly positive in the + /// signed view — i.e. `min f > 0` for [`Extremum::Min`], `max f < 0` for + /// [`Extremum::Max`]. + /// + /// `tol` is the wrong stopping rule for a sign question, because it is an + /// **absolute** width. On an expression whose extremum is closer to zero + /// than `tol`, the tolerance test fires while the enclosure still straddles + /// zero, and the predicate can only answer `Undecided` — the search stopped + /// on a criterion unrelated to the question being asked. Under this goal + /// the pass instead refines until the sign is settled, or until the work + /// budget or the width floor stops it. + DecideSign, +} + /// Moore–Skelboe branch-and-bound for one end of the range. /// /// Returns `(bound, subdivisions, budget_exhausted)` where `bound` is a @@ -286,6 +319,7 @@ fn extremum_search( boxes0: &[FBox], opts: &BoundOptions, which: Extremum, + goal: SearchGoal, ) -> Result<(Float, usize, bool)> { let prec = opts.prec; let floor = max_dim_width(boxes0, prec) * 2f64.powi(-SINGULARITY_BISECTION_LIMIT); @@ -309,10 +343,24 @@ fn extremum_search( // Active list of (lower bound of signed range, box). `best_ub` is the best // proven upper bound on the signed extremum, from any box seen so far. let mut active: Vec<(Float, Vec)> = Vec::new(); + // Smallest key among boxes that have been bisected down to `floor` and so + // can never be refined again. Such a box must leave `active`: pushing it + // back would make it the argmin again on the very next iteration with + // nothing changed, and the loop would spin forever without consuming any + // budget. Only the key is kept, which is all the final bound needs. + let mut retired: Option = None; let mut best_ub: Option = None; let mut subdivisions = 0usize; let mut exhausted = false; + // Running minimum of two optional keys. + let keep_min = |cur: Option, k: Float| -> Option { + Some(match cur { + Some(c) if c <= k => c, + _ => k, + }) + }; + let seed = evaluate_box(expr, pool, boxes0, opts, floor, prec)?; match seed { BoxOutcome::Range(r) => { @@ -329,16 +377,34 @@ fn extremum_search( while let Some(idx) = argmin_key(&active) { let (key, b) = active.swap_remove(idx); + // Rigorous lower bound on the signed extremum as things stand: `key` is + // the smallest key still active, and no retired box holds anything + // smaller than `retired`. + let overall = match &retired { + Some(r) if r < &key => r.clone(), + _ => key.clone(), + }; + + if goal == SearchGoal::DecideSign && overall > 0 { + // `min f > 0` (or, in the signed view of a Max pass, `max f < 0`). + // The sign question is settled and no amount of further refinement + // can unsettle it. + active.push((key, b)); + break; + } + // Prune: this box cannot contain the extremum. if let Some(ub_best) = &best_ub { if &key > ub_best { continue; } // Converged: the uncertainty in the extremum is within tol. - let gap = Float::with_val(prec, ub_best - &key); - if gap <= tol_f { - active.push((key, b)); - break; + if goal == SearchGoal::Tolerance { + let gap = Float::with_val(prec, ub_best - &overall); + if gap <= tol_f { + active.push((key, b)); + break; + } } } @@ -348,16 +414,8 @@ fn extremum_search( break; } if max_dim_width(&b, prec) <= floor { - // Cannot refine further; keep it as-is so its bound still counts. - active.push((key, b)); - // Everything else is either pruned or equally unrefinable. - let stuck = active - .iter() - .all(|(_, bx)| max_dim_width(bx, prec) <= floor); - if stuck { - exhausted = true; - break; - } + // Cannot refine further: retire it, so the loop makes progress. + retired = keep_min(retired, key); continue; } @@ -380,22 +438,26 @@ fn extremum_search( } } - // The rigorous bound on the signed extremum is the smallest lower bound - // still active (pruned boxes provably cannot beat `best_ub`). + // Every point of the original box lies in an active box, a retired box, or + // a pruned one — and a pruned box provably cannot beat `best_ub`. So the + // smallest key over active ∪ retired is a rigorous bound on the signed + // extremum, whenever the loop happened to stop. let mut bound = active .iter() .map(|(k, _)| k.clone()) - .fold(None::, |acc, k| { - Some(match acc { - Some(cur) if cur <= k => cur, - _ => k, - }) - }) + .chain(retired.clone()) + .fold(None::, keep_min) .or_else(|| best_ub.clone()) .ok_or_else(|| ValidatedError::InvalidInput { what: "no enclosure was produced".into(), })?; + // Retiring every remaining box means the search ran out of width, not that + // it converged — the old `stuck` case, reported the same way. + if active.is_empty() && retired.is_some() { + exhausted = true; + } + if let Some(ub_best) = &best_ub { if &bound > ub_best { bound = ub_best.clone(); @@ -1327,9 +1389,52 @@ pub fn verified_sign( return Ok(Verdict::False); } + // The search above stopped on `tol`, which is an absolute width and so says + // nothing about the sign. Re-run it with the sign as the stopping rule: it + // keeps refining while the running bound still straddles zero, which is + // what decides an inequality whose margin is narrower than `tol`. + let strictly_signed = match predicate { + // `min f > 0` proves `f > 0` everywhere, hence also `f >= 0`. + SignPredicate::Positive | SignPredicate::NonNegative => { + sign_targeted_bound(expr, pool, boxes, opts, Extremum::Min)? > 0 + } + SignPredicate::Negative | SignPredicate::NonPositive => { + sign_targeted_bound(expr, pool, boxes, opts, Extremum::Max)? < 0 + } + }; + if strictly_signed { + return Ok(Verdict::True); + } + + // Still undecided, which is what a margin that *vanishes* at an endpoint + // always looks like to a subdivision search. Try the series argument, which + // is the only one of the two that can reach a tight endpoint at all. + if let Some(v) = endpoint_series_verdict(expr, pool, boxes, predicate, opts) { + return Ok(v); + } + Ok(Verdict::Undecided) } +/// Rigorous one-sided bound computed with [`SearchGoal::DecideSign`]: a lower +/// bound on `min f` for [`Extremum::Min`], an upper bound on `max f` for +/// [`Extremum::Max`]. +fn sign_targeted_bound( + expr: ExprId, + pool: &ExprPool, + boxes: &[(ExprId, f64, f64)], + opts: &BoundOptions, + which: Extremum, +) -> Result { + let prec = opts.prec; + let boxes0: Vec = boxes + .iter() + .map(|(v, lo, hi)| (*v, Float::with_val(prec, lo), Float::with_val(prec, hi))) + .collect(); + let (bound, _, _) = extremum_search(expr, pool, &boxes0, opts, which, SearchGoal::DecideSign)?; + Ok(bound) +} + /// Rigorously evaluate `expr` at a handful of points of the box (corners and /// centre) and report whether any of them *proves* the predicate false there. /// @@ -1394,6 +1499,374 @@ fn violates_at_some_sample( Ok(false) } +// --------------------------------------------------------------------------- +// Endpoint series proof — inequalities that are tight where the box ends +// --------------------------------------------------------------------------- + +/// Highest derivative order the endpoint expansion will search for the first +/// coefficient that is not proven to vanish. +const MAX_VANISHING_ORDER: usize = 16; + +/// Taylor terms kept past the leading one before the Lagrange remainder takes +/// over. Three gives the halving loop room to work while keeping the number of +/// symbolic derivatives — and so their size — small. +const SERIES_TAIL_TERMS: usize = 3; + +/// Halvings of the candidate sub-interval tried while looking for one on which +/// the series argument closes. +const SERIES_HALVINGS: u32 = 240; + +/// A rigorously bounded Taylor expansion of `g` at one endpoint of the box. +/// +/// Write `p` for the endpoint and `s = ±1` for the direction pointing into the +/// box. This describes `h(t) = g(p + s·t)` on `t ∈ [0, span]`: +/// +/// ```text +/// h(t) = Σ_{k, + /// Enclosure of the Lagrange constant `sup|g^{(m)}| / m!` over the whole + /// candidate interval. + rem: ArbBall, + /// Width of the box, rounded **down**: the largest collar the certificate + /// covers, and the starting point of the halving sequence. + span: Float, +} + +/// Build the endpoint expansion of `g` at `p`, looking a distance `span` into +/// the box in direction `inward` (`+1` for a left endpoint, `-1` for a right +/// one). +/// +/// # Why the result is rigorous +/// +/// Two separate facts are established, by two different mechanisms: +/// +/// * **The vanishing coefficients really are zero.** `c_k = s^k g^{(k)}(p)/k!` +/// is accepted as zero only when [`vanishes_exactly`] — substitution followed +/// by `simplify`, checked to land on the literal integer `0` — says so. A +/// numeric enclosure can never prove a value is zero, and none is used for +/// that here. As elsewhere in this module the symbolic verdict is +/// cross-checked against ball arithmetic, and a disagreement (which would +/// mean a simplifier bug) abandons the expansion rather than building on it. +/// * **The remainder really is bounded.** `bound_on_fboxes` is run on every +/// derivative `g, g', …, g^{(m)}` over the closed candidate interval. Those +/// calls succeed only if each is analytic throughout — [`super::taylor`] +/// refuses otherwise — which is exactly the hypothesis of Taylor's theorem +/// with Lagrange remainder on that interval. `rem` is then an outward-rounded +/// bound on `sup|g^{(m)}|/m!` there, so `|R(t)| ≤ rem·t^m` holds for every +/// `t` in the interval, and a fortiori on any sub-interval `[0, δ]` of it. +/// This is a *proven* remainder, not a truncation assumed to be small. +/// +/// Returns `None` — never an error — whenever any step declines. The series +/// argument is an extra attempt layered on top of the subdivision search, so +/// its failure must leave the caller's verdict at `Undecided`, not turn it into +/// a refusal. +fn expand_at_endpoint( + g: ExprId, + pool: &ExprPool, + var: ExprId, + ilo: &Float, + ihi: &Float, + inward: i32, + opts: &BoundOptions, +) -> Option { + let prec = opts.prec; + // Expand at whichever end `inward` points away from. + let p = if inward > 0 { ilo } else { ihi }; + // Rounded **down**, so `p ± span` can never leave the box and every `δ` the + // halving loop considers names a collar the analyticity certificate below + // actually covers. + let span = Float::with_val_round(prec, ihi - ilo, Round::Down).0; + if span <= 0 { + return None; + } + + // Successive symbolic derivatives, simplified at each step to keep them + // from growing. + let mut derivs: Vec = vec![g]; + let extend = |derivs: &mut Vec| -> bool { + let last = *derivs.last().expect("derivs is never empty"); + match diff(last, var, pool) { + Ok(d) => { + derivs.push(simplify(d.value, pool).value); + true + } + Err(_) => false, + } + }; + + // Locate `j`, the first coefficient not proven to vanish at `p`. + let mut found = None; + for k in 0..=MAX_VANISHING_ORDER { + if k > 0 && !extend(&mut derivs) { + return None; + } + if vanishes_exactly(derivs[k], pool, var, p) { + // `simplify` claims an exact zero; ball arithmetic must agree. + if !enclosure_admits_zero(derivs[k], pool, var, p, opts.order, prec) { + return None; + } + continue; + } + found = Some(k); + break; + } + let j = found?; + + let m = j + 1 + SERIES_TAIL_TERMS; + while derivs.len() <= m { + if !extend(&mut derivs) { + return None; + } + } + + // Analyticity certificate for Taylor's theorem, plus the sup bound feeding + // the Lagrange remainder. Both are taken over the *whole* box, so they hold + // on every collar the halving loop can propose. + let interval = vec![(var, ilo.clone(), ihi.clone())]; + let mut sup_m = None; + for (k, &d) in derivs.iter().enumerate().take(m + 1) { + let r = bound_on_fboxes(d, pool, &interval, opts).ok()?; + if !is_finite(r.enclosure()) { + return None; + } + if k == m { + sup_m = Some(mag(r.enclosure())); + } + } + + // Exact factorials (20! < 2^62, so these are exact at working precision). + let mut fact = vec![Float::with_val(prec, 1)]; + for k in 1..=m { + let next = Float::with_val(prec, &fact[k - 1] * k as u32); + fact.push(next); + } + + let point = vec![(var, p.clone(), p.clone())]; + let at_point = |k: usize| taylor_range(derivs[k], pool, &point, opts.order, prec).ok(); + + // c_j = s^j · g^{(j)}(p) / j! + let mut cj = (at_point(j)? / from_float(&fact[j], prec))?; + if inward < 0 && j % 2 == 1 { + cj = -cj; + } + if !is_finite(&cj) { + return None; + } + + // |c_k| for j < k < m — the sign is irrelevant, so `s` does not enter. + let mut tail = Vec::with_capacity(m - j - 1); + for (k, fk) in fact.iter().enumerate().take(m).skip(j + 1) { + let b = (at_point(k)?.abs_ball() / from_float(fk, prec))?; + if !is_finite(&b) { + return None; + } + tail.push(b); + } + + let rem = (from_float(&sup_m?, prec) / from_float(&fact[m], prec))?; + if !is_finite(&rem) { + return None; + } + + Some(EndpointSeries { + j, + m, + cj, + tail, + rem, + span, + }) +} + +impl EndpointSeries { + /// Upper bound on the tail `Σ_{j Float { + let d = from_float(delta, prec); + let mut acc = ArbBall::from_f64(0.0, prec); + for (i, c) in self.tail.iter().enumerate() { + acc = acc + c.clone() * d.powi((i + 1) as i64); + } + acc = acc + self.rem.clone() * d.powi((self.m - self.j) as i64); + ub(&acc) + } + + /// The largest `δ` of the halving sequence starting at `span` on which + /// `h(t)` is proven to keep the sign of `c_j` throughout `t ∈ [0, δ]`, + /// together with that sign (`true` for positive). + /// + /// # The argument + /// + /// With `c_0 … c_{j-1}` proven zero, + /// + /// ```text + /// h(t) = c_j t^j + Σ_{j tail_bound(δ)` gives `h ≥ 0` on `[0, δ]`, and + /// symmetrically a proven `c_j < −tail_bound(δ)` gives `h ≤ 0` there, with + /// `h < 0` for `t > 0`. + fn reach(&self, prec: u32) -> Option<(Float, bool)> { + let (lo, hi) = (lb(&self.cj), ub(&self.cj)); + // Only a *proven* sign for the leading coefficient is usable. + let (positive, margin) = if lo > 0 { + (true, lo) + } else if hi < 0 { + (false, Float::with_val(prec, -hi)) + } else { + return None; + }; + + let mut delta = self.span.clone(); + for _ in 0..SERIES_HALVINGS { + if margin > self.tail_bound(&delta, prec) { + return Some((delta, positive)); + } + delta = Float::with_val(prec, &delta / 2u32); + if delta <= 0 { + break; + } + } + None + } +} + +/// Try to settle a sign predicate whose margin vanishes at an endpoint of the +/// box, by splitting it into series-proved collars at the ends and an ordinary +/// branch-and-bound in the middle. +/// +/// Subdivision alone provably cannot do this: where the margin goes to zero, +/// every enclosure of the range straddles zero no matter how fine the boxes +/// get, so the honest verdict from that machinery is always `Undecided`. The +/// standard remedy, and the one used here, is a truncated Taylor expansion at +/// the endpoint with a rigorous remainder — see [`expand_at_endpoint`] for why +/// the remainder is proven rather than assumed, and [`EndpointSeries::reach`] +/// for the positivity argument. +/// +/// # Composition +/// +/// The three pieces are `[a, a+δₗ]`, `[a+δₗ, b−δᵣ]` and `[b−δᵣ, b]`. They are +/// closed and share their endpoints, so their union is exactly `[a, b]` with no +/// gap — in particular the join points themselves are covered twice, never +/// zero times. A collar is only claimed when its own expansion closed, and the +/// middle is only claimed when the subdivision search proves a strict sign +/// there; `True` is returned only if every piece that exists is proven. If the +/// collars already meet or overlap there is no middle piece to prove. +/// +/// Restricted to one dimension, since the expansion is in a single variable. +fn endpoint_series_verdict( + expr: ExprId, + pool: &ExprPool, + boxes: &[(ExprId, f64, f64)], + predicate: SignPredicate, + opts: &BoundOptions, +) -> Option { + if boxes.len() != 1 { + return None; + } + let (var, lo, hi) = boxes[0]; + if lo >= hi || !lo.is_finite() || !hi.is_finite() { + return None; + } + let prec = opts.prec; + + // Reduce every predicate to a statement about `g ≥ 0` / `g > 0`. + let g = match predicate { + SignPredicate::NonNegative | SignPredicate::Positive => expr, + SignPredicate::NonPositive | SignPredicate::Negative => { + pool.mul(vec![pool.integer(-1_i32), expr]) + } + }; + let strict = matches!(predicate, SignPredicate::Positive | SignPredicate::Negative); + + let a = Float::with_val(prec, lo); + let b = Float::with_val(prec, hi); + + let left = expand_at_endpoint(g, pool, var, &a, &b, 1, opts); + let right = expand_at_endpoint(g, pool, var, &a, &b, -1, opts); + + // A leading coefficient proven *negative* at an endpoint disproves the + // predicate outright: `h(t) < 0` for every `t ∈ (0, δ]`, and those points + // are in the box. + for s in [&left, &right].into_iter().flatten() { + if let Some((_, false)) = s.reach(prec) { + return Some(Verdict::False); + } + } + + // `j ≥ 1` means `g` was *proven* to be exactly zero at that endpoint, which + // is a point of the box — so a strict inequality fails there. + if strict + && [&left, &right] + .iter() + .any(|e| e.as_ref().is_some_and(|s| s.j >= 1)) + { + return Some(Verdict::False); + } + + let collar = |e: &Option| -> Option { + e.as_ref() + .and_then(|s| s.reach(prec)) + .filter(|(_, positive)| *positive) + .map(|(d, _)| d) + }; + let dl = collar(&left); + let dr = collar(&right); + if dl.is_none() && dr.is_none() { + return None; + } + + // The join points are rounded **into** the collars — `a+δₗ` down, `b−δᵣ` up + // — so the middle piece can only overlap what the series already proved, + // never leave a sliver between them uncovered. Rounding the other way would + // open a gap narrower than an ulp, and a gap is a hole in the proof however + // narrow it is. + let mlo = dl.map_or_else( + || a.clone(), + |d| Float::with_val_round(prec, &a + &d, Round::Down).0, + ); + let mhi = dr.map_or_else( + || b.clone(), + |d| Float::with_val_round(prec, &b - &d, Round::Up).0, + ); + + // The collars already cover the box; nothing left to prove. + if mlo >= mhi { + return Some(Verdict::True); + } + + let middle = vec![(var, mlo, mhi)]; + let bound = extremum_search( + g, + pool, + &middle, + opts, + Extremum::Min, + SearchGoal::DecideSign, + ) + .ok()?; + (bound.0 > 0).then_some(Verdict::True) +} + #[cfg(test)] mod tests { use super::*; @@ -2020,4 +2493,176 @@ mod tests { ); } } + + // ── endpoint-tight inequalities ───────────────────────────────────── + + /// `x·(2 + cos x) − 3·sin x` — Cusa–Huygens with the denominator cleared. + /// Non-negative on `[0, π/2)` and tight as `x → 0`, where it vanishes to + /// fifth order. + fn cusa_huygens(pool: &ExprPool, x: ExprId) -> ExprId { + let two_plus_cos = pool.add(vec![pool.integer(2_i32), pool.func("cos", vec![x])]); + sub( + pool, + pool.mul(vec![x, two_plus_cos]), + pool.mul(vec![pool.integer(3_i32), pool.func("sin", vec![x])]), + ) + } + + #[test] + fn tight_at_the_endpoint_is_certified_by_the_series_split() { + // Subdivision alone can never do this: at x = 0 the margin is zero, so + // every enclosure of the range straddles zero however fine the boxes. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let v = verified_sign( + cusa_huygens(&pool, x), + &pool, + &[(x, 0.0, 1.5)], + SignPredicate::NonNegative, + &opts(), + ) + .unwrap(); + assert_eq!(v, Verdict::True); + } + + #[test] + fn reversing_a_tight_inequality_is_refuted_not_certified() { + // The control for the above: a sign error in the series argument would + // surface here as a `True`. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let reversed = pool.mul(vec![pool.integer(-1_i32), cusa_huygens(&pool, x)]); + let v = verified_sign( + reversed, + &pool, + &[(x, 0.0, 1.5)], + SignPredicate::NonNegative, + &opts(), + ) + .unwrap(); + assert_eq!(v, Verdict::False); + } + + #[test] + fn a_violation_only_next_to_the_endpoint_is_never_certified() { + // `x^3 − x^2/1000` is negative exactly on (0, 1/1000): invisible to + // endpoint and centre sampling, and far narrower than `tol`. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = sub( + &pool, + pool.pow(x, pool.integer(3_i32)), + pool.mul(vec![ + pool.rational(1_i32, 1000_i32), + pool.pow(x, pool.integer(2_i32)), + ]), + ); + let v = verified_sign( + e, + &pool, + &[(x, 0.0, 1.5)], + SignPredicate::NonNegative, + &opts(), + ) + .unwrap(); + assert_ne!(v, Verdict::True); + } + + #[test] + fn tightness_in_the_interior_is_still_out_of_reach() { + // `(x − 7/10)^2 · (x + 1)` is non-negative and touches zero at an + // interior point. The endpoint expansion does not apply, and the + // honest answer stays Undecided rather than becoming a wrong True. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let shifted = sub(&pool, x, pool.rational(7_i32, 10_i32)); + let e = pool.mul(vec![ + shifted, + shifted, + pool.add(vec![x, pool.integer(1_i32)]), + ]); + let v = verified_sign( + e, + &pool, + &[(x, 0.0, 1.5)], + SignPredicate::NonNegative, + &opts(), + ) + .unwrap(); + assert_eq!(v, Verdict::Undecided); + } + + #[test] + fn strict_positivity_is_false_where_the_function_vanishes_exactly() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = pool.mul(vec![x, x]); + assert_eq!( + verified_sign(e, &pool, &[(x, 0.0, 1.0)], SignPredicate::Positive, &opts()).unwrap(), + Verdict::False + ); + assert_eq!( + verified_sign( + e, + &pool, + &[(x, 0.0, 1.0)], + SignPredicate::NonNegative, + &opts() + ) + .unwrap(), + Verdict::True + ); + } + + // ── termination ───────────────────────────────────────────────────── + + #[test] + fn an_unreachable_tolerance_terminates_instead_of_spinning() { + // Regression: a sub-box bisected down to the width floor used to be + // pushed back onto the active list and immediately re-selected, so the + // loop spun without ever consuming its subdivision budget. It only + // showed up once `tol` was out of reach — which a large rational + // coefficient arranges, because `tol` is an *absolute* width. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = sub( + &pool, + pool.mul(vec![ + pool.integer(1_000_000_000_000_i64), + pool.func("sin", vec![x]), + ]), + pool.mul(vec![pool.integer(636_619_772_368_i64), x]), + ); + let tight = BoundOptions { + tol: 1e-40, + ..opts() + }; + let r = bound_on_box(e, &pool, &[(x, 0.0, 1.5)], &tight).unwrap(); + assert!(r.lower() <= 0.0 && r.upper() >= 0.0, "{r:?}"); + } + + #[test] + fn a_sharp_rational_constant_is_certified_rather_than_timing_out() { + // `10^12·sin x − 636619772368·x >= 0` on [0, 3/2]: Jordan's inequality + // with 2/π rationalised to twelve digits, tight at x = 0. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = sub( + &pool, + pool.mul(vec![ + pool.integer(1_000_000_000_000_i64), + pool.func("sin", vec![x]), + ]), + pool.mul(vec![pool.integer(636_619_772_368_i64), x]), + ); + let v = verified_sign( + e, + &pool, + &[(x, 0.0, 1.5)], + SignPredicate::NonNegative, + &opts(), + ) + .unwrap(); + assert_eq!(v, Verdict::True); + } } diff --git a/alkahest-py/src/lib.rs b/alkahest-py/src/lib.rs index 73a3ce37..699a9a41 100644 --- a/alkahest-py/src/lib.rs +++ b/alkahest-py/src/lib.rs @@ -9349,6 +9349,124 @@ fn parse_box( (pool, out) } +/// Whether the validated-bounds subsystem can bound an expression at all. +/// +/// Returned by :func:`alkahest.bounds_supported`. Truthy exactly when every +/// construct in the expression has a rigorous Taylor-model rule, so it drops +/// into ``if ak.bounds_supported(f):``. +/// +/// Attributes +/// ---------- +/// supported : bool +/// The verdict. ``bool(self)`` is the same value. +/// blocker : str or None +/// The evaluator's own description of the **first** construct it has no +/// rule for (``"function `bessel_j0`"``), or ``None`` when supported. +/// functions : list of str +/// Every function in the expression with no Taylor-model rule, sorted. +/// Empty when ``supported`` is ``True`` — and possibly empty when it is +/// ``False``, if the blocker is a node kind rather than a function. +/// detail : str +/// One-sentence human explanation. +#[pyclass(name = "BoundsSupport", frozen)] +struct PyBoundsSupport { + #[pyo3(get)] + supported: bool, + #[pyo3(get)] + blocker: Option, + #[pyo3(get)] + functions: Vec, + #[pyo3(get)] + detail: String, +} + +#[pymethods] +impl PyBoundsSupport { + fn __bool__(&self) -> bool { + self.supported + } + + /// The verdict as a plain dict, for logging and JSON. + fn as_dict(&self, py: Python<'_>) -> PyResult { + let d = PyDict::new_bound(py); + d.set_item("supported", self.supported)?; + d.set_item("blocker", self.blocker.clone())?; + d.set_item("functions", self.functions.clone())?; + d.set_item("detail", &self.detail)?; + Ok(d.into_py(py)) + } + + fn __repr__(&self) -> String { + match &self.blocker { + None => "BoundsSupport(supported=True)".to_string(), + Some(what) => format!("BoundsSupport(supported=False, blocker={what:?})"), + } + } +} + +/// `alkahest.bounds_supported(expr) -> BoundsSupport` +/// +/// Can the validated-bounds subsystem bound this expression at all? +/// +/// :func:`bound_on_box`, :func:`verified_integral`, :func:`verified_no_roots` +/// and :func:`verified_sign` all evaluate through the same Taylor-model +/// evaluator, and it refuses any construct it has no rigorous rule for with +/// ``E-VALIDATED-001``. This answers that question *without running the +/// bound*, so a planning loop can pick a certifiable route instead of +/// discovering the boundary by hitting it. +/// +/// The answer is produced by running the real evaluator on a probe box, so it +/// cannot drift from what :func:`bound_on_box` does. The same information per +/// primitive is the ``taylor_model`` flag in +/// ``capabilities()["primitives"]``. +/// +/// **What it does not promise.** ``True`` means no ``E-VALIDATED-001``. A +/// supported function can still be refused on a *particular box* for a +/// domain violation (``E-VALIDATED-003``, e.g. ``log`` on ``[-2, -1]``) or a +/// non-finite enclosure (``E-VALIDATED-004``) — those depend on the box, not +/// on the expression, so no box-free predicate can rule them out. +/// +/// Note this is a different question from :func:`alkahest.certifiable`, which +/// asks whether an operation emits a **Lean** certificate. A rigorous +/// enclosure is not a Lean proof term, and the validated subsystem is not in +/// the certificate ledger; conflating the two under one predicate would make +/// a ``True`` mean two different kinds of evidence. +/// +/// Examples +/// -------- +/// >>> import alkahest as ak +/// >>> p = ak.ExprPool(); x = p.symbol("x") +/// >>> bool(ak.bounds_supported(ak.sin(x) * ak.exp(x))) +/// True +/// >>> answer = ak.bounds_supported(ak.bessel_j0(x)) +/// >>> bool(answer), answer.functions +/// (False, ['bessel_j0']) +#[pyfunction] +#[pyo3(name = "bounds_supported")] +fn py_bounds_supported(py: Python<'_>, expr: PyRef) -> PyResult { + guard_expr_depth(py, &expr)?; + let pool_py = expr.pool.clone_ref(py); + let pool = pool_py.borrow(py); + let blocker = alkahest_core::taylor_model_refusal(expr.id, &pool.inner); + let functions = alkahest_core::taylor_model_blockers(expr.id, &pool.inner); + let detail = match &blocker { + None => "every construct has a rigorous Taylor-model rule; the validated-bounds \ + entry points will not refuse this expression with E-VALIDATED-001 (a \ + particular box may still hit a domain violation)" + .to_string(), + Some(what) => format!( + "the validated-bounds subsystem has no rigorous Taylor-model rule for \ + {what}; bound_on_box would refuse with E-VALIDATED-001" + ), + }; + Ok(PyBoundsSupport { + supported: blocker.is_none(), + blocker, + functions, + detail, + }) +} + /// `alkahest.bound_on_box(expr, box, *, order=6, prec=128, tol=1e-9, max_subdivisions=2048)` /// /// Rigorous enclosure of the **range** of `expr` over an axis-aligned box, @@ -9503,6 +9621,21 @@ fn py_verified_no_roots( /// ``"nonpositive"``. A ``"false"`` verdict is itself certified — either the /// enclosure proves the predicate fails everywhere, or a rigorously evaluated /// point witnesses the failure. +/// +/// An inequality that is **tight at an endpoint** of the box is still decided. +/// Subdivision alone cannot do it — where the margin vanishes, every enclosure +/// of the range straddles zero — so the box is split: a collar at the endpoint +/// is handled by a truncated Taylor expansion with a proven Lagrange remainder, +/// the rest by branch-and-bound. That covers the classical sharp trigonometric +/// inequalities (Cusa–Huygens, Mitrinović–Adamović, Wilker, Huygens, Jordan) on +/// boxes reaching ``x = 0``. Tightness in the *interior* is not covered and +/// stays ``"undecided"``. +/// +/// ``tol`` sets the tolerance of the *enclosure*, not of the verdict: it is an +/// absolute width, so it does not bound how close to zero the answer may be. +/// Once the enclosure has been computed, the search is re-run with the sign +/// itself as the stopping rule, which is what decides an inequality whose +/// margin is narrower than ``tol``. #[allow(clippy::too_many_arguments)] #[pyfunction] #[pyo3(name = "verified_sign", signature = (expr, r#box, predicate, *, order = 6, prec = 128, tol = 1e-9, max_subdivisions = 2048))] @@ -9756,6 +9889,7 @@ impl PyPrimitiveRegistry { ("numeric_ball", caps.contains(Capabilities::NUMERIC_BALL)), ("lower_llvm", caps.contains(Capabilities::LOWER_LLVM)), ("lean_theorem", caps.contains(Capabilities::LEAN_THEOREM)), + ("taylor_model", caps.contains(Capabilities::TAYLOR_MODEL)), ] .into_iter() .map(|(k, v)| (k.to_string(), v)) @@ -9811,6 +9945,15 @@ impl PyPrimitiveRegistry { "lean_theorem", caps.contains(Capabilities::LEAN_THEOREM).into_py(py), ), + // Not implied by `numeric_ball`: ball arithmetic is + // pointwise, a Taylor model needs a rule with a + // rigorous remainder. Derived by running the + // validated evaluator, so it cannot drift from what + // `bound_on_box` actually accepts. + ( + "taylor_model", + caps.contains(Capabilities::TAYLOR_MODEL).into_py(py), + ), ] .into_iter() .map(|(k, v)| (k.to_string(), v)) @@ -12351,6 +12494,8 @@ fn alkahest(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(py_cad_project, m)?)?; // P1 item 9 — rigorous global bounds (Taylor models / validated numerics) m.add_class::()?; + m.add_class::()?; + m.add_function(wrap_pyfunction!(py_bounds_supported, m)?)?; m.add_function(wrap_pyfunction!(py_bound_on_box, m)?)?; m.add_function(wrap_pyfunction!(py_verified_integral, m)?)?; m.add_function(wrap_pyfunction!(py_verified_no_roots, m)?)?; diff --git a/alkahest-skill/alkahest.md b/alkahest-skill/alkahest.md index 2e7bb02f..d5eceb40 100644 --- a/alkahest-skill/alkahest.md +++ b/alkahest-skill/alkahest.md @@ -481,6 +481,26 @@ Escalation when `decide` refuses: `sos_decompose` / `prove_nonneg` for a positiv certificate, `alkahest.smt` (z3's `nlsat` is complete over the reals), or `bound_on_box` / `verified_sign` if a rigorous statement over a box is enough. +**Check that route before you build the workload: `ak.bounds_supported(expr)`.** The +validated-bounds entry points (`bound_on_box`, `verified_integral`, +`verified_no_roots`, `verified_sign`) reach the elementary fragment only — `sin`, `cos`, +`tan`, `exp`, `log`, `sqrt`, `abs`, the inverse-trig and hyperbolic functions — and +refuse everything else with `E-VALIDATED-001`. Every special function is outside it: +`erf`, `bessel_j0/j1`, `digamma`, `lambert_w`, `gamma`, the elliptic integrals, +`floor`/`ceil`, and the two-argument `atan2`. **`capabilities()["primitives"][i]` +carries this as `taylor_model`; do not read `numeric_ball` as the coverage flag** — it +is pointwise ball arithmetic, it is `True` for `erf` and `bessel_j0`, and it says +nothing about whether a bound can be certified. `bounds_supported` answers for a whole +expression without running anything, and names the blocking functions: + +```python +answer = ak.bounds_supported(ak.bessel_j0(x) * x) +bool(answer), answer.functions # (False, ['bessel_j0']) +``` + +A `True` means "not `E-VALIDATED-001`"; a bad box can still refuse with +`E-VALIDATED-003` (domain violation) or `-004` (non-finite enclosure). + **`E-SOS-002` from that escalation is `unknown`, not "not SOS".** The SOS search covers an LP-representable subcone of the PSD cone at one `basis_degree` (one `level` for Handelman), so a refusal is compatible with `p` being SOS outside that subcone, SOS at a diff --git a/docs/features.md b/docs/features.md index 2ad879d2..e0ab0888 100644 --- a/docs/features.md +++ b/docs/features.md @@ -61,7 +61,7 @@ Current stable feature surface. - Coefficient asymptotics of rational generating functions (`experimental.coefficient_asymptotics`): singularity analysis with the leading constant by Richardson extrapolation; declines when the dominant singularity is not unique (equal-modulus poles make the coefficients oscillate) - Asymptotics of sums (`experimental.euler_maclaurin`): Euler–Maclaurin expansion of `Σ_{k=a}^{n} f(k)` with Bernoulli corrections, numerically gated, returning an `AsymptoticReport` that marks each hypothesis checked or assumed (the additive constant — γ for the harmonic numbers — is fitted, not proved, and labelled as such) -- Validated numerics (`bound_on_box`, `verified_integral`, `verified_no_roots`, `verified_sign`): Taylor models over a box with Moore–Skelboe branch-and-bound; rigorous range enclosures, definite-integral enclosures and three-valued (`true`/`false`/`undecided`) predicates. Sound before tight — a wide bound is returned rather than a wrong one, and unbounded cases refuse (`E-VALIDATED-*`) +- Validated numerics (`bound_on_box`, `verified_integral`, `verified_no_roots`, `verified_sign`): Taylor models over a box with Moore–Skelboe branch-and-bound; rigorous range enclosures, definite-integral enclosures and three-valued (`true`/`false`/`undecided`) predicates. Sound before tight — a wide bound is returned rather than a wrong one, and unbounded cases refuse (`E-VALIDATED-*`). Coverage is the elementary fragment and is queryable before you commit to a route: `bounds_supported(expr)`, and `taylor_model` per primitive in `capabilities()["primitives"]` (not `numeric_ball`, which is pointwise ball arithmetic and reaches further) ## Discrete mathematics diff --git a/docs/mdbook/src/validated-bounds.md b/docs/mdbook/src/validated-bounds.md index 42a730e0..2730c7a0 100644 --- a/docs/mdbook/src/validated-bounds.md +++ b/docs/mdbook/src/validated-bounds.md @@ -115,11 +115,89 @@ The last row is the honest limit. A double root never changes sign, so no witness pair exists; `"undecided"` is the answer, and it is not upgraded to `"false"` on the strength of an enclosure that merely touches zero. +### Inequalities that are tight at an endpoint + +The interesting inequalities are usually the sharp ones, and sharp means the +margin goes to zero somewhere. Subdivision alone cannot certify those: where the +margin vanishes, every enclosure of the range straddles zero however fine the +boxes get. + +Two separate things are done about it. `tol` is an **absolute** width, so it is +the wrong stopping rule for a sign question — an expression whose minimum is +`10⁻¹³` meets a `1e-9` tolerance while its enclosure still straddles zero. +`verified_sign` therefore re-runs the search with the sign itself as the goal, +refining while the bound straddles zero rather than to a fixed width. And where +the margin genuinely reaches zero, the box is **split**: a collar `[a, a+δ]` at +the endpoint is handled by a truncated Taylor expansion with a proven Lagrange +remainder, the rest by ordinary branch-and-bound. The pieces are closed and +share the join point, so their union is the original box. + +```python +x = pool.symbol("x") +# Cusa–Huygens, denominator cleared: x(2 + cos x) − 3 sin x ≥ 0, tight at x = 0 +f = x * (pool.integer(2) + ak.cos(x)) - pool.integer(3) * ak.sin(x) +ak.verified_sign(f, [(x, 0.0, 1.5)], "nonnegative") # "true" +``` + +Mitrinović–Adamović, Wilker, Huygens and Jordan's inequality behave the same +way. The remainder is *proven*, not assumed: a Taylor coefficient counts as zero +only when substitution and `simplify` land on a literal integer `0` — no numeric +enclosure can prove a value is zero — and the tail is bounded by +`sup|g⁽ᵐ⁾|/m!` enclosed over the whole collar, with analyticity certified by +requiring every derivative up to `g⁽ᵐ⁾` to enclose successfully there. + +The limits are worth knowing: + +| Case | Verdict | Why | +|---|---|---| +| tight at an endpoint of the box | `"true"` | the expansion applies there | +| leading coefficient proven negative | `"false"` | `g < 0` just inside the endpoint | +| `"positive"` where `g` provably vanishes | `"false"` | a strict claim fails at that point | +| tight in the **interior** | `"undecided"` | the expansion does not apply | + +`(x − 7/10)²(x + 1)` on `[0, 3/2]` is non-negative and touches zero in the +middle; it stays `"undecided"` rather than being upgraded on the strength of an +enclosure that merely touches zero. + +## Which functions are covered — ask before you build the workload + +Taylor models reach the **elementary fragment**: `exp`, `log`, `sqrt`, `sin`, +`cos`, `tan`, `asin`, `acos`, `atan`, `sinh`, `cosh`, `tanh`, `abs`, plus +arithmetic and integer/rational powers. Every special function is outside it — +`erf`, `erfc`, `bessel_j0`, `bessel_j1`, `digamma`, `lambert_w`, `gamma`, the +elliptic integrals, `floor`, `ceil`, `acosh`, `asinh` — and so is any +two-argument function such as `atan2`. + +That boundary is queryable, so a search loop can choose a certifiable route +instead of discovering it by hitting `E-VALIDATED-001`: + +```python +ak.bounds_supported(ak.sin(x) * ak.exp(x)) # truthy +answer = ak.bounds_supported(ak.bessel_j0(x)) +bool(answer), answer.functions # (False, ['bessel_j0']) +answer.blocker # "function `bessel_j0`" + +# Per primitive, in the agent contract: +{row["name"] for row in ak.capabilities()["primitives"] if row["taylor_model"]} +``` + +**`numeric_ball` is not this flag.** It reports pointwise ball arithmetic, +which `erf`, `bessel_j0`, `digamma` and `floor` all have; a Taylor model +additionally needs a rule with a rigorous Lagrange remainder, which they do +not. Both bits are honest — they answer different questions. `taylor_model` +and `bounds_supported` are derived by *running* the Taylor evaluator, not from +a maintained list, so neither can drift from what `bound_on_box` accepts. + +A `True` answer means "will not be refused with `E-VALIDATED-001`". It is not +a promise of success: a covered function can still hit a domain violation or +an infinite enclosure on a *particular* box, which is a property of the box +and not of the expression. + ## Refusals | Code | Meaning | |---|---| -| `E-VALIDATED-001` | No rigorous Taylor model rule for some primitive in the expression | +| `E-VALIDATED-001` | No rigorous Taylor model rule for some primitive in the expression (ask `bounds_supported` first — see above) | | `E-VALIDATED-002` | A free symbol has no interval in the box | | `E-VALIDATED-003` | A singularity or branch cut inside the box (e.g. `1/x` over a box containing 0) | | `E-VALIDATED-004` | An enclosure overflowed to infinity | diff --git a/python/alkahest/__init__.py b/python/alkahest/__init__.py index bd54fb78..09104787 100644 --- a/python/alkahest/__init__.py +++ b/python/alkahest/__init__.py @@ -99,6 +99,8 @@ ArbBall, # Explicit positive/nonzero refinement for conservative simplification Assumptions, + # Verdict from `bounds_supported` — see `Enclosure` below + BoundsSupport, CompileCache, # Phase 21: JIT compiled evaluation CompiledFn, @@ -163,6 +165,8 @@ bessel_j0, bessel_j1, bound_on_box, + # Is the validated-bounds subsystem even reachable for this expression? + bounds_supported, cad_lift, cad_project, # Rational-function cancel/together @@ -1926,6 +1930,20 @@ def capabilities() -> dict: ledger (:func:`certificate_coverage`), which is also where each primitive's ``lean_theorem`` bit comes from — so the capability bits and the coverage table cannot disagree. + + Each ``primitives`` row also carries ``taylor_model``: whether the + validated-bounds subsystem (:func:`bound_on_box`, + :func:`verified_integral`, :func:`verified_no_roots`, + :func:`verified_sign`) has a rigorous Taylor-model rule for it, i.e. + whether those entry points can bound it at all. **It is a different + question from ``numeric_ball``**, which is pointwise ball arithmetic + and is ``True`` for `erf`, `bessel_j0`, `digamma` and `floor` — + none of which can be bounded over a box. Reading the latter as + coverage is what made the boundary invisible before 3.9.0. The bit is + derived by running the Taylor evaluator, so it cannot drift from what + `bound_on_box` accepts; :func:`bounds_supported` asks the same + question for a whole expression. + Symbolic linear algebra (``Matrix.rref``, ``rank``, ``nullspace``, ``jordan_form``, ``minimal_polynomial``, LU/QR/Cholesky, etc.) is available on :class:`Matrix`; unsupported inputs raise :class:`LinearAlgebraError` @@ -2090,6 +2108,8 @@ def wrapper(*args, **kwargs): "Assumptions", # P1 search plumbing items 4–5 — budgets + batch/streaming fan-out "BatchItem", + # Verdict from `bounds_supported` + "BoundsSupport", "Budget", "BudgetExceededError", "CadError", @@ -2216,6 +2236,7 @@ def wrapper(*args, **kwargs): "bessel_j0", "bessel_j1", "bound_on_box", + "bounds_supported", # P1 search plumbing item 4 "budget_seed", "cad_lift", diff --git a/python/alkahest/_types.pyi b/python/alkahest/_types.pyi index 130523eb..400db3f7 100644 --- a/python/alkahest/_types.pyi +++ b/python/alkahest/_types.pyi @@ -48,6 +48,19 @@ class Certifiability: def __bool__(self) -> bool: ... def as_dict(self) -> dict[str, object]: ... +class BoundsSupport: + """Verdict from :func:`bounds_supported` — truthy iff every construct in + the expression has a rigorous Taylor-model rule.""" + + supported: bool + blocker: str | None + functions: list[str] + detail: str + + def __bool__(self) -> bool: ... + def as_dict(self) -> dict[str, object]: ... + +def bounds_supported(expr: Expr) -> BoundsSupport: ... def certifiable( op: str | object, *args: object, diff --git a/tests/test_agent_contract.py b/tests/test_agent_contract.py index a1301ddf..5accae9d 100644 --- a/tests/test_agent_contract.py +++ b/tests/test_agent_contract.py @@ -65,6 +65,11 @@ def test_cranelift_jit_enables_session_jit_flag(): "numeric_ball", "lower_llvm", "lean_theorem", + # Validated-bounds coverage. Distinct from `numeric_ball`, which is + # pointwise ball arithmetic and does *not* imply a Taylor-model rule — + # see tests/test_taylor_model_coverage.py, which cross-checks this bit + # against `bound_on_box` for every primitive. + "taylor_model", } == primitives[0].keys() diff --git a/tests/test_taylor_model_coverage.py b/tests/test_taylor_model_coverage.py new file mode 100644 index 00000000..2c78ca87 --- /dev/null +++ b/tests/test_taylor_model_coverage.py @@ -0,0 +1,235 @@ +"""The `taylor_model` capability flag must equal what `bound_on_box` does. + +`capabilities()["primitives"][i]["numeric_ball"]` is *not* the flag that +governs the validated-bounds subsystem, and used to be the only per-function +coverage bit exposed. Ball arithmetic is pointwise; a Taylor model needs a +rule with a rigorous remainder, written per function in +`alkahest-core/src/validated/taylor.rs`. Ten primitives — `erf`, `erfc`, +`bessel_j0`, `bessel_j1`, `digamma`, `lambert_w`, `acosh`, `asinh`, `floor`, +`ceil` — have real Arb ball arithmetic and no Taylor-model rule, so +`bound_on_box` refuses them with `E-VALIDATED-001` while `numeric_ball` says +``True``. That boundary was correct at runtime and invisible beforehand. + +The whole point of these tests is that the *new* flag cannot repeat that +mistake. `taylor_model` is derived by running the real evaluator (see +`alkahest-core/src/primitive/taylor_support.rs`), never from a list, and +``test_taylor_model_flag_agrees_with_bound_on_box`` re-derives it here the +only other way there is — by calling `bound_on_box` on every registered +primitive — and fails if the two ever disagree. +""" + +from __future__ import annotations + +import alkahest as ak +import pytest + +# Small budget: these tests ask *which error* comes back, never how tight the +# bound is, so there is no reason to pay for convergence. +_OPTS = {"order": 2, "prec": 64, "tol": 1e-3, "max_subdivisions": 4} +_MAX_ARITY = 3 +_UNSUPPORTED = "E-VALIDATED-001" + + +def _primitive_rows(): + return ak.capabilities()["primitives"] + + +def _bound_on_box_accepts(name: str, arity: int) -> bool: + """Does `bound_on_box` get past dispatch for ``name`` at this arity? + + ``True`` means "not refused with ``E-VALIDATED-001``". Every other + outcome — success, a domain violation, a non-finite enclosure — means the + evaluator *had* a rule and something about this particular box stopped it, + which is a different question and not what the flag claims. + """ + pool = ak.ExprPool() + args = [pool.symbol(f"x{i}") for i in range(arity)] + call = pool.func(name, args) + box = [(a, 0.25, 0.5) for a in args] + try: + ak.bound_on_box(call, box, **_OPTS) + except ak.ValidatedError as exc: + return exc.code != _UNSUPPORTED + except Exception: + # Any other failure still means the evaluator dispatched on the name. + return True + return True + + +@pytest.mark.parametrize("row", _primitive_rows(), ids=lambda row: row["name"]) +def test_taylor_model_flag_agrees_with_bound_on_box(row): + """The flag and the subsystem it describes, cross-checked per primitive. + + This is the guard that stops the two drifting apart. If it fails, either + a Taylor-model rule was added/removed without the flag following (it + cannot be — the flag is derived) or the derivation itself is asking the + wrong question. + """ + name = row["name"] + observed = any(_bound_on_box_accepts(name, n) for n in range(1, _MAX_ARITY + 1)) + assert row["taylor_model"] is observed, ( + f"capabilities() says taylor_model={row['taylor_model']} for `{name}`, " + f"but bound_on_box {'accepts' if observed else 'refuses ' + _UNSUPPORTED + ' for'} it" + ) + + +def test_the_flag_is_not_a_restatement_of_numeric_ball(): + """The bug this fixes: `numeric_ball` was read as the coverage flag. + + `numeric_ball` is *not* wrong for these ten — they really do have Arb ball + arithmetic (`alkahest-core/src/primitive/taylor_support.rs` pins that in + Rust). It answers a different question, which is why reading it as the + validated-bounds coverage bit cost a whole workload. + """ + rows = {row["name"]: row for row in _primitive_rows()} + ball_only = { + name for name, row in rows.items() if row["numeric_ball"] and not row["taylor_model"] + } + assert ball_only == { + "acosh", + "asinh", + "bessel_j0", + "bessel_j1", + "ceil", + "digamma", + "erf", + "erfc", + "floor", + "lambert_w", + } + + +def test_supported_set_is_the_elementary_fragment(): + """A pin on the boundary as it stands, so a change to it is deliberate. + + Widening this set is a feature (add the rule, then add the name here); + narrowing it silently would be a regression an agent's plan depends on. + """ + supported = {row["name"] for row in _primitive_rows() if row["taylor_model"]} + assert supported == { + "abs", + "acos", + "asin", + "atan", + "cos", + "cosh", + "exp", + "log", + "sin", + "sinh", + "sqrt", + "tan", + "tanh", + } + + +def test_registry_capabilities_and_the_table_agree(): + reg = ak.PrimitiveRegistry.default_registry() + for row in _primitive_rows(): + assert reg.capabilities(row["name"])["taylor_model"] is row["taylor_model"] + + +# --------------------------------------------------------------------------- +# bounds_supported — the same question for a whole expression +# --------------------------------------------------------------------------- + + +def test_bounds_supported_matches_bound_on_box_on_composite_expressions(): + pool = ak.ExprPool() + x = pool.symbol("x") + y = pool.symbol("y") + cases = [ + ak.sin(x) * ak.exp(x) + x * x, + x - x, + ak.sqrt(x + 1) / (x + 2), + ak.bessel_j0(x), + x * ak.erf(x), + ak.sin(x) + ak.digamma(x), + ak.atan2(x, y), + ak.gamma(x), + ak.floor(x) + ak.ceil(x), + ak.tanh(x * y) - ak.log(x + 3), + ] + box = [(x, 0.25, 0.5), (y, 0.25, 0.5)] + for expr in cases: + answer = ak.bounds_supported(expr) + try: + ak.bound_on_box(expr, box, **_OPTS) + refused = False + except ak.ValidatedError as exc: + refused = exc.code == _UNSUPPORTED + assert bool(answer) is (not refused), f"{expr}: {answer.detail}" + + +def test_bounds_supported_names_every_blocking_function(): + pool = ak.ExprPool() + x = pool.symbol("x") + answer = ak.bounds_supported(ak.bessel_j0(x) + ak.erf(x) + ak.sin(x)) + + assert not answer + assert answer.supported is False + assert answer.functions == ["bessel_j0", "erf"] + assert "bessel_j0" in answer.blocker + assert _UNSUPPORTED in answer.detail + + +def test_bounds_supported_is_truthy_and_quiet_when_everything_is_covered(): + pool = ak.ExprPool() + x = pool.symbol("x") + answer = ak.bounds_supported(ak.sin(x) ** 2 + ak.cosh(x)) + + assert answer + assert answer.blocker is None + assert answer.functions == [] + assert "E-VALIDATED-001" in answer.detail + + +def test_a_domain_violation_is_not_reported_as_unsupported(): + """`log` has a rule; `[-2, -1]` is just a bad box. + + Reporting the refusal as "unsupported" would send a planner off a route + that works everywhere else, which is the failure this predicate exists to + prevent — in the other direction. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + f = ak.log(x) + + assert ak.bounds_supported(f) + with pytest.raises(ak.ValidatedError) as excinfo: + ak.bound_on_box(f, [(x, -2.0, -1.0)], **_OPTS) + assert excinfo.value.code != _UNSUPPORTED + + +def test_arity_is_part_of_the_question(): + """`atan2` is refused for its arity, not for being `atan2`.""" + pool = ak.ExprPool() + x = pool.symbol("x") + answer = ak.bounds_supported(pool.func("atan2", [x, x])) + + assert not answer + assert "2 arguments" in answer.blocker + + +def test_constant_expressions_are_classified_too(): + pool = ak.ExprPool() + assert ak.bounds_supported(ak.sin(pool.integer(2))) + assert not ak.bounds_supported(ak.erf(pool.integer(2))) + + +def test_bounds_support_surface(): + pool = ak.ExprPool() + x = pool.symbol("x") + + for name in ("bounds_supported", "BoundsSupport"): + assert name in ak.__all__, name + + answer = ak.bounds_supported(ak.erf(x)) + assert isinstance(answer, ak.BoundsSupport) + assert "BoundsSupport(" in repr(answer) + assert answer.as_dict() == { + "supported": False, + "blocker": answer.blocker, + "functions": ["erf"], + "detail": answer.detail, + } diff --git a/tests/test_validated_bounds.py b/tests/test_validated_bounds.py index 3fbf9f30..dc69ceec 100644 --- a/tests/test_validated_bounds.py +++ b/tests/test_validated_bounds.py @@ -425,3 +425,187 @@ def test_enclosure_repr_and_helpers(): assert r.width >= 1.0 assert r.subdivisions >= 0 assert isinstance(r.budget_exhausted, bool) + + +# --------------------------------------------------------------------------- +# Inequalities that are tight where the box ends +# +# The classical trigonometric inequalities are asymptotically tight as x -> 0. +# Subdivision alone provably cannot certify them there: wherever the margin goes +# to zero, every enclosure of the range straddles zero however fine the boxes +# get. `verified_sign` therefore splits -- a truncated Taylor expansion with a +# proven Lagrange remainder on a collar at the endpoint, branch-and-bound on the +# rest -- and the two pieces share their join point, so nothing is left out. +# +# The soundness half of this section matters more than the reach half: a `true` +# is a certificate, so every "must not be true" case below is load-bearing. +# --------------------------------------------------------------------------- + + +_CLASSICAL_NAMES = ("cusa_huygens", "huygens", "mitrinovic_adamovic", "wilker") + + +def _classical(pool, x): + """The four classical inequalities, cleared of denominators. + + Each is stated as `f(x) >= 0` on `(0, pi/2)`, tight as `x -> 0`: + + * Cusa-Huygens `sin x / x < (2 + cos x) / 3` + * Mitrinovic-Adamovic `(sin x / x)^3 > cos x` + * Wilker `(sin x / x)^2 + tan x / x > 2` + * Huygens `2 sin x / x + tan x / x > 3` + + Wilker and Huygens are multiplied through by `x^2 cos x` and `x cos x`, + which are positive on the open interval, so the sign is unchanged. + """ + return { + "cusa_huygens": x * (pool.integer(2) + ak.cos(x)) - pool.integer(3) * ak.sin(x), + "mitrinovic_adamovic": ak.sin(x) ** 3 - x**3 * ak.cos(x), + "wilker": ak.sin(x) ** 2 * ak.cos(x) + x * ak.sin(x) - pool.integer(2) * x**2 * ak.cos(x), + "huygens": pool.integer(2) * ak.sin(x) * ak.cos(x) + + ak.sin(x) + - pool.integer(3) * x * ak.cos(x), + } + + +@pytest.mark.parametrize("name", _CLASSICAL_NAMES) +@pytest.mark.parametrize("lo", [0.0, 0.01, 0.1]) +def test_classical_trig_inequalities_are_certified_up_to_the_tight_endpoint(name, lo): + """All four, including at `x = 0` itself where the margin vanishes.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + assert ak.verified_sign(_classical(pool, x)[name], _box(x, lo, 1.5), "nonnegative") == "true" + + +def test_jordan_inequality_with_an_exactly_rational_bound(): + """`sin x >= (2/pi) x`, stated exactly. + + `pi` is a plain symbol here, so the constant is rationalised instead: + `636619772368/10^12 > 2/pi`, which makes `D sin x - N x >= 0` *stronger* + than Jordan's inequality on the same interval. It is tight at `x = 0`. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + f = ak.sin(x) * pool.integer(10**12) - pool.integer(636619772368) * x + + assert ak.verified_sign(f, _box(x, 0.0, 1.5), "nonnegative") == "true" + + +def test_jordan_with_too_large_a_constant_is_refuted_at_the_far_endpoint(): + """The same constant on `[0, pi/2]`, where `N/D > 2/pi` makes it false.""" + pool = ak.ExprPool() + x = pool.symbol("x") + f = ak.sin(x) * pool.integer(10**12) - pool.integer(636619772368) * x + + assert ak.verified_sign(f, _box(x, 0.0, math.pi / 2), "nonnegative") == "false" + + +@pytest.mark.parametrize("name", _CLASSICAL_NAMES) +def test_reversing_a_tight_inequality_never_certifies_it(name): + """The reverse of each is false, and tight at the same endpoint. + + This is the control that makes the endpoint machinery falsifiable: a sign + error in the series argument would show up here as a `true`. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + reversed_f = pool.integer(-1) * _classical(pool, x)[name] + + assert ak.verified_sign(reversed_f, _box(x, 0.0, 1.5), "nonnegative") != "true" + + +def test_false_only_in_a_tiny_neighbourhood_of_the_endpoint_is_not_certified(): + """`x^3 - x^2/1000` is negative exactly on `(0, 1/1000)`. + + The violation is invisible to endpoint and centre sampling and is far + narrower than the default tolerance, so nothing but the endpoint expansion + can see it at all. It must never come back `true`. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + f = x**3 - pool.rational(1, 1000) * x**2 + + assert ak.verified_sign(f, _box(x, 0.0, 1.5), "nonnegative") != "true" + + +def test_a_strict_inequality_is_false_where_the_function_vanishes_exactly(): + """`x^2 > 0` fails at `x = 0`; `x^2 >= 0` holds.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + assert ak.verified_sign(x**2, _box(x, 0.0, 1.0), "positive") == "false" + assert ak.verified_sign(x**2, _box(x, 0.0, 1.0), "nonnegative") == "true" + + +def test_tightness_away_from_an_endpoint_stays_undecided(): + """A margin that vanishes in the *interior* is still out of reach. + + `(x - 7/10)^2 (x + 1)` is non-negative on `[0, 3/2]` and touches zero at + `x = 7/10`. The endpoint expansion does not apply there, and the honest + answer remains `undecided` -- it must not be upgraded to a wrong `true`. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + f = (x - pool.rational(7, 10)) ** 2 * (x + pool.integer(1)) + + assert ak.verified_sign(f, _box(x, 0.0, 1.5), "nonnegative") == "undecided" + + +def test_a_shallow_interior_dip_is_never_certified_true(): + """`(x - 7/10)^2 - 10^-6` dips below zero only near `x = 7/10`.""" + pool = ak.ExprPool() + x = pool.symbol("x") + f = (x - pool.rational(7, 10)) ** 2 - pool.rational(1, 10**6) + + assert ak.verified_sign(f, _box(x, 0.0, 1.5), "nonnegative") != "true" + + +# --------------------------------------------------------------------------- +# Termination: a tolerance the search cannot reach +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "digits", + [3, 6, 9, 12], +) +def test_verified_sign_terminates_for_any_size_of_rational_constant(digits): + """Rationalising a constant to more digits must not change the cost class. + + Before the fix, `N/D` at 12 digits ran for over 300 s while 9 digits took + 0.08 s: a sub-box bisected down to the width floor was pushed back onto the + active list and immediately re-selected, so the loop spun without ever + consuming its subdivision budget. Three extra digits were enough to push the + tolerance out of reach and trigger it. + """ + import time + + d = 10**digits + n = int(0.636619772368 * d) + pool = ak.ExprPool() + x = pool.symbol("x") + f = ak.sin(x) * pool.integer(d) - pool.integer(n) * x + + start = time.monotonic() + verdict = ak.verified_sign(f, _box(x, 0.0, 1.5), "nonnegative") + elapsed = time.monotonic() - start + + assert verdict in ("true", "false", "undecided") + assert elapsed < 60.0, f"{digits} digits took {elapsed:.1f}s" + + +def test_bound_on_box_terminates_when_the_tolerance_is_unreachable(): + """`tol` below what the width floor can deliver must stop, not spin.""" + import time + + pool = ak.ExprPool() + x = pool.symbol("x") + f = ak.sin(x) * pool.integer(10**12) - pool.integer(636619772368) * x + + start = time.monotonic() + r = ak.bound_on_box(f, _box(x, 0.0, 1.5), tol=1e-40) + elapsed = time.monotonic() - start + + assert elapsed < 60.0, f"took {elapsed:.1f}s" + assert r.lower <= 0.0 <= r.upper