diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8da17c72..1e348371 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -89,6 +89,12 @@ jobs: - name: cargo clippy (parallel feature) run: cargo clippy --all-targets --features parallel -- -D warnings + # `cranelift` ships in the default PyPI wheel (PR #299) but had no clippy + # step, so three lints accumulated in code every wheel user runs. It needs + # no system LLVM — the backend is pure Rust — so it lints anywhere. + - name: cargo clippy (cranelift feature) + run: cargo clippy --all-targets --features cranelift -- -D warnings + - name: cargo clippy (jit feature) run: cargo clippy --all-targets --features jit -- -D warnings diff --git a/CHANGELOG.md b/CHANGELOG.md index f58c947f..63dd55db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -109,6 +109,33 @@ Both are detailed under *Behaviour changes to plan for*. ### Fixed +- **Interval evaluation refused `bessel_j0` / `bessel_j1`, and its Bessel ball + kernel was not an enclosure.** `evaluate(bessel_j0(x), {x: ArbBall(...)}, + mode="interval")` came back `status="unsupported"` with `E-EVAL-010` even + though both functions have rigorous ball kernels and `capabilities()` reported + `numeric_ball: True` for both. The evaluator carried its own hand-written + match over function names — the third such list in this area — and these two + had simply never been added to it. It no longer has one: every `Func` node is + now dispatched through the primitive registry, so the set of functions + interval evaluation accepts **is** the set the registry advertises a + `numeric_ball` kernel for, by construction rather than by agreement. Auditing + the two sets against each other also turned up the gap in the other direction: + `atanh` had a ball kernel and did *not* advertise it, because the capability + probe tested ball kernels at `1.0` only and `atanh`'s domain is the open + interval `(-1, 1)` — it declined the sole probe point and lost a bit it had + earned, while keeping `numeric_f64` because that probe already tried `0.5`. + Ball kernels are now probed at the same points as `f64` ones. Wrong-arity + calls are declined rather than silently bounding the first argument + (`sin(x, y)` was `sin(x)`), which the generic dispatch made reachable. + + **`ArbBall::bessel_jn` was separately unsound and is rewritten.** It hulled + `Jₙ(lo)` and `Jₙ(hi)`, which is only an enclosure for a *monotone* function — + `J₀` on `[-1, 1]` has equal endpoints (`≈ 0.7652`), so the hull collapsed to a + point that excluded `J₀(0) = 1`, the function's own maximum. It is now a + midpoint evaluation plus a mean-value bound with `L = 1`, which is rigorous at + every order: `|Jₙ| ≤ 1` for all real `x`, and `J₀′ = −J₁`, + `Jₙ′ = (Jₙ₋₁ − Jₙ₊₁)/2`, so `|Jₙ′| ≤ 1`. A randomised sweep samples the true + function inside 200 random intervals and checks nothing escapes the enclosure. - **`verified_no_roots` could not prove a root *exists* past an even root count.** The `"false"` direction fired only when a sign change was visible at the box's own endpoints, so any even number of roots defeated it however @@ -124,9 +151,31 @@ Both are detailed under *Behaviour changes to plan for*. several variables — `x − y` on `[-1, 1]²` is now `"false"` where it used to be `"undecided"`. Continuity, which the argument needs, is exactly what the full-box enclosure succeeding already certifies. **Nothing was weakened to buy - 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. + this**: a root that never produces a sign change and never lands on a point the + search can *prove* is a root — a double root like `(x−1/3)²` on `[0, 2]`, or + `(x²−1)²` on `[-2, 2]` — still answers `"undecided"`, because no witness pair + exists and none is invented. +- **`verified_no_roots` could not see a root sitting exactly on a box + endpoint.** `x` on `[0, 1]` and `x−1` on `[0, 1]` both came back + `"undecided"`, as did `sin(x)` on `[0, 1]` and `log(x)` on `[1, 2]`. The + sign-change search above provably cannot settle them — `x` is non-negative + everywhere on `[0, 1]`, so the negative witness it looks for does not exist to + be found, however long it runs. But no search is needed: the box is **closed**, + so a point of it at which the expression is *proven* to be zero is a root in + the box, full stop. `"false"` is now returned when either of two independent + proofs is in hand — the existing sign change, or a point whose value is pinned + to zero. That also settles a root of even multiplicity that lands on a point + the search visits (`(x−1)²` on `[0, 1]` and on `[0, 2]`, and the multivariate + `(x−½)² + (y−½)²` on `[0, 1]²`), which a sign change can never reach. + **The certificate is not weakened**: a value is pinned to zero only by a + degenerate `[0, 0]` enclosure — an enclosure is a superset of the value, so + `[0, 0]` forces it — or by substituting the point's exact rational coordinates + and simplifying to the literal `0`, cross-checked against the enclosure exactly + as the removable-singularity path already is. An enclosure that merely + *contains* zero proves nothing and is not used: `exp(x) − 1 + 10⁻⁴⁰` on + `[0, 1]` has no root at all, but its value at `x = 0` is far below the width of + any enclosure computable there, and it stays `"undecided"` rather than being + claimed either way. - **`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 @@ -418,6 +467,50 @@ Both are detailed under *Behaviour changes to plan for*. taking one point as a single sequence — answers `f(1.0, 2.0)` and `f(1.0)` with that convention and a pointer to `numpy_eval` for batches. Exception types are unchanged. +- **The Cranelift backend was never linted, and had stopped being clean.** + `ci.yml` ran `cargo clippy -- -D warnings` for default / `egraph` / + `parallel` / `jit` / `groebner-cuda` but not for `cranelift`, so three + warnings accumulated in code that ships in the **default PyPI wheel** — and + `CONTRIBUTING.md` tells contributors to run `--all-features`, which therefore + failed on a clean checkout. A `cargo clippy (cranelift feature)` step now runs + alongside the others; it needs no system LLVM, because the backend is pure + Rust. The three lints are fixed rather than suppressed: `emit_eval_body` took + nine positional arguments and now takes an `EvalTarget` (root node, inputs, + pool) and an `InputLayout` enum (`Scalar { ptr }` / `Batch { ptr, point_idx, + n_points }`) — which also makes the half-specified batch layout a point index + with no point count, previously two independent `Option`s — unrepresentable; + and the two `return`s in `compile_jit_only`'s cranelift arm collapse to the + one-line `return` the other backends' arms already use. +- **`docs/mdbook/src/representations.md` documented a method that does not + exist, and five outputs that are not what the code prints.** The page showed + `p.leading_coeff()` on `UniPoly`, which had no leading-coefficient accessor at + all (it now has one — see *Added*); `sparse_interp_univariate(..., T=3)` and + `sparse_interp(..., T=2, D=5)`, whose parameters are `term_bound` and + `degree_bound`; and `p.to_symbolic(pool)`, under the claim that "all + specialized types can be converted back to a generic `Expr`" — no polynomial + type exposes a symbolic conversion to Python, so that whole section was + fiction. Every code block on the page is now executed against the built + extension, and the `# output` comments match what `print` actually emits + (`MultiPoly` prints in ascending exponent-vector order, an `ArbBall` prints as + `midpoint ± radius` rather than as an interval, and + `sparse_interp_univariate` returns a list of `(coefficient, exponent)` pairs + rather than a `MultiPolyFp`). The round-trip section now shows the conversion + that does exist, `GbPoly.to_expr()`, and says plainly that `UniPoly`, + `MultiPoly` and `RationalFunction` have none. +- **Four `examples/` scripts had rotted against the 3.9.0 API**, two of them + *silently* — they exited 0 while printing the wrong thing. + `phase3_polynomials.py` died at line 30 on `ExprPool.pow`, which no longer + exists on the Python side (use `a ** b`), and then on the trailing `pool` + argument that `UniPoly` / `MultiPoly` / `RationalFunction.from_symbolic` no + longer take, and on `UniPoly.pow(3)`. `agent_workflow.py` called `round()` on + the exact symbolic solutions `solve` now returns, and passes `numeric=True`. + The two silent ones: `risch_integration.py` claimed a non-elementary refusal + for `∫√(x³+1)dx`, which is genus 1 and now returns `EllipticF`, so the guard + moved to a genus-2 integrand that still raises `E-INT-004`; and + `lean_certificates.py` printed an empty Lean export, because `to_lean` withholds + a vacuous `e = e := rfl` when the simplifier found no rewrite. All 17 example + scripts now run end to end. `examples/` is not covered by CI's ruff or pytest, + which is why this went unnoticed. ### Added @@ -426,9 +519,10 @@ Both are detailed under *Behaviour changes to plan for*. 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 + Taylor model needs a rule with a rigorous Lagrange remainder, and eleven primitives have the first without the second — `erf`, `erfc`, `bessel_j0`, - `bessel_j1`, `digamma`, `lambert_w`, `acosh`, `asinh`, `floor`, `ceil`. So + `bessel_j1`, `digamma`, `lambert_w`, `acosh`, `asinh`, `atanh`, `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 @@ -453,7 +547,7 @@ Both are detailed under *Behaviour changes to plan for*. 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 + `numeric_ball` itself is *accurate* and stays as it is: those eleven 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 @@ -521,6 +615,14 @@ Both are detailed under *Behaviour changes to plan for*. differentiation rounds — and the equations it appended, plus the higher jets (`ddx/dt/dt`, …) they introduced, are now visible in `equations()` and `derivatives()`. +- **`UniPoly.leading_coeff`** — the leading (highest-degree) coefficient, `0` + for the zero polynomial so it pairs with `degree == -1`. A **property**, per + the accessor convention: it is a single FLINT coefficient read. It is also + *exact*, returned as a Python `int` of any size, which `coefficients()` is + not — that one is `i64` and truncates silently, so `3x² + 1` scaled by `2¹⁰⁰` + reports a leading coefficient of `1` through the list and the true value + through this accessor. Documented in `representations.md`, which had been + showing a `leading_coeff()` *method* that never existed. ### Performance @@ -572,6 +674,23 @@ Both are detailed under *Behaviour changes to plan for*. (`tests/silent_errors/`), still at **0 silent errors**. The new cases cover the PSLQ precision verdict as a *word* rather than a truthy value, so an `unknown` verdict cannot silently collapse into a pass. +- **`test_run_with_wall_fallback_bounds_a_cooperative_callee` no longer fails + because the box is busy.** Its bound read `elapsed_ms < 20 * 300` against a + call that measurably costs 2.7–5.3 s of work when it is behaving — a 13% + margin — so it went red repeatedly during saturated parallel runs and always + passed in isolation, which is the pattern that teaches people to ignore red. + The bound is now on **process CPU time**, which tracks the work the callee + actually did rather than the time the harness spent waiting for it. Wall time + here is `wall_ms` — a real-time timer, which does not stretch — plus the join + of a still-running worker, and only that second term was being measured + against load. Idle vs. 24 spinners on 12 cores: wall 5.3 s → 24.9 s (4× over + the old bound, a guaranteed failure), CPU 5.3 s → 8.8 s against a 60 s + ceiling. The + property is unchanged and still enforced from both ends — a callee that stops + seeing the budget burns CPU without limit and trips the assertion, and one + that stops coming back at all trips the `timeout(120)` marker — plus the test + now also pins that the call ended on the fallback's own join rather than on + some other check that raises the same code. ## 3.8.0 — 2026-08-12 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 0908791d..214317ed 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -46,13 +46,21 @@ See [`TESTING.md`](TESTING.md) for the full testing strategy (fuzzing, oracle cr ```bash cargo fmt --all -cargo clippy --all --all-features -- -D warnings +# One pass per feature set, which is what CI runs. `--all-features` is not the +# same thing and is not a substitute: it needs LLVM (`jit`) and the CUDA toolkit +# (`cuda`), so on a machine without them it fails to build rather than reporting +# lints, and the feature sets that *do* build stop being checked. +for f in "" egraph parallel cranelift jit groebner-cuda; do + cargo clippy --all-targets ${f:+--features $f} -- -D warnings +done uv run ruff check python/ tests/ --fix uv run ruff format python/ tests/ uv run ty check python/alkahest/ ``` -CI enforces all of the above on every PR. +CI enforces all of the above on every PR. Ruff and `ty` are scoped to `python/` +and `tests/` — `examples/`, `benchmarks/` and `scripts/` are not linted or run +by CI, so changes to a public API have to be swept through them by hand. ## Adding a new mathematical primitive diff --git a/alkahest-core/src/ball/mod.rs b/alkahest-core/src/ball/mod.rs index 1aee6dc3..3e1a40fa 100644 --- a/alkahest-core/src/ball/mod.rs +++ b/alkahest-core/src/ball/mod.rs @@ -64,9 +64,11 @@ use crate::kernel::expr::PredicateKind; use crate::kernel::{ExprData, ExprId, ExprPool}; +use crate::primitive::PrimitiveRegistry; use rug::{ops::Pow, Float}; use std::collections::HashMap; use std::fmt; +use std::sync::OnceLock; // --------------------------------------------------------------------------- // Precision constant @@ -765,20 +767,35 @@ impl ArbBall { } /// Bessel function of the first kind Jₙ(x) for integer order `n`. + /// + /// # Why this is midpoint + Lipschitz rather than an endpoint hull + /// + /// Jₙ oscillates, so `hull(Jₙ(lo), Jₙ(hi))` is **not** an enclosure of its + /// range: on `[-1, 1]` the two endpoints agree (`J₀(±1) ≈ 0.7652`) and the + /// hull collapses to a point that excludes `J₀(0) = 1`. An endpoint hull is + /// only valid for a monotone function, which every other kernel here that + /// uses one (`exp`, `log`, `sqrt`, `tanh`, the inverse trig family, + /// `digamma` between its poles, `floor`, `ceil`) is on its stated domain. + /// + /// The mean value theorem gives a sound enclosure instead: + /// `|Jₙ(x) − Jₙ(m)| ≤ L·|x − m|` with `L = sup|Jₙ′|`. For every integer `n` + /// and every real `x`, `|Jₙ(x)| ≤ 1`; with `J₀′ = −J₁` and + /// `Jₙ′ = (Jₙ₋₁ − Jₙ₊₁)/2` for `n ≥ 1` that yields `|Jₙ′| ≤ 1`, so `L = 1` + /// is rigorous at every order (and `J₋ₙ = (−1)ⁿ Jₙ` covers negative `n`). + /// The bound is loose — the true suprema are ≈ 0.582 for `J₀` and ≈ 0.582 + /// for `J₁` — but it is a *bound*, which the hull was not. pub fn bessel_jn(&self, n: i32) -> Self { let prec = self.prec; - let mut flo = Float::with_val(prec, self.lo().to_f64()); - flo.jn_mut(n); - let mut fhi = Float::with_val(prec, self.hi().to_f64()); - fhi.jn_mut(n); - let sum = Float::with_val(prec, &flo + &fhi); - let diff = Float::with_val(prec, &fhi - &flo); + // Evaluate at the midpoint at full working precision (MPFR's `jn` is + // correctly rounded, so the error is under half an ulp and is absorbed + // by `add_rounding_error` below). + let mut mid = Float::with_val(prec, &self.mid); + mid.jn_mut(n); let mut b = ArbBall { - mid: sum / 2_f64, - rad: diff / 2_f64, + mid, + rad: self.rad.clone(), prec, }; - // Endpoints are rounded to `prec`; see `exp`. b.add_rounding_error(); b } @@ -884,6 +901,29 @@ impl fmt::Display for AcbBall { // IntervalEval — expression evaluator using ArbBall // --------------------------------------------------------------------------- +/// Lazily-initialised, process-wide [`PrimitiveRegistry`] used to evaluate +/// `ExprData::Func` nodes. +/// +/// A singleton, because building it is far too much work to repeat per `Func` +/// node; the same pattern is used by the JIT's tree-walking interpreter +/// (`crate::jit::registry`). +/// +/// `dispatch_registry`, not `default_registry`: this path calls `numeric_ball` +/// on the primitive and treats `None` as unsupported, so it never reads a +/// capability bit — and probing 41 primitives across six argument shapes cost +/// ~1.2 ms of one-time work on whichever call touched the registry first, +/// against 0.02 ms steady state. That is invisible to a wall-clock test and +/// very visible to an instruction-counting one, which is how CodSpeed caught it +/// as a 24x regression on `test_ball_sin_cos_eps1e2`. +/// +/// The anti-drift guarantee is unchanged, and in fact stronger: the set of +/// functions accepted here is the set whose `numeric_ball` kernel returns +/// `Some`, which is ground truth rather than a probed summary of it. +fn registry() -> &'static PrimitiveRegistry { + static REGISTRY: OnceLock = OnceLock::new(); + REGISTRY.get_or_init(PrimitiveRegistry::dispatch_registry) +} + /// Evaluates a symbolic expression using rigorous ball arithmetic. /// /// Each variable can be bound to an `ArbBall` interval. The result is an @@ -1020,40 +1060,24 @@ impl IntervalEval { } Some(b.pow_f(&e)) } - ExprData::Func { name, args } if args.len() == 1 => { - let x = self.eval_node(args[0], pool)?; - match name.as_str() { - "sin" => Some(x.sin()), - "cos" => Some(x.cos()), - "exp" => Some(x.exp()), - "log" => x.log(), - "sqrt" => x.sqrt(), - // Every arm below is an already-implemented, outward-rounded - // `ArbBall` kernel that was simply not reachable from an - // expression. Leaving them unwired made `eval` answer `None` - // for `tan(x)` and friends, which callers that read `None` - // as "cannot bound this" — the zero test in - // `crate::matrix`, `crate::validated::bounds` — had to treat - // as a refusal. - "tan" => x.tan(), - "sinh" => Some(x.sinh()), - "cosh" => Some(x.cosh()), - "tanh" => Some(x.tanh()), - "asin" => x.asin(), - "acos" => x.acos(), - "atan" => Some(x.atan()), - "asinh" => Some(x.asinh()), - "acosh" => x.acosh(), - "atanh" => x.atanh(), - "erf" => Some(x.erf()), - "erfc" => Some(x.erfc()), - "abs" => Some(x.abs_ball()), - "floor" => Some(x.floor_ball()), - "ceil" => Some(x.ceil_ball()), - "digamma" => x.digamma(), - "lambert_w" => x.lambert_w0(), - _ => None, + // Every `Func` node is dispatched through the primitive registry's + // `numeric_ball` slot, so the set of functions interval evaluation + // accepts *is* the set the registry advertises as + // `Capabilities::NUMERIC_BALL` — by construction, not by + // agreement. The hand-written match this replaces had drifted: + // `bessel_j0`/`bessel_j1` had outward-rounded ball kernels and a + // `numeric_ball: true` capability bit, and were still refused here + // with `E-EVAL-010` for no reason anyone had recorded. + // + // Arity is the primitive's own business: each kernel declines an + // argument list it does not have a rule for (see `builtins::unary`), + // which is what makes dispatching the whole list safe. + ExprData::Func { name, args } if !args.is_empty() => { + let mut vals = Vec::with_capacity(args.len()); + for &a in &args { + vals.push(self.eval_node(a, pool)?); } + registry().numeric_ball(&name, &vals) } ExprData::Piecewise { branches, default } => { for (c, v) in branches { @@ -1250,6 +1274,80 @@ mod tests { assert!(r.rad_f64() < 1e-30, "rad={}", r.rad_f64()); } + /// The dispatch is *derived* from the registry, so this cannot drift the + /// way the hand-written match did. Pin it anyway: the property that used to + /// be violated (a primitive advertising `numeric_ball` that interval + /// evaluation refuses) is the whole point, and a future refactor that + /// reintroduces a list here fails this test rather than a user's proof. + #[test] + fn interval_eval_accepts_every_primitive_that_advertises_numeric_ball() { + use crate::primitive::{Capabilities, PrimitiveRegistry}; + let reg = PrimitiveRegistry::default_registry(); + let pool = p(); + let x = pool.symbol("x", Domain::Real); + // No single probe point is in every kernel's domain — `acosh` needs + // `x ≥ 1` and `atanh` needs `|x| < 1` — so a primitive counts as + // reachable if *some* probe evaluates. A domain refusal is a different + // answer from "I have no rule for this name", and only the second is + // the coverage gap under test. + let probes = [1.5_f64, 0.5]; + let refused: Vec = reg + .iter() + .filter(|(_, caps)| caps.contains(Capabilities::NUMERIC_BALL)) + .filter(|(name, _)| { + let call = pool.func(*name, vec![x]); + probes.iter().all(|&probe| { + let mut ev = IntervalEval::new(128); + ev.bind(x, ArbBall::from_f64(probe, 128)); + ev.eval(call, &pool).is_none() + }) + }) + .map(|(name, _)| name.to_string()) + .collect(); + assert!( + refused.is_empty(), + "these primitives advertise numeric_ball but interval evaluation \ + refuses them: {refused:?}" + ); + } + + /// The two the audit was opened for, pinned by name and by value. + #[test] + fn interval_eval_evaluates_bessel() { + let pool = p(); + let x = pool.symbol("x", Domain::Real); + let mut ev = IntervalEval::new(128); + ev.bind(x, ArbBall::from_f64(1.0, 128)); + let j0 = ev.eval(pool.func("bessel_j0", vec![x]), &pool).unwrap(); + let j1 = ev.eval(pool.func("bessel_j1", vec![x]), &pool).unwrap(); + // The enclosure at an exact point is far tighter than an `f64` literal, + // so compare midpoints rather than asking it to `contain` one. + assert!( + (j0.mid_f64() - 0.765_197_686_557_966_5).abs() < 1e-15, + "J0(1) = {j0}" + ); + assert!( + (j1.mid_f64() - 0.440_050_585_744_933_5).abs() < 1e-15, + "J1(1) = {j1}" + ); + assert!(j0.rad_f64() < 1e-30 && j1.rad_f64() < 1e-30); + } + + /// A unary kernel handed the wrong number of arguments must decline, not + /// quietly bound the first one: a rigorous enclosure of the wrong function + /// is the worst answer this module can give. + #[test] + fn interval_eval_refuses_a_unary_primitive_at_the_wrong_arity() { + let pool = p(); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let mut ev = IntervalEval::new(128); + ev.bind(x, ArbBall::from_f64(1.0, 128)); + ev.bind(y, ArbBall::from_f64(2.0, 128)); + assert!(ev.eval(pool.func("sin", vec![x, y]), &pool).is_none()); + assert!(ev.eval(pool.func("sin", vec![]), &pool).is_none()); + } + #[test] fn acb_modulus() { // |3 + 4i| = 5 @@ -1302,6 +1400,60 @@ mod rounding_soundness_tests { ); } + /// Jₙ oscillates, so the endpoint hull `bessel_jn` used to take was not an + /// enclosure: on `[-1, 1]` both endpoints give `J₀(±1) ≈ 0.7652`, the hull + /// collapsed to that point, and `J₀(0) = 1` — the maximum of the function — + /// fell outside the "rigorous" ball. + #[test] + fn bessel_encloses_an_interior_extremum() { + let b = ArbBall::from_midpoint_radius(0.0, 1.0, PREC); // [-1, 1] + let j0 = b.bessel_jn(0); + assert!(j0.contains(1.0), "J0 on [-1,1] misses its maximum: {j0}"); + assert!(j0.contains(0.7651976865579666), "{j0}"); + + // Same failure one period out, where the endpoints straddle a zero + // rather than a peak. + let b = ArbBall::from_midpoint_radius(3.0, 1.0, PREC); // [2, 4] + let j1 = b.bessel_jn(1); + for x in [2.0_f64, 2.5, 3.0, 3.5, 4.0] { + let mut v = Float::with_val(PREC, x); + v.jn_mut(1); + assert!(j1.contains(v.to_f64()), "J1({x}) outside {j1}"); + } + } + + /// Randomised enclosure check: sampling the true function inside the ball + /// must never escape the reported enclosure. + #[test] + fn bessel_enclosure_holds_on_random_intervals() { + let mut seed = 0x2545_F491_4F6C_DD1D_u64; + let mut next = move || { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + (seed >> 11) as f64 / (1u64 << 53) as f64 + }; + for _ in 0..200 { + let centre = (next() - 0.5) * 40.0; + let radius = next() * 5.0; + for n in [0_i32, 1] { + let ball = ArbBall::from_midpoint_radius(centre, radius, PREC); + let out = ball.bessel_jn(n); + for k in 0..=20 { + let x = centre - radius + 2.0 * radius * (k as f64 / 20.0); + let mut v = Float::with_val(PREC, x); + v.jn_mut(n); + assert!( + out.contains(v.to_f64()), + "J{n}({x}) = {v} escapes {out} for [{}, {}]", + centre - radius, + centre + radius + ); + } + } + } + } + /// The radius must stay at the working-precision scale, not balloon. /// /// Soundness is trivially achievable by making every ball enormous; this diff --git a/alkahest-core/src/jit/cranelift_backend.rs b/alkahest-core/src/jit/cranelift_backend.rs index 3185cd47..3477f029 100644 --- a/alkahest-core/src/jit/cranelift_backend.rs +++ b/alkahest-core/src/jit/cranelift_backend.rs @@ -224,52 +224,75 @@ fn codegen_node( // Load input variables into the values map (scalar or batch layout) // --------------------------------------------------------------------------- +/// What is being compiled: the root node, its input variables, and the pool +/// they all live in. These three always travel together, so they are passed +/// as one value rather than as three positional arguments. +#[derive(Clone, Copy)] +struct EvalTarget<'a> { + expr: ExprId, + inputs: &'a [ExprId], + pool: &'a ExprPool, +} + +/// Where the emitted body reads variable `i` from. +/// +/// The scalar entry point reads it at `ptr[i]`; the bulk entry point reads it +/// at `ptr[i * n_points + point_idx]`. Modelling this as an enum instead of +/// three loose `Option`s makes the half-specified batch layout (a point index +/// with no point count, or vice versa) unrepresentable. +#[derive(Clone, Copy)] +enum InputLayout { + /// One point, variables contiguous from `ptr`. + Scalar { ptr: cranelift_codegen::ir::Value }, + /// `n_points` points, variable-major; emitted inside the bulk loop body. + Batch { + ptr: cranelift_codegen::ir::Value, + point_idx: cranelift_codegen::ir::Value, + n_points: cranelift_codegen::ir::Value, + }, +} + fn load_input_vars( builder: &mut FunctionBuilder, - inputs_ptr: cranelift_codegen::ir::Value, inputs: &[ExprId], values: &mut HashMap, - point_idx: Option, - n_points: Option, + layout: InputLayout, ) { for (i, &var) in inputs.iter().enumerate() { - let val = if let (Some(idx), Some(n_pts)) = (point_idx, n_points) { - let var_i = builder.ins().iconst(types::I64, i as i64); - let stride = builder.ins().imul(var_i, n_pts); - let elem = builder.ins().iadd(stride, idx); - let byte_off = builder.ins().imul_imm(elem, 8); - let addr = builder.ins().iadd(inputs_ptr, byte_off); - builder.ins().load(types::F64, MemFlags::trusted(), addr, 0) - } else { - let byte_offset = (i * std::mem::size_of::()) as i32; - builder - .ins() - .load(types::F64, MemFlags::trusted(), inputs_ptr, byte_offset) + let val = match layout { + InputLayout::Batch { + ptr, + point_idx, + n_points, + } => { + let var_i = builder.ins().iconst(types::I64, i as i64); + let stride = builder.ins().imul(var_i, n_points); + let elem = builder.ins().iadd(stride, point_idx); + let byte_off = builder.ins().imul_imm(elem, 8); + let addr = builder.ins().iadd(ptr, byte_off); + builder.ins().load(types::F64, MemFlags::trusted(), addr, 0) + } + InputLayout::Scalar { ptr } => { + let byte_offset = (i * std::mem::size_of::()) as i32; + builder + .ins() + .load(types::F64, MemFlags::trusted(), ptr, byte_offset) + } }; values.insert(var, val); } } fn emit_eval_body( - expr: ExprId, - inputs: &[ExprId], - pool: &ExprPool, + target: EvalTarget<'_>, builder: &mut FunctionBuilder, module: &mut JITModule, math: &MathFuncIds, - inputs_ptr: cranelift_codegen::ir::Value, - point_idx: Option, - n_points: Option, + layout: InputLayout, ) -> Result { + let EvalTarget { expr, inputs, pool } = target; let mut values: HashMap = HashMap::new(); - load_input_vars( - builder, - inputs_ptr, - inputs, - &mut values, - point_idx, - n_points, - ); + load_input_vars(builder, inputs, &mut values, layout); let topo = topo_sort(expr, pool); for &node in &topo { if values.contains_key(&node) { @@ -297,6 +320,10 @@ pub fn compile_cranelift( inputs: &[ExprId], pool: &ExprPool, ) -> Result { + // Both entry points below emit the same DAG against the same pool; only the + // input layout differs. + let target = EvalTarget { expr, inputs, pool }; + // ------------------------------------------------------------------ // 1. ISA — native host architecture, speed optimised // ------------------------------------------------------------------ @@ -409,15 +436,11 @@ pub fn compile_cranelift( builder.seal_block(block); let inputs_ptr = builder.block_params(block)[0]; let result = emit_eval_body( - expr, - inputs, - pool, + target, &mut builder, &mut module, &math, - inputs_ptr, - None, - None, + InputLayout::Scalar { ptr: inputs_ptr }, )?; builder.ins().return_(&[result]); builder.finalize(); @@ -471,15 +494,15 @@ pub fn compile_cranelift( builder.switch_to_block(loop_body); let result = emit_eval_body( - expr, - inputs, - pool, + target, &mut builder, &mut module, &math, - bulk_inputs_ptr, - Some(loop_idx), - Some(bulk_n_points), + InputLayout::Batch { + ptr: bulk_inputs_ptr, + point_idx: loop_idx, + n_points: bulk_n_points, + }, )?; let out_byte_off = builder.ins().imul_imm(loop_idx, 8); let out_addr = builder.ins().iadd(bulk_outputs_ptr, out_byte_off); diff --git a/alkahest-core/src/jit/mod.rs b/alkahest-core/src/jit/mod.rs index 6e345baa..e75f58f6 100644 --- a/alkahest-core/src/jit/mod.rs +++ b/alkahest-core/src/jit/mod.rs @@ -544,10 +544,7 @@ pub fn compile_jit_only( pool: &ExprPool, ) -> Result { #[cfg(feature = "cranelift")] - match cranelift_backend::compile_cranelift(expr, inputs, pool) { - Ok(f) => return Ok(f), - Err(e) => return Err(e), - } + return cranelift_backend::compile_cranelift(expr, inputs, pool); #[cfg(all(feature = "jit", not(feature = "cranelift")))] return compile_llvm(expr, inputs, pool); diff --git a/alkahest-core/src/poly/unipoly.rs b/alkahest-core/src/poly/unipoly.rs index c3a7a431..25c65afd 100644 --- a/alkahest-core/src/poly/unipoly.rs +++ b/alkahest-core/src/poly/unipoly.rs @@ -388,6 +388,12 @@ impl UniPoly { .collect() } + /// Leading (highest-degree) coefficient, exact. `0` for the zero + /// polynomial, matching `degree() == -1` there. + pub fn leading_coeff(&self) -> rug::Integer { + self.coeffs.leading_coeff_fmpz().to_rug() + } + pub fn degree(&self) -> i64 { self.coeffs.degree() } @@ -652,6 +658,26 @@ mod tests { assert_eq!(poly.coefficients_i64(), vec![1, 2, 1]); } + #[test] + fn leading_coeff_is_exact_and_zero_for_the_zero_polynomial() { + let (p, x) = pool_and_var(); + // 3*x^2 + 1 + let xsq = p.pow(x, p.integer(2_i32)); + let expr = p.add(vec![p.mul(vec![p.integer(3_i32), xsq]), p.integer(1_i32)]); + let poly = UniPoly::from_symbolic(expr, x, &p).unwrap(); + assert_eq!(poly.leading_coeff(), rug::Integer::from(3)); + + // Beyond i64, where `coefficients_i64` truncates. + let big: rug::Integer = rug::Integer::from(1) << 100; + let big_expr = p.mul(vec![p.integer(big.clone()), xsq]); + let big_poly = UniPoly::from_symbolic(big_expr, x, &p).unwrap(); + assert_eq!(big_poly.leading_coeff(), big); + + let zero = UniPoly::from_symbolic(p.integer(0_i32), x, &p).unwrap(); + assert_eq!(zero.degree(), -1); + assert_eq!(zero.leading_coeff(), rug::Integer::ZERO); + } + #[test] fn from_symbolic_constant() { let (p, x) = pool_and_var(); diff --git a/alkahest-core/src/primitive/mod.rs b/alkahest-core/src/primitive/mod.rs index bcc1329a..1c629da3 100644 --- a/alkahest-core/src/primitive/mod.rs +++ b/alkahest-core/src/primitive/mod.rs @@ -271,6 +271,25 @@ impl PrimitiveRegistry { self.map.insert(name, Entry { primitive: p, caps }); } + /// Register without probing capabilities. + /// + /// For registries used purely to *dispatch* — where the caller invokes a + /// slot such as `numeric_ball` and takes `None` as "unsupported" — the + /// probed bits are never read, and probing 41 primitives across six + /// argument shapes costs ~1.2 ms of one-time work on whatever call happens + /// to touch the registry first. `capabilities()` on such a registry reports + /// `empty()`; use [`Self::default_registry`] if you need the bits. + pub fn register_unprobed(&mut self, p: Box) { + let name = p.name(); + self.map.insert( + name, + Entry { + primitive: p, + caps: Capabilities::empty(), + }, + ); + } + /// Look up a primitive by name. pub fn get(&self, name: &str) -> Option<&dyn Primitive> { self.map.get(name).map(|e| &*e.primitive) @@ -339,52 +358,75 @@ impl PrimitiveRegistry { /// Return a registry pre-populated with Alkahest's built-in primitives. pub fn default_registry() -> Self { + Self::build(true) + } + + /// Every built-in primitive, registered without capability probing. + /// + /// Same dispatch behaviour as [`Self::default_registry`] — identical names, + /// identical primitives — but `capabilities()` reports `empty()`. Use it + /// when you only call a primitive slot and treat `None` as unsupported. + pub fn dispatch_registry() -> Self { + Self::build(false) + } + + fn build(probe: bool) -> Self { let mut reg = Self::new(); - reg.register(Box::new(builtins::SinPrimitive)); - reg.register(Box::new(builtins::CosPrimitive)); - reg.register(Box::new(builtins::ExpPrimitive)); - reg.register(Box::new(builtins::LogPrimitive)); - reg.register(Box::new(builtins::SqrtPrimitive)); + reg.register_unprobed(Box::new(builtins::SinPrimitive)); + reg.register_unprobed(Box::new(builtins::CosPrimitive)); + reg.register_unprobed(Box::new(builtins::ExpPrimitive)); + reg.register_unprobed(Box::new(builtins::LogPrimitive)); + reg.register_unprobed(Box::new(builtins::SqrtPrimitive)); // V1-12: expanded registry - reg.register(Box::new(builtins::TanPrimitive)); - reg.register(Box::new(builtins::SinhPrimitive)); - reg.register(Box::new(builtins::CoshPrimitive)); - reg.register(Box::new(builtins::TanhPrimitive)); - reg.register(Box::new(builtins::AsinPrimitive)); - reg.register(Box::new(builtins::AcosPrimitive)); - reg.register(Box::new(builtins::AtanPrimitive)); - reg.register(Box::new(builtins::AsinhPrimitive)); - reg.register(Box::new(builtins::AcoshPrimitive)); - reg.register(Box::new(builtins::AtanhPrimitive)); - reg.register(Box::new(builtins::ErfPrimitive)); - reg.register(Box::new(builtins::ErfcPrimitive)); + reg.register_unprobed(Box::new(builtins::TanPrimitive)); + reg.register_unprobed(Box::new(builtins::SinhPrimitive)); + reg.register_unprobed(Box::new(builtins::CoshPrimitive)); + reg.register_unprobed(Box::new(builtins::TanhPrimitive)); + reg.register_unprobed(Box::new(builtins::AsinPrimitive)); + reg.register_unprobed(Box::new(builtins::AcosPrimitive)); + reg.register_unprobed(Box::new(builtins::AtanPrimitive)); + reg.register_unprobed(Box::new(builtins::AsinhPrimitive)); + reg.register_unprobed(Box::new(builtins::AcoshPrimitive)); + reg.register_unprobed(Box::new(builtins::AtanhPrimitive)); + reg.register_unprobed(Box::new(builtins::ErfPrimitive)); + reg.register_unprobed(Box::new(builtins::ErfcPrimitive)); // Elliptic special functions (parameter convention m = k²). - reg.register(Box::new(builtins::EllipticKPrimitive)); - reg.register(Box::new(builtins::EllipticEPrimitive)); - reg.register(Box::new(builtins::EllipticFPrimitive)); - reg.register(Box::new(builtins::EllipticPiPrimitive)); - reg.register(Box::new(builtins::AbsPrimitive)); - reg.register(Box::new(builtins::SignPrimitive)); - reg.register(Box::new(builtins::HeavisidePrimitive)); - reg.register(Box::new(builtins::DiracDeltaPrimitive)); - reg.register(Box::new(builtins::FloorPrimitive)); - reg.register(Box::new(builtins::CeilPrimitive)); - reg.register(Box::new(builtins::RoundPrimitive)); - reg.register(Box::new(builtins::Atan2Primitive)); - reg.register(Box::new(builtins::GammaPrimitive)); - reg.register(Box::new(builtins::LambertWPrimitive)); - reg.register(Box::new(builtins::DigammaPrimitive)); - reg.register(Box::new(builtins::BesselJ0Primitive)); - reg.register(Box::new(builtins::BesselJ1Primitive)); - reg.register(Box::new(builtins::MinPrimitive)); - reg.register(Box::new(builtins::MaxPrimitive)); - reg.register(Box::new(builtins::ConjugatePrimitive)); - reg.register(Box::new(builtins::RePrimitive)); - reg.register(Box::new(builtins::ImPrimitive)); - reg.register(Box::new(builtins::ArgPrimitive)); + reg.register_unprobed(Box::new(builtins::EllipticKPrimitive)); + reg.register_unprobed(Box::new(builtins::EllipticEPrimitive)); + reg.register_unprobed(Box::new(builtins::EllipticFPrimitive)); + reg.register_unprobed(Box::new(builtins::EllipticPiPrimitive)); + reg.register_unprobed(Box::new(builtins::AbsPrimitive)); + reg.register_unprobed(Box::new(builtins::SignPrimitive)); + reg.register_unprobed(Box::new(builtins::HeavisidePrimitive)); + reg.register_unprobed(Box::new(builtins::DiracDeltaPrimitive)); + reg.register_unprobed(Box::new(builtins::FloorPrimitive)); + reg.register_unprobed(Box::new(builtins::CeilPrimitive)); + reg.register_unprobed(Box::new(builtins::RoundPrimitive)); + reg.register_unprobed(Box::new(builtins::Atan2Primitive)); + reg.register_unprobed(Box::new(builtins::GammaPrimitive)); + reg.register_unprobed(Box::new(builtins::LambertWPrimitive)); + reg.register_unprobed(Box::new(builtins::DigammaPrimitive)); + reg.register_unprobed(Box::new(builtins::BesselJ0Primitive)); + reg.register_unprobed(Box::new(builtins::BesselJ1Primitive)); + reg.register_unprobed(Box::new(builtins::MinPrimitive)); + reg.register_unprobed(Box::new(builtins::MaxPrimitive)); + reg.register_unprobed(Box::new(builtins::ConjugatePrimitive)); + reg.register_unprobed(Box::new(builtins::RePrimitive)); + reg.register_unprobed(Box::new(builtins::ImPrimitive)); + reg.register_unprobed(Box::new(builtins::ArgPrimitive)); + if probe { + reg.probe_all(); + } reg } + /// Fill in every entry's capability bits, probing each primitive. + fn probe_all(&mut self) { + for entry in self.map.values_mut() { + entry.caps = probe_caps(&*entry.primitive); + } + } + /// Returns true if a primitive with this name is registered. pub fn is_registered(&self, name: &str) -> bool { self.map.contains_key(name) @@ -434,10 +476,18 @@ fn probe_caps(p: &dyn Primitive) -> Capabilities { } } - let ball1 = [ArbBall::from_f64(1.0, 128)]; - let ball2 = [ArbBall::from_f64(1.0, 128), ArbBall::from_f64(2.0, 128)]; - if p.numeric_ball(&ball1).is_some() || p.numeric_ball(&ball2).is_some() { - caps |= Capabilities::NUMERIC_BALL; + // The *same* probe sets as `numeric_f64`, deliberately: a primitive whose + // two kernels cover the same domain must not be reported as having one and + // not the other. Probing balls at `1.0` alone did exactly that to `atanh`, + // whose domain is the open interval `(-1, 1)` — it declined the probe point + // and lost a `numeric_ball` bit it had earned, while keeping `numeric_f64` + // because that probe already tried `0.5`. + for args in probe_f64_sets { + let balls: Vec = args.iter().map(|&v| ArbBall::from_f64(v, 128)).collect(); + if p.numeric_ball(&balls).is_some() { + caps |= Capabilities::NUMERIC_BALL; + break; + } } // diff_forward / diff_reverse / simplify: probe with a fresh pool @@ -514,6 +564,22 @@ pub mod builtins { use crate::kernel::expr::ExprData; use crate::kernel::{ExprId, ExprPool}; + /// The single argument of a unary primitive, or `None` if the caller + /// supplied a different number of them. + /// + /// `numeric_ball` is dispatched generically — [`crate::ball::IntervalEval`] + /// hands the registry whatever argument list the `Func` node carries — so a + /// unary kernel that reached for `args[0]` unconditionally would answer + /// `sin(x, y)` with `sin(x)`: a *confidently wrong* rigorous enclosure, the + /// one failure mode this subsystem must not have. Declining the arity is + /// the honest answer, and it costs the caller only an `E-EVAL-010`. + fn unary(args: &[ArbBall]) -> Option<&ArbBall> { + match args { + [only] => Some(only), + _ => None, + } + } + macro_rules! symbolic_complex_primitive { ($type:ident, $name:literal) => { pub struct $type; @@ -568,7 +634,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].sin()) + Some(unary(args)?.sin()) } fn lean_theorem(&self) -> Option<&'static str> { @@ -614,7 +680,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].cos()) + Some(unary(args)?.cos()) } fn lean_theorem(&self) -> Option<&'static str> { @@ -658,7 +724,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].exp()) + Some(unary(args)?.exp()) } fn lean_theorem(&self) -> Option<&'static str> { @@ -703,7 +769,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].log() + unary(args)?.log() } fn lean_theorem(&self) -> Option<&'static str> { @@ -758,7 +824,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].sqrt() + unary(args)?.sqrt() } fn lean_theorem(&self) -> Option<&'static str> { @@ -812,7 +878,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].tan() + unary(args)?.tan() } fn lean_theorem(&self) -> Option<&'static str> { @@ -861,7 +927,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].sinh()) + Some(unary(args)?.sinh()) } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -904,7 +970,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].cosh()) + Some(unary(args)?.cosh()) } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -956,7 +1022,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].tanh()) + Some(unary(args)?.tanh()) } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -1010,7 +1076,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].asin() + unary(args)?.asin() } } @@ -1063,7 +1129,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].acos() + unary(args)?.acos() } } @@ -1110,7 +1176,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].atan()) + Some(unary(args)?.atan()) } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -1162,7 +1228,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].asinh()) + Some(unary(args)?.asinh()) } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -1214,7 +1280,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].acosh() + unary(args)?.acosh() } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -1264,7 +1330,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - args[0].atanh() + unary(args)?.atanh() } // NOTE: no `lean_theorem` override — see the `tan` primitive above @@ -1314,7 +1380,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].erf()) + Some(unary(args)?.erf()) } } @@ -1360,7 +1426,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].erfc()) + Some(unary(args)?.erfc()) } } @@ -1918,7 +1984,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].abs_ball()) + Some(unary(args)?.abs_ball()) } } @@ -2030,7 +2096,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].floor_ball()) + Some(unary(args)?.floor_ball()) } } @@ -2052,7 +2118,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - Some(args[0].ceil_ball()) + Some(unary(args)?.ceil_ball()) } } @@ -2185,11 +2251,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - if args.len() == 1 { - args[0].lambert_w0() - } else { - None - } + unary(args)?.lambert_w0() } } @@ -2223,11 +2285,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - if args.len() == 1 { - args[0].digamma() - } else { - None - } + unary(args)?.digamma() } } @@ -2271,11 +2329,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - if args.len() == 1 { - Some(args[0].bessel_jn(0)) - } else { - None - } + Some(unary(args)?.bessel_jn(0)) } } @@ -2323,11 +2377,7 @@ pub mod builtins { } fn numeric_ball(&self, args: &[ArbBall]) -> Option { - if args.len() == 1 { - Some(args[0].bessel_jn(1)) - } else { - None - } + Some(unary(args)?.bessel_jn(1)) } } diff --git a/alkahest-core/src/primitive/taylor_support.rs b/alkahest-core/src/primitive/taylor_support.rs index fdfcdff7..ee791b10 100644 --- a/alkahest-core/src/primitive/taylor_support.rs +++ b/alkahest-core/src/primitive/taylor_support.rs @@ -287,9 +287,13 @@ mod tests { continue; } differing += 1; - let arg = [ArbBall::from_f64(1.0, 128)]; + // Two probe points, because no single one is inside every domain: + // `acosh` needs `x ≥ 1` and `atanh` needs `|x| < 1`. Declining an + // out-of-domain argument is not the failure under test. assert!( - reg.numeric_ball(name, &arg).is_some(), + [1.0_f64, 0.5].into_iter().any(|v| reg + .numeric_ball(name, &[ArbBall::from_f64(v, 128)]) + .is_some()), "`{name}` advertises numeric_ball but has none" ); } diff --git a/alkahest-core/src/validated/bounds.rs b/alkahest-core/src/validated/bounds.rs index 8beaafbc..ae883141 100644 --- a/alkahest-core/src/validated/bounds.rs +++ b/alkahest-core/src/validated/bounds.rs @@ -661,17 +661,35 @@ fn split_quotient(expr: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId)> { /// singularity test has to go through the symbolic path; anything short of an /// exact zero is treated as "not removable", which is the safe direction. fn vanishes_exactly(expr: ExprId, pool: &ExprPool, var: ExprId, at: &Float) -> bool { - let Some(rational) = at.to_rational() else { - return false; - }; - let (n, d) = rational.into_numer_denom(); - let point = if d == 1 { - pool.integer(n) - } else { - pool.rational(n, d) - }; + vanishes_exactly_at(expr, pool, &[(var, at.clone(), at.clone())]) +} + +/// `expr` with **every** variable of the degenerate box `point` replaced by its +/// exact rational coordinate, simplified. +/// +/// Returns `true` only when the result is the literal integer (or rational) +/// zero — a symbolic proof that `expr` vanishes at that point, of the same kind +/// [`vanishes_exactly`] provides in one dimension. A `Float` coordinate is a +/// binary rational and converts exactly, so nothing is rounded on the way in; +/// a coordinate whose interval is not degenerate names a set rather than a +/// point and is declined. +fn vanishes_exactly_at(expr: ExprId, pool: &ExprPool, point: &[FBox]) -> bool { let mut mapping = HashMap::new(); - mapping.insert(var, point); + for (var, lo, hi) in point { + if lo != hi { + return false; + } + let Some(rational) = lo.to_rational() else { + return false; + }; + let (n, d) = rational.into_numer_denom(); + let value = if d == 1 { + pool.integer(n) + } else { + pool.rational(n, d) + }; + mapping.insert(*var, value); + } let substituted = subs(expr, &mapping, pool); let reduced = simplify(substituted, pool).value; match pool.get(reduced) { @@ -1120,23 +1138,85 @@ impl SignWitnesses { } } -/// Rigorously evaluate `expr` at one point of the box and report its sign, or -/// `None` when the enclosure straddles zero or the evaluation refused. +/// True only for the degenerate enclosure `[0, 0]`. +/// +/// An enclosure is a superset of the value it describes, so `[0, 0]` — and +/// *only* `[0, 0]` — pins that value to zero. An enclosure that merely +/// *contains* zero proves nothing about whether the value is zero, which is why +/// this is a two-sided test on the outward-rounded [`lb`]/[`ub`] rather than +/// [`contains_zero`]. +fn is_exact_zero(b: &ArbBall) -> bool { + lb(b) == 0 && ub(b) == 0 +} + +/// What a rigorous evaluation at a single point of the box established there. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum PointOutcome { + /// The value at that point is **proven** to be exactly zero, so the point + /// is a root — and being a point of the (closed) box, it is a root *in* the + /// box. + Zero, + /// The value is proven strictly positive (`true`) or strictly negative + /// (`false`). + Sign(bool), + /// Nothing was proven: the enclosure straddles zero without pinning it, or + /// the evaluation refused. + Unknown, +} + +/// Rigorously evaluate `expr` at one point of the box and report what that +/// proves. /// /// The point is passed as a degenerate box, so the Taylor model collapses to a /// ball evaluation and the sign test goes through [`lb`]/[`ub`], which round /// outward. A reported sign is therefore a proof, never a rounding artefact. -fn point_sign( +/// +/// [`PointOutcome::Zero`] is likewise only ever a proof, by one of two routes: +/// +/// * the enclosure is the degenerate interval `[0, 0]`, which forces the value +/// it encloses to be zero — this settles the exactly-representable cases +/// (`x` at `0`, `x − 1` at `1`, `(x − 1)²` at `1`, `sin x` at `0`), where the +/// arithmetic is exact and no rounding term is ever added because every +/// intermediate midpoint is zero; or +/// * substituting the point's exact rational coordinates and simplifying lands +/// on the literal `0` — the same symbolic argument the removable-singularity +/// path already relies on, and the only kind of argument that can prove a +/// transcendental combination such as `exp(x) − 1` vanishes at `0`, where the +/// enclosure is `[-ε, ε]` and can never be tightened to a point. +/// +/// The symbolic route is cross-checked against the enclosure exactly as +/// [`enclosure_admits_zero`] does elsewhere: a `simplify` that claims an exact +/// zero which outward-rounded ball arithmetic contradicts would be a simplifier +/// bug, and the safe response is to prove nothing rather than to certify on top +/// of it. That check is free here — reaching it means [`determined_sign`] +/// already declined, which is precisely `contains_zero`. +/// +/// `symbolic` selects whether the second route is attempted. It substitutes and +/// simplifies, which interns a fresh copy of `expr` in the pool, so callers +/// that run this once per bisection leave it off; the box's own distinguished +/// points — endpoints, corners, centre — are where a root "sitting on the box" +/// can be, and they are a fixed, small set. +fn point_outcome( expr: ExprId, pool: &ExprPool, point: &[FBox], order: usize, prec: u32, -) -> Option { - match taylor_range(expr, pool, point, order, prec) { - Ok(r) => determined_sign(&r), - Err(_) => None, + symbolic: bool, +) -> PointOutcome { + let Ok(r) = taylor_range(expr, pool, point, order, prec) else { + return PointOutcome::Unknown; + }; + if let Some(s) = determined_sign(&r) { + return PointOutcome::Sign(s); + } + if is_exact_zero(&r) { + return PointOutcome::Zero; + } + if symbolic && contains_zero(&r) && vanishes_exactly_at(expr, pool, point) { + return PointOutcome::Zero; } + PointOutcome::Unknown } /// Degenerate boxes for the centre, the per-axis endpoints and (in low @@ -1186,10 +1266,22 @@ fn seed_points(boxes0: &[FBox], prec: u32) -> Vec> { .collect() } -/// Search for a proof that `expr` has a root in the box, by finding one point -/// where it is provably positive and one where it is provably negative. +/// Search for a proof that `expr` has a root in the box. /// -/// # Why this is a proof +/// Two independent proofs are looked for, and either one suffices: +/// +/// 1. **A point proven to be a root.** The box is closed, so a point of it at +/// which `expr` is *proven* zero — by a degenerate `[0, 0]` enclosure or by +/// exact symbolic substitution, see [`point_outcome`] — is a root in the +/// box, full stop. No continuity argument, no sign change, and no relation +/// between that point and any other is needed. This is what settles a root +/// sitting exactly on an endpoint (`x` on `[0, 1]`), where the function has +/// one sign throughout the interior and the search below provably cannot +/// find a witness pair, and what settles a root of even multiplicity that a +/// sign change cannot see either (`(x − 1)²` on `[0, 1]`). +/// 2. **A sign change**, described next. +/// +/// # Why the sign-change route is a proof /// /// The caller must already have obtained a successful [`bound_on_box`] over /// the same box. That is the continuity certificate: the branch-and-bound @@ -1228,9 +1320,15 @@ fn root_exists_witness( let mut w = SignWitnesses::default(); for point in seed_points(boxes0, prec) { - w.record(point_sign(expr, pool, &point, order, prec)); - if w.both() { - return true; + match point_outcome(expr, pool, &point, order, prec, true) { + PointOutcome::Zero => return true, + PointOutcome::Sign(s) => { + w.record(Some(s)); + if w.both() { + return true; + } + } + PointOutcome::Unknown => {} } } @@ -1255,6 +1353,12 @@ fn root_exists_witness( } continue; } + // An enclosure of the *range* over a whole sub-box that is the + // degenerate `[0, 0]` proves `expr` vanishes identically there, and + // the sub-box is non-empty — so every one of its points is a root. + if is_exact_zero(&r) { + return true; + } } let centre: Vec = b @@ -1264,9 +1368,15 @@ fn root_exists_witness( (*v, m.clone(), m) }) .collect(); - w.record(point_sign(expr, pool, ¢re, order, prec)); - if w.both() { - return true; + match point_outcome(expr, pool, ¢re, order, prec, false) { + PointOutcome::Zero => return true, + PointOutcome::Sign(s) => { + w.record(Some(s)); + if w.both() { + return true; + } + } + PointOutcome::Unknown => {} } if max_dim_width(&b, prec) <= floor { @@ -1286,20 +1396,29 @@ fn root_exists_witness( /// - [`Verdict::True`] when the rigorous range enclosure of `expr` over the /// whole box does not contain zero — `expr` is certified to have no root /// anywhere in the box. -/// - [`Verdict::False`] when the box is proven free of poles/branch cuts (the -/// full-box enclosure succeeded, so `expr` is continuous on the box) *and* -/// two points of the box are found at which `expr` is rigorously proven to -/// have opposite signs — a root is then certified to exist by the -/// intermediate value theorem along the segment joining them, which stays in -/// the box because a box is convex. The points are looked for by subdividing -/// the box, so an even number of roots no longer defeats the test: `x² − 2` -/// on `[-2, 2]` has two roots and two positive endpoints, and is settled at -/// the first bisection. -/// - [`Verdict::Undecided`] otherwise: the enclosure straddles zero and the -/// search found no pair of opposite-signed points within the budget. This is -/// the honest answer for a root that never produces a sign change at all — -/// a double root such as `(x − 1)²` — and it is never collapsed into either -/// of the other two verdicts. +/// - [`Verdict::False`] when a root is *certified to exist* in the box, by +/// either of two independent arguments (see the private `root_exists_witness`): +/// - **a point of the box proven to be a root** — its value pinned to zero by +/// a degenerate `[0, 0]` enclosure or by exact symbolic substitution. The +/// box is closed, so this covers a root sitting exactly on an endpoint +/// (`x` on `[0, 1]`, `x − 1` on `[0, 1]`) and a root of even multiplicity +/// that produces no sign change at all (`(x − 1)²` on `[0, 1]`); or +/// - **a sign change**: the box is proven free of poles/branch cuts (the +/// full-box enclosure succeeded, so `expr` is continuous on the box) *and* +/// two points of the box are found at which `expr` is rigorously proven to +/// have opposite signs — a root then exists by the intermediate value +/// theorem along the segment joining them, which stays in the box because +/// a box is convex. The points are looked for by subdividing the box, so +/// an even number of roots does not defeat the test: `x² − 2` on `[-2, 2]` +/// has two roots and two positive endpoints, and is settled at the first +/// bisection. +/// - [`Verdict::Undecided`] otherwise: the enclosure straddles zero, no point +/// of the box was *proven* to be a root, and the search found no pair of +/// opposite-signed points within the budget. An enclosure that merely +/// *contains* zero is never enough — a value known only to lie in `[-ε, ε]` +/// is not a proven root — so a function that grazes zero to within the +/// working precision without provably touching it stays `Undecided`, and it +/// is never collapsed into either of the other two verdicts. /// /// Propagates a [`ValidatedError`] (refuses) exactly when /// [`bound_on_box`] would: unsupported primitives, unbound symbols, or a @@ -2290,10 +2409,40 @@ mod tests { #[test] fn no_roots_undecided_for_a_multivariate_tangential_zero() { - // (x-1/2)² + (y-1/2)² is zero at exactly one point of [0,1]² and + // (x-1/3)² + (y-1/3)² is zero at exactly one point of [0,1]² and // positive everywhere else, so no sign-change witness can exist. The // enclosure straddles zero, so `True` is unavailable too: `Undecided` // is the only honest answer and must not collapse either way. + // + // The zero sits at a point no search this module performs can name — + // seed points and bisection midpoints are all dyadic — so the + // point-root proof cannot reach it either. That is deliberate: it is + // what keeps this test about the *absence* of a proof rather than + // about arithmetic luck. When the tangential zero does land on a point + // the search visits, the honest verdict changes; see + // `no_roots_false_for_a_tangential_zero_at_a_point_the_search_visits`. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let third = pool.rational(1_i32, 3_i32); + let dx = sub(&pool, x, third); + let dy = sub(&pool, y, third); + let e = pool.add(vec![pool.mul(vec![dx, dx]), pool.mul(vec![dy, dy])]); + let cheap = BoundOptions { + max_subdivisions: 64, + ..opts() + }; + let v = verified_no_roots(e, &pool, &[(x, 0.0, 1.0), (y, 0.0, 1.0)], &cheap).unwrap(); + assert_eq!(v, Verdict::Undecided); + } + + /// The same shape, with its single zero at the centre of the box — a point + /// the seed sweep evaluates. Substituting `x = y = 1/2` and simplifying + /// lands on the literal `0`, which *proves* the value there is zero, and + /// the centre is a point of the box. `False` is then a certificate, not the + /// guess the sign-change search would have had to make. + #[test] + fn no_roots_false_for_a_tangential_zero_at_a_point_the_search_visits() { let pool = ExprPool::new(); let x = pool.symbol("x", Domain::Real); let y = pool.symbol("y", Domain::Real); @@ -2306,7 +2455,7 @@ mod tests { ..opts() }; let v = verified_no_roots(e, &pool, &[(x, 0.0, 1.0), (y, 0.0, 1.0)], &cheap).unwrap(); - assert_eq!(v, Verdict::Undecided); + assert_eq!(v, Verdict::False); } /// `x² - 2` on every box from the issue-13 table, plus the product that @@ -2335,6 +2484,123 @@ mod tests { } } + /// A root sitting exactly *on* an endpoint of the box. + /// + /// Subdivision provably cannot settle these: `x` on `[0, 1]` is + /// non-negative throughout, so no negative witness exists anywhere in the + /// box and the IVT search must come back empty however long it runs. The + /// proof is not a search at all — the box is closed, `x = 0` is a point of + /// it, and `0` is exactly zero there. + #[test] + fn no_roots_false_for_a_root_on_a_box_endpoint() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + for (expr, lo, hi, label) in [ + (x, 0.0, 1.0, "x on [0,1] (root at the left endpoint)"), + ( + sub(&pool, x, pool.integer(1_i32)), + 0.0, + 1.0, + "x-1 on [0,1] (root at the right endpoint)", + ), + ( + pool.func("sin", vec![x]), + 0.0, + 1.0, + "sin on [0,1] (root at the left endpoint)", + ), + ( + pool.func("log", vec![x]), + 1.0, + 2.0, + "log on [1,2] (root at the left endpoint)", + ), + ( + sub(&pool, pool.func("exp", vec![x]), pool.integer(1_i32)), + 0.0, + 1.0, + "exp-1 on [0,1] (root at the left endpoint)", + ), + ] { + let v = verified_no_roots(expr, &pool, &[(x, lo, hi)], &opts()).unwrap(); + assert_eq!(v, Verdict::False, "{label}"); + } + } + + /// The other direction, and the one that keeps `False` a certificate: an + /// enclosure at an endpoint that merely *straddles* zero proves nothing. + /// + /// `exp(x) − 1 + 10⁻⁴⁰` is strictly positive on `[0, 1]` — it has no root + /// at all — but its value at `x = 0` is `10⁻⁴⁰`, far below the `2⁻¹²⁸`-scale + /// width of the enclosure there, so the enclosure contains zero. `True` is + /// therefore out of reach, and `False` must **not** be claimed: containing + /// zero is not being zero. `Undecided` is the only sound answer. + #[test] + fn no_roots_undecided_when_the_endpoint_enclosure_only_straddles_zero() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + // 10^40, written out so the test does not depend on an integer-power + // helper. + let ten_to_40 = + rug::Integer::from_str_radix("10000000000000000000000000000000000000000", 10).unwrap(); + let tiny = pool.rational(rug::Integer::from(1), ten_to_40); + let e = pool.add(vec![pool.func("exp", vec![x]), pool.integer(-1_i32), tiny]); + let v = verified_no_roots(e, &pool, &[(x, 0.0, 1.0)], &opts()).unwrap(); + assert_eq!(v, Verdict::Undecided); + } + + /// Randomised sweep over expressions whose roots are known exactly. + /// + /// Every box endpoint and every root is a multiple of 1/4, so roots land on + /// endpoints, on bisection midpoints and strictly between them — the three + /// cases the two proof routes divide up. The assertion is one-sided in each + /// direction and is the whole contract: `True` may only be returned when + /// there is genuinely no root in the closed box, `False` only when there + /// genuinely is one. `Undecided` is always permitted. + #[test] + fn no_roots_verdict_is_never_wrong_on_a_randomised_sweep() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let mut seed = 0x9E37_79B9_7F4A_7C15_u64; + let mut next = move |n: i64| -> i64 { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + (seed >> 33) as i64 % n + }; + let cheap = BoundOptions { + max_subdivisions: 64, + tol: 1e-6, + ..opts() + }; + + for _ in 0..120 { + // Roots and box endpoints on the same quarter-integer grid. + let (r1, r2) = (next(17) - 8, next(17) - 8); + let lo_q = next(17) - 8; + let hi_q = lo_q + 1 + next(8); + let (lo, hi) = (lo_q as f64 / 4.0, hi_q as f64 / 4.0); + + let d1 = sub(&pool, x, pool.rational(r1, 4_i64)); + let d2 = sub(&pool, x, pool.rational(r2, 4_i64)); + for (expr, roots) in [(d1, vec![r1]), (pool.mul(vec![d1, d2]), vec![r1, r2])] { + let has_root = roots.iter().any(|r| *r >= lo_q && *r <= hi_q); + let v = verified_no_roots(expr, &pool, &[(x, lo, hi)], &cheap).unwrap(); + match v { + Verdict::True => assert!( + !has_root, + "certified root-free but roots {roots:?}/4 are in [{lo}, {hi}]" + ), + Verdict::False => assert!( + has_root, + "certified a root but roots {roots:?}/4 are outside [{lo}, {hi}]" + ), + Verdict::Undecided => {} + } + } + } + } + #[test] fn no_roots_stays_true_where_it_was_true() { // The witness search must never be reached when the enclosure already @@ -2354,18 +2620,39 @@ mod tests { #[test] fn no_roots_undecided_for_a_double_root_that_cannot_be_witnessed() { - // (x-1)² has a genuine root at x = 1 inside [0,2], but it never + // (x-1/3)² has a genuine root at x = 1/3 inside [0,2], but it never // changes sign, so no IVT witness exists and none may be invented: - // turning this into `False` would be a lucky guess, not a proof. - // `True` is also unavailable (the enclosure contains zero). + // turning this into `False` on the strength of the search having got + // *close* would be a lucky guess, not a proof. `True` is also + // unavailable (the enclosure contains zero). The root is at a + // non-dyadic point, so no seed point or bisection midpoint ever lands + // on it and the point-root proof cannot fire either. let pool = ExprPool::new(); let x = pool.symbol("x", Domain::Real); - let d = sub(&pool, x, pool.integer(1_i32)); + let d = sub(&pool, x, pool.rational(1_i32, 3_i32)); let e = pool.mul(vec![d, d]); let v = verified_no_roots(e, &pool, &[(x, 0.0, 2.0)], &opts()).unwrap(); assert_eq!(v, Verdict::Undecided); } + /// The same double root, this time at a point the seed sweep visits. + /// A root of even multiplicity produces no sign change anywhere, so the + /// IVT search provably cannot settle it — but `(1-1)² = 0` is an exact + /// symbolic identity at a point of the box, which settles it outright. + #[test] + fn no_roots_false_for_a_double_root_proven_at_a_point() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let d = sub(&pool, x, pool.integer(1_i32)); + let e = pool.mul(vec![d, d]); + // x = 1 is the centre of [0, 2] and the right endpoint of [0, 1]: + // interior and boundary alike are points of the closed box. + for (lo, hi) in [(0.0, 2.0), (0.0, 1.0), (1.0, 3.0)] { + let v = verified_no_roots(e, &pool, &[(x, lo, hi)], &opts()).unwrap(); + assert_eq!(v, Verdict::False, "(x-1)^2 on [{lo}, {hi}]"); + } + } + #[test] fn no_roots_undecided_for_a_quartic_double_root() { // (x²-1)² — two double roots at ±1 in [-2,2], again with no sign diff --git a/alkahest-py/src/lib.rs b/alkahest-py/src/lib.rs index 699a9a41..69c4bab3 100644 --- a/alkahest-py/src/lib.rs +++ b/alkahest-py/src/lib.rs @@ -3247,8 +3247,38 @@ impl PyUniPoly { }) } - fn coefficients(&self) -> Vec { - self.inner.coefficients_i64() + /// Coefficients in ascending degree order, as exact Python `int`s. + /// + /// Lossless for coefficients of any size. This used to go through + /// `coefficients_i64()`, which does not merely saturate — it returned `0` + /// for anything past `i64`, so `2**100 * x**2 + 1` came back as + /// `[1, 0, 0]`: a quadratic reading as the constant `1`, with no exception + /// and no flag. That is exactly the silent-error class this project gates + /// against, and it is reachable from ordinary use, since `factor_z`, + /// resultants and pseudo-division all grow coefficients past 64 bits. + fn coefficients(&self, py: Python<'_>) -> PyResult> { + let int_cls = py.get_type_bound::(); + self.inner + .coefficients() + .into_iter() + .map(|c| Ok(int_cls.call1((c.to_string(),))?.into_py(py))) + .collect() + } + + /// Leading (highest-degree) coefficient, as an exact Python `int`. + /// + /// `0` for the zero polynomial, matching `degree == -1` there. Lossless + /// for coefficients of any size, as `coefficients()` also now is — so both + /// are safe on the output of `factor_z`, resultants and pseudo-division, + /// where the coefficients grow past 64 bits. + /// + /// A property, not a method: it is a single FLINT coefficient read. + #[getter] + fn leading_coeff(&self, py: Python<'_>) -> PyResult { + let int_cls = py.get_type_bound::(); + Ok(int_cls + .call1((self.inner.leading_coeff().to_string(),))? + .into_py(py)) } /// Degree of the polynomial (`-1` for the zero polynomial). diff --git a/docs/mdbook/src/representations.md b/docs/mdbook/src/representations.md index 101d9daa..14bb8e9a 100644 --- a/docs/mdbook/src/representations.md +++ b/docs/mdbook/src/representations.md @@ -44,20 +44,25 @@ p = UniPoly.from_symbolic(x**3 + pool.integer(-2) * x + pool.integer(1), x) print(p.degree) # 3 print(p.coefficients()) # [1, -2, 0, 1] (constant first) -print(p.leading_coeff()) # 1 +print(p.leading_coeff) # 1 (a property, and an exact Python int) # Arithmetic — all FLINT-backed, exact q = UniPoly.from_symbolic(x + pool.integer(-1), x) -print(p * q) # x^4 - x^3 - 2x^2 + 3x - 1 -print(p.gcd(q)) # x - 1 -print(p // q) # x^2 + x - 1 +print(p * q) # x^4-x^3-2*x^2+3*x-1 +print(p.gcd(q)) # x-1 +print(p // q) # x^2+x-1 print(p % q) # 0 # Powers r = UniPoly.from_symbolic(x + pool.integer(1), x) -print(r ** 3) # x^3 + 3x^2 + 3x + 1 +print(r ** 3) # x^3+3*x^2+3*x+1 ``` +`degree`, `is_zero` and `leading_coeff` are properties (zero-argument O(1) +accessors); `coefficients()` is a method because it allocates a list. Note that +`coefficients()` is `i64` and truncates coefficients that do not fit, whereas +`leading_coeff` is exact at any size. + `UniPoly` is the right choice when you are doing heavy univariate polynomial arithmetic (GCD chains, resultants, factorization) because FLINT applies highly optimized algorithms with exact arithmetic. ## MultiPoly @@ -76,10 +81,13 @@ print(mp.integer_content()) # 1 # Arithmetic mp2 = MultiPoly.from_symbolic(x * y, [x, y]) -print(mp + mp2) # x^2*y + x*y^2 + x*y - 1 -print(mp * mp2) # x^3*y^2 + x^2*y^3 - x*y +print(mp + mp2) # -1 + xy + xy^2 + x^2y +print(mp * mp2) # -xy + x^2y^3 + x^3y^2 ``` +Terms print in the polynomial's internal (ascending exponent-vector) order, not +in descending degree. + Variable order matters for the exponent-vector key. Pass variables in a consistent order when constructing `MultiPoly` objects that will be combined. ## MultiPolyFp @@ -91,25 +99,27 @@ from alkahest import sparse_interp_univariate, sparse_interp, gcd_sparse, MultiP p = 32749 # prime -# Recover a sparse univariate from 2T black-box evaluations (Ben-Or/Tiwari) -f = sparse_interp_univariate(lambda v: (v**5 + 3*v**3 + 7) % p, T=3, prime=p) -print(f) # x^5 + 3*x^3 + 7 (as MultiPolyFp) +# Recover a sparse univariate from 2·term_bound black-box evaluations +# (Ben-Or/Tiwari). Returns a list of (coefficient, exponent) pairs. +f = sparse_interp_univariate(lambda v: (v**5 + 3*v**3 + 7) % p, term_bound=3, prime=p) +print(f) # [(7, 0), (3, 3), (1, 5)] i.e. x^5 + 3*x^3 + 7 -# Recover a sparse multivariate via Zippel's algorithm +# Recover a sparse multivariate via Zippel's algorithm — this one is a MultiPolyFp, +# printed over positional variables x0, x1, ... in the order given by `vars` f2 = sparse_interp( lambda vals: (vals[0]**3 * vals[1]**2 + vals[0] * vals[1]**4) % p, - vars=[x, y], T=2, D=5, prime=p, + vars=[x, y], term_bound=2, degree_bound=5, prime=p, ) -print(f2) # x^3*y^2 + x*y^4 +print(f2) # 1*x0*x1^4 + 1*x0^3*x1^2 (mod 32749) # Sparse modular GCD over ℤ[x₁,...,xₙ] — substrate for exact GCD algorithms a = MultiPoly.from_symbolic((x + y) * (x - y), [x, y]) b = MultiPoly.from_symbolic((x + y) * (x + pool.integer(1)), [x, y]) h = gcd_sparse(a, b, term_bound=4, degree_bound=4) -print(h) # x + y +print(h) # y + x ``` -`sparse_interp_univariate` uses Berlekamp–Massey + BSGS discrete-log + Vandermonde solve and requires exactly `2T` oracle calls. `sparse_interp` uses Zippel's variable-by-variable algorithm with batched Vandermonde lifting. +`sparse_interp_univariate` uses Berlekamp–Massey + BSGS discrete-log + Vandermonde solve and requires exactly `2 * term_bound` oracle calls. `sparse_interp` uses Zippel's variable-by-variable algorithm with batched Vandermonde lifting. ## RationalFunction @@ -122,12 +132,12 @@ from alkahest import RationalFunction numer = x**2 + pool.integer(-1) denom = x + pool.integer(-1) rf = RationalFunction.from_symbolic(numer, denom, [x]) -print(rf) # x + 1 +print(rf) # 1 + x # Arithmetic preserves the rational form rf_x = RationalFunction.from_symbolic(x, pool.integer(1), [x]) rf_inv = RationalFunction.from_symbolic(pool.integer(1), x, [x]) -print(rf_x + rf_inv) # (x^2 + 1) / x +print(rf_x + rf_inv) # (1 + x^2) / (x) ``` GCD normalization runs at construction, so every `RationalFunction` is in lowest terms. @@ -137,14 +147,15 @@ GCD normalization runs at construction, so every `RationalFunction` is in lowest A real interval `[midpoint ± radius]` backed by FLINT's Arb library. Arithmetic on `ArbBall` values produces guaranteed enclosures of the true result. ```python -from alkahest import ArbBall, interval_eval, sin +from alkahest import ArbBall, ExprPool, interval_eval, sin # ArbBall(midpoint, radius, precision_bits=53) a = ArbBall(2.0, 0.5) # [1.5, 2.5] b = ArbBall(3.0, 0.0) # exactly 3 -print(a + b) # [4.5, 5.5] -print(a * b) # [4.5, 7.5] +# A ball prints as midpoint ± radius, not as an interval +print(a + b) # ArbBall(5.000000 ± 5.00e-1) i.e. [4.5, 5.5] +print(a * b) # ArbBall(6.000000 ± 1.50e0) i.e. [4.5, 7.5] # Evaluate a symbolic expression rigorously pool = ExprPool() @@ -162,10 +173,18 @@ See [Ball arithmetic](./ball-arithmetic.md) for more detail. ## Converting back to Expr -All specialized types can be converted back to a generic `Expr` for further symbolic manipulation: +Conversion back to `Expr` is per-type, not universal. `GbPoly` — the Gröbner +representation — round-trips: ```python -p = UniPoly.from_symbolic(x**2 + pool.integer(1), x) -expr_again = p.to_symbolic(pool) -dr = diff(expr_again, x) +from alkahest import expr_to_gbpoly, diff + +g = expr_to_gbpoly(x**2 + pool.integer(1), [x]) +expr_again = g.to_expr() # (x^2 + 1) +dr = diff(expr_again, x) # DerivedResult(value=(x * 2)) ``` + +`UniPoly`, `MultiPoly` and `RationalFunction` do **not** currently expose a +symbolic conversion. Read their coefficients (`UniPoly.coefficients()`, +`MultiPolyFp.terms`) or keep the original `Expr` alongside the polynomial — +`from_symbolic` does not consume it. diff --git a/examples/agent_workflow.py b/examples/agent_workflow.py index acfea79c..e22baab8 100644 --- a/examples/agent_workflow.py +++ b/examples/agent_workflow.py @@ -144,9 +144,12 @@ print("=" * 62) # Circle–line intersection: x² + y² = 1, y = x +# `numeric=True` evaluates the exact roots to floats; drop it to get the +# symbolic (radical) solutions back instead. solutions = solve( [x**2 + y**2 + pool.integer(-1), y + pool.integer(-1)*x], - [x, y] + [x, y], + numeric=True, ) print(f"Circle ∩ line (x²+y²=1, y=x):") for s in solutions: @@ -165,7 +168,8 @@ # Linear system linear_soln = solve( [x + y + pool.integer(-3), x + pool.integer(-1)*y + pool.integer(-1)], - [x, y] + [x, y], + numeric=True, ) print(f"\nx+y=3, x-y=1 → {[{str(k): round(v,4) for k,v in s.items()} for s in linear_soln]}") diff --git a/examples/lean_certificates.py b/examples/lean_certificates.py index 851b405d..dce79a86 100644 --- a/examples/lean_certificates.py +++ b/examples/lean_certificates.py @@ -16,7 +16,10 @@ def main(): x = pool.symbol("x") print("=== Lean export from expressions ===") - lean_expr = ak.to_lean(x**2 + pool.integer(1)) + # `to_lean(expr)` simplifies first and certifies the resulting derivation + # log; it returns "" when the simplifier had nothing to rewrite (an + # already-normal expression has no theorem to state), so give it one. + lean_expr = ak.to_lean(x**2 + pool.integer(1) + pool.integer(0)) print(lean_expr[:120] + ("…" if len(lean_expr) > 120 else "")) print("\n=== Lean certificate on diff (deriv goals) ===") diff --git a/examples/phase3_polynomials.md b/examples/phase3_polynomials.md index 5974ca8d..695f4168 100644 --- a/examples/phase3_polynomials.md +++ b/examples/phase3_polynomials.md @@ -41,10 +41,11 @@ x = pool.symbol("x") two = pool.integer(2) # x^2 + 2x + 1 -expr = pool.add([pool.pow(x, two), pool.mul([two, x]), pool.integer(1)]) +expr = pool.add([x ** two, pool.mul([two, x]), pool.integer(1)]) ``` -Available builders: `symbol`, `integer`, `add`, `mul`, `pow`. +Available builders: `symbol`, `integer`, `add`, `mul`; powers are written with +the `**` operator on an `Expr`. --- @@ -61,11 +62,10 @@ x = pool.symbol("x") # x^2 + 2x + 1 p = ak.UniPoly.from_symbolic( - pool.add([pool.pow(x, pool.integer(2)), + pool.add([x ** 2, pool.mul([pool.integer(2), x]), pool.integer(1)]), x, - pool, ) print(p) # x^2+2*x+1 @@ -78,20 +78,20 @@ print(p.coefficients()) # [1, 2, 1] All four operations are implemented via FLINT: ```python -xp1 = ak.UniPoly.from_symbolic(pool.add([x, pool.integer(1)]), x, pool) -xm1 = ak.UniPoly.from_symbolic(pool.add([x, pool.integer(-1)]), x, pool) +xp1 = ak.UniPoly.from_symbolic(pool.add([x, pool.integer(1)]), x) +xm1 = ak.UniPoly.from_symbolic(pool.add([x, pool.integer(-1)]), x) print(xp1 + xm1) # 2*x print(xp1 - xm1) # 2 print(xp1 * xm1) # x^2-1 -print(xp1.pow(3)) # x^3+3*x^2+3*x+1 +print(xp1 ** 3) # x^3+3*x^2+3*x+1 ``` ### GCD ```python x2m1 = ak.UniPoly.from_symbolic( - pool.add([pool.pow(x, pool.integer(2)), pool.integer(-1)]), x, pool + pool.add([x ** 2, pool.integer(-1)]), x ) g = x2m1.gcd(xm1) print(g) # x-1 @@ -109,8 +109,9 @@ leading coefficient. `rug::Integer` coefficient. The variable ordering is fixed at construction time by the `vars` list. -> **Display note:** variable names in the string representation use -> positional labels (`x0`, `x1`, …) matching the index in `vars`. +> **Display note:** the string representation uses the symbol names from the +> pool, and prints terms in the monomial order of the representation rather +> than the order they were written in. ### Construction @@ -121,13 +122,12 @@ y = pool.symbol("y") # x^2 + xy + y^2 a = ak.MultiPoly.from_symbolic( - pool.add([pool.pow(x, pool.integer(2)), + pool.add([x ** 2, pool.mul([x, y]), - pool.pow(y, pool.integer(2))]), - [x, y], # variable ordering: x=x0, y=x1 - pool, + y ** 2]), + [x, y], # variable ordering ) -print(a) # x1^2 + x0x1 + x0^2 +print(a) # y^2 + xy + x^2 print(a.total_degree) # 2 ``` @@ -140,7 +140,7 @@ the *primitive part*. # 6x + 4 → content = 2, primitive part = 3x + 2 b = ak.MultiPoly.from_symbolic( pool.add([pool.mul([pool.integer(6), x]), pool.integer(4)]), - [x, y], pool, + [x, y], ) print(b.integer_content()) # 2 ``` @@ -148,11 +148,11 @@ print(b.integer_content()) # 2 ### Arithmetic ```python -c = ak.MultiPoly.from_symbolic(pool.add([x, y]), [x, y], pool) -d = ak.MultiPoly.from_symbolic(pool.add([x, pool.integer(-1)]), [x, y], pool) +c = ak.MultiPoly.from_symbolic(pool.add([x, y]), [x, y]) +d = ak.MultiPoly.from_symbolic(pool.add([x, pool.integer(-1)]), [x, y]) -print(c + d) # 2x + y - 1 -print(c * d) # x^2 + xy - x - y +print(c + d) # -1 + y + 2x +print(c * d) # -y - x + xy + x^2 ``` --- @@ -178,19 +178,18 @@ rf = ak.RationalFunction.from_symbolic( pool.mul([pool.integer(6), x]), # numerator pool.integer(4), # denominator [x, y], - pool, ) -print(rf) # (3x0) / (2) -print(rf.numer()) # 3x0 +print(rf) # (3x) / (2) +print(rf.numer()) # 3x print(rf.denom()) # 2 # Denominator 1 is elided -rf2 = ak.RationalFunction.from_symbolic(x, pool.integer(1), [x, y], pool) -print(rf2) # x0 +rf2 = ak.RationalFunction.from_symbolic(x, pool.integer(1), [x, y]) +print(rf2) # x # Sign normalisation: denom leading coefficient is always positive -rf3 = ak.RationalFunction.from_symbolic(x, pool.integer(-2), [x, y], pool) -print(rf3) # (-x0) / (2) +rf3 = ak.RationalFunction.from_symbolic(x, pool.integer(-2), [x, y]) +print(rf3) # (-x) / (2) ``` --- @@ -213,7 +212,7 @@ x = pool.symbol("x") y = pool.symbol("y") try: - ak.UniPoly.from_symbolic(pool.add([x, y]), x, pool) + ak.UniPoly.from_symbolic(pool.add([x, y]), x) except ValueError as e: print(e) # unexpected free symbol 'y' ... ``` diff --git a/examples/phase3_polynomials.py b/examples/phase3_polynomials.py index 3eeee6c9..d759db53 100644 --- a/examples/phase3_polynomials.py +++ b/examples/phase3_polynomials.py @@ -27,23 +27,23 @@ x = pool.symbol("x") # Build x^2 + 2x + 1 -xsq = pool.pow(x, pool.integer(2)) +xsq = x ** 2 two_x = pool.mul([pool.integer(2), x]) one = pool.integer(1) expr = pool.add([xsq, two_x, one]) -p = ak.UniPoly.from_symbolic(expr, x, pool) +p = ak.UniPoly.from_symbolic(expr, x) print(f"p = {p}") # x^2+2*x+1 print(f" degree : {p.degree}") # 2 print(f" coeffs : {p.coefficients()}") # [1, 2, 1] (ascending degree) # q = x + 1 xp1 = ak.UniPoly.from_symbolic( - pool.add([x, pool.integer(1)]), x, pool + pool.add([x, pool.integer(1)]), x ) # r = x - 1 xm1 = ak.UniPoly.from_symbolic( - pool.add([x, pool.integer(-1)]), x, pool + pool.add([x, pool.integer(-1)]), x ) print(f"\nq = x+1 = {xp1}") @@ -54,12 +54,12 @@ print(f"q * r = {xp1 * xm1}") # x^2-1 (difference of squares) # Powers -print(f"\n(x+1)^3 = {xp1.pow(3)}") # x^3+3*x^2+3*x+1 -print(f"(x-1)^3 = {xm1.pow(3)}") # x^3-3*x^2+3*x-1 +print(f"\n(x+1)^3 = {xp1 ** 3}") # x^3+3*x^2+3*x+1 +print(f"(x-1)^3 = {xm1 ** 3}") # x^3-3*x^2+3*x-1 # GCD: gcd(x^2-1, x-1) → x-1 (up to leading coefficient sign) x2m1 = ak.UniPoly.from_symbolic( - pool.add([pool.pow(x, pool.integer(2)), pool.integer(-1)]), x, pool + pool.add([x ** 2, pool.integer(-1)]), x ) g = x2m1.gcd(xm1) print(f"\ngcd(x^2-1, x-1) = {g} (degree {g.degree})") @@ -69,8 +69,8 @@ # --------------------------------------------------------------------------- print("\n=== MultiPoly (multivariate polynomial over ℤ) ===\n") -# MultiPoly uses positional variable notation in display: the first element of -# the vars list is x0, the second is x1, etc. +# MultiPoly displays the symbol names from the pool, with terms in the +# monomial order of the representation (not necessarily as written below). pool2 = ak.ExprPool() x2 = pool2.symbol("x") @@ -78,25 +78,25 @@ # x^2 + x*y + y^2 expr_a = pool2.add([ - pool2.pow(x2, pool2.integer(2)), + x2 ** 2, pool2.mul([x2, y2]), - pool2.pow(y2, pool2.integer(2)), + y2 ** 2, ]) -a = ak.MultiPoly.from_symbolic(expr_a, [x2, y2], pool2) +a = ak.MultiPoly.from_symbolic(expr_a, [x2, y2]) print(f"a = x^2 + xy + y^2 = {a}") print(f" total_degree : {a.total_degree}") # 2 print(f" integer_content : {a.integer_content()}") # 1 # 6x + 4 → content = 2 expr_b = pool2.add([pool2.mul([pool2.integer(6), x2]), pool2.integer(4)]) -b = ak.MultiPoly.from_symbolic(expr_b, [x2, y2], pool2) +b = ak.MultiPoly.from_symbolic(expr_b, [x2, y2]) print(f"\nb = 6x+4 = {b}") print(f" integer_content : {b.integer_content()}") # 2 # Arithmetic -c = ak.MultiPoly.from_symbolic(pool2.add([x2, y2]), [x2, y2], pool2) +c = ak.MultiPoly.from_symbolic(pool2.add([x2, y2]), [x2, y2]) d = ak.MultiPoly.from_symbolic( - pool2.add([x2, pool2.integer(-1)]), [x2, y2], pool2 + pool2.add([x2, pool2.integer(-1)]), [x2, y2] ) print(f"\nc = x+y = {c}") @@ -118,21 +118,20 @@ pool3.mul([pool3.integer(6), x3]), pool3.integer(4), [x3, y3], - pool3, ) print(f"(6x)/(4) → {rf1}") -print(f" numer : {rf1.numer()}") # 3x0 +print(f" numer : {rf1.numer()}") # 3x print(f" denom : {rf1.denom()}") # 2 # x / 1 → denominator 1 is elided in display rf2 = ak.RationalFunction.from_symbolic( - x3, pool3.integer(1), [x3, y3], pool3 + x3, pool3.integer(1), [x3, y3] ) print(f"\nx/1 → {rf2}") # x / (-2) → sign normalised so denom has positive leading coefficient rf3 = ak.RationalFunction.from_symbolic( - x3, pool3.integer(-2), [x3, y3], pool3 + x3, pool3.integer(-2), [x3, y3] ) print(f"x/(-2) → {rf3}") @@ -147,14 +146,14 @@ # y is a free symbol when building a UniPoly in x try: - ak.UniPoly.from_symbolic(pool4.add([x4, y4]), x4, pool4) + ak.UniPoly.from_symbolic(pool4.add([x4, y4]), x4) except ValueError as e: print(f"Free symbol : {e}") # Zero denominator try: ak.RationalFunction.from_symbolic( - x4, pool4.integer(0), [x4], pool4 + x4, pool4.integer(0), [x4] ) except ValueError as e: print(f"Zero denom : {e}") @@ -162,7 +161,7 @@ # Negative exponent try: ak.UniPoly.from_symbolic( - pool4.pow(x4, pool4.integer(-1)), x4, pool4 + x4 ** pool4.integer(-1), x4 ) except ValueError as e: print(f"Negative exp : {e}") diff --git a/examples/risch_integration.py b/examples/risch_integration.py index 00809535..db5bbfe6 100644 --- a/examples/risch_integration.py +++ b/examples/risch_integration.py @@ -140,13 +140,17 @@ def main(): pool.integer(3) * x + pool.integer(2) * sqrt(p_q4), x) # ----------------------------------------------------------------------- - section("5. NonElementary guard — elliptic integrals") + section("5. NonElementary guard — hyperelliptic integrals") # ----------------------------------------------------------------------- + # Genus 1 (deg P = 3) is no longer a decline: it comes back in terms of + # the elliptic-integral primitives, e.g. + # ∫ sqrt(x³+1) dx → (2/5)·x·sqrt(x³+1) + (3/5)·3^(-1/4)·EllipticF(…) + # The guard now fires for genus ≥ 2, where no such reduction exists. pool = ExprPool(); x = pool.symbol("x") - p_ell = x ** 3 + pool.integer(1) + p_ell = x ** 5 + x + pool.integer(1) s_ell = sqrt(p_ell) - print(f"\n[∫ sqrt(x³+1) dx — should raise NonElementary]") + print(f"\n[∫ sqrt(x⁵+x+1) dx — should raise NonElementary]") display(pool, "f", s_ell) try: integrate(s_ell, x) diff --git a/tests/test_api.py b/tests/test_api.py index f2491688..f3356677 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -317,6 +317,39 @@ def test_is_zero(self): p = self._make([0]) assert p.is_zero + def test_leading_coeff(self): + # 3x^2 + 2x + 1 -> 3 + assert self._make([1, 2, 3]).leading_coeff == 3 + # The zero polynomial has no leading term; 0 pairs with degree == -1. + zero = self._make([0]) + assert zero.leading_coeff == 0 + assert zero.degree == -1 + + def test_coefficients_are_exact_beyond_i64(self): + """Neither accessor may truncate a coefficient past 64 bits. + + `coefficients()` used to go through `coefficients_i64()`, which does not + saturate — it returned `0`. So `2**100 * x**2 + 1` came back as + `[1, 0, 0]`: a quadratic reading as the constant `1`, with no exception + and no flag. Reachable from ordinary use, since `factor_z`, resultants + and pseudo-division all grow coefficients past 64 bits. + """ + pool, x = self.pool, self.x + big = 2**100 + 1 + p = UniPoly.from_symbolic(pool.integer(big) * x**2 + pool.integer(1), x) + + assert p.leading_coeff == big + assert isinstance(p.leading_coeff, int) + assert p.coefficients() == [1, 0, big] + assert all(isinstance(c, int) for c in p.coefficients()) + # The two accessors must agree — the bug was that they did not. + assert p.coefficients()[-1] == p.leading_coeff + + def test_coefficients_unchanged_for_small_values(self): + """The exact path must not perturb the ordinary case.""" + assert self._make([-5, 3, 1]).coefficients() == [-5, 3, 1] + assert self._make([0]).coefficients() in ([], [0]) + def test_add(self): p = self._make([1, 2]) q = self._make([3, 4]) diff --git a/tests/test_budget.py b/tests/test_budget.py index 574ca827..ebabbb02 100644 --- a/tests/test_budget.py +++ b/tests/test_budget.py @@ -18,6 +18,14 @@ #: call is a bug, not a slow machine. HEAVY_TIMEOUT = 120 +#: Ceiling on the **CPU** time a bounded cooperative callee may burn, in +#: milliseconds. See +#: ``test_run_with_wall_fallback_bounds_a_cooperative_callee`` for why the bound +#: is measured in CPU time rather than wall-clock time. Measured cost of that +#: case on a 12-core box: 2.7-5.3 s idle, 8.8 s under 2x oversubscription, so +#: this leaves ~7x headroom against a saturated machine. +COOPERATIVE_CALLEE_CPU_BOUND_MS = 60_000 + @pytest.fixture def pool() -> ak.ExprPool: @@ -551,15 +559,39 @@ def test_run_with_wall_fallback_bounds_a_cooperative_callee(pool, x): cooperative checkpoint stops on its own budget rather than only on the global cancel flag. Unbudgeted this integrand does not come back at all (see ``test_wall_budget_stops_a_hard_trig_integral``); the loose bound is - the property under test.""" + the property under test. + + The bound is on **CPU** time, not wall-clock time. It used to read + ``elapsed_ms < 20 * 300`` against a call that measurably costs 2.7-5.3 s + when it works -- a 13% margin -- so the test went red whenever the machine + was busy, which is a fact about the box and not about the property. Wall + time here is ``wall_ms`` (a real-time timer, so it does not stretch) plus + the join of a worker that is still running, and contention inflates that + second term without the callee doing any more work. Measured on one box, + idle vs. 24 spinners on 12 cores: wall 5.3 s -> 24.9 s (4x over the old + bound, a guaranteed failure), CPU 5.3 s -> 8.8 s. CPU time is not perfectly + load-free -- contention costs some real cycles -- but it tracks the work + done rather than the waiting, which is the thing under test, and it leaves + real headroom under ``COOPERATIVE_CALLEE_CPU_BOUND_MS``. + + The property is still enforced from both ends: a callee that stops seeing + the budget burns CPU without limit and trips the assertion, and one that + stops coming back at all trips the ``timeout`` marker above. + """ s = ak.sin(x) hard = ak.cos(x) * s**60 / (s**31 + s + 1) - started = time.perf_counter() + cpu_started = time.process_time() with pytest.raises(ak.BudgetExceededError) as excinfo: ak.run_with_wall_fallback(ak.integrate, hard, x, budget=ak.Budget(wall_ms=300)) - elapsed_ms = (time.perf_counter() - started) * 1000.0 + cpu_ms = (time.process_time() - cpu_started) * 1000.0 assert excinfo.value.code == "E-BUDGET-001" - assert elapsed_ms < 20 * 300 + # The call ended on the wall-clock fallback's own join, not on some other + # budget check that happens to raise the same code. + assert "returned control after" in str(excinfo.value) + assert cpu_ms < COOPERATIVE_CALLEE_CPU_BOUND_MS, ( + f"the callee burned {cpu_ms:.0f} ms of CPU against a 300 ms budget: the " + "wall-clock fallback is no longer bounding a cooperative callee" + ) @pytest.mark.timeout(HEAVY_TIMEOUT) diff --git a/tests/test_evaluate.py b/tests/test_evaluate.py index 93dfd4b4..57e105ce 100644 --- a/tests/test_evaluate.py +++ b/tests/test_evaluate.py @@ -48,3 +48,63 @@ def test_unsupported_evaluation_returns_stable_status(): assert result.status == "unsupported" assert result.value is None assert result.reason == "E-EVAL-001" + + +def test_interval_evaluation_covers_every_primitive_that_advertises_numeric_ball(): + """The interval evaluator dispatches through the primitive registry, so the + set of functions it accepts *is* the set the registry reports a + ``numeric_ball`` kernel for. It used to carry its own hand-written list, + which had silently dropped ``bessel_j0`` / ``bessel_j1``: both had rigorous + ball kernels and both were refused with ``E-EVAL-010``. + + No single probe point is inside every kernel's domain (``acosh`` needs + x >= 1, ``atanh`` needs |x| < 1), so a primitive counts as covered if some + probe evaluates — a domain refusal is a different answer from "no rule for + this name", and only the second is a coverage gap. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + registry = ak.PrimitiveRegistry.default_registry() + + refused = [] + for row in registry.coverage_report(): + if not row["numeric_ball"]: + continue + call = pool.func(row["name"], [x]) + if all( + evaluate(call, {x: ak.ArbBall(probe, 0.0)}, mode="interval").status != "ok" + for probe in (1.5, 0.5) + ): + refused.append(row["name"]) + + assert not refused, f"advertised numeric_ball but interval evaluation refuses: {refused}" + + +def _at(build, xi): + """`build(y)` at the exact point `xi`, as a tight point enclosure.""" + pool = ak.ExprPool() + y = pool.symbol("y") + result = evaluate(build(y), {y: ak.ArbBall(xi, 0.0)}, mode="interval") + assert result.status == "ok" + return result.value.mid + + +def test_interval_evaluation_of_bessel(): + """The two the coverage audit was opened for.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + for build, expected in ((ak.bessel_j0, 0.7651976865579665), (ak.bessel_j1, 0.4400505857449335)): + result = evaluate(build(x), {x: ak.ArbBall(1.0, 0.0)}, mode="interval") + + assert result.status == "ok" + assert abs(result.value.mid - expected) < 1e-15 + + # J0 and J1 oscillate, so the enclosure has to come from a Lipschitz + # bound around the midpoint, not from hulling the two endpoints; on a + # wide interval an endpoint hull misses the function's own extrema. + wide = evaluate(build(x), {x: ak.ArbBall(0.0, 1.0)}, mode="interval") + assert wide.status == "ok" + for i in range(21): + xi = -1.0 + 2.0 * i / 20 + assert wide.value.contains(_at(build, xi)), f"J({xi}) escaped {wide.value}" diff --git a/tests/test_taylor_model_coverage.py b/tests/test_taylor_model_coverage.py index 2c78ca87..f76e4bba 100644 --- a/tests/test_taylor_model_coverage.py +++ b/tests/test_taylor_model_coverage.py @@ -76,10 +76,15 @@ def test_taylor_model_flag_agrees_with_bound_on_box(row): 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 + `numeric_ball` is *not* wrong for these eleven — 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. + + `atanh` joined the set when the capability probe stopped testing ball + kernels at `1.0` only: its domain is the open interval `(-1, 1)`, so it + declined the sole probe point and lost a bit it had earned, while keeping + `numeric_f64` because that probe already tried `0.5`. """ rows = {row["name"]: row for row in _primitive_rows()} ball_only = { @@ -88,6 +93,7 @@ def test_the_flag_is_not_a_restatement_of_numeric_ball(): assert ball_only == { "acosh", "asinh", + "atanh", "bessel_j0", "bessel_j1", "ceil", diff --git a/tests/test_validated_bounds.py b/tests/test_validated_bounds.py index dc69ceec..bacd32dd 100644 --- a/tests/test_validated_bounds.py +++ b/tests/test_validated_bounds.py @@ -309,10 +309,10 @@ def test_no_roots_true_cases_stay_true(build, lo, hi): ("build", "lo", "hi", "why"), [ ( - lambda pool, x: (x - pool.integer(1)) * (x - pool.integer(1)), + lambda pool, x: (x - pool.rational(1, 3)) * (x - pool.rational(1, 3)), 0.0, 2.0, - "double root at x=1: no sign change anywhere", + "double root at x=1/3: no sign change anywhere", ), ( lambda pool, x: (x * x - pool.integer(1)) * (x * x - pool.integer(1)), @@ -324,14 +324,67 @@ def test_no_roots_true_cases_stay_true(build, lo, hi): ) def test_a_root_that_cannot_be_witnessed_stays_undecided(build, lo, hi, why): """These expressions *do* have roots in the box, but they never change - sign, so no intermediate-value witness exists. `"undecided"` is the honest - answer; reporting `"false"` here would be a guess dressed as a proof.""" + sign, so no intermediate-value witness exists, and no point the search + evaluates is *proven* to be a root either — the search only ever visits + endpoints, corners and repeated midpoints, all of which are dyadic. + `"undecided"` is the honest answer; reporting `"false"` here would be a + guess dressed as a proof.""" pool = ak.ExprPool() x = pool.symbol("x") assert ak.verified_no_roots(build(pool, x), [(x, lo, hi)]) == "undecided", why +@pytest.mark.parametrize( + ("build", "lo", "hi", "why"), + [ + (lambda pool, x: x, 0.0, 1.0, "root at the left endpoint"), + (lambda pool, x: x - pool.integer(1), 0.0, 1.0, "root at the right endpoint"), + (lambda pool, x: ak.sin(x), 0.0, 1.0, "sin vanishes at the left endpoint"), + (lambda pool, x: ak.log(x), 1.0, 2.0, "log vanishes at the left endpoint"), + (lambda pool, x: ak.exp(x) - pool.integer(1), 0.0, 1.0, "exp-1 at the left endpoint"), + ( + lambda pool, x: (x - pool.integer(1)) * (x - pool.integer(1)), + 0.0, + 1.0, + "double root at the right endpoint", + ), + ( + lambda pool, x: (x - pool.integer(1)) * (x - pool.integer(1)), + 0.0, + 2.0, + "double root at the centre, which is a point of the box too", + ), + ], +) +def test_a_root_on_the_box_boundary_is_proven(build, lo, hi, why): + """The box is closed, so a point of it at which the function is *proven* + zero is a root in the box — no continuity argument and no sign change + required. Subdivision alone can never settle these: `x` on [0,1] is + non-negative throughout, so a negative witness does not exist to be found. + """ + pool = ak.ExprPool() + x = pool.symbol("x") + + assert ak.verified_no_roots(build(pool, x), [(x, lo, hi)]) == "false", why + + +def test_an_endpoint_enclosure_that_only_straddles_zero_stays_undecided(): + """The soundness side of the same coin. `exp(x) - 1 + 10**-40` is strictly + positive on [0,1] — it has no root at all — but its value at x=0 is far + below the width of any enclosure computable there, so the enclosure + contains zero. Containing zero is not being zero: `"false"` would be a + false certificate, and `"true"` is out of reach, so `"undecided"` is the + only sound answer.""" + pool = ak.ExprPool() + x = pool.symbol("x") + # 10**-40 as an exact power; `pool.rational` takes C longs, and 10**40 + # does not fit in one. + f = ak.exp(x) - pool.integer(1) + pool.integer(10) ** pool.integer(-40) + + assert ak.verified_no_roots(f, _box(x, 0.0, 1.0)) == "undecided" + + def test_sign_positive_verified_true(): pool = ak.ExprPool() x = pool.symbol("x")