diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c103969..a107f8c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,137 +1,20 @@ # Changelog -## 3.8.0 — 2026-08-12 +## 3.9.0 — 2026-08-14 -### Silent errors fixed — do results you already computed need rechecking? +Everything in this section landed **after `v3.8.0` was tagged and published**, +so none of it is in the released 3.8.0 wheel. It is the result of acting on a +trial autoresearch run: fifteen logged issues, of which this release fixes +fourteen (issue #6 was expectation-setting, not a defect, and issue #10 is +reduced but not closed — see below). -A *silent error* is a confident, plausible, mathematically wrong answer with no -exception, no `NaN` and no verification flag. Eleven were found and fixed this -release. **Eight of them shipped in 3.7 or earlier**, so if you have results -from an affected call, re-run them. The other three were in code added during -this release cycle and never reached a published wheel. - -| Affected call | Wrong answer it gave | First shipped in | Recheck? | -|---|---|---|---| -| `decide(Forall(x, φ))` where the counterexample is a rational root whose denominator is **not a power of two** | `True` for a **false** universal theorem, e.g. `∀x. (3x+2)² > 0` (false at `x = −2/3`) | ≤ 3.7 | **Yes** — any `decide` verdict | -| `decide(Exists(x, φ))` with an `=` atom | `(True, witness)` where the witness does **not** satisfy the sentence, e.g. `∃x. 3x−2 = 0 → x = 1/2` | ≤ 3.7 | **Yes** — any cited witness | -| `Matrix.nullspace()` on a 2×2 with a symbolic determinant | A confident wrong kernel basis; `[[x,0],[0,1]]` returned `(0, x)`, for which `M·v ≠ 0` | 3.7 | **Yes** — verify `M·v = 0` numerically | -| `simplify` / `simplify_egraph` on a product containing `0⁻¹` | `1`, or `0`, depending on the engine — for an expression with no value at all. Reachable from `diff(2/(x − x), x)` | ≤ 3.7 | Yes, if any input could reduce to `0⁻¹` | -| `decide` on a two-variable sentence true only at an irrational point | `False` for a satisfiable `∃x∃y`, and `True` for its false `∀x∀y` dual | this cycle (2-var `decide` is new) | No published release affected | -| `batch_map(..., parallel=True)` under `context(budget=…)` | Ran **unbudgeted**, so candidates a sequential sweep reported as `E-BUDGET-001` came back as `E-INT-001` — a *mathematical* verdict | this cycle (batch APIs are new) | No published release affected | -| `product_definite` on a term with any non-integer coefficient | Off by `c^(hi−lo+1)`: `Π_{k=1}^{5} ½` returned `1` instead of `1/32`, `Π (2k−1)/(2k)` at `n = 6` returned `14.4375` instead of `0.2255859375` | ≤ 3.7 | **Yes** — any `product_definite` / `product_indefinite` result | -| `sum_definite` where the summand has a pole strictly *between* the bounds | A clean finite number for a sum with an undefined term: `Σ_{k=1}^{10} 1/((k−3)(k−2))` returned `−5/8` | ≤ 3.7 | **Yes** — any `sum_definite` over a range containing a denominator root | -| `euler_maclaurin` when `corrections` is too small for the summand | A fabricated additive constant — the missing term frozen at the fitting point. `Σ k⁹` at the default `corrections = 2` acquired `34359738368 = 512⁴/2` in a Faulhaber polynomial whose constant term is `0` | this cycle (Euler–Maclaurin is new) | No published release affected | -| `rsolve` on a **forward-shift** spelling with a non-zero right-hand side | The solution of a *different* equation: `f(n+1) − f(n) = n²` with `f(0) = 0` returned `Σ_{j=1}^{n} j²` instead of `Σ_{j=0}^{n−1} j²` | ≤ 3.7 | **Yes** — any `rsolve` written with `f(n+i)`, `i > 0`, and an inhomogeneous term | -| `rsolve` / `solve_linear_recurrence_homogeneous` on an order-2 recurrence with a **repeated** characteristic root | `C₀·rⁿ + C₁·rⁿ` — a one-parameter family presented as the general solution of a second-order equation, losing the `n·rⁿ` branch | ≤ 3.7 | **Yes** — check the discriminant of `r² + b r + c` | - -Also fixed, and not a silent error but worse for an unattended loop: a Rust -panic escaped `interval_eval` as `pyo3_runtime.PanicException`, which inherits -from `BaseException` and therefore slips past `except Exception`. Shipped in -3.7 — a loop that survived everything else died on `x^(3/2)` over a negative -ball. - -The deterministic silent-error gate (`tests/silent_errors/`, Tier-1 CI) now -scores **0 silent errors out of 241 scored cases** across evaluation, -integration, limits, linear algebra, number theory, real QE, series, -simplification, solving, and sums/products. Every trap added this cycle was -re-run against a build with the fix reverted and confirmed to score -`silent_error` there, and every trap is paired with a **control** — its nearest -convergent neighbour — so a subsystem cannot pass the gate by refusing -everything. That is a statement about the corpus, not a guarantee about the -library. +**Upgrading from 3.8.0 is not transparent.** Two changes need a read before you +upgrade: `relation_confidence` now returns a *tri-state* verdict, and seventeen +zero-argument accessors became properties, which is a hard break with no alias. +Both are detailed under *Behaviour changes to plan for*. ### Behaviour changes to plan for -Fixing a silent error means some calls that used to return now refuse. Every one -of these is a call whose previous answer was not justified: - -- **`decide` raises `CadError` (`E-CAD-001`) where it used to answer**, whenever - the formula has a non-strict atom (`=`, `≠`, `≤`, `≥`) and a boundary root has - not been shown rational. This includes mixed-alternation sentences that route - through De Morgan — `∀x∃y. p > 0` becomes `¬∃x∀y. p ≤ 0`, and the negation - makes a strict body non-strict. `decide` is **not** a complete decision - procedure in this implementation; treat `E-CAD-001` as *undecided*, never as - *false*. -- **`rank`, `rref`, `nullspace`, `eigenvects`, `jordan_form` raise - `E-LINALG-010`, and `inverse` raises the new `E-MAT-004`**, when an entry's or - the determinant's vanishing can be decided neither way. Previously "could not - prove non-zero" was silently read as "zero". -- **`simplify` leaves `0 · 0⁻¹` unevaluated** instead of returning `1` (or `0`). - A result containing `(0 * 0^-1)` is Alkahest declining to give an - indeterminate form a value, not a simplifier failure. -- **`sum_definite` raises `SumError` (`E-SUM-003`) when the summand is undefined - at an integer inside `[lo, hi]`**, not only when the pole lands on `lo` or - `hi+1`. The refusal names the offending index. Sums whose poles lie outside - the range are unaffected: `Σ_{k=4}^{10} 1/((k−3)(k−2))` still returns `7/8`. -- **`euler_maclaurin` may return a shorter expansion, with no additive - constant.** The constant is now fitted at a point outside the gate's check - points and re-fitted at a second one; if the two disagree it is not a constant - and none is claimed. The report says which way that went in `derivation`, and - the `"fitted numerically"` hypothesis is only listed when a fitted constant is - actually part of the answer. Genuine constants (`γ`, `ζ(2)`, `½log 2π`, …) are - unaffected — they agree across fitting points to 13+ digits. -- **`product_definite(term, k, lo, hi)` with `lo > hi` returns `1` even for a - zero term.** The empty product takes no factors; it previously returned `0` - for `Π_{k=1}^{0} 0` while returning `1` for `Π_{k=1}^{0} k`. -- **`capabilities()["contract_version"]` is `3`, and `features` lost two keys: - `groebner_cuda` and `numpy`.** Indexing either now raises `KeyError`; use - `features.get(name, False)` if you need to span versions. Both were removed - rather than wired up because neither was *falsifiable* — no observation a - Python caller could make distinguished `True` from `False`: - - `groebner_cuda` reported that the CUDA Macaulay-matrix kernel had been - compiled in. The string `groebner_cuda` occurred exactly once anywhere in - `alkahest-py` — the capability line itself. There was no binding, no - `*gpu*` name in the public or the private module, and `GroebnerBasis` - exposes only CPU methods. The kernel is unchanged and still reachable from - Rust as `alkahest_cas::poly::groebner::compute_groebner_basis_gpu`; if - dispatch ever prefers it, the binding lands first and a bit follows it. - - `numpy` mapped to a Cargo feature gating the `numpy` crate, which - `alkahest-py` never used an item from. The feature and the dependency are - both gone. `ak.numpy_eval` and `ak.numpy_eval_par` are unaffected — they go - through the buffer protocol and always worked with the bit `False`, which - is its value on every wheel ever published. - - An unreachable `True` makes a caller trust something it should not, which is - the same class of defect as a silent wrong answer; a bit that correlates with - nothing is better removed than left to be misread. - `tests/test_agent_contract.py::test_every_advertised_feature_has_an_entry_point` - now walks `features` and fails on any key without a named, reachable entry - point, so the next one cannot ship. -- **Rust, `--features groebner-cuda`: `compute_groebner_basis_gpu` and - `reduce_batch` return `(polys, GpuBackendReport)` instead of `polys`.** Both - fall back to CPU row reduction — when `device_id` is `None`, and when the - driver fails — and the basis is identical either way, so a caller previously - had no way to tell a GPU run from a CPU one. `GpuBackendReport::ran_on_gpu()` - is true only when at least one mod-p reduction ran on a device and none fell - back; `reductions_on_gpu`, `reductions_on_cpu` and `first_gpu_error` carry - the detail. A compile error on upgrade is the intended failure mode for code - that was recording these results as GPU results. Nothing at the Python - surface changes: the feature has no binding. -- **`residue(f, z, point)` refuses a non-constant `point` with - `AlkahestError` / `E-RESIDUE-005`** instead of leaking - `AttributeError: 'Expr' object has no attribute 'numerator'` from the - argument parser. `AttributeError` is not an `AlkahestError`, so - `except ak.AlkahestError` missed it entirely. The existing `E-RESIDUE-001..4` - refusals are now `AlkahestError`s carrying `.code` and `.remediation` too, - rather than bare `ValueError`s with the code glued into the message; - `AlkahestError` subclasses `ValueError`, so `except ValueError` still works. -- **`series` refuses instead of running forever, with `SeriesError` / - `E-SERIES-003`.** `series(sqrt(t**-2 + t**-1), t, 0, 32)` never returned: - coefficients are formed by repeated differentiation without re-simplifying, so - a nested radical's derivatives grow by a constant factor per coefficient and - the cost doubles per order. It now honours an active `Budget` (raising - `BudgetExceededError`) and, with none, an internal work ceiling. It never - returns a *shorter* series: `O(h^order)` on fewer coefficients than were asked - for is a false statement about the remainder, which is worse than the refusal. - Ordinary expansions are unaffected — the heaviest in the suites intern a few - thousand nodes against a ceiling of 50 000. -- **`simplify_expanded` records a derivation step when its expansion bound stops - it** (`expand_pow_limit_reached`, a no-op step naming the power it declined), - and the bound itself is now a budget on the number of distributed products - rather than a flat exponent cap. `(x+y)**6` and `(x+y+1)**7` now expand where - the exponent-only cap refused them while permitting a twenty-term sum to the - fourth power; anything above the budget comes back unexpanded *and says so* - instead of looking like an expression that was already expanded. - **`relation_confidence` returns `credible: None` — *unknown* — for inputs whose precision it cannot establish, where it used to return `True`.** It judged only `float` inputs, on the premise that "decimal strings and ints are @@ -175,74 +58,124 @@ of these is a call whose previous answer was not justified: emitted. Nothing was unsound; the question being answered had changed. Code that relied on pool symbols being `Real` inside an integer context should pass `domain="real"` explicitly. - -### Known limits — documented, not fixed - -These are properties of the design as it stands. They are called out here -because 3.8 is aimed at long unattended loops, and each of them is a way such a -loop fails. - -- **`ExprPool` never reclaims.** The arena is append-only: no `clear`, no - refcount, no GC, and the storage cannot shrink. The only way to free interned - nodes is to **drop the whole pool** — and every `Expr`, `Matrix`, `Series` and - `DerivedResult` holds a *strong* reference to its pool, so retaining one - interesting result retains everything. Growth on a shared pool is linear and - unbounded (~200 bytes/node; measured ~2 KB per `integrate` call over 20 000 - calls, 0 B/call with a fresh pool per iteration) while per-call latency stays - **flat**, so the failure mode is a clean OOM with no slowdown to warn you - first. `ExprPool` also exposes no `__len__` or `stats()`, so the growth is not - observable from Python. The supported pattern is **one pool per problem**, - documented in [`budgets.md`](docs/mdbook/src/budgets.md#exprpool-never-reclaims). -- **`run_with_wall_fallback` does not bound wall time for an uncooperative - callee.** It joins its worker before the exception propagates, so it returns - when the callee returns: `run_with_wall_fallback(time.sleep, 3.0, - budget=Budget(wall_ms=50))` raises `E-BUDGET-001` after 3000 ms, and the - message reports the real elapsed time so this shows up in a log rather than - being inferred later. Python cannot kill a thread, and abandoning one would - leak a live thread that still allocates into the pool and can only be stopped - through the process-wide cancel flag. Only an **OS-level bound** (subprocess, - process watchdog) is a hard deadline. -- **`wall_ms` granularity is one primitive operation, and FLINT calls cannot be - interrupted.** After the checkpoint work above the overshoot is a small - additive term (1.0–1.2×), but past a certain degree a single operation is a - FLINT factorisation or resultant — one foreign-function call, ~2 s on a - degree-62 integrand, which no cooperative mechanism can stop part-way. -- **`Matrix.eigenvals()` grows the pool on identical input** (~1.9 KB/call, - measured over 20 000 calls on the same 2×2 integer matrix): it interns a fresh - `__eigen_lambda_N` gensym per call. Every other Python-facing entry point - measured is flat on repeated input. Cache eigenvalue results. -- **`Matrix.eigenvals()` can emit casus-irreducibilis cube roots** — correct - under Alkahest's real cube-root convention (and honestly refused by - `eval_expr` with `E-EVAL-009`, with `interval_eval` returning an unbounded - ball) but evaluated on the **principal** branch by SymPy, NumPy and most other - tools, which return a confident number that is not an eigenvalue. 14 of 720 - random integer matrices produced one. An honest refusal here becomes somebody - else's silent error the moment the expression crosses the boundary, so - evaluate inside Alkahest before exporting, or export a verified numeric - enclosure instead. See [`interop.md`](docs/mdbook/src/interop.md). -- **The LLVM JIT leaks an LLVM `Context` per compile** (`Box::leak`, on the - error paths as well as the success path). Feature-gated behind `jit`, so - default PyPI wheels (Cranelift) are unaffected; do not compile in a loop under - a `+jit` / `+full` wheel. -- **No sanitizer covers any Python-facing path.** The PR-gating ASan job runs - with `detect_leaks=0`, the nightly LSan shard cannot reach a `cdylib` with no - `#[test]` functions, and `pytest` is never run under a sanitizer. The - behavioural substitute is the fresh-pool sweep described in - [`TESTING.md`](TESTING.md#3-memory-safety--sanitizers). -- **There is no `cuda_device_count()`.** `CudaCompiledFn.call_batch_on(ordinal, - …)` selects a device, but the valid range can only be discovered by trying an - ordinal and catching `CudaError` (`E-CUDA-003`); the loop that does it is in - [`gpu.md`](docs/mdbook/src/gpu.md#discovering-the-valid-device-ordinals). Not - added yet on purpose: `cuda` implies LLVM 15 with NVPTX, so such a binding - cannot be compiled on an ordinary dev box, no CI job builds the Python - extension with either CUDA feature, and exercising it needs a device — it - would ship with no verification of any kind, which is the provenance of the - capability overclaims fixed above. It belongs in the same change as the - missing `maturin develop --features cuda` + `pytest tests/test_cuda.py` - nightly step. +- **Seventeen zero-argument scalar accessors became properties — drop the + `()`.** There was no rule a caller could predict: `Enclosure.width` and + `RegularChain.n_vars` were properties while `DAE.n_equations()` and + `MultiPoly.total_degree()` were methods, so the same shape of question was + asked two different ways depending on the class. The convention now is: **a + zero-argument, O(1), non-allocating accessor returning a scalar or a flag is a + property; anything that returns a collection, allocates, or does real work is + a method.** No compatibility alias is provided, deliberately — an accessor + that answers to both forms leaves the ambiguity in place. Migration is + mechanical: + + | Before | After | + |---|---| + | `UniPoly.degree()` | `UniPoly.degree` | + | `UniPoly.is_zero()` | `UniPoly.is_zero` | + | `MultiPoly.is_zero()` | `MultiPoly.is_zero` | + | `MultiPoly.total_degree()` | `MultiPoly.total_degree` | + | `MultiPolyFp.is_zero()` | `MultiPolyFp.is_zero` | + | `MultiPolyFp.total_degree()` | `MultiPolyFp.total_degree` | + | `RationalFunction.is_zero()` | `RationalFunction.is_zero` | + | `GbPoly.is_zero()` | `GbPoly.is_zero` | + | `GbPoly.n_vars()` | `GbPoly.n_vars` | + | `ODE.order()` | `ODE.order` | + | `DAE.n_equations()` | `DAE.n_equations` | + | `DAE.n_variables()` | `DAE.n_variables` | + | `HybridODE.n_events()` | `HybridODE.n_events` | + | `Component.n_equations()` | `Component.n_equations` | + | `Component.n_ports()` | `Component.n_ports` | + | `OdeTrajectory.t_final()` | `OdeTrajectory.t_final` | + | `ArbBall.is_exact()` | `ArbBall.is_exact` | + + **Calling the old form now raises `TypeError: 'int' object is not callable` + (or `'bool'`, `'float'`), which is loud. The reverse mistake is not.** + Writing `if dae.n_equations:` against a *pre*-3.8.0 build silently reads a + bound method, which is always truthy, and `f"{dae.n_equations}"` formats as + ``; so grep for these names rather than waiting for a + traceback. Accessors that were *already* properties (`Enclosure.lower`, + `.upper`, `.width`, `.subdivisions`, `Matrix.rows`, `.cols`, + `RegularChain.n_vars`, `RosenfeldGroebnerResult.consistent`, `.truncated`, + `ArbBall.mid`, `.rad`, `.lo`, `.hi`, …) are unchanged; they were already + correct under the rule, as were the collection-returning methods that sit + beside them (`RegularChain.polys()`, `RosenfeldGroebnerResult.final_basis()`). + Three zero-argument scalar calls stay methods because they do real work rather + than read a field: `Matrix.rank()` (Gaussian elimination), `ODE.is_autonomous()` + (walks every RHS expression) and `PositivityCertificate.verify()` (re-runs the + exact SOS identity check). `tests/test_accessor_convention.py` pins the + converted set and scans `alkahest-py/src/lib.rs` for new offenders, so the + inconsistency cannot creep back. ### Fixed +- **`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 + obvious they were: `x²−2` on `[-2, 0]` was `"false"`, but the same expression + on `[-2, 2]` — which contains *both* roots — was `"undecided"`, as were + `[-10, 10]` and `(x²−2)(x²+1)` on `[-2, 2]`. The machinery to settle it was + already running; one bisection of `[-2, 2]` makes each half answer + immediately. The `"false"` direction now searches for its two witness points + by subdividing the box, and the intermediate-value argument is stated over the + box rather than over an interval: a box is convex, so two points at which the + expression is *rigorously proven* to have opposite signs certify a root on the + segment joining them, which stays inside the box. That also lifts the test to + 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. +- **`verified_integral` refused removable singularities.** Taylor-model + quadrature raised `E-VALIDATED-003` on any sub-interval where the reciprocal's + enclosure contained zero, which put `∫₀¹ ln(1+x)/x dx = π²/12` out of reach + even though nothing about that integral is singular — only the expression as + written is, and the integrand extends continuously to 1 at `x = 0`. An + integrand that splits as `N(x)/D(x)` with `N(p) = D(p) = 0` at a point `p` of + the offending sub-interval is now enclosed through **Cauchy's mean value + theorem** instead: `N(x)/D(x) = N′(ξ)/D′(ξ)` for some `ξ` in the sub-interval, + so the piece is bounded by an enclosure of `N′/D′`, which is regular. The + value returned is the integral of the continuous extension. + `∫₀¹ ln(1+x)/x`, `∫_{-1}^{1} sin(x)/x`, `∫₀¹ (eˣ−1)/x` and + `∫₀¹ (1−cos x)/x` now come back as enclosures that bracket `π²/12`, + `2·Si(1)`, `Σ 1/(n·n!)` and `Cin(1)` to better than `10⁻⁹` wide. + Three guards keep a genuine pole out: the two zeros are established + **symbolically** (substitute the exact rational `p`, simplify, require a + literal zero — no numeric enclosure can prove a value is exactly zero, so none + is asked to), `D′` must be certified non-vanishing on the sub-interval, and + `N` and `D` must each enclose successfully over it, which is what certifies + they are analytic and hence that the symbolic derivatives are the real ones. + `1/x`, `sin(x)/x²` and `(x−p)²/(x−p)³` on boxes containing the pole are all + still refused. +- **A `verified_integral` refusal now says whether the *integral* fails to exist + or only the enclosure of the *integrand* does.** An integrable singularity + that is not removable — `∫₀¹ −ln x dx = 1`, `∫₀¹ (ln x)² dx = 2`, + `∫₀¹ dx/√(1−x²) = π/2`, `∫₀¹ xˣ dx`, `∫₀¹ ln(x)·ln(1−x) dx` — is still + refused, because no rigorous bound on the singular tail can be derived from + the expression alone today. But the `E-VALIDATED-003` message now names the + location (left endpoint, right endpoint or interior, with the approximate + coordinate), reports the underlying cause, and states explicitly that an + integrable singularity still has a finite integral that this routine cannot + certify. Widening an enclosure to make those cases "pass" would have broken + the contract that makes the module worth using, so they refuse honestly + instead. +- **The ODE/DAE guide documented an API that does not exist.** + [`docs/mdbook/src/ode-dae.md`](docs/mdbook/src/ode-dae.md) showed keyword + constructors — `ODE(state=…, derivatives=…, independent=…)`, + `DAE(equations=…, variables=…, independent=…)` — and a one-argument + `lower_to_first_order(higher_order_ode)`, none of which were ever real: the + actual calls are `DAE.new(equations, variables, derivatives, time_var)` with + derivatives as *separate symbols* like `pool.symbol("dx/dt")`, and + `lower_to_first_order(var, rhs, order, time_var)`. It also printed + `reduced.differentiated` on the `pantelides` result, which the returned `DAE` + does not have. Following the page failed on its first line, and it was the + only documentation there was, because `DAE.new`, `ODE.new`, + `rosenfeld_groebner` and the `RosenfeldGroebnerResult` accessors had no + docstrings at all. The page is rewritten around calls that run — every block + is executed by `tests/test_docs_ode_dae.py` — those docstrings now exist, + `reduced.index` works as documented (`DAE.index` is exposed), and + `.differentiated` is gone rather than faked. - **`zeilberger` no longer refuses a constant base just because it is not already a literal.** `(-one)**(n+k)`, with `one = pool.integer(1)`, builds the node `Mul(1, -1)` — the pool does no arithmetic at construction — and the @@ -436,6 +369,309 @@ loop fails. 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. + +### Added + +- **Gröbner results can be read back — `GbPoly.to_expr`, iteration over a + `GroebnerBasis`, and `expr_to_gbpoly`.** Everything that returned a basis + returned a handle nobody could open. `GbPoly` exposed only `is_zero` and + `n_vars`; `GroebnerBasis` exposed only its own constructors plus `reduce` and + `contains`, and `reduce` took a `GbPoly` that no exported function could + build — `expr_to_gbpoly`, named in `compute_raw`'s own docstring, was never + registered on the module. So `rosenfeld_groebner(...).final_basis()`, + `triangularize(...)`, `primary_decomposition(...)` and a parametric `solve` + all handed back objects whose only readable property was how many generators + they had. Differential elimination was write-only: the input–output equations + it computes could not be looked at, which is the whole of structural + identifiability. + + Now: `alkahest.expr_to_gbpoly(expr, vars)` converts in, `GbPoly.to_expr()` + converts back out, and `GbPoly.terms()` gives `(exponent tuple, exact + int/Fraction)` pairs. A `GroebnerBasis` is a sequence — `len()`, indexing, + iteration — with `polynomials()`, `to_exprs()`, `variables()` and an `order` + property. `reduce()` now accepts an `Expr` as well as a `GbPoly`, so the + membership and reduction API is reachable from expressions alone. + `GroebnerBasis.eliminate(vars)` is bound too — the mdbook and Sphinx pages + had documented it for releases, but it existed only in Rust, so the + implicitization example on the solving page could not run. + + The part that needed a real fix rather than an accessor is the **variable + ordering**. A `GbPoly` stores exponent vectors, not names, so a basis without + its variable list cannot be read at all — and `rosenfeld_groebner` discovered + its jet variables internally (`t`, `x`, `dx/dt`, `ddx/dt/dt`, …) and threw + them away. Every object that hands out a `GbPoly` now carries that list: + `RosenfeldGroebnerResult.variables()`, `GroebnerBasis.variables()`, + `RegularChain.variables()`, and for a parametric `solve` the solve variables + *followed by the free parameters*, which is the order the exponent vectors + were actually built in. Asking for an `Expr` with too few variables named + raises `ValueError` rather than silently misreading exponent slots. + + New in Rust: `alkahest_cas::gbpoly_to_expr`, `GroebnerBasis::order()`, + `MonomialOrder::as_str()`, `solver::collect_parameters`, and + `rosenfeld_groebner_ranked` / `dae_index_reduce_ranked`, which return the + `DifferentialRanking` alongside the result (the existing entry points are + unchanged wrappers, so no struct gained a field). + +- **`DAE` can be read: `equations()`, `variables()`, `derivatives()`, + `time_var`, `index`.** It previously exposed `n_equations` and `n_variables` + and nothing else, so a prolonged system that reported six variables for a + two-variable input gave no way to find out what any of them were. + `pantelides` sets `index` on the DAE it returns — the number of + differentiation rounds — and the equations it appended, plus the higher jets + (`ddx/dt/dt`, …) they introduced, are now visible in `equations()` and + `derivatives()`. + +### Performance + +- **`zeilberger`'s exact `Q(n)(k)` post-processing no longer swells its own + coefficients.** With the search fixed (below), what was left was entirely + after it: on `Σ_k C(n,k)³` the search reached `(order 2, degree 3)` in 0.22 s + and the run then spent ~29 s normalising the certificate and re-verifying it. + The cause was `PolyK::gcd` — a textbook Euclidean remainder sequence over the + *field* `Q(n)`, whose coefficients are rational functions in `n`: every + division step adds numerator and denominator degrees and no step ever removes + content, the classic intermediate-expression-swell blowup, and + `RatK::normalize` ran it on every normalisation. The gcd now leaves the field + and runs **Brown's subresultant PRS in `Z[n][k]`** (Collins 1967, Brown 1971; + Knuth TAOCP 2 § 4.6.1), with both cofactors divided out in the same integral + domain, and `Q[n]` gcds (`rn_mul` / `rn_add` / `rn_inv`, which cancel + crosswise now rather than reducing the full cross-multiplied product) go + through the same subresultant sequence over `Z[n]`. At the shipped defaults, + measured before and after on one machine: `Σ (−1)^k C(n,k)³` **1.6 s → + 0.11 s** (15×), `Σ_k C(n,k)³` **56 s → 0.07 s** (800×), + `Σ_k C(n,k)²C(n+k,k)²` **16.5 s → 0.05 s** (330×, and still Apéry's + recurrence coefficient for coefficient). Two OEIS targets that timed out past 300 s at certificate + degree ≥ 3 are now decided: **A357510** `Σ k·C(n,k)²·C(n+k,k)²` and + **A357512** `Σ k⁵·C(n,k)²·C(n+k,k)²` both yield a verified order-3 recurrence + in under a second. This is a change of algorithm, not of contract: a monic + gcd is unique, so every certificate is the same one as before, and every + certificate is still checked as an exact `Q(n)(k)` identity before it is + returned — nothing here is probabilistic and no verification was weakened. +- **`zeilberger`'s `max_order` / `max_degree` are now upper bounds instead of + starting points.** The search used to sweep certificate degrees + `d = 0..=max_degree` at order 1 before ever trying order 2, and a single + degree probe gets ~3× more expensive per degree step (measured on + `Σ (−1)^k C(n,k)³`: 0.7 ms at `d = 0`, 0.6 s at `d = 7`, 84 s at `d = 12`). + Every order ≥ 2 identity — Dixon, Franel, Apéry — therefore ran for minutes + or never at the shipped defaults while being seconds away at `max_degree=4`, + i.e. **raising the bound made easy inputs slower rather than admitting harder + ones**. The `(order, degree)` grid is now visited by iterative deepening, + cheapest estimated probe first (one extra order is priced at three extra + degrees, which is what the measurements say), and the first *verified* + relation is returned. `Σ (−1)^k C(n,k)³` at the defaults goes from >400 s + (killed) to **0.67 s**; `Σ_k C(n,k)³` from >400 s (killed) to ~31 s. Nothing + is skipped — the plan still visits every pair inside the bounds, so an + exhausted search costs what it always did — and verification is unchanged: a + candidate that fails the exact `Q(n)(k)` check is still discarded, never + returned. + +### Testing + +- The deterministic silent-error gate grew from **213 to 241 scored cases** + (`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. + +## 3.8.0 — 2026-08-12 + +### Silent errors fixed — do results you already computed need rechecking? + +A *silent error* is a confident, plausible, mathematically wrong answer with no +exception, no `NaN` and no verification flag. Eleven were found and fixed this +release. **Eight of them shipped in 3.7 or earlier**, so if you have results +from an affected call, re-run them. The other three were in code added during +this release cycle and never reached a published wheel. + +| Affected call | Wrong answer it gave | First shipped in | Recheck? | +|---|---|---|---| +| `decide(Forall(x, φ))` where the counterexample is a rational root whose denominator is **not a power of two** | `True` for a **false** universal theorem, e.g. `∀x. (3x+2)² > 0` (false at `x = −2/3`) | ≤ 3.7 | **Yes** — any `decide` verdict | +| `decide(Exists(x, φ))` with an `=` atom | `(True, witness)` where the witness does **not** satisfy the sentence, e.g. `∃x. 3x−2 = 0 → x = 1/2` | ≤ 3.7 | **Yes** — any cited witness | +| `Matrix.nullspace()` on a 2×2 with a symbolic determinant | A confident wrong kernel basis; `[[x,0],[0,1]]` returned `(0, x)`, for which `M·v ≠ 0` | 3.7 | **Yes** — verify `M·v = 0` numerically | +| `simplify` / `simplify_egraph` on a product containing `0⁻¹` | `1`, or `0`, depending on the engine — for an expression with no value at all. Reachable from `diff(2/(x − x), x)` | ≤ 3.7 | Yes, if any input could reduce to `0⁻¹` | +| `decide` on a two-variable sentence true only at an irrational point | `False` for a satisfiable `∃x∃y`, and `True` for its false `∀x∀y` dual | this cycle (2-var `decide` is new) | No published release affected | +| `batch_map(..., parallel=True)` under `context(budget=…)` | Ran **unbudgeted**, so candidates a sequential sweep reported as `E-BUDGET-001` came back as `E-INT-001` — a *mathematical* verdict | this cycle (batch APIs are new) | No published release affected | +| `product_definite` on a term with any non-integer coefficient | Off by `c^(hi−lo+1)`: `Π_{k=1}^{5} ½` returned `1` instead of `1/32`, `Π (2k−1)/(2k)` at `n = 6` returned `14.4375` instead of `0.2255859375` | ≤ 3.7 | **Yes** — any `product_definite` / `product_indefinite` result | +| `sum_definite` where the summand has a pole strictly *between* the bounds | A clean finite number for a sum with an undefined term: `Σ_{k=1}^{10} 1/((k−3)(k−2))` returned `−5/8` | ≤ 3.7 | **Yes** — any `sum_definite` over a range containing a denominator root | +| `euler_maclaurin` when `corrections` is too small for the summand | A fabricated additive constant — the missing term frozen at the fitting point. `Σ k⁹` at the default `corrections = 2` acquired `34359738368 = 512⁴/2` in a Faulhaber polynomial whose constant term is `0` | this cycle (Euler–Maclaurin is new) | No published release affected | +| `rsolve` on a **forward-shift** spelling with a non-zero right-hand side | The solution of a *different* equation: `f(n+1) − f(n) = n²` with `f(0) = 0` returned `Σ_{j=1}^{n} j²` instead of `Σ_{j=0}^{n−1} j²` | ≤ 3.7 | **Yes** — any `rsolve` written with `f(n+i)`, `i > 0`, and an inhomogeneous term | +| `rsolve` / `solve_linear_recurrence_homogeneous` on an order-2 recurrence with a **repeated** characteristic root | `C₀·rⁿ + C₁·rⁿ` — a one-parameter family presented as the general solution of a second-order equation, losing the `n·rⁿ` branch | ≤ 3.7 | **Yes** — check the discriminant of `r² + b r + c` | + +Also fixed, and not a silent error but worse for an unattended loop: a Rust +panic escaped `interval_eval` as `pyo3_runtime.PanicException`, which inherits +from `BaseException` and therefore slips past `except Exception`. Shipped in +3.7 — a loop that survived everything else died on `x^(3/2)` over a negative +ball. + +The deterministic silent-error gate (`tests/silent_errors/`, Tier-1 CI) now +scores **0 silent errors out of 213 scored cases** across evaluation, +integration, limits, linear algebra, number theory, real QE, series, +simplification, solving, and sums/products. Every trap added this cycle was +re-run against a build with the fix reverted and confirmed to score +`silent_error` there, and every trap is paired with a **control** — its nearest +convergent neighbour — so a subsystem cannot pass the gate by refusing +everything. That is a statement about the corpus, not a guarantee about the +library. + +### Behaviour changes to plan for + +Fixing a silent error means some calls that used to return now refuse. Every one +of these is a call whose previous answer was not justified: + +- **`decide` raises `CadError` (`E-CAD-001`) where it used to answer**, whenever + the formula has a non-strict atom (`=`, `≠`, `≤`, `≥`) and a boundary root has + not been shown rational. This includes mixed-alternation sentences that route + through De Morgan — `∀x∃y. p > 0` becomes `¬∃x∀y. p ≤ 0`, and the negation + makes a strict body non-strict. `decide` is **not** a complete decision + procedure in this implementation; treat `E-CAD-001` as *undecided*, never as + *false*. +- **`rank`, `rref`, `nullspace`, `eigenvects`, `jordan_form` raise + `E-LINALG-010`, and `inverse` raises the new `E-MAT-004`**, when an entry's or + the determinant's vanishing can be decided neither way. Previously "could not + prove non-zero" was silently read as "zero". +- **`simplify` leaves `0 · 0⁻¹` unevaluated** instead of returning `1` (or `0`). + A result containing `(0 * 0^-1)` is Alkahest declining to give an + indeterminate form a value, not a simplifier failure. +- **`sum_definite` raises `SumError` (`E-SUM-003`) when the summand is undefined + at an integer inside `[lo, hi]`**, not only when the pole lands on `lo` or + `hi+1`. The refusal names the offending index. Sums whose poles lie outside + the range are unaffected: `Σ_{k=4}^{10} 1/((k−3)(k−2))` still returns `7/8`. +- **`euler_maclaurin` may return a shorter expansion, with no additive + constant.** The constant is now fitted at a point outside the gate's check + points and re-fitted at a second one; if the two disagree it is not a constant + and none is claimed. The report says which way that went in `derivation`, and + the `"fitted numerically"` hypothesis is only listed when a fitted constant is + actually part of the answer. Genuine constants (`γ`, `ζ(2)`, `½log 2π`, …) are + unaffected — they agree across fitting points to 13+ digits. +- **`product_definite(term, k, lo, hi)` with `lo > hi` returns `1` even for a + zero term.** The empty product takes no factors; it previously returned `0` + for `Π_{k=1}^{0} 0` while returning `1` for `Π_{k=1}^{0} k`. +- **`capabilities()["contract_version"]` is `3`, and `features` lost two keys: + `groebner_cuda` and `numpy`.** Indexing either now raises `KeyError`; use + `features.get(name, False)` if you need to span versions. Both were removed + rather than wired up because neither was *falsifiable* — no observation a + Python caller could make distinguished `True` from `False`: + - `groebner_cuda` reported that the CUDA Macaulay-matrix kernel had been + compiled in. The string `groebner_cuda` occurred exactly once anywhere in + `alkahest-py` — the capability line itself. There was no binding, no + `*gpu*` name in the public or the private module, and `GroebnerBasis` + exposes only CPU methods. The kernel is unchanged and still reachable from + Rust as `alkahest_cas::poly::groebner::compute_groebner_basis_gpu`; if + dispatch ever prefers it, the binding lands first and a bit follows it. + - `numpy` mapped to a Cargo feature gating the `numpy` crate, which + `alkahest-py` never used an item from. The feature and the dependency are + both gone. `ak.numpy_eval` and `ak.numpy_eval_par` are unaffected — they go + through the buffer protocol and always worked with the bit `False`, which + is its value on every wheel ever published. + + An unreachable `True` makes a caller trust something it should not, which is + the same class of defect as a silent wrong answer; a bit that correlates with + nothing is better removed than left to be misread. + `tests/test_agent_contract.py::test_every_advertised_feature_has_an_entry_point` + now walks `features` and fails on any key without a named, reachable entry + point, so the next one cannot ship. +- **Rust, `--features groebner-cuda`: `compute_groebner_basis_gpu` and + `reduce_batch` return `(polys, GpuBackendReport)` instead of `polys`.** Both + fall back to CPU row reduction — when `device_id` is `None`, and when the + driver fails — and the basis is identical either way, so a caller previously + had no way to tell a GPU run from a CPU one. `GpuBackendReport::ran_on_gpu()` + is true only when at least one mod-p reduction ran on a device and none fell + back; `reductions_on_gpu`, `reductions_on_cpu` and `first_gpu_error` carry + the detail. A compile error on upgrade is the intended failure mode for code + that was recording these results as GPU results. Nothing at the Python + surface changes: the feature has no binding. +- **`residue(f, z, point)` refuses a non-constant `point` with + `AlkahestError` / `E-RESIDUE-005`** instead of leaking + `AttributeError: 'Expr' object has no attribute 'numerator'` from the + argument parser. `AttributeError` is not an `AlkahestError`, so + `except ak.AlkahestError` missed it entirely. The existing `E-RESIDUE-001..4` + refusals are now `AlkahestError`s carrying `.code` and `.remediation` too, + rather than bare `ValueError`s with the code glued into the message; + `AlkahestError` subclasses `ValueError`, so `except ValueError` still works. +- **`series` refuses instead of running forever, with `SeriesError` / + `E-SERIES-003`.** `series(sqrt(t**-2 + t**-1), t, 0, 32)` never returned: + coefficients are formed by repeated differentiation without re-simplifying, so + a nested radical's derivatives grow by a constant factor per coefficient and + the cost doubles per order. It now honours an active `Budget` (raising + `BudgetExceededError`) and, with none, an internal work ceiling. It never + returns a *shorter* series: `O(h^order)` on fewer coefficients than were asked + for is a false statement about the remainder, which is worse than the refusal. + Ordinary expansions are unaffected — the heaviest in the suites intern a few + thousand nodes against a ceiling of 50 000. +- **`simplify_expanded` records a derivation step when its expansion bound stops + it** (`expand_pow_limit_reached`, a no-op step naming the power it declined), + and the bound itself is now a budget on the number of distributed products + rather than a flat exponent cap. `(x+y)**6` and `(x+y+1)**7` now expand where + the exponent-only cap refused them while permitting a twenty-term sum to the + fourth power; anything above the budget comes back unexpanded *and says so* + instead of looking like an expression that was already expanded. + +### Known limits — documented, not fixed + +These are properties of the design as it stands. They are called out here +because 3.8 is aimed at long unattended loops, and each of them is a way such a +loop fails. + +- **`ExprPool` never reclaims.** The arena is append-only: no `clear`, no + refcount, no GC, and the storage cannot shrink. The only way to free interned + nodes is to **drop the whole pool** — and every `Expr`, `Matrix`, `Series` and + `DerivedResult` holds a *strong* reference to its pool, so retaining one + interesting result retains everything. Growth on a shared pool is linear and + unbounded (~200 bytes/node; measured ~2 KB per `integrate` call over 20 000 + calls, 0 B/call with a fresh pool per iteration) while per-call latency stays + **flat**, so the failure mode is a clean OOM with no slowdown to warn you + first. `ExprPool` also exposes no `__len__` or `stats()`, so the growth is not + observable from Python. The supported pattern is **one pool per problem**, + documented in [`budgets.md`](docs/mdbook/src/budgets.md#exprpool-never-reclaims). +- **`run_with_wall_fallback` does not bound wall time for an uncooperative + callee.** It joins its worker before the exception propagates, so it returns + when the callee returns: `run_with_wall_fallback(time.sleep, 3.0, + budget=Budget(wall_ms=50))` raises `E-BUDGET-001` after 3000 ms, and the + message reports the real elapsed time so this shows up in a log rather than + being inferred later. Python cannot kill a thread, and abandoning one would + leak a live thread that still allocates into the pool and can only be stopped + through the process-wide cancel flag. Only an **OS-level bound** (subprocess, + process watchdog) is a hard deadline. +- **`wall_ms` granularity is one primitive operation, and FLINT calls cannot be + interrupted.** After the checkpoint work above the overshoot is a small + additive term (1.0–1.2×), but past a certain degree a single operation is a + FLINT factorisation or resultant — one foreign-function call, ~2 s on a + degree-62 integrand, which no cooperative mechanism can stop part-way. +- **`Matrix.eigenvals()` grows the pool on identical input** (~1.9 KB/call, + measured over 20 000 calls on the same 2×2 integer matrix): it interns a fresh + `__eigen_lambda_N` gensym per call. Every other Python-facing entry point + measured is flat on repeated input. Cache eigenvalue results. +- **`Matrix.eigenvals()` can emit casus-irreducibilis cube roots** — correct + under Alkahest's real cube-root convention (and honestly refused by + `eval_expr` with `E-EVAL-009`, with `interval_eval` returning an unbounded + ball) but evaluated on the **principal** branch by SymPy, NumPy and most other + tools, which return a confident number that is not an eigenvalue. 14 of 720 + random integer matrices produced one. An honest refusal here becomes somebody + else's silent error the moment the expression crosses the boundary, so + evaluate inside Alkahest before exporting, or export a verified numeric + enclosure instead. See [`interop.md`](docs/mdbook/src/interop.md). +- **The LLVM JIT leaks an LLVM `Context` per compile** (`Box::leak`, on the + error paths as well as the success path). Feature-gated behind `jit`, so + default PyPI wheels (Cranelift) are unaffected; do not compile in a loop under + a `+jit` / `+full` wheel. +- **No sanitizer covers any Python-facing path.** The PR-gating ASan job runs + with `detect_leaks=0`, the nightly LSan shard cannot reach a `cdylib` with no + `#[test]` functions, and `pytest` is never run under a sanitizer. The + behavioural substitute is the fresh-pool sweep described in + [`TESTING.md`](TESTING.md#3-memory-safety--sanitizers). +- **There is no `cuda_device_count()`.** `CudaCompiledFn.call_batch_on(ordinal, + …)` selects a device, but the valid range can only be discovered by trying an + ordinal and catching `CudaError` (`E-CUDA-003`); the loop that does it is in + [`gpu.md`](docs/mdbook/src/gpu.md#discovering-the-valid-device-ordinals). Not + added yet on purpose: `cuda` implies LLVM 15 with NVPTX, so such a binding + cannot be compiled on an ordinary dev box, no CI job builds the Python + extension with either CUDA feature, and exercising it needs a device — it + would ship with no verification of any kind, which is the provenance of the + capability overclaims fixed above. It belongs in the same change as the + missing `maturin develop --features cuda` + `pytest tests/test_cuda.py` + nightly step. + +### Fixed + - **`cargo test --features groebner-cuda` could not pass on a machine with no NVIDIA driver**, contradicting the header comment of `alkahest-core/tests/groebner_cuda.rs`. `cudarc` *panics* rather than @@ -990,47 +1226,6 @@ loop fails. ### Performance -- **`zeilberger`'s exact `Q(n)(k)` post-processing no longer swells its own - coefficients.** With the search fixed (below), what was left was entirely - after it: on `Σ_k C(n,k)³` the search reached `(order 2, degree 3)` in 0.22 s - and the run then spent ~29 s normalising the certificate and re-verifying it. - The cause was `PolyK::gcd` — a textbook Euclidean remainder sequence over the - *field* `Q(n)`, whose coefficients are rational functions in `n`: every - division step adds numerator and denominator degrees and no step ever removes - content, the classic intermediate-expression-swell blowup, and - `RatK::normalize` ran it on every normalisation. The gcd now leaves the field - and runs **Brown's subresultant PRS in `Z[n][k]`** (Collins 1967, Brown 1971; - Knuth TAOCP 2 § 4.6.1), with both cofactors divided out in the same integral - domain, and `Q[n]` gcds (`rn_mul` / `rn_add` / `rn_inv`, which cancel - crosswise now rather than reducing the full cross-multiplied product) go - through the same subresultant sequence over `Z[n]`. At the shipped defaults, - measured before and after on one machine: `Σ (−1)^k C(n,k)³` **1.6 s → - 0.11 s** (15×), `Σ_k C(n,k)³` **56 s → 0.07 s** (800×), - `Σ_k C(n,k)²C(n+k,k)²` **16.5 s → 0.05 s** (330×, and still Apéry's - recurrence coefficient for coefficient). Two OEIS targets that timed out past 300 s at certificate - degree ≥ 3 are now decided: **A357510** `Σ k·C(n,k)²·C(n+k,k)²` and - **A357512** `Σ k⁵·C(n,k)²·C(n+k,k)²` both yield a verified order-3 recurrence - in under a second. This is a change of algorithm, not of contract: a monic - gcd is unique, so every certificate is the same one as before, and every - certificate is still checked as an exact `Q(n)(k)` identity before it is - returned — nothing here is probabilistic and no verification was weakened. -- **`zeilberger`'s `max_order` / `max_degree` are now upper bounds instead of - starting points.** The search used to sweep certificate degrees - `d = 0..=max_degree` at order 1 before ever trying order 2, and a single - degree probe gets ~3× more expensive per degree step (measured on - `Σ (−1)^k C(n,k)³`: 0.7 ms at `d = 0`, 0.6 s at `d = 7`, 84 s at `d = 12`). - Every order ≥ 2 identity — Dixon, Franel, Apéry — therefore ran for minutes - or never at the shipped defaults while being seconds away at `max_degree=4`, - i.e. **raising the bound made easy inputs slower rather than admitting harder - ones**. The `(order, degree)` grid is now visited by iterative deepening, - cheapest estimated probe first (one extra order is priced at three extra - degrees, which is what the measurements say), and the first *verified* - relation is returned. `Σ (−1)^k C(n,k)³` at the defaults goes from >400 s - (killed) to **0.67 s**; `Σ_k C(n,k)³` from >400 s (killed) to ~31 s. Nothing - is skipped — the plan still visits every pair inside the bounds, so an - exhausted search costs what it always did — and verification is unchanged: a - candidate that fails the exact `Q(n)(k)` check is still discarded, never - returned. - **`numpy_eval` / `numpy_eval_par` no longer round-trip through a Python list of floats.** The previous implementation converted every NumPy array to a flat Python list via `.tolist()` before crossing into Rust diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index a6a55bc0..0908791d 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -80,6 +80,31 @@ Rules live in `alkahest-core/src/simplify/`. Each rule is a `RewriteRule` with a - Add the rule to the appropriate rule set (`arithmetic_rules`, `trig_rules`, `log_exp_rules_safe`, etc.). - Add a proptest case verifying the rule is idempotent: `simplify(simplify(expr)) == simplify(expr)`. +## Accessors: property or method? + +One rule, applied to every `#[pymethods]` entry in `alkahest-py/src/lib.rs`: + +> **A zero-argument, O(1), non-allocating accessor that returns a scalar or a flag is a `#[getter]` (a Python property). Anything that returns a collection, allocates, or does real work is a method.** + +```rust +#[getter] // property: reads a field, cannot fail +fn n_equations(&self) -> usize { self.inner.n_equations() } + +fn polys(&self) -> Vec { … } // method: allocates a collection +fn rank(&self, py: Python<'_>) -> PyResult { … } // method: real work, can fail +``` + +A property and a method sitting side by side on the same class is not in itself a problem — `RegularChain.n_vars` (property) next to `RegularChain.polys()` (method) is exactly what the rule asks for. What the rule rules out is the *same* kind of question being asked two different ways on two different classes, which is what made the surface unpredictable before 3.8.0. + +Why the split falls there: + +- A property that can raise, block, or take a noticeable amount of time is a trap — the caller reads `x.rank` as a field access. Real work stays behind parentheses. +- A method that returns a scalar is the more dangerous mistake in the other direction: `if x.n_equations:` on a bound method is always `True` and `f"{x.n_equations}"` prints ``. Neither raises. Converting these was the whole point of the 3.8.0 sweep. + +`tests/test_accessor_convention.py` enforces this. It pins the converted accessors at runtime and statically scans `alkahest-py/src/lib.rs` for zero-argument scalar-returning methods; if a new one is genuinely doing real work, add it to `REAL_WORK_EXEMPTIONS` there with a one-line reason. A handful of pre-3.8.0 getters return small collections (`AsymptoticReport.terms`, `CertifiedSolution.coordinates`, `PositivityCertificate.log`, …); they are grandfathered, not precedent. + +Changing an existing accessor's form is a **breaking change**: record it in `CHANGELOG.md` under the release's "Behaviour changes to plan for" with a before/after line, and update every caller in `tests/`, `examples/`, `benchmarks/`, `docs/mdbook/`, `alkahest-skill/alkahest.md` and the `.pyi` stubs. + ## Pull requests - Keep PRs focused on one item from `ROADMAP.md` or one issue. diff --git a/Cargo.lock b/Cargo.lock index 4062bad3..336f098c 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -13,7 +13,7 @@ dependencies = [ [[package]] name = "alkahest-cas" -version = "3.8.0" +version = "3.9.0" dependencies = [ "bitflags 2.12.1", "boxcar", @@ -38,7 +38,7 @@ dependencies = [ [[package]] name = "alkahest-mlir" -version = "3.8.0" +version = "3.9.0" dependencies = [ "alkahest-cas", "proptest", @@ -46,7 +46,7 @@ dependencies = [ [[package]] name = "alkahest-py" -version = "3.8.0" +version = "3.9.0" dependencies = [ "alkahest-cas", "pyo3", diff --git a/Cargo.toml b/Cargo.toml index 010f5967..f4e29ae7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,7 +4,7 @@ exclude = ["fuzz"] resolver = "2" [workspace.package] -version = "3.8.0" +version = "3.9.0" edition = "2021" authors = ["Alkahest Contributors"] license = "Apache-2.0" diff --git a/alkahest-core/src/diffalg/mod.rs b/alkahest-core/src/diffalg/mod.rs index a33be0a7..b152ae55 100644 --- a/alkahest-core/src/diffalg/mod.rs +++ b/alkahest-core/src/diffalg/mod.rs @@ -240,6 +240,22 @@ pub fn rosenfeld_groebner_with_options( order: MonomialOrder, max_prolong_rounds: usize, ) -> Result { + rosenfeld_groebner_ranked(dae, pool, order, max_prolong_rounds).map(|(result, _)| result) +} + +/// [`rosenfeld_groebner_with_options`] plus the jet [`DifferentialRanking`] that +/// indexes the exponent vectors of [`RosenfeldGroebnerResult::final_basis`]. +/// +/// The elimination result is unreadable without this: a [`GbPoly`] stores only +/// exponent vectors, so `ranking.vars[i]` is what exponent slot `i` refers to. +/// Pair it with [`crate::solver::gbpoly_to_expr`] to recover the input–output +/// equations as [`ExprId`]s. +pub fn rosenfeld_groebner_ranked( + dae: &DAE, + pool: &ExprPool, + order: MonomialOrder, + max_prolong_rounds: usize, +) -> Result<(RosenfeldGroebnerResult, DifferentialRanking), DiffAlgError> { if dae.equations.is_empty() { return Err(DiffAlgError::EmptySystem); } @@ -256,14 +272,17 @@ pub fn rosenfeld_groebner_with_options( for round in 0..max_prolong_rounds { let gb = GroebnerBasis::compute(active.clone(), order); if is_unit_ideal_gb(&gb) { - return Ok(RosenfeldGroebnerResult { - consistent: false, - chains: vec![], - working_dae: work, - final_basis: None, - prolongation_rounds, - truncated: false, - }); + return Ok(( + RosenfeldGroebnerResult { + consistent: false, + chains: vec![], + working_dae: work, + final_basis: None, + prolongation_rounds, + truncated: false, + }, + DifferentialRanking { vars }, + )); } let mut next_prolong = Vec::with_capacity(prolong_exprs.len()); @@ -305,14 +324,17 @@ pub fn rosenfeld_groebner_with_options( } else { vec![] }; - return Ok(RosenfeldGroebnerResult { - consistent, - chains, - working_dae: work, - final_basis: if consistent { Some(final_basis) } else { None }, - prolongation_rounds, - truncated: false, - }); + return Ok(( + RosenfeldGroebnerResult { + consistent, + chains, + working_dae: work, + final_basis: if consistent { Some(final_basis) } else { None }, + prolongation_rounds, + truncated: false, + }, + DifferentialRanking { vars }, + )); } active.extend(to_add); @@ -328,33 +350,39 @@ pub fn rosenfeld_groebner_with_options( } else { vec![] }; - return Ok(RosenfeldGroebnerResult { - consistent, - chains, - working_dae: work, - final_basis: if consistent { Some(final_basis) } else { None }, - prolongation_rounds, - truncated: true, - }); + return Ok(( + RosenfeldGroebnerResult { + consistent, + chains, + working_dae: work, + final_basis: if consistent { Some(final_basis) } else { None }, + prolongation_rounds, + truncated: true, + }, + DifferentialRanking { vars }, + )); } } let final_basis = GroebnerBasis::compute(active, order); let consistent = !is_unit_ideal_gb(&final_basis); - Ok(RosenfeldGroebnerResult { - consistent, - chains: if consistent { - vec![RegularDifferentialChain { - basis: final_basis.clone(), - }] - } else { - vec![] + Ok(( + RosenfeldGroebnerResult { + consistent, + chains: if consistent { + vec![RegularDifferentialChain { + basis: final_basis.clone(), + }] + } else { + vec![] + }, + working_dae: work, + final_basis: if consistent { Some(final_basis) } else { None }, + prolongation_rounds, + truncated: true, }, - working_dae: work, - final_basis: if consistent { Some(final_basis) } else { None }, - prolongation_rounds, - truncated: true, - }) + DifferentialRanking { vars }, + )) } /// Calls [`rosenfeld_groebner_with_options`] with the default maximum prolongation rounds. @@ -372,16 +400,29 @@ pub fn dae_index_reduce( pool: &ExprPool, order: MonomialOrder, ) -> Result { + dae_index_reduce_ranked(dae, pool, order).map(|(r, _)| r) +} + +/// [`dae_index_reduce`] plus the jet [`DifferentialRanking`] for the Gröbner +/// fallback — `None` when Pantelides succeeded and no basis was built. +pub fn dae_index_reduce_ranked( + dae: &DAE, + pool: &ExprPool, + order: MonomialOrder, +) -> Result<(DaeIndexReduction, Option), DaeError> { match pantelides(dae, pool) { - Ok(p) => Ok(DaeIndexReduction::Pantelides(p)), + Ok(p) => Ok((DaeIndexReduction::Pantelides(p), None)), Err(DaeError::IndexTooHigh) => { - let r = rosenfeld_groebner(dae, pool, order).map_err(|e| match e { - DiffAlgError::DiffError(s) | DiffAlgError::NotPolynomial(s) => { - DaeError::DiffError(s) - } - DiffAlgError::EmptySystem => DaeError::StructurallyInconsistent, - })?; - Ok(DaeIndexReduction::Rosenfeld(r)) + let (r, ranking) = + rosenfeld_groebner_ranked(dae, pool, order, DEFAULT_MAX_PROLONG_ROUNDS).map_err( + |e| match e { + DiffAlgError::DiffError(s) | DiffAlgError::NotPolynomial(s) => { + DaeError::DiffError(s) + } + DiffAlgError::EmptySystem => DaeError::StructurallyInconsistent, + }, + )?; + Ok((DaeIndexReduction::Rosenfeld(r), Some(ranking))) } Err(e) => Err(e), } diff --git a/alkahest-core/src/lib.rs b/alkahest-core/src/lib.rs index 3eabedc2..5155ed2a 100644 --- a/alkahest-core/src/lib.rs +++ b/alkahest-core/src/lib.rs @@ -204,9 +204,10 @@ pub use lean::{ // V2-1 — Modular / CRT framework #[cfg(feature = "groebner")] pub use diffalg::{ - dae_index_reduce, rosenfeld_groebner, rosenfeld_groebner_algebraic, - rosenfeld_groebner_with_options, DaeIndexReduction, DiffAlgError, DifferentialIdeal, - DifferentialRanking, DifferentialRing, RegularDifferentialChain, RosenfeldGroebnerResult, + dae_index_reduce, dae_index_reduce_ranked, rosenfeld_groebner, rosenfeld_groebner_algebraic, + rosenfeld_groebner_ranked, rosenfeld_groebner_with_options, DaeIndexReduction, DiffAlgError, + DifferentialIdeal, DifferentialRanking, DifferentialRing, RegularDifferentialChain, + RosenfeldGroebnerResult, }; #[cfg(feature = "groebner")] pub use ideal::{ @@ -224,10 +225,10 @@ pub use number_theory::{ pub use primitive::{Capabilities, CoverageReport, CoverageRow, Primitive, PrimitiveRegistry}; #[cfg(feature = "groebner")] pub use solver::{ - diophantine, expr_to_gbpoly, extract_regular_chain_from_basis, main_variable_recursive, - solve_numerical, solve_polynomial_system, solve_transcendental, triangularize, CertifiedPoint, - DiophantineError, DiophantineSolution, HomotopyError, HomotopyOpts, RegularChain, Solution, - SolutionSet, SolverError, TranscendentalOutcome, + diophantine, expr_to_gbpoly, extract_regular_chain_from_basis, gbpoly_to_expr, + main_variable_recursive, solve_numerical, solve_polynomial_system, solve_transcendental, + triangularize, CertifiedPoint, DiophantineError, DiophantineSolution, HomotopyError, + HomotopyOpts, RegularChain, Solution, SolutionSet, SolverError, TranscendentalOutcome, }; pub fn version() -> &'static str { @@ -252,9 +253,10 @@ pub mod stable { pub use crate::diff::{diff, diff_forward, grad, DiffError}; #[cfg(feature = "groebner")] pub use crate::diffalg::{ - dae_index_reduce, rosenfeld_groebner, rosenfeld_groebner_algebraic, - rosenfeld_groebner_with_options, DaeIndexReduction, DiffAlgError, DifferentialIdeal, - DifferentialRanking, DifferentialRing, RegularDifferentialChain, RosenfeldGroebnerResult, + dae_index_reduce, dae_index_reduce_ranked, rosenfeld_groebner, + rosenfeld_groebner_algebraic, rosenfeld_groebner_ranked, rosenfeld_groebner_with_options, + DaeIndexReduction, DiffAlgError, DifferentialIdeal, DifferentialRanking, DifferentialRing, + RegularDifferentialChain, RosenfeldGroebnerResult, }; pub use crate::errors::AlkahestError; pub use crate::eval::{ @@ -319,10 +321,10 @@ pub mod stable { }; #[cfg(feature = "groebner")] pub use crate::solver::{ - diophantine, expr_to_gbpoly, extract_regular_chain_from_basis, main_variable_recursive, - solve_numerical, solve_polynomial_system, triangularize, CertifiedPoint, DiophantineError, - DiophantineSolution, HomotopyError, HomotopyOpts, RegularChain, Solution, SolutionSet, - SolverError, + diophantine, expr_to_gbpoly, extract_regular_chain_from_basis, gbpoly_to_expr, + main_variable_recursive, solve_numerical, solve_polynomial_system, triangularize, + CertifiedPoint, DiophantineError, DiophantineSolution, HomotopyError, HomotopyOpts, + RegularChain, Solution, SolutionSet, SolverError, }; pub use crate::stablehlo::emit_stablehlo; pub use crate::sum::{ diff --git a/alkahest-core/src/poly/groebner/mod.rs b/alkahest-core/src/poly/groebner/mod.rs index d4347a4b..8d7422fc 100644 --- a/alkahest-core/src/poly/groebner/mod.rs +++ b/alkahest-core/src/poly/groebner/mod.rs @@ -103,6 +103,14 @@ impl GroebnerBasis { &self.generators } + /// The monomial order the generators were reduced under. + /// + /// Needed to interpret leading terms of [`Self::generators`] and to build a + /// compatible polynomial before calling [`Self::reduce`]. + pub fn order(&self) -> MonomialOrder { + self.order + } + /// Reduce a polynomial by this basis. Returns the remainder. pub fn reduce(&self, p: &GbPoly) -> GbPoly { reduce(p, &self.generators, self.order) diff --git a/alkahest-core/src/poly/groebner/monomial_order.rs b/alkahest-core/src/poly/groebner/monomial_order.rs index 1dfb52bb..a5bffdc9 100644 --- a/alkahest-core/src/poly/groebner/monomial_order.rs +++ b/alkahest-core/src/poly/groebner/monomial_order.rs @@ -74,6 +74,15 @@ impl MonomialOrder { _ => None, } } + + /// The canonical name, round-tripping through [`Self::from_str`]. + pub fn as_str(self) -> &'static str { + match self { + MonomialOrder::Lex => "lex", + MonomialOrder::GrLex => "grlex", + MonomialOrder::GRevLex => "grevlex", + } + } } #[cfg(test)] diff --git a/alkahest-core/src/solver/mod.rs b/alkahest-core/src/solver/mod.rs index f273934d..200015d2 100644 --- a/alkahest-core/src/solver/mod.rs +++ b/alkahest-core/src/solver/mod.rs @@ -277,6 +277,54 @@ fn expr_to_gbpoly_rec( } } +// --------------------------------------------------------------------------- +// GbPoly → Expr conversion +// --------------------------------------------------------------------------- + +/// Rebuild an `Expr` from a [`GbPoly`] — the inverse of [`expr_to_gbpoly`]. +/// +/// `vars` must be the same variable list, in the same order, that produced the +/// polynomial's exponent vectors: exponent slot `i` names `vars[i]`. +/// +/// Returns `None` when `vars` is too short to name every variable the +/// polynomial actually uses; silently mis-naming exponent slots would be worse +/// than refusing. The zero polynomial converts to the integer `0`. +pub fn gbpoly_to_expr(poly: &GbPoly, vars: &[ExprId], pool: &ExprPool) -> Option { + let mut terms: Vec = Vec::with_capacity(poly.terms.len()); + for (exp, coeff) in &poly.terms { + if *coeff == 0 { + continue; + } + let mut factors: Vec = Vec::new(); + for (i, &e) in exp.iter().enumerate() { + if e == 0 { + continue; + } + let v = *vars.get(i)?; + factors.push(if e == 1 { + v + } else { + pool.pow(v, pool.integer(e)) + }); + } + // Keep an explicit coefficient factor unless it is a bare `1` in front + // of at least one variable. + if factors.is_empty() || *coeff != 1 { + factors.insert(0, rational_to_expr(coeff, pool)); + } + terms.push(if factors.len() == 1 { + factors[0] + } else { + pool.mul(factors) + }); + } + Some(match terms.len() { + 0 => pool.integer(0), + 1 => terms[0], + _ => pool.add(terms), + }) +} + // --------------------------------------------------------------------------- // ExprId builders // --------------------------------------------------------------------------- @@ -878,7 +926,12 @@ fn refine_solutions( /// Free symbols in `equations` that are not among the declared solve `vars`, /// in stable [`ExprId`] order (via [`collect_free_vars`]'s `BTreeSet`). -fn collect_parameters(equations: &[ExprId], vars: &[ExprId], pool: &ExprPool) -> Vec { +/// +/// [`solve_polynomial_system`] appends these after `vars`, and the resulting +/// concatenation is the exponent-vector ordering of any +/// [`SolutionSet::Parametric`] basis it returns — so a caller that wants to +/// read that basis back with [`gbpoly_to_expr`] needs this list. +pub fn collect_parameters(equations: &[ExprId], vars: &[ExprId], pool: &ExprPool) -> Vec { let declared: BTreeSet = vars.iter().copied().collect(); let mut params = BTreeSet::new(); for &eq in equations { @@ -1021,6 +1074,99 @@ mod tests { }) } + /// `expr_to_gbpoly` ∘ `gbpoly_to_expr` is the identity on the canonical + /// side: an `Expr` rebuilt from a polynomial converts back to the same + /// polynomial. + #[test] + fn gbpoly_expr_round_trip() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let vars = vec![x, y]; + + // 3/2·x²y − y + 7 + let expr = pool.add(vec![ + pool.mul(vec![ + pool.rational(3_i32, 2_i32), + pool.pow(x, pool.integer(2_i32)), + y, + ]), + pool.mul(vec![pool.integer(-1_i32), y]), + pool.integer(7_i32), + ]); + + let p = expr_to_gbpoly(expr, &vars, &pool).unwrap(); + let back = gbpoly_to_expr(&p, &vars, &pool).expect("named every variable"); + let p2 = expr_to_gbpoly(back, &vars, &pool).unwrap(); + + assert_eq!(p.n_vars, p2.n_vars); + assert_eq!(p.terms, p2.terms); + } + + #[test] + fn gbpoly_to_expr_zero_and_constant() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + + let zero = GbPoly::zero(1); + assert_eq!( + gbpoly_to_expr(&zero, &[x], &pool), + Some(pool.integer(0_i32)) + ); + + let five = GbPoly::constant(Rational::from(5), 1); + assert_eq!( + gbpoly_to_expr(&five, &[x], &pool), + Some(pool.integer(5_i32)) + ); + } + + /// A short `vars` list must refuse rather than silently rename exponent + /// slots — a wrong-but-plausible polynomial is the worst outcome here. + #[test] + fn gbpoly_to_expr_refuses_a_short_variable_list() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + + let p = expr_to_gbpoly(pool.mul(vec![x, y]), &[x, y], &pool).unwrap(); + + assert_eq!(gbpoly_to_expr(&p, &[x], &pool), None); + assert!(gbpoly_to_expr(&p, &[x, y], &pool).is_some()); + } + + /// The Gröbner basis of an ideal must survive being read out as `Expr` and + /// fed back in — otherwise elimination results cannot be reused. + #[test] + fn basis_generators_round_trip_through_expr() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let vars = vec![x, y]; + let neg_one = pool.integer(-1_i32); + + // x² + y² − 1, x − y + let circle = pool.add(vec![ + pool.pow(x, pool.integer(2_i32)), + pool.pow(y, pool.integer(2_i32)), + neg_one, + ]); + let line = pool.add(vec![x, pool.mul(vec![neg_one, y])]); + + let gens = vec![ + expr_to_gbpoly(circle, &vars, &pool).unwrap(), + expr_to_gbpoly(line, &vars, &pool).unwrap(), + ]; + let gb = GroebnerBasis::compute_lex(gens); + assert_eq!(gb.order(), MonomialOrder::Lex); + + for g in gb.generators() { + let e = gbpoly_to_expr(g, &vars, &pool).expect("named every variable"); + let reparsed = expr_to_gbpoly(e, &vars, &pool).unwrap(); + assert!(gb.contains(&reparsed), "generator left the ideal"); + } + } + #[test] fn linear_system() { // x + y - 1 = 0, x - y = 0 → x = 1/2, y = 1/2 diff --git a/alkahest-core/src/validated/bounds.rs b/alkahest-core/src/validated/bounds.rs index c372601b..9a67e9a9 100644 --- a/alkahest-core/src/validated/bounds.rs +++ b/alkahest-core/src/validated/bounds.rs @@ -16,6 +16,11 @@ //! value theorem argument on a box already proven to be free of poles and //! branch cuts). `Undecided` means the enclosure was not informative //! enough to prove either case — it is never conflated with `False`. +//! * [`verified_integral`] integrates the **continuous extension** of the +//! integrand across a *removable* singularity `N(x)/D(x)` with +//! `N(p) = D(p) = 0` and `D'` non-vanishing — the enclosure there comes +//! from Cauchy's mean value theorem, not from ignoring the singular point. +//! Genuine (non-removable) singularities are still refused. //! * A sub-box that cannot be bounded rigorously (branch cut, pole, or //! domain violation persisting after the box has been bisected far below //! the scale of the original box) causes the whole call to **refuse** @@ -25,9 +30,13 @@ use super::taylor::{taylor_range, TaylorContext, MAX_ORDER}; use super::{contains_zero, from_bounds, from_float, is_finite, lb, ub, width, ValidatedError}; use crate::ball::ArbBall; -use crate::kernel::{ExprId, ExprPool}; +use crate::diff::diff; +use crate::kernel::subs::subs; +use crate::kernel::{ExprData, ExprId, ExprPool}; +use crate::simplify::simplify; use rug::float::Round; use rug::Float; +use std::collections::{HashMap, VecDeque}; type Result = std::result::Result; @@ -116,6 +125,16 @@ fn max_dim_width(boxes: &[FBox], prec: u32) -> f64 { .fold(0.0_f64, f64::max) } +/// Midpoint of `[lo, hi]` at working precision. +/// +/// Round-to-nearest of a value that lies between two `prec`-representable +/// numbers cannot leave `[lo, hi]`, so the result is always a point *of the +/// interval* — which is what every use here relies on (it is evaluated as a +/// witness point, or used as a split boundary). +fn midpoint(lo: &Float, hi: &Float, prec: u32) -> Float { + Float::with_val(prec, Float::with_val(prec, lo + hi) / 2u32) +} + /// Split the widest dimension of `boxes` at its midpoint into two boxes that /// exactly cover the original (they share the midpoint as a boundary, so /// their union is the original box with no gap). @@ -130,7 +149,7 @@ fn split_widest(boxes: &[FBox], prec: u32) -> (Vec, Vec) { } } let (v, lo, hi) = &boxes[best]; - let mid = Float::with_val(prec, Float::with_val(prec, lo + hi) / 2u32); + let mid = midpoint(lo, hi, prec); let mut b1 = boxes.to_vec(); let mut b2 = boxes.to_vec(); b1[best] = (*v, lo.clone(), mid.clone()); @@ -197,7 +216,28 @@ pub fn bound_on_box( boxes: &[(ExprId, f64, f64)], opts: &BoundOptions, ) -> Result { - if boxes.is_empty() { + let prec = opts.prec; + let boxes0: Vec = boxes + .iter() + .map(|(v, lo, hi)| (*v, Float::with_val(prec, lo), Float::with_val(prec, hi))) + .collect(); + bound_on_fboxes(expr, pool, &boxes0, opts) +} + +/// [`bound_on_box`] over a box whose endpoints are already exact working-precision +/// `Float`s. +/// +/// The public entry point rounds `f64` endpoints into `Float`s; sub-interval +/// endpoints produced by repeated bisection need more than 53 bits, so the +/// internal callers use this form rather than round-tripping through `f64` +/// (which would silently bound a *different* box). +fn bound_on_fboxes( + expr: ExprId, + pool: &ExprPool, + boxes0: &[FBox], + opts: &BoundOptions, +) -> Result { + if boxes0.is_empty() { return Err(ValidatedError::InvalidInput { what: "the box must constrain at least one variable".into(), }); @@ -208,15 +248,11 @@ pub fn bound_on_box( }); } let prec = opts.prec; - let boxes0: Vec = boxes - .iter() - .map(|(v, lo, hi)| (*v, Float::with_val(prec, lo), Float::with_val(prec, hi))) - .collect(); let (lo_bound, used_lo, exhausted_lo) = - extremum_search(expr, pool, &boxes0, opts, Extremum::Min)?; + extremum_search(expr, pool, boxes0, opts, Extremum::Min)?; let (hi_bound, used_hi, exhausted_hi) = - extremum_search(expr, pool, &boxes0, opts, Extremum::Max)?; + extremum_search(expr, pool, boxes0, opts, Extremum::Max)?; let enclosure = from_bounds(&lo_bound, &hi_bound, prec); if !is_finite(&enclosure) { @@ -510,6 +546,305 @@ fn local_integral( Ok(piece) } +// --------------------------------------------------------------------------- +// Removable singularities +// --------------------------------------------------------------------------- + +/// Split a product into `(numerator, denominator)` at its negative integer +/// powers, so that `expr = numerator / denominator` as real functions wherever +/// the denominator does not vanish. +/// +/// Returns `None` when there is no negative power to divide by — such an +/// expression cannot have a `0/0` removable singularity in this shape. +fn split_quotient(expr: ExprId, pool: &ExprPool) -> Option<(ExprId, ExprId)> { + let factors: Vec = match pool.get(expr) { + ExprData::Mul(args) => args.to_vec(), + ExprData::Pow { .. } => vec![expr], + _ => return None, + }; + let mut num: Vec = Vec::new(); + let mut den: Vec = Vec::new(); + for f in factors { + match pool.get(f) { + ExprData::Pow { base, exp } => match pool.get(exp) { + ExprData::Integer(n) if n.0 < 0 => { + let k: rug::Integer = -n.0.clone(); + den.push(pool.pow(base, pool.integer(k))); + } + _ => num.push(f), + }, + _ => num.push(f), + } + } + if den.is_empty() { + return None; + } + let n = match num.len() { + 0 => pool.integer(1_i32), + 1 => num[0], + _ => pool.mul(num), + }; + let d = if den.len() == 1 { + den[0] + } else { + pool.mul(den) + }; + Some((n, d)) +} + +/// `expr` with `var` replaced by the *exact* rational `at`, simplified. +/// +/// Returns `true` only when the result is the literal integer `0`. A numeric +/// enclosure can never prove a value is exactly zero, so the removable- +/// 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) + }; + let mut mapping = HashMap::new(); + mapping.insert(var, point); + let substituted = subs(expr, &mapping, pool); + let reduced = simplify(substituted, pool).value; + match pool.get(reduced) { + ExprData::Integer(n) => n.0.cmp0() == std::cmp::Ordering::Equal, + ExprData::Rational(r) => r.0.cmp0() == std::cmp::Ordering::Equal, + _ => false, + } +} + +/// Whether the rigorous enclosure of `expr` at the single point `at` contains +/// zero. +/// +/// This is *not* a proof that the value is zero — no numeric enclosure can be +/// one — and it is not used as such. It is a cross-check on +/// [`vanishes_exactly`]: `simplify` reporting an exact zero that outward-rounded +/// ball arithmetic contradicts would mean a bug in the simplifier, and the safe +/// response is to decline the removable branch rather than certify anything on +/// top of it. +fn enclosure_admits_zero( + expr: ExprId, + pool: &ExprPool, + var: ExprId, + at: &Float, + order: usize, + prec: u32, +) -> bool { + let point = vec![(var, at.clone(), at.clone())]; + match taylor_range(expr, pool, &point, order, prec) { + Ok(r) => contains_zero(&r), + Err(_) => false, + } +} + +/// Midpoint of the enclosure of `expr` at the point `at`, or `None` if it +/// could not be evaluated. +fn point_value( + expr: ExprId, + pool: &ExprPool, + var: ExprId, + at: &Float, + order: usize, + prec: u32, +) -> Option { + let point = vec![(var, at.clone(), at.clone())]; + let r = taylor_range(expr, pool, &point, order, prec).ok()?; + r.mid.is_finite().then(|| r.mid.clone()) +} + +/// The ingredients of the L'Hôpital enclosure for an integrand written as a +/// quotient: `N`, `D`, `D'`, and the mean-value quotient `N'/D'`. +struct RemovableQuotient { + num: ExprId, + den: ExprId, + dden: ExprId, + /// `N' · (D')⁻¹`. + ratio: ExprId, +} + +impl RemovableQuotient { + /// Structural analysis of the integrand, done once per + /// [`verified_integral`] call. This is only *shape* recognition — nothing + /// here asserts that a removable singularity exists. + fn detect(expr: ExprId, pool: &ExprPool, var: ExprId) -> Option { + let (num, den) = split_quotient(expr, pool)?; + let dnum = diff(num, var, pool).ok()?.value; + let dden = diff(den, var, pool).ok()?.value; + let ratio = pool.mul(vec![dnum, pool.pow(dden, pool.integer(-1_i32))]); + Some(RemovableQuotient { + num, + den, + dden, + ratio, + }) + } + + /// Newton iterates of `D` inside `[lo, hi]`, as *candidate* locations for + /// the singular point. + /// + /// This is pure heuristic search and carries no soundness weight: whatever + /// points come out are still put through the exact symbolic vanishing test + /// before anything is certified. Its only job is to name the singular point + /// exactly when it is a working-precision number that the dyadic bisection + /// grid would never land on — for a denominator like `x - 1/4` a single + /// Newton step from anywhere produces exactly `1/4`. + fn newton_candidates( + &self, + pool: &ExprPool, + var: ExprId, + lo: &Float, + hi: &Float, + opts: &IntegralOptions, + ) -> Vec { + const STEPS: usize = 6; + let (order, prec) = (opts.order, opts.prec); + let mut out = Vec::new(); + let mut z = midpoint(lo, hi, prec); + for _ in 0..STEPS { + let Some(f) = point_value(self.den, pool, var, &z, order, prec) else { + break; + }; + let Some(df) = point_value(self.dden, pool, var, &z, order, prec) else { + break; + }; + if df.is_zero() || !df.is_finite() { + break; + } + let next = Float::with_val(prec, &z - Float::with_val(prec, &f / &df)); + if !next.is_finite() || next < *lo || next > *hi || next == z { + break; + } + z = next; + out.push(z.clone()); + } + out + } + + /// Rigorous enclosure of `∫_lo^hi N/D dx` when `N/D` has a **removable** + /// singularity somewhere in `[lo, hi]`; `None` when that cannot be + /// established, in which case the caller must keep refusing. + /// + /// # Why this is an enclosure + /// + /// Write `J = [lo, hi]` and let `p ∈ J` be a point at which `N` and `D` + /// both vanish *exactly* (checked symbolically, not numerically). The + /// checks below establish, in order: + /// + /// 1. `N` and `D` are analytic on the whole of `J` — a successful + /// [`bound_on_fboxes`] covers `J` by sub-boxes on each of which the + /// Taylor model was built without a domain violation, and every rule in + /// [`super::taylor`] refuses unless its argument stays strictly inside + /// the analytic interior of the primitive's domain. Analytic implies + /// differentiable, so the symbolic derivatives `N'`, `D'` really are the + /// derivatives of `N`, `D` on `J`. + /// 2. `D'` has no zero on `J` (its enclosure excludes zero). Hence `D` is + /// strictly monotone on `J`, so `p` is its *only* zero there and + /// `D(x) ≠ 0` for every other `x ∈ J`. + /// 3. `R` is a rigorous enclosure of the range of `N'/D'` over `J`. + /// + /// Cauchy's mean value theorem then gives, for every `x ∈ J \ {p}`, some + /// `ξ` strictly between `p` and `x` with + /// `(N(x) − N(p))·D'(ξ) = (D(x) − D(p))·N'(ξ)`; since `N(p) = D(p) = 0`, + /// `D(x) ≠ 0` and `D'(ξ) ≠ 0`, this is `N(x)/D(x) = N'(ξ)/D'(ξ) ∈ R`. + /// So the integrand is bounded by `R` on `J` minus a single point, and + /// `∫_J N/D dx ∈ (hi − lo)·R`. + /// + /// Note what is being integrated: the integrand is *undefined* at `p`, and + /// the value returned is the integral of its continuous extension (which + /// is the same number for every extension, `{p}` being a null set). + fn piece( + &self, + pool: &ExprPool, + var: ExprId, + lo: &Float, + hi: &Float, + opts: &IntegralOptions, + ) -> Option { + let prec = opts.prec; + if lo >= hi { + return None; + } + let bopts = BoundOptions { + order: opts.order, + prec, + tol: opts.tol, + max_subdivisions: opts.max_subdivisions, + }; + // The singular point has to be *named* exactly, since only a symbolic + // test can prove a value is zero. Candidates are the two endpoints, the + // midpoint, and whatever Newton's method on the denominator turns up — + // the last of those is what reaches a singularity sitting at a point + // the dyadic bisection grid never visits. A singularity at a point that + // no candidate names exactly is simply refused. + let mut candidates = vec![lo.clone(), hi.clone(), midpoint(lo, hi, prec)]; + candidates.extend(self.newton_candidates(pool, var, lo, hi, opts)); + candidates.iter().find(|p| { + vanishes_exactly(self.den, pool, var, p) + && vanishes_exactly(self.num, pool, var, p) + && enclosure_admits_zero(self.den, pool, var, p, opts.order, prec) + && enclosure_admits_zero(self.num, pool, var, p, opts.order, prec) + })?; + + let j = vec![(var, lo.clone(), hi.clone())]; + // (1) N and D analytic on J. + bound_on_fboxes(self.num, pool, &j, &bopts).ok()?; + bound_on_fboxes(self.den, pool, &j, &bopts).ok()?; + // (2) D' non-vanishing on J. + let dd = bound_on_fboxes(self.dden, pool, &j, &bopts).ok()?; + if contains_zero(dd.enclosure()) { + return None; + } + // (3) the mean-value quotient. + let r = bound_on_fboxes(self.ratio, pool, &j, &bopts).ok()?; + + let w = Float::with_val(prec, hi - lo); + let piece = from_float(&w, prec) * r.enclosure().clone(); + if !is_finite(&piece) { + return None; + } + Some(piece) + } +} + +/// Turn a refusal that survived bisection into a message that says *what* is +/// singular and *where*, and that distinguishes "no rigorous enclosure of the +/// integrand exists here" from "the integral does not exist". +fn describe_singularity( + cause: ValidatedError, + lo: &Float, + hi: &Float, + a: &Float, + b: &Float, +) -> ValidatedError { + let ValidatedError::DomainViolation { what } = &cause else { + return cause; + }; + let where_ = if lo <= a { + "the left endpoint" + } else if hi >= b { + "the right endpoint" + } else { + "an interior point" + }; + let at = midpoint(lo, hi, lo.prec()).to_f64(); + ValidatedError::DomainViolation { + what: format!( + "the integrand is singular at {where_} x ≈ {at:e} ({what}). \ + This is a statement about the *integrand*, not about the integral: \ + an integrable singularity still has a finite integral, which this \ + routine cannot certify. Removable singularities written as N(x)/D(x) \ + with N(p) = D(p) = 0 exactly and D'(p) ≠ 0 are handled automatically \ + (their continuous extension is integrated); this one is not of that form" + ), + } +} + /// Rigorous enclosure of `∫_a^b f(x) dx`, via adaptive Taylor-model /// quadrature. /// @@ -522,16 +857,32 @@ fn local_integral( /// (interval addition is exact for disjoint, adjacent pieces), so the total /// is always a sound outer bound of the true integral. /// +/// # Removable singularities +/// +/// An integrand written as `N(x)/D(x)` with `N(p) = D(p) = 0` — `ln(1+x)/x` on +/// `[0, 1]`, `sin(x)/x` on `[-1, 1]` — has no rigorous Taylor model at `p` +/// (the reciprocal's enclosure contains zero) even though nothing about the +/// integral is singular. Such a sub-interval is enclosed through Cauchy's mean +/// value theorem instead: `N(x)/D(x) = N'(ξ)/D'(ξ)` for some `ξ` in the +/// sub-interval, so an enclosure `R` of `N'/D'` there — which is perfectly +/// regular — gives `∫_J N/D dx ∈ |J| · R`. The vanishing of `N` and `D` +/// at `p` is checked *symbolically* — a numeric enclosure cannot prove a value +/// is exactly zero — and `D'` must be certified non-vanishing on the +/// sub-interval, so a genuine pole is never mistaken for a removable one. The +/// number returned is the integral of the continuous extension. +/// /// Refuses (does not guess) when: /// - `a` or `b` is non-finite (infinite-limit improper integrals are not /// supported — there is no box to Taylor-expand over), /// - `a > b`, -/// - the integrand has a genuine singularity in `[a, b]` (e.g. `1/sqrt(x)` -/// on `[0, 1]`) — subdivision is tried first in case the domain violation -/// is only a boundary artefact of a coarse box, but a true interior -/// singularity persists down to the bisection floor and the call refuses -/// with the underlying [`ValidatedError`] rather than silently skipping -/// the offending piece. +/// - the integrand has a singularity in `[a, b]` that is not removable in the +/// above sense (e.g. `1/sqrt(x)` on `[0, 1]`, or the *integrable* endpoint +/// singularity of `-log(x)` on `[0, 1]`) — subdivision is tried first in +/// case the domain violation is only a boundary artefact of a coarse box, +/// but a persistent one refuses with a [`ValidatedError`] that names the +/// location and distinguishes "the integrand is singular here" from "the +/// integral does not exist", rather than silently skipping the offending +/// piece. pub fn verified_integral( expr: ExprId, pool: &ExprPool, @@ -570,19 +921,40 @@ pub fn verified_integral( let floor = total_width.to_f64_round(Round::Up) * 2f64.powi(-SINGULARITY_BISECTION_LIMIT); let tol_total = Float::with_val(prec, opts.tol); - let mut stack: Vec<(Float, Float)> = vec![(a_f, b_f)]; + let mut stack: Vec<(Float, Float)> = vec![(a_f.clone(), b_f.clone())]; let mut total: Option = None; let mut subdivisions = 0usize; let mut exhausted = false; + // Structural `N/D` analysis of the integrand, built at most once and only + // if a sub-interval actually refuses. + let mut removable: Option> = None; while let Some((lo, hi)) = stack.pop() { - match local_integral(expr, pool, var, &lo, &hi, opts.order, prec) { + let piece_w = Float::with_val(prec, &hi - &lo); + let piece_tol = Float::with_val( + prec, + &tol_total * Float::with_val(prec, &piece_w / &total_width), + ); + let outcome = match local_integral(expr, pool, var, &lo, &hi, opts.order, prec) { + Ok(piece) => Ok(piece), + Err(e) if is_recoverable_domain_issue(&e) => { + // The Taylor model refuses here. Before bisecting (and before + // eventually giving up), see whether the refusal is only the + // `0/0` shape of a removable singularity, which has a rigorous + // enclosure of its own. + let q = removable + .get_or_insert_with(|| RemovableQuotient::detect(expr, pool, var)) + .as_ref(); + match q.and_then(|q| q.piece(pool, var, &lo, &hi, opts)) { + Some(piece) => Ok(piece), + None => Err(e), + } + } + Err(e) => return Err(e), + }; + + match outcome { Ok(piece) => { - let piece_w = Float::with_val(prec, &hi - &lo); - let piece_tol = Float::with_val( - prec, - &tol_total * Float::with_val(prec, &piece_w / &total_width), - ); let w = width(&piece); if w <= piece_tol || subdivisions >= opts.max_subdivisions { if w > piece_tol { @@ -594,22 +966,21 @@ pub fn verified_integral( }); } else { subdivisions += 1; - let mid = Float::with_val(prec, Float::with_val(prec, &lo + &hi) / 2u32); + let mid = midpoint(&lo, &hi, prec); stack.push((mid.clone(), hi)); stack.push((lo, mid)); } } - Err(e) if is_recoverable_domain_issue(&e) => { - let piece_w = Float::with_val(prec, &hi - &lo).to_f64_round(Round::Up); - if subdivisions >= opts.max_subdivisions || piece_w <= floor { - return Err(e); + Err(e) => { + let w = piece_w.to_f64_round(Round::Up); + if subdivisions >= opts.max_subdivisions || w <= floor { + return Err(describe_singularity(e, &lo, &hi, &a_f, &b_f)); } subdivisions += 1; - let mid = Float::with_val(prec, Float::with_val(prec, &lo + &hi) / 2u32); + let mid = midpoint(&lo, &hi, prec); stack.push((mid.clone(), hi)); stack.push((lo, mid)); } - Err(e) => return Err(e), } } @@ -661,19 +1032,212 @@ fn determined_sign(b: &ArbBall) -> Option { } } +/// Proven sign witnesses collected anywhere in one box. +/// +/// Only *proven* signs are recorded: a value whose enclosure straddles zero +/// contributes nothing. Once both a positive and a negative witness exist the +/// intermediate value theorem applies along the segment joining them, which +/// stays inside the (convex) box. +#[derive(Default)] +struct SignWitnesses { + positive: bool, + negative: bool, +} + +impl SignWitnesses { + fn record(&mut self, s: Option) { + match s { + Some(true) => self.positive = true, + Some(false) => self.negative = true, + None => {} + } + } + + fn both(&self) -> bool { + self.positive && self.negative + } +} + +/// Rigorously evaluate `expr` at one point of the box and report its sign, or +/// `None` when the enclosure straddles zero or the evaluation refused. +/// +/// 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( + 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, + } +} + +/// Degenerate boxes for the centre, the per-axis endpoints and (in low +/// dimension) the corners of `boxes0` — the cheap witnesses that settle most +/// boxes before any subdivision happens. +fn seed_points(boxes0: &[FBox], prec: u32) -> Vec> { + let n = boxes0.len(); + let centre: Vec = boxes0 + .iter() + .map(|(_, lo, hi)| midpoint(lo, hi, prec)) + .collect(); + let mut points: Vec> = vec![centre.clone()]; + for (i, (_, lo, hi)) in boxes0.iter().enumerate() { + for end in [lo, hi] { + let mut p = centre.clone(); + p[i] = end.clone(); + points.push(p); + } + } + // 2ⁿ corners is only affordable in low dimension. + if n <= 3 { + for mask in 0..(1usize << n) { + points.push( + boxes0 + .iter() + .enumerate() + .map(|(i, (_, lo, hi))| { + if mask & (1 << i) == 0 { + lo.clone() + } else { + hi.clone() + } + }) + .collect(), + ); + } + } + points + .into_iter() + .map(|p| { + boxes0 + .iter() + .zip(p) + .map(|((v, _, _), c)| (*v, c.clone(), c)) + .collect() + }) + .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. +/// +/// # Why this 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 +/// search only returns `Ok` once every part of the box lies in some sub-box on +/// which a Taylor model was built, and [`super::taylor`] refuses to build one +/// unless every elementary step stayed strictly inside the analytic interior +/// of its domain (`log`/`sqrt` of a strictly positive enclosure, reciprocals of +/// an enclosure bounded away from zero, `tan` where `cos` cannot vanish, …). +/// A box on which a `Refine` outcome is never resolved keeps a `-∞` key, is +/// never pruned and never satisfies the convergence test, so it can only end in +/// a refusal — `Ok` really does mean "analytic on a finite closed cover of the +/// box", hence continuous on the box. +/// +/// The box is a product of intervals and therefore convex, so for any two of +/// its points `p`, `q` the segment `[p, q]` stays inside it. If `f(p) > 0` and +/// `f(q) < 0` — each proven by an outward-rounded enclosure at a degenerate box +/// — then `t ↦ f(p + t(q − p))` is continuous on `[0, 1]` and changes sign, so +/// it vanishes somewhere: a root of `f` inside the box. This is exactly the +/// intermediate value theorem, and it needs no relationship between `p` and `q` +/// beyond both lying in the box — which is what lets an even number of roots +/// (`x² − 2` on `[-2, 2]`, whose *endpoints* are both positive) be settled. +/// +/// The search itself is the same Moore–Skelboe subdivision used elsewhere: +/// sub-boxes whose enclosure has a determined sign are recorded as witnesses +/// and dropped (they cannot contain a root); the rest are sampled at their +/// centre and split. It is a *search*, so failing to find a witness pair proves +/// nothing — the caller must answer `Undecided`, never `True`. +fn root_exists_witness( + expr: ExprId, + pool: &ExprPool, + boxes0: &[FBox], + opts: &BoundOptions, +) -> bool { + let (prec, order) = (opts.prec, opts.order); + let floor = max_dim_width(boxes0, prec) * 2f64.powi(-SINGULARITY_BISECTION_LIMIT); + 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; + } + } + + let mut queue: VecDeque> = VecDeque::new(); + queue.push_back(boxes0.to_vec()); + let mut budget = opts.max_subdivisions; + + while let Some(b) = queue.pop_front() { + if budget == 0 { + break; + } + budget -= 1; + + // A sub-box whose enclosure has a determined sign is non-empty, so it + // contains a point of that sign — a witness — and it cannot contain a + // root, so there is nothing left to look for inside it. + if let Ok(r) = taylor_range(expr, pool, &b, order, prec) { + if let Some(s) = determined_sign(&r) { + w.record(Some(s)); + if w.both() { + return true; + } + continue; + } + } + + let centre: Vec = b + .iter() + .map(|(v, lo, hi)| { + let m = midpoint(lo, hi, prec); + (*v, m.clone(), m) + }) + .collect(); + w.record(point_sign(expr, pool, ¢re, order, prec)); + if w.both() { + return true; + } + + if max_dim_width(&b, prec) <= floor { + continue; + } + let (b1, b2) = split_widest(&b, prec); + queue.push_back(b1); + queue.push_back(b2); + } + + false +} + /// Verified check for the absence of roots of `expr` on `boxes`. /// /// Returns: /// - [`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`] in the univariate case, when the box is proven free -/// of poles/branch cuts (the full-box enclosure succeeded) *and* the -/// rigorously evaluated endpoint values have determined, opposite signs — -/// a root is certified to exist by the intermediate value theorem. -/// - [`Verdict::Undecided`] otherwise: the enclosure straddles zero and no -/// sign-change proof was available (multivariate, or same-signed / -/// indeterminate endpoints). +/// - [`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. /// /// Propagates a [`ValidatedError`] (refuses) exactly when /// [`bound_on_box`] would: unsupported primitives, unbound symbols, or a @@ -689,26 +1253,16 @@ pub fn verified_no_roots( return Ok(Verdict::True); } - // Root-existence (Verdict::False) via IVT: only sound in 1-D, and only - // once we already know (from the successful full-box call above) that - // `expr` is finite and pole/branch-cut free across the whole box, so it - // is continuous on it. - if boxes.len() == 1 { - let (v, lo, hi) = boxes[0]; - if lo < hi { - let flo = bound_on_box(expr, pool, &[(v, lo, lo)], opts); - let fhi = bound_on_box(expr, pool, &[(v, hi, hi)], opts); - if let (Ok(flo), Ok(fhi)) = (flo, fhi) { - if let (Some(slo), Some(shi)) = ( - determined_sign(flo.enclosure()), - determined_sign(fhi.enclosure()), - ) { - if slo != shi { - return Ok(Verdict::False); - } - } - } - } + // The call above is the continuity certificate the IVT argument needs; see + // `root_exists_witness` for why, and for why a witness pair anywhere in the + // box (not just at its endpoints) is enough. + let prec = opts.prec; + let boxes0: Vec = boxes + .iter() + .map(|(v, lo, hi)| (*v, Float::with_val(prec, lo), Float::with_val(prec, hi))) + .collect(); + if root_exists_witness(expr, pool, &boxes0, opts) { + return Ok(Verdict::False); } Ok(Verdict::Undecided) @@ -1059,6 +1613,152 @@ mod tests { )); } + // ── verified_integral: removable singularities ───────────────────── + + /// Assert that `r` really brackets `exact` and is no wider than `tol`. + /// An enclosure that is merely *returned* is worth nothing. + fn assert_brackets(r: &IntegralResult, exact: f64, tol: f64, label: &str) { + assert!( + r.lower() <= exact && exact <= r.upper(), + "{label}: {exact} not in [{}, {}]", + r.lower(), + r.upper() + ); + assert!( + r.upper() - r.lower() <= tol, + "{label}: enclosure [{}, {}] is wider than {tol}", + r.lower(), + r.upper() + ); + } + + #[test] + fn integral_of_log1p_over_x_is_pi_squared_over_twelve() { + // ∫₀¹ ln(1+x)/x dx = π²/12. The reciprocal's enclosure contains zero + // at x = 0, but the singularity is removable, so the continuous + // extension is integrated instead of refusing. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = div( + &pool, + pool.func("log", vec![pool.add(vec![pool.integer(1_i32), x])]), + x, + ); + let r = verified_integral(e, &pool, x, 0.0, 1.0, &iopts()).unwrap(); + let exact = std::f64::consts::PI * std::f64::consts::PI / 12.0; + assert_brackets(&r, exact, 1e-6, "ln(1+x)/x"); + } + + #[test] + fn integral_of_sinc_is_twice_si_of_one() { + // ∫_{-1}^{1} sin(x)/x dx = 2·Si(1) = 1.892166140734366…, and the + // removable point x = 0 is *interior*: it only becomes a sub-interval + // endpoint after the first bisection. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = div(&pool, pool.func("sin", vec![x]), x); + let r = verified_integral(e, &pool, x, -1.0, 1.0, &iopts()).unwrap(); + assert_brackets(&r, 1.892_166_140_734_366_4, 1e-6, "sin(x)/x"); + } + + #[test] + fn integral_of_expm1_over_x() { + // ∫₀¹ (exp(x)-1)/x dx = Σ_{n≥1} 1/(n·n!) = 1.317902151454403… + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let num = sub(&pool, pool.func("exp", vec![x]), pool.integer(1_i32)); + let e = div(&pool, num, x); + let r = verified_integral(e, &pool, x, 0.0, 1.0, &iopts()).unwrap(); + assert_brackets(&r, 1.317_902_151_454_403_9, 1e-6, "(exp(x)-1)/x"); + } + + #[test] + fn removable_extension_agrees_with_the_regular_quadrature_away_from_zero() { + // Cross-check: ∫₀¹ ln(1+x)/x must equal ∫₀^{1/2} + ∫_{1/2}^1, and the + // second piece needs no special handling at all. If the removable + // branch were biased, the two would disagree. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = div( + &pool, + pool.func("log", vec![pool.add(vec![pool.integer(1_i32), x])]), + x, + ); + let whole = verified_integral(e, &pool, x, 0.0, 1.0, &iopts()).unwrap(); + let left = verified_integral(e, &pool, x, 0.0, 0.5, &iopts()).unwrap(); + let right = verified_integral(e, &pool, x, 0.5, 1.0, &iopts()).unwrap(); + assert!( + whole.lower() <= left.upper() + right.upper() + && whole.upper() >= left.lower() + right.lower(), + "split disagrees: whole [{}, {}] vs pieces [{}, {}] + [{}, {}]", + whole.lower(), + whole.upper(), + left.lower(), + left.upper(), + right.lower(), + right.upper() + ); + } + + #[test] + fn genuine_pole_is_not_mistaken_for_a_removable_one() { + // 1/x on [-1,1]: the numerator does not vanish, so the removable + // branch must decline and the call must still refuse. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = div(&pool, pool.integer(1_i32), x); + let err = verified_integral(e, &pool, x, -1.0, 1.0, &iopts()).unwrap_err(); + assert_eq!(crate::errors::AlkahestError::code(&err), "E-VALIDATED-003"); + } + + #[test] + fn double_pole_with_a_simple_numerator_zero_is_refused() { + // sin(x)/x² ~ 1/x near 0: the numerator vanishes to order 1 but the + // denominator to order 2, so the singularity is *not* removable and + // the integral does not converge. D' = 2x vanishes at 0, which is + // exactly the check that stops the L'Hôpital argument here. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = div( + &pool, + pool.func("sin", vec![x]), + pool.pow(x, pool.integer(2_i32)), + ); + let err = verified_integral(e, &pool, x, -1.0, 1.0, &iopts()).unwrap_err(); + assert_eq!(crate::errors::AlkahestError::code(&err), "E-VALIDATED-003"); + } + + #[test] + fn integrable_endpoint_singularity_is_refused_with_a_message_that_says_so() { + // ∫₀¹ -log x dx = 1 exists, but no rigorous enclosure of the + // *integrand* does. Option (c) of the issue: the refusal must not read + // as "the integral does not exist". + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = pool.mul(vec![pool.integer(-1_i32), pool.func("log", vec![x])]); + let err = verified_integral(e, &pool, x, 0.0, 1.0, &iopts()).unwrap_err(); + assert_eq!(crate::errors::AlkahestError::code(&err), "E-VALIDATED-003"); + let msg = err.to_string(); + assert!(msg.contains("the left endpoint"), "{msg}"); + assert!(msg.contains("integrable singularity"), "{msg}"); + } + + #[test] + fn interior_singularity_is_reported_as_interior() { + // 1/(x - 1/2) on [0,1]: the refusal should localise the pole away from + // both endpoints. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let e = div( + &pool, + pool.integer(1_i32), + sub(&pool, x, pool.rational(1_i32, 2_i32)), + ); + let err = verified_integral(e, &pool, x, 0.0, 1.0, &iopts()).unwrap_err(); + let msg = err.to_string(); + assert!(msg.contains("an interior point"), "{msg}"); + } + #[test] fn integral_refuses_on_infinite_bounds() { let pool = ExprPool::new(); @@ -1101,14 +1801,107 @@ mod tests { } #[test] - fn no_roots_undecided_when_enclosure_straddles_zero_without_sign_change() { + fn no_roots_false_across_a_convex_multivariate_box() { let pool = ExprPool::new(); let x = pool.symbol("x", Domain::Real); let y = pool.symbol("y", Domain::Real); - // x - y over a 2-D box that straddles zero: no 1-D IVT argument - // applies, so this must not be misreported as True or False. + // x - y on [-1,1]²: it is +2 at (1,-1) and -2 at (-1,1), and the + // segment between them stays in the box, so a root is certified. This + // is the honest verdict — the whole plane x = y of roots is inside the + // box — and it is why the witness search does not need the two points + // to be endpoints of a 1-D interval. let e = sub(&pool, x, y); let v = verified_no_roots(e, &pool, &[(x, -1.0, 1.0), (y, -1.0, 1.0)], &opts()).unwrap(); + assert_eq!(v, Verdict::False); + } + + #[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 + // 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. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let half = pool.rational(1_i32, 2_i32); + let dx = sub(&pool, x, half); + let dy = sub(&pool, y, half); + 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); + } + + /// `x² - 2` on every box from the issue-13 table, plus the product that + /// hides the same two roots behind a strictly positive factor. Before the + /// witness search only the boxes with a *top-level* endpoint sign change + /// could be settled; an even number of roots defeated the rest. + #[test] + fn no_roots_false_regardless_of_the_root_count_in_the_box() { + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let x2 = pool.pow(x, pool.integer(2_i32)); + let f = sub(&pool, x2, pool.integer(2_i32)); + // (x²-2)(x²+1) — same two roots, second factor never vanishes. + let g = pool.mul(vec![f, pool.add(vec![x2, pool.integer(1_i32)])]); + + for (expr, lo, hi, label) in [ + (f, -2.0, 0.0, "x^2-2 on [-2,0] (1 root)"), + (f, 0.0, 2.0, "x^2-2 on [0,2] (1 root)"), + (f, 1.3, 1.5, "x^2-2 on [1.3,1.5] (1 root)"), + (f, -2.0, 2.0, "x^2-2 on [-2,2] (2 roots)"), + (f, -10.0, 10.0, "x^2-2 on [-10,10] (2 roots)"), + (g, -2.0, 2.0, "(x^2-2)(x^2+1) on [-2,2] (2 roots)"), + ] { + let v = verified_no_roots(expr, &pool, &[(x, lo, hi)], &opts()).unwrap(); + assert_eq!(v, Verdict::False, "{label}"); + } + } + + #[test] + fn no_roots_stays_true_where_it_was_true() { + // The witness search must never be reached when the enclosure already + // proves absence — these are the `True` rows that have to stay `True`. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let x2 = pool.pow(x, pool.integer(2_i32)); + for (expr, lo, hi) in [ + (pool.add(vec![x2, pool.integer(1_i32)]), -10.0, 10.0), + (pool.add(vec![x2, pool.integer(2_i32)]), -2.0, 2.0), + (pool.func("exp", vec![x]), -5.0, 5.0), + ] { + let v = verified_no_roots(expr, &pool, &[(x, lo, hi)], &opts()).unwrap(); + assert_eq!(v, Verdict::True); + } + } + + #[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 + // 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). + 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]); + let v = verified_no_roots(e, &pool, &[(x, 0.0, 2.0)], &opts()).unwrap(); + assert_eq!(v, Verdict::Undecided); + } + + #[test] + fn no_roots_undecided_for_a_quartic_double_root() { + // (x²-1)² — two double roots at ±1 in [-2,2], again with no sign + // change anywhere. + let pool = ExprPool::new(); + let x = pool.symbol("x", Domain::Real); + let d = sub(&pool, pool.pow(x, pool.integer(2_i32)), pool.integer(1_i32)); + let e = pool.mul(vec![d, d]); + let v = verified_no_roots(e, &pool, &[(x, -2.0, 2.0)], &opts()).unwrap(); assert_eq!(v, Verdict::Undecided); } diff --git a/alkahest-py/src/lib.rs b/alkahest-py/src/lib.rs index af15d58b..73a3ce37 100644 --- a/alkahest-py/src/lib.rs +++ b/alkahest-py/src/lib.rs @@ -3251,10 +3251,14 @@ impl PyUniPoly { self.inner.coefficients_i64() } + /// Degree of the polynomial (`-1` for the zero polynomial). + #[getter] fn degree(&self) -> i64 { self.inner.degree() } + /// True if this is the zero polynomial. + #[getter] fn is_zero(&self) -> bool { self.inner.is_zero() } @@ -3354,10 +3358,14 @@ impl PyMultiPoly { .map_err(conv_error_to_py) } + /// True if this is the zero polynomial. + #[getter] fn is_zero(&self) -> bool { self.inner.is_zero() } + /// Highest total degree over all terms (`0` for the zero polynomial). + #[getter] fn total_degree(&self) -> u32 { self.inner.total_degree() } @@ -3519,6 +3527,8 @@ impl PyRationalFunction { .map_err(conv_error_to_py) } + /// True if the numerator is the zero polynomial. + #[getter] fn is_zero(&self) -> bool { self.inner.is_zero() } @@ -4021,6 +4031,8 @@ impl PyOdeTrajectory { PyList::new_bound(py, rows).into_py(py) } + /// The final time point, or ``None`` for an empty trajectory. + #[getter] fn t_final(&self) -> Option { self.inner.t_final() } @@ -6170,6 +6182,26 @@ impl PyODE { PyODE::new(py, state_vars, rhs, time_var) } + /// Build an explicit first-order system ``d(state_vars)/dt = rhs``. + /// + /// Parameters + /// ---------- + /// state_vars : list[Expr] + /// Symbolic state variables. + /// rhs : list[Expr] + /// Right-hand sides, one per state variable and in the same order. + /// time_var : Expr + /// The independent variable. + /// + /// All three arguments are **positional**; there are no keyword forms. + /// ``ODE(state_vars, rhs, time_var)`` is accepted as a synonym. + /// For an implicit or constrained system use :class:`DAE` instead, and for + /// a scalar higher-order equation use :func:`lower_to_first_order`. + /// + /// Example:: + /// + /// # Harmonic oscillator as a first-order system: x' = v, v' = -x + /// ode = alkahest.ODE.new([x, v], [v, -x], t) #[staticmethod] fn new( py: Python<'_>, @@ -6197,6 +6229,8 @@ impl PyODE { } } + /// Number of state variables (the order of the system). + #[getter] fn order(&self) -> usize { self.inner.order() } @@ -6255,7 +6289,34 @@ impl PyODE { } } -/// `alkahest.lower_to_first_order(var, rhs, order, time_var)` — lower a scalar ODE to first-order. +/// `alkahest.lower_to_first_order(var, rhs, order, time_var)` — lower a scalar +/// higher-order ODE to an equivalent first-order system. +/// +/// Takes the *pieces* of the equation, not an :class:`ODE` object: the scalar +/// unknown, the right-hand side, the order of the derivative on the left, and +/// the independent variable. All four are required and positional. +/// +/// Parameters +/// ---------- +/// var : Expr +/// The scalar unknown, e.g. ``x``. +/// rhs : Expr +/// Right-hand side of ``d^order(var)/dt^order = rhs``. +/// order : int +/// Order of the equation (``2`` for ``x''``). +/// time_var : Expr +/// The independent variable. +/// +/// Returns +/// ------- +/// ODE +/// A first-order system whose state is ``[var, var', …, var^(order-1)]``. +/// +/// Example:: +/// +/// # x'' = -4x → [x' = x', (x')' = -4x] +/// ode = alkahest.lower_to_first_order(x, -4 * x, 2, t) +/// ode.order # 2 #[pyfunction] #[pyo3(name = "lower_to_first_order")] fn py_lower_to_first_order( @@ -6387,6 +6448,22 @@ fn py_adjoint_system( /// The independent time variable. /// /// Use :func:`pantelides` to reduce the differential index before simulation. +/// +/// Attributes +/// ---------- +/// n_equations : int +/// n_variables : int +/// time_var : Expr +/// index : int | None +/// Differentiation index, set by :func:`pantelides` on the DAE it +/// returns; ``None`` on a DAE that has not been index-reduced. +/// +/// Example:: +/// +/// p = alkahest.ExprPool() +/// t, x, dx = p.symbol("t"), p.symbol("x"), p.symbol("dx/dt") +/// dae = alkahest.DAE.new([dx - x], [x], [dx], t) # x' = x +/// dae.equations() # [-x + dx/dt] #[pyclass(name = "DAE")] struct PyDAE { inner: DAE, @@ -6395,6 +6472,32 @@ struct PyDAE { #[pymethods] impl PyDAE { + /// Build a DAE from implicit equations. + /// + /// Parameters + /// ---------- + /// equations : list[Expr] + /// Implicit equations, each meaning ``g = 0``. Write ``x' = f`` as + /// ``dx - f``, using a *separate symbol* for the derivative. + /// variables : list[Expr] + /// Dependent variables (state and algebraic). + /// derivatives : list[Expr] + /// Symbol standing for the time derivative of ``variables[i]``, e.g. + /// ``pool.symbol("dx/dt")``. Alkahest does not parse the name; the + /// positional pairing is what makes it a derivative. + /// time_var : Expr + /// The independent variable. + /// + /// Example:: + /// + /// # Pendulum in Cartesian coordinates (index 3): + /// # x' = u, u' = -lam*x, x**2 + y**2 - 1 = 0 + /// dae = alkahest.DAE.new( + /// [dx - u, du + lam * x, x**2 + y**2 - one], + /// [x, u, lam], + /// [dx, du, dlam], + /// t, + /// ) #[staticmethod] fn new( py: Python<'_>, @@ -6414,20 +6517,108 @@ impl PyDAE { } } + /// Number of equations in the system. + #[getter] fn n_equations(&self) -> usize { self.inner.n_equations() } + + /// Number of dependent variables in the system. + #[getter] fn n_variables(&self) -> usize { self.inner.n_variables() } + /// The independent (time) variable. + #[getter] + fn time_var(&self, py: Python<'_>) -> PyExpr { + PyExpr { + id: self.inner.time_var, + pool: self.pool.clone_ref(py), + } + } + + /// Differentiation index, or ``None`` if it has not been computed. + /// + /// :func:`pantelides` sets it on the DAE it returns: `0` means the input + /// already had a perfect structural matching, `k` means `k` rounds of + /// differentiation were needed. + #[getter] + fn index(&self) -> Option { + self.inner.index + } + + /// The implicit equations, each meaning ``g = 0``. + /// + /// After :func:`pantelides` or :func:`rosenfeld_groebner` this is how you + /// see which equations were added by differentiation. + fn equations(&self, py: Python<'_>) -> Vec { + self.inner + .equations + .iter() + .map(|&id| PyExpr { + id, + pool: self.pool.clone_ref(py), + }) + .collect() + } + + /// The dependent variables, parallel to :meth:`derivatives`. + fn variables(&self, py: Python<'_>) -> Vec { + self.inner + .variables + .iter() + .map(|&id| PyExpr { + id, + pool: self.pool.clone_ref(py), + }) + .collect() + } + + /// The derivative symbols: ``derivatives()[i]`` is ``d(variables()[i])/dt``. + /// + /// Index reduction appends higher jets here (``d2x/dt2``, …), so this is + /// how the extra variables of a prolonged system get names. + fn derivatives(&self, py: Python<'_>) -> Vec { + self.inner + .derivatives + .iter() + .map(|&id| PyExpr { + id, + pool: self.pool.clone_ref(py), + }) + .collect() + } + fn __repr__(&self, py: Python<'_>) -> String { let pool = self.pool.borrow(py); format!("DAE(\n{}\n)", self.inner.display(&pool.inner)) } } -/// `alkahest.pantelides(dae)` — apply the Pantelides index-reduction algorithm. +/// `alkahest.pantelides(dae)` — structural index reduction of a DAE. +/// +/// Repeatedly differentiates the equations that a maximum bipartite matching +/// leaves unmatched, until every equation is matched to a variable, and +/// returns the **reduced** :class:`DAE` — not a report object. +/// +/// The reduced system has more equations and more derivative jets than the +/// input; read them with :meth:`DAE.equations` and :meth:`DAE.derivatives`. +/// :attr:`DAE.index` on the result is the number of differentiation rounds it +/// took (``0`` = the input already matched). +/// +/// Raises `ValueError` (`E-DAE-002`) if the index exceeds 10 — try +/// :func:`dae_index_reduce`, which falls back to :func:`rosenfeld_groebner`. +/// +/// This is a *structural* algorithm: it looks at which variables occur in +/// which equations, not at whether the coefficients actually make the system +/// solvable. +/// +/// Example:: +/// +/// reduced = alkahest.pantelides(dae) +/// reduced.index # differentiation rounds used +/// reduced.equations() # original equations plus the differentiated ones #[pyfunction] #[pyo3(name = "pantelides")] fn py_pantelides(py: Python<'_>, dae: PyRef) -> PyResult { @@ -6523,6 +6714,8 @@ impl PyHybridODE { } } + /// Number of registered discrete events. + #[getter] fn n_events(&self) -> usize { self.inner.events.len() } @@ -6601,11 +6794,13 @@ impl PyComponent { } /// Number of constitutive equations contributed by this component. + #[getter] fn n_equations(&self) -> usize { self.inner.equations.len() } /// Number of external connection ports. + #[getter] fn n_ports(&self) -> usize { self.inner.ports.len() } @@ -7375,6 +7570,8 @@ impl PyArbBall { self.inner.contains(v) } + /// True if the ball has zero radius (a single exact value). + #[getter] fn is_exact(&self) -> bool { self.inner.is_exact() } @@ -9204,6 +9401,18 @@ fn py_bound_on_box( /// unlike a floating-point quadrature, the answer is a *theorem*: the true /// integral is guaranteed to lie in the returned interval. Refuses on /// singular or improper integrands rather than guessing. +/// +/// A **removable** singularity is not a refusal: an integrand written as +/// ``N(x)/D(x)`` with ``N(p) = D(p) = 0`` and ``D'(p) != 0`` — ``log(1+x)/x`` +/// on ``[0, 1]``, ``sin(x)/x`` on ``[-1, 1]`` — is enclosed via Cauchy's mean +/// value theorem, and the value returned is the integral of the continuous +/// extension. The two zeros are checked *symbolically*, so a genuine pole is +/// never mistaken for a removable one. +/// +/// A singularity that is integrable but not removable (``-log(x)`` on +/// ``[0, 1]``, ``1/sqrt(1-x*x)`` on ``[0, 1]``) is still refused: the integral +/// exists, but no rigorous enclosure of the *integrand* does. The +/// :class:`ValidatedError` message says which of the two situations it is. #[allow(clippy::too_many_arguments)] #[pyfunction] #[pyo3(name = "verified_integral", signature = (expr, var, a, b, *, order = 8, prec = 128, tol = 1e-9, max_subdivisions = 4096))] @@ -9244,6 +9453,19 @@ fn py_verified_integral( /// ``"false"`` proves it does, and ``"undecided"`` means neither could be /// established within the budget. The third is never collapsed into the /// other two. +/// +/// ``"false"`` is certified by the intermediate value theorem: the box is +/// first proven free of poles and branch cuts (so `expr` is continuous on it), +/// then the box is subdivided until two points are found at which `expr` is +/// rigorously proven to have opposite signs. A box is convex, so the segment +/// joining them stays inside it and `expr` must vanish somewhere along it. +/// Because the two points need not be the box's own endpoints, an even number +/// of roots no longer defeats the test — ``x*x - 2`` on ``[-2, 2]`` is +/// ``"false"``. +/// +/// A root that never produces a sign change — a double root such as +/// ``(x-1)**2`` — stays ``"undecided"``: no witness exists, and none is +/// invented. #[pyfunction] #[pyo3(name = "verified_no_roots", signature = (expr, r#box, *, order = 6, prec = 128, tol = 1e-9, max_subdivisions = 2048))] fn py_verified_no_roots( @@ -9805,32 +10027,237 @@ fn py_cuda_device_count() -> usize { #[cfg(feature = "groebner")] use alkahest_core::{ - dae_index_reduce, expr_to_gbpoly, primary_decomposition, radical as core_ideal_radical, - rosenfeld_groebner_with_options, DaeIndexReduction, GbPoly, GroebnerBasis, MonomialOrder, + dae_index_reduce_ranked, expr_to_gbpoly, gbpoly_to_expr, primary_decomposition, + radical as core_ideal_radical, rosenfeld_groebner_ranked, DaeIndexReduction, GbPoly, + GroebnerBasis, MonomialOrder, }; +/// A sparse multivariate polynomial over ℚ, as used by the Gröbner machinery. +/// +/// A `GbPoly` stores exponent *vectors*, not variable names, so reading one +/// back needs the variable list its exponent slots refer to. Polynomials +/// handed out by Alkahest carry that list with them (see :meth:`variables`), +/// so :meth:`to_expr` normally takes no arguments; pass `vars` explicitly only +/// for a polynomial built against a different list. +/// +/// Attributes +/// ---------- +/// is_zero : bool +/// n_vars : int +/// n_terms : int #[cfg(feature = "groebner")] #[pyclass(name = "GbPoly")] struct PyGbPoly { inner: GbPoly, + /// Pool that `var_ids` belong to, when the variable context is known. + pool: Option>, + /// Variables in the order used for exponent vectors. + var_ids: Vec, +} + +#[cfg(feature = "groebner")] +impl PyGbPoly { + /// Wrap a core polynomial together with the variable context that names + /// its exponent slots. + fn with_ctx( + py: Python<'_>, + inner: GbPoly, + pool: Option<&Py>, + var_ids: &[ExprId], + ) -> PyGbPoly { + PyGbPoly { + inner, + pool: pool.map(|p| p.clone_ref(py)), + var_ids: var_ids.to_vec(), + } + } +} + +/// Resolve the `(pool, var_ids)` pair for an `Expr` conversion: an explicit +/// `vars` argument wins, otherwise the stored context, otherwise an error. +#[cfg(feature = "groebner")] +fn resolve_gb_ctx( + py: Python<'_>, + stored_pool: Option<&Py>, + stored_vars: &[ExprId], + vars: Option>>, + what: &str, +) -> PyResult<(Py, Vec)> { + match vars { + Some(v) if !v.is_empty() => Ok((v[0].pool.clone_ref(py), v.iter().map(|e| e.id).collect())), + _ => match stored_pool { + Some(p) => Ok((p.clone_ref(py), stored_vars.to_vec())), + None => Err(pyo3::exceptions::PyValueError::new_err(format!( + "{what} carries no variable context; pass vars=[...] naming exponent slots 0, 1, …" + ))), + }, + } } #[cfg(feature = "groebner")] #[pymethods] impl PyGbPoly { + /// True if this is the zero polynomial. + #[getter] fn is_zero(&self) -> bool { self.inner.is_zero() } + /// Number of variables in the ambient ring. + #[getter] fn n_vars(&self) -> usize { self.inner.n_vars } + /// Number of non-zero terms. + #[getter] + fn n_terms(&self) -> usize { + self.inner.terms.len() + } + + /// The variables naming this polynomial's exponent slots, in order. + /// + /// Empty when the polynomial carries no variable context — see + /// :meth:`to_expr`. + fn variables(&self, py: Python<'_>) -> Vec { + match &self.pool { + None => vec![], + Some(pool) => self + .var_ids + .iter() + .map(|&id| PyExpr { + id, + pool: pool.clone_ref(py), + }) + .collect(), + } + } + + /// The terms as ``(exponents, coefficient)`` pairs. + /// + /// `exponents` is a tuple of `int` parallel to :meth:`variables`; + /// `coefficient` is an exact Python `int` or `fractions.Fraction`. + /// Terms come in ascending exponent-vector order. + /// + /// Example:: + /// + /// p = alkahest.expr_to_gbpoly(x**2 - 3*y, [x, y]) + /// p.terms() # [((0, 1), -3), ((2, 0), 1)] + fn terms(&self, py: Python<'_>) -> PyResult> { + let mut out = Vec::with_capacity(self.inner.terms.len()); + for (exp, coeff) in &self.inner.terms { + let exps = pyo3::types::PyTuple::new_bound(py, exp.iter().map(|&e| e as u64)); + out.push((exps.into_py(py), rational_to_py(py, coeff)?)); + } + Ok(out) + } + + /// Convert back to an :class:`Expr`. + /// + /// Parameters + /// ---------- + /// vars : list[Expr], optional + /// Variables naming exponent slots 0, 1, …. Defaults to + /// :meth:`variables`, which is what you want for a polynomial that came + /// out of a :class:`GroebnerBasis`, a :class:`RegularChain` or + /// :func:`rosenfeld_groebner`. + /// + /// Raises + /// ------ + /// ValueError + /// If no variable context is available, or `vars` names fewer variables + /// than the polynomial actually uses. + /// + /// Example:: + /// + /// gb = alkahest.GroebnerBasis.compute([x**2 - y, x - y], [x, y]) + /// [g.to_expr() for g in gb] + #[pyo3(signature = (vars=None))] + fn to_expr(&self, py: Python<'_>, vars: Option>>) -> PyResult { + let (pool_py, var_ids) = + resolve_gb_ctx(py, self.pool.as_ref(), &self.var_ids, vars, "GbPoly")?; + let id = { + let pool = pool_py.borrow(py); + gbpoly_to_expr(&self.inner, &var_ids, &pool.inner) + }; + match id { + Some(id) => Ok(PyExpr { id, pool: pool_py }), + None => Err(pyo3::exceptions::PyValueError::new_err(format!( + "GbPoly is over {} variables but only {} were named; \ + pass the full vars list it was built with", + self.inner.n_vars, + var_ids.len() + ))), + } + } + fn __repr__(&self) -> String { format!("GbPoly(n_terms={})", self.inner.terms.len()) } } +/// `alkahest.expr_to_gbpoly(expr, vars)` — convert a polynomial :class:`Expr` +/// into the :class:`GbPoly` representation. +/// +/// The inverse of :meth:`GbPoly.to_expr`. Exponent slot `i` of the result +/// refers to ``vars[i]``, and the polynomial remembers `vars`, so it can be fed +/// straight to :meth:`GroebnerBasis.reduce`, :meth:`GroebnerBasis.contains` or +/// :meth:`GroebnerBasis.compute_raw`. +/// +/// Raises `ValueError` if *expr* is not polynomial in *vars* — a free symbol +/// outside *vars*, a negative or symbolic exponent, or a transcendental call. +/// +/// Example:: +/// +/// p = alkahest.expr_to_gbpoly(x**2 + y**2 - pool.integer(1), [x, y]) +/// gb = alkahest.GroebnerBasis.compute_raw([p]) +#[cfg(feature = "groebner")] +#[pyfunction] +#[pyo3(name = "expr_to_gbpoly", signature = (expr, vars))] +fn py_expr_to_gbpoly( + py: Python<'_>, + expr: PyRef, + vars: Vec>, +) -> PyResult { + if vars.is_empty() { + return Err(pyo3::exceptions::PyValueError::new_err( + "expr_to_gbpoly requires at least one variable", + )); + } + let pool_py = expr.pool.clone_ref(py); + let var_ids: Vec = vars.iter().map(|v| v.id).collect(); + let inner = { + let pool = pool_py.borrow(py); + expr_to_gbpoly(expr.id, &var_ids, &pool.inner) + .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))? + }; + Ok(PyGbPoly { + inner, + pool: Some(pool_py), + var_ids, + }) +} + +/// A computed Gröbner basis for a polynomial ideal. +/// +/// The basis is a **sequence**: it supports `len()`, integer indexing and +/// iteration, yielding its generators as :class:`GbPoly`. Read the generators +/// back as expressions with :meth:`to_exprs`, or one at a time with +/// :meth:`GbPoly.to_expr` — that is how an elimination ideal (a ``"lex"`` basis, +/// or the differential elimination performed by :func:`rosenfeld_groebner`) is +/// turned into readable relations. +/// +/// Attributes +/// ---------- +/// order : str +/// Monomial order the generators are reduced under: ``"lex"``, ``"grlex"`` +/// or ``"grevlex"``. +/// +/// Example:: +/// +/// gb = alkahest.GroebnerBasis.compute([x**2 + y**2 - one, x - y], [x, y]) +/// len(gb) # number of generators +/// [g.to_expr() for g in gb] # generators as Expr #[cfg(feature = "groebner")] #[pyclass(name = "GroebnerBasis")] struct PyGroebnerBasis { @@ -9842,6 +10269,35 @@ struct PyGroebnerBasis { var_ids: Vec, } +#[cfg(feature = "groebner")] +impl PyGroebnerBasis { + /// A generator wrapped with this basis's variable context. + fn wrap(&self, py: Python<'_>, p: GbPoly) -> PyGbPoly { + PyGbPoly::with_ctx(py, p, self.pool.as_ref(), &self.var_ids) + } + + /// Accept a `GbPoly` as-is, or convert an `Expr` using this basis's + /// variable ordering. + fn coerce_to_gbpoly(&self, py: Python<'_>, p: &Bound<'_, PyAny>) -> PyResult { + if let Ok(gbp) = p.downcast::() { + return Ok(gbp.borrow().inner.clone()); + } + if let Ok(expr) = p.downcast::() { + let pool_py = self.pool.as_ref().ok_or_else(|| { + pyo3::exceptions::PyValueError::new_err( + "GroebnerBasis has no variable context; use GroebnerBasis.compute() to build one that accepts Expr, or pass a GbPoly from expr_to_gbpoly()", + ) + })?; + let pool = pool_py.borrow(py); + return expr_to_gbpoly(expr.borrow().id, &self.var_ids, &pool.inner) + .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string())); + } + Err(pyo3::exceptions::PyTypeError::new_err( + "expected a GbPoly or an Expr", + )) + } +} + #[cfg(feature = "groebner")] #[pymethods] impl PyGroebnerBasis { @@ -9939,8 +10395,13 @@ impl PyGroebnerBasis { } /// Low-level entry point that accepts already-converted ``GbPoly`` objects, - /// bypassing the ``expr_to_gbpoly`` conversion. Useful when the polynomial - /// representation is already known (e.g., from ``MultiPoly`` reconstruction). + /// bypassing the :func:`expr_to_gbpoly` conversion. Useful when the + /// polynomial representation is already known (e.g., from ``MultiPoly`` + /// reconstruction). + /// + /// The variable context of the *first* input polynomial, if it has one, is + /// carried onto the resulting basis, so a basis built from + /// :func:`expr_to_gbpoly` output stays readable via :meth:`to_exprs`. /// /// Parameters /// ---------- @@ -9952,6 +10413,7 @@ impl PyGroebnerBasis { #[staticmethod] #[pyo3(signature = (gb_polys, order=None))] fn compute_raw( + py: Python<'_>, gb_polys: Vec>, order: Option<&str>, ) -> PyResult { @@ -9960,6 +10422,8 @@ impl PyGroebnerBasis { "GroebnerBasis.compute_raw requires at least one GbPoly", )); } + let pool = gb_polys[0].pool.as_ref().map(|p| p.clone_ref(py)); + let var_ids = gb_polys[0].var_ids.clone(); let raw: Vec = gb_polys.iter().map(|p| p.inner.clone()).collect(); let parsed_order = order .and_then(MonomialOrder::from_str) @@ -9970,44 +10434,177 @@ impl PyGroebnerBasis { }; Ok(PyGroebnerBasis { inner, - pool: None, - var_ids: vec![], + pool, + var_ids, }) } - fn reduce(&self, p: PyRef) -> PyGbPoly { - PyGbPoly { - inner: self.inner.reduce(&p.inner), + /// The monomial order the generators are reduced under. + #[getter] + fn order(&self) -> &'static str { + self.inner.order().as_str() + } + + /// The elimination ideal `I ∩ k[remaining vars]`, as a `GroebnerBasis`. + /// + /// Drops every generator whose support mentions one of *vars*. Under a + /// ``"lex"`` basis with the eliminated variables ordered **first**, what is + /// left is a Gröbner basis for the elimination ideal — the relations among + /// the remaining variables alone. Read them with :meth:`to_exprs`. + /// + /// The basis must know its variable ordering; *vars* must be among + /// :meth:`variables`. + /// + /// This is the implicitization move: parametrize a curve or surface, then + /// eliminate the parameters. + /// + /// Example:: + /// + /// # Implicitize (t, t**2): eliminate t from {x - t, y - t**2} + /// gb = alkahest.GroebnerBasis.compute([x - t, y - t**2], [t, x, y]) + /// gb.eliminate([t]).to_exprs() # [((y * -1) + x^2)] i.e. y = x**2 + fn eliminate(&self, py: Python<'_>, vars: Vec>) -> PyResult { + if self.pool.is_none() { + return Err(pyo3::exceptions::PyValueError::new_err( + "GroebnerBasis has no variable context; use GroebnerBasis.compute() to build one that can eliminate", + )); + } + let mut indices = Vec::with_capacity(vars.len()); + for v in &vars { + match self.var_ids.iter().position(|&id| id == v.id) { + Some(i) => indices.push(i), + None => { + return Err(pyo3::exceptions::PyValueError::new_err( + "eliminate() was given a variable this basis is not written over; \ + see GroebnerBasis.variables()", + )) + } + } + } + Ok(PyGroebnerBasis { + inner: self.inner.eliminate(&indices), + pool: self.pool.as_ref().map(|p| p.clone_ref(py)), + var_ids: self.var_ids.clone(), + }) + } + + /// The variables naming exponent slots 0, 1, … of the generators. + /// + /// Empty when the basis carries no variable context (a basis built by + /// :meth:`compute_raw` from context-free polynomials). + fn variables(&self, py: Python<'_>) -> Vec { + match &self.pool { + None => vec![], + Some(pool) => self + .var_ids + .iter() + .map(|&id| PyExpr { + id, + pool: pool.clone_ref(py), + }) + .collect(), + } + } + + /// The basis generators as :class:`GbPoly`, interreduced and monic. + /// + /// Equivalent to ``list(basis)``. + fn polynomials(&self, py: Python<'_>) -> Vec { + self.inner + .generators() + .iter() + .map(|p| self.wrap(py, p.clone())) + .collect() + } + + /// The basis generators as :class:`Expr`, each meaning ``g = 0``. + /// + /// This is the read path for elimination: with a ``"lex"`` basis, the + /// generators free of the eliminated variables are the eliminated + /// relations. + /// + /// Parameters + /// ---------- + /// vars : list[Expr], optional + /// Override the stored variable ordering — see :meth:`GbPoly.to_expr`. + /// + /// Example:: + /// + /// gb = alkahest.GroebnerBasis.compute([x**2 + y**2 - one, x - y], [x, y]) + /// gb.to_exprs() + #[pyo3(signature = (vars=None))] + fn to_exprs(&self, py: Python<'_>, vars: Option>>) -> PyResult> { + let (pool_py, var_ids) = + resolve_gb_ctx(py, self.pool.as_ref(), &self.var_ids, vars, "GroebnerBasis")?; + let ids: Option> = { + let pool = pool_py.borrow(py); + self.inner + .generators() + .iter() + .map(|g| gbpoly_to_expr(g, &var_ids, &pool.inner)) + .collect() + }; + match ids { + Some(ids) => Ok(ids + .into_iter() + .map(|id| PyExpr { + id, + pool: pool_py.clone_ref(py), + }) + .collect()), + None => Err(pyo3::exceptions::PyValueError::new_err(format!( + "basis is over more variables than the {} named; \ + pass the full vars list it was built with", + var_ids.len() + ))), } } + /// Reduce a polynomial modulo this basis and return the remainder. + /// + /// Accepts a :class:`GbPoly` or an :class:`Expr`; passing an ``Expr`` + /// requires the basis to know its variable ordering (see + /// :meth:`variables`). The remainder is a :class:`GbPoly` — call + /// :meth:`GbPoly.to_expr` on it to read it back. It is zero exactly when + /// :meth:`contains` is true. + fn reduce(&self, py: Python<'_>, p: &Bound<'_, PyAny>) -> PyResult { + let poly = self.coerce_to_gbpoly(py, p)?; + Ok(self.wrap(py, self.inner.reduce(&poly))) + } + /// Test membership. Accepts either a ``GbPoly`` or an ``Expr``; when /// passing an ``Expr`` the basis must have been created via ``compute()`` /// so that the variable order is known. fn contains(&self, py: Python<'_>, p: &Bound<'_, PyAny>) -> PyResult { - if let Ok(gbp) = p.downcast::() { - return Ok(self.inner.contains(&gbp.borrow().inner)); - } - if let Ok(expr) = p.downcast::() { - let pool_py = self.pool.as_ref().ok_or_else(|| { - pyo3::exceptions::PyValueError::new_err( - "GroebnerBasis has no variable context; use GroebnerBasis.compute() to build one that accepts Expr", - ) - })?; - let pool = pool_py.borrow(py); - let gbp = expr_to_gbpoly(expr.borrow().id, &self.var_ids, &pool.inner) - .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?; - return Ok(self.inner.contains(&gbp)); - } - Err(pyo3::exceptions::PyTypeError::new_err( - "contains() expects a GbPoly or an Expr", - )) + let poly = self.coerce_to_gbpoly(py, p)?; + Ok(self.inner.contains(&poly)) } fn __len__(&self) -> usize { self.inner.len() } + /// Generator `i` as a :class:`GbPoly`; negative indices count from the end. + fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult { + let n = self.inner.len() as isize; + let i = if index < 0 { index + n } else { index }; + if i < 0 || i >= n { + return Err(pyo3::exceptions::PyIndexError::new_err( + "GroebnerBasis index out of range", + )); + } + Ok(self.wrap(py, self.inner.generators()[i as usize].clone())) + } + + /// Iterate over the generators as :class:`GbPoly`. + fn __iter__(&self, py: Python<'_>) -> PyResult { + let list = pyo3::types::PyList::empty_bound(py); + for p in self.inner.generators() { + list.append(Py::new(py, self.wrap(py, p.clone()))?)?; + } + Ok(list.as_any().iter()?.into_py(py)) + } + fn __repr__(&self) -> String { format!("GroebnerBasis(n_generators={})", self.inner.len()) } @@ -10020,7 +10617,23 @@ fn py_monomial_order_for_dae(order: Option<&str>) -> MonomialOrder { .unwrap_or(MonomialOrder::GRevLex) } -/// V2-13 — Rosenfeld–Gröbner-style differential elimination result. +/// V2-13 — result of Rosenfeld–Gröbner-style differential elimination. +/// +/// Returned by :func:`rosenfeld_groebner`. The eliminated relations are in +/// :meth:`final_basis`; read them with ``result.final_basis().to_exprs()``. +/// +/// Attributes +/// ---------- +/// consistent : bool +/// ``False`` iff the unit ideal was reached — the system has no common +/// jet solution over ℚ, i.e. the equations are contradictory. +/// truncated : bool +/// ``True`` if prolongation stopped at ``max_prolong_rounds`` rather than +/// because the differential chain saturated. A truncated basis is a +/// *sound* set of consequences of the system but need not be complete, so +/// "not in the basis" does not mean "not a consequence". +/// prolongation_rounds : int +/// Number of prolongation rounds that contributed new relations. #[cfg(feature = "groebner")] #[pyclass(name = "RosenfeldGroebnerResult")] struct PyRosenfeldGroebnerResult { @@ -10033,11 +10646,15 @@ struct PyRosenfeldGroebnerResult { working_dae: DAE, final_basis: Option, pool: Py, + /// Jet variables indexing the exponent vectors of `final_basis`. + var_ids: Vec, } #[cfg(feature = "groebner")] #[pymethods] impl PyRosenfeldGroebnerResult { + /// The prolonged :class:`DAE`: the input system plus the derivative jets + /// introduced while differentiating it. fn working_dae(&self, py: Python<'_>) -> PyDAE { PyDAE { inner: self.working_dae.clone(), @@ -10045,6 +10662,33 @@ impl PyRosenfeldGroebnerResult { } } + /// The jet variables indexing the basis, in exponent-slot order. + /// + /// These are the symbols the elimination actually ran over — the time + /// variable, the declared states and derivatives, and every higher jet + /// (``d2x/dt2``, …) introduced by prolongation. + fn variables(&self, py: Python<'_>) -> Vec { + self.var_ids + .iter() + .map(|&id| PyExpr { + id, + pool: self.pool.clone_ref(py), + }) + .collect() + } + + /// The saturated Gröbner basis, or ``None`` when the system is + /// inconsistent. + /// + /// The returned :class:`GroebnerBasis` knows its variable ordering, so + /// ``final_basis().to_exprs()`` gives the eliminated relations as + /// :class:`Expr`. + /// + /// Example:: + /// + /// r = alkahest.rosenfeld_groebner(dae, max_prolong_rounds=1) + /// for eq in r.final_basis().to_exprs(): + /// print(eq, "= 0") fn final_basis(&self, py: Python<'_>) -> PyResult>> { match &self.final_basis { None => Ok(None), @@ -10053,7 +10697,7 @@ impl PyRosenfeldGroebnerResult { PyGroebnerBasis { inner: gb.clone(), pool: Some(self.pool.clone_ref(py)), - var_ids: vec![], + var_ids: self.var_ids.clone(), }, )?)), } @@ -10072,6 +10716,8 @@ impl PyRosenfeldGroebnerResult { struct PyDaeIndexReduction { inner: DaeIndexReduction, pool: Py, + /// Jet variables for the Gröbner fallback; empty when Pantelides won. + var_ids: Vec, } #[cfg(feature = "groebner")] @@ -10110,6 +10756,7 @@ impl PyDaeIndexReduction { working_dae: r.working_dae.clone(), final_basis: r.final_basis.clone(), pool: self.pool.clone_ref(py), + var_ids: self.var_ids.clone(), }, ) .ok(), @@ -10128,6 +10775,39 @@ impl PyDaeIndexReduction { } } +/// `alkahest.rosenfeld_groebner(dae, order=None, max_prolong_rounds=None)` — +/// Rosenfeld–Gröbner-style differential elimination. +/// +/// Prolongs the system (differentiates each equation, introducing higher jets +/// as new indeterminates) and computes a Gröbner basis after each round, until +/// differentiating adds nothing new to the ideal or the round budget runs out. +/// The basis is the set of *algebraic consequences* of the differential +/// system — the input–output relations elimination is after. +/// +/// Parameters +/// ---------- +/// dae : DAE +/// The system, polynomial in its variables and derivative symbols. +/// order : str, optional +/// Monomial order — ``"grevlex"`` (default), ``"grlex"`` or ``"lex"``. +/// Use ``"lex"`` when you want elimination-ordered generators. +/// max_prolong_rounds : int, optional +/// Prolongation budget (default 8). Nonlinear jets often do not saturate +/// in finitely many algebraic steps, so hitting the budget is normal and +/// sets :attr:`RosenfeldGroebnerResult.truncated`. +/// +/// Returns +/// ------- +/// RosenfeldGroebnerResult +/// Read the relations with ``result.final_basis().to_exprs()``. +/// +/// Example:: +/// +/// t, x, dx = p.symbol("t"), p.symbol("x"), p.symbol("dx/dt") +/// dae = alkahest.DAE.new([dx - x], [x], [dx], t) +/// r = alkahest.rosenfeld_groebner(dae, max_prolong_rounds=1) +/// r.consistent # True +/// r.final_basis().to_exprs() # the eliminated relations, as Expr #[cfg(feature = "groebner")] #[pyfunction] #[pyo3(name = "rosenfeld_groebner", signature = (dae, order=None, max_prolong_rounds=None))] @@ -10140,14 +10820,14 @@ fn py_rosenfeld_groebner( let pool_py = dae.pool.clone_ref(py); let r = { let pool = pool_py.borrow(py); - rosenfeld_groebner_with_options( + rosenfeld_groebner_ranked( &dae.inner, &pool.inner, py_monomial_order_for_dae(order), max_prolong_rounds.unwrap_or(8), ) }; - let r = r.map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?; + let (r, ranking) = r.map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?; Ok(PyRosenfeldGroebnerResult { consistent: r.consistent, truncated: r.truncated, @@ -10155,6 +10835,7 @@ fn py_rosenfeld_groebner( working_dae: r.working_dae, final_basis: r.final_basis, pool: pool_py, + var_ids: ranking.vars, }) } @@ -10167,14 +10848,16 @@ fn py_dae_index_reduce( order: Option<&str>, ) -> PyResult { let pool_py = dae.pool.clone_ref(py); - let inner = { + let out = { let pool = pool_py.borrow(py); - dae_index_reduce(&dae.inner, &pool.inner, py_monomial_order_for_dae(order)) + dae_index_reduce_ranked(&dae.inner, &pool.inner, py_monomial_order_for_dae(order)) }; - let inner = inner.map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?; + let (inner, ranking) = + out.map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?; Ok(PyDaeIndexReduction { inner, pool: pool_py, + var_ids: ranking.map(|r| r.vars).unwrap_or_default(), }) } @@ -10750,7 +11433,7 @@ fn py_solve( let solutions: Vec> = values.into_iter().map(|v| vec![v]).collect(); let result: Result = Ok(SolutionSet::Finite(solutions)); - return finite_solutions_to_py(py, result, &pool_py, &var_ids, numeric); + return finite_solutions_to_py(py, result, &pool_py, &var_ids, &var_ids, numeric); } } @@ -10794,7 +11477,20 @@ fn py_solve( } } - finite_solutions_to_py(py, result, &pool_py, &var_ids, numeric) + // A `Parametric` basis is indexed by the solve variables *followed by* the + // free parameters — the same concatenation `solve_polynomial_system` builds + // its exponent vectors from. Without it the returned basis cannot be read. + let basis_var_ids: Vec = { + let pool = pool_py.borrow(py); + let mut all = var_ids.clone(); + all.extend(alkahest_core::solver::collect_parameters( + &eq_ids, + &var_ids, + &pool.inner, + )); + all + }; + finite_solutions_to_py(py, result, &pool_py, &var_ids, &basis_var_ids, numeric) } #[cfg(feature = "groebner")] @@ -10860,6 +11556,7 @@ fn finite_solutions_to_py( result: Result, pool_py: &Py, var_ids: &[ExprId], + basis_var_ids: &[ExprId], numeric: bool, ) -> PyResult { match result { @@ -10870,8 +11567,8 @@ fn finite_solutions_to_py( Ok(SolutionSet::NoSolution) => Ok(pyo3::types::PyList::empty_bound(py).into()), Ok(SolutionSet::Parametric(gb)) => Ok(PyGroebnerBasis { inner: gb, - pool: None, - var_ids: vec![], + pool: Some(pool_py.clone_ref(py)), + var_ids: basis_var_ids.to_vec(), } .into_py(py)), Ok(SolutionSet::Finite(solutions)) => { @@ -10926,9 +11623,17 @@ fn finite_solutions_to_py( // --------------------------------------------------------------------------- #[cfg(feature = "groebner")] +/// One triangular component of a decomposition, as returned by +/// :func:`triangularize`. +/// +/// The component's polynomials are :class:`GbPoly`; each carries the variable +/// ordering :func:`triangularize` was called with, so ``p.to_expr()`` reads it +/// back as an :class:`Expr`. #[pyclass(name = "RegularChain")] struct PyRegularChain { inner: RegularChain, + pool: Option>, + var_ids: Vec, } #[cfg(feature = "groebner")] @@ -10940,14 +11645,60 @@ impl PyRegularChain { } /// Gröbner-style polynomial tiles (``GbPoly``), ascending by main variable. - fn polys(&self) -> Vec { + /// + /// Each tile knows the variables it is written over, so + /// ``[p.to_expr() for p in chain.polys()]`` gives the triangular system as + /// :class:`Expr` (see :meth:`to_exprs`). + fn polys(&self, py: Python<'_>) -> Vec { self.inner .polys .iter() - .map(|p| PyGbPoly { inner: p.clone() }) + .map(|p| PyGbPoly::with_ctx(py, p.clone(), self.pool.as_ref(), &self.var_ids)) .collect() } + /// The variables the chain is written over, in exponent-slot order. + fn variables(&self, py: Python<'_>) -> Vec { + match &self.pool { + None => vec![], + Some(pool) => self + .var_ids + .iter() + .map(|&id| PyExpr { + id, + pool: pool.clone_ref(py), + }) + .collect(), + } + } + + /// The triangular system as :class:`Expr`, each meaning ``p = 0``. + fn to_exprs(&self, py: Python<'_>) -> PyResult> { + let (pool_py, var_ids) = + resolve_gb_ctx(py, self.pool.as_ref(), &self.var_ids, None, "RegularChain")?; + let ids: Option> = { + let pool = pool_py.borrow(py); + self.inner + .polys + .iter() + .map(|p| gbpoly_to_expr(p, &var_ids, &pool.inner)) + .collect() + }; + ids.map(|ids| { + ids.into_iter() + .map(|id| PyExpr { + id, + pool: pool_py.clone_ref(py), + }) + .collect() + }) + .ok_or_else(|| { + pyo3::exceptions::PyValueError::new_err( + "regular chain is over more variables than were named", + ) + }) + } + fn __len__(&self) -> usize { self.inner.len() } @@ -10981,7 +11732,7 @@ fn py_triangularize( let result = { let pool = pool_py.borrow(py); - triangularize(eq_ids, var_ids, &pool.inner) + triangularize(eq_ids, var_ids.clone(), &pool.inner) }; match result { @@ -11003,7 +11754,14 @@ fn py_triangularize( Ok(chains) => { let list = pyo3::types::PyList::empty_bound(py); for chain in chains { - list.append(PyRegularChain { inner: chain }.into_py(py))?; + list.append( + PyRegularChain { + inner: chain, + pool: Some(pool_py.clone_ref(py)), + var_ids: var_ids.clone(), + } + .into_py(py), + )?; } Ok(list.into()) } @@ -11024,10 +11782,14 @@ struct PyMultiPolyFp { #[pymethods] impl PyMultiPolyFp { + /// True if this is the zero polynomial. + #[getter] fn is_zero(&self) -> bool { self.inner.is_zero() } + /// Highest total degree over all terms (`0` for the zero polynomial). + #[getter] fn total_degree(&self) -> u32 { self.inner.total_degree() } @@ -11631,6 +12393,7 @@ fn alkahest(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(py_ideal_radical, m)?)?; m.add_function(wrap_pyfunction!(py_rosenfeld_groebner, m)?)?; m.add_function(wrap_pyfunction!(py_dae_index_reduce, m)?)?; + m.add_function(wrap_pyfunction!(py_expr_to_gbpoly, m)?)?; } // V2-2 — Resultants and subresultant PRS m.add_function(wrap_pyfunction!(py_resultant, m)?)?; diff --git a/alkahest-skill/alkahest.md b/alkahest-skill/alkahest.md index ab7b5f7e..2e7bb02f 100644 --- a/alkahest-skill/alkahest.md +++ b/alkahest-skill/alkahest.md @@ -518,7 +518,7 @@ from alkahest import UniPoly, MultiPoly, RationalFunction # Dense univariate polynomial p = UniPoly.from_symbolic(x**3 + pool.integer(-2)*x + pool.integer(1), x) -p.degree() # 3 +p.degree # 3 p.coefficients() # [1, -2, 0, 1] (constant first) # `UniPoly` has no numeric-eval method. Its full surface is: # from_symbolic, from_coefficients, coefficients, degree, is_zero, gcd, factor_z @@ -532,7 +532,7 @@ a.gcd(b) # x - 1 # Sparse multivariate polynomial (over ℤ) mp = MultiPoly.from_symbolic(x**2 * y + x * y**2, [x, y]) -mp.total_degree() # 3 +mp.total_degree # 3 # Rational function (GCD-normalized automatically) rf = RationalFunction.from_symbolic(x**2 + pool.integer(-1), x + pool.integer(-1), [x]) @@ -866,7 +866,7 @@ y = pool.symbol("y") k = pool.symbol("k") ode = ODE.new([y], [pool.integer(-1)*k*y], t) -ode.order() +ode.order ode.is_autonomous() ode.state_vars() ode.rhs() diff --git a/docs/mdbook/src/ode-dae.md b/docs/mdbook/src/ode-dae.md index 1c00c781..248e05b3 100644 --- a/docs/mdbook/src/ode-dae.md +++ b/docs/mdbook/src/ode-dae.md @@ -2,88 +2,158 @@ Alkahest provides symbolic infrastructure for ordinary differential equations (ODEs) and differential-algebraic equations (DAEs), including structural analysis and automatic index reduction. +Every snippet on this page is executed by `tests/test_docs_ode_dae.py`, so it runs against the version of Alkahest it ships with. + ## ODE -`ODE` represents an ordinary differential equation system. Build one from symbolic expressions: +`ODE` represents an explicit first-order system `d(state_vars)/dt = rhs`. The constructor takes three **positional** arguments — state variables, right-hand sides, and the independent variable — and `ODE(...)` and `ODE.new(...)` are the same thing: ```python -from alkahest import ExprPool, ODE, lower_to_first_order, sin +from alkahest import ExprPool, ODE pool = ExprPool() t = pool.symbol("t") x = pool.symbol("x") v = pool.symbol("v") -# Simple harmonic oscillator: x'' + ω²x = 0 -# Represented as a first-order system: [x' = v, v' = -ω²x] -omega = pool.integer(1) -ode = ODE( - state=[x, v], - derivatives=[v, pool.integer(-1) * omega**pool.integer(2) * x], - independent=t, -) +# Simple harmonic oscillator x'' + x = 0, written as a first-order system: +# x' = v, v' = -x +ode = ODE.new([x, v], [v, pool.integer(-1) * x], t) + +ode.order # 2 — number of state variables +ode.state_vars() # [x, v] +ode.rhs() # [v, (x * -1)] +ode.is_autonomous() # True — t does not appear on the right +``` + +Initial conditions are attached one at a time, and `with_ic` returns a new `ODE` rather than mutating: + +```python +ode_with_ic = ode.with_ic(x, pool.integer(1)).with_ic(v, pool.integer(0)) ``` ### Lowering to first order -Higher-order ODEs are automatically lowered to first-order form: +`lower_to_first_order` takes the *pieces* of a scalar higher-order equation — the unknown, the right-hand side, the order, and the independent variable — not an existing `ODE`: ```python from alkahest import lower_to_first_order -first_order_system = lower_to_first_order(higher_order_ode) +# x'' = -4x → [x' = x_1, x_1' = -4x] +ode = lower_to_first_order(x, pool.integer(-4) * x, 2, t) + +ode.state_vars() # [x, x_1] +ode.rhs() # [x_1, (x * -4)] ``` +The auxiliary states are named `x_1`, `x_2`, … for the successive derivatives. + ## DAE -`DAE` represents a differential-algebraic system where some equations are algebraic constraints rather than differential equations. +`DAE` represents a system of implicit equations mixing differential and algebraic constraints, `F(t, variables, derivatives) = 0`. Build one with `DAE.new(equations, variables, derivatives, time_var)` — four positional arguments: + +- `equations` — each expression means `g = 0`. Write `x' = f` as `dx - f`. +- `variables` — the dependent variables. +- `derivatives` — a **separate symbol** standing for the time derivative of `variables[i]`, e.g. `pool.symbol("dx/dt")`. Alkahest does not parse the name; the positional pairing is what makes it a derivative. +- `time_var` — the independent variable. + +A variable with **no** entry in `derivatives` is purely algebraic. That is how you declare a Lagrange multiplier, and it is what makes the system high-index: ```python -from alkahest import DAE, pantelides +from alkahest import DAE, ExprPool pool = ExprPool() t = pool.symbol("t") -x = pool.symbol("x") # differential variable -y = pool.symbol("y") # algebraic variable (constrained) -lam = pool.symbol("lam") # Lagrange multiplier - -# Pendulum: differential equations + constraint -dae = DAE( - equations=[...], # system of equations - variables=[x, y, lam], - independent=t, +x, y, u, w = (pool.symbol(n) for n in ("x", "y", "u", "w")) +lam = pool.symbol("lam") # Lagrange multiplier — algebraic +dx, dy, du, dw = (pool.symbol(n) for n in ("dx/dt", "dy/dt", "du/dt", "dw/dt")) +one, two = pool.integer(1), pool.integer(2) + +# Cartesian pendulum (index 3): x' = u, y' = w, u' = -lam*x, w' = -lam*y - 1, +# subject to x**2 + y**2 = 1. +dae = DAE.new( + [dx - u, dy - w, du + lam * x, dw + lam * y + one, x**two + y**two - one], + [x, y, u, w, lam], # five variables … + [dx, dy, du, dw], # … but only four have derivatives + t, ) + +dae.n_equations # 5 +dae.n_variables # 5 +dae.equations() # the five expressions, each meaning "= 0" +dae.derivatives() # [dx/dt, dy/dt, du/dt, dw/dt] +dae.time_var # t ``` ## Pantelides algorithm -The Pantelides algorithm performs structural index reduction on DAEs. It identifies which equations need to be differentiated to make the system structurally regular: +The Pantelides algorithm performs *structural* index reduction: it repeatedly differentiates the equations that a maximum bipartite matching leaves unmatched, until every equation is matched to a variable. It returns the **reduced `DAE`** — not a separate report object: ```python from alkahest import pantelides reduced = pantelides(dae) -print(reduced.index) # structural index of the reduced system -print(reduced.differentiated) # which equations were differentiated + +reduced.index # 1 — differentiation rounds used +reduced.n_equations # 6 — one more than the input +reduced.equations()[-1] +# ((x * dx/dt * 2) + (y * dy/dt * 2)) — the differentiated constraint ``` -Index reduction converts a high-index DAE (index > 1) into an index-1 system that ODE solvers can handle. The result includes the differentiated equations as symbolic expressions. +`reduced.index` is the number of rounds, so `0` means the input already had a perfect structural matching and nothing was differentiated. The equations appended by differentiation are visible in `reduced.equations()`, and the higher jets they introduce (`ddx/dt/dt`, …) in `reduced.derivatives()`. + +`pantelides` raises `ValueError` (`E-DAE-002`) above index 10. `dae_index_reduce(dae)` runs Pantelides first and falls back to `rosenfeld_groebner` when it hits that cap. + +Because the algorithm looks only at *which* variables occur in *which* equations, it can reduce a system whose coefficients make it unsolvable; structural regularity is not numerical regularity. + +## Differential elimination + +`rosenfeld_groebner` prolongs the system — differentiating each equation and treating the new jets as fresh indeterminates — and computes a Gröbner basis after each round. The basis holds the algebraic consequences of the differential system, which is what elimination-based work (input–output equations, structural identifiability) needs: + +```python +from alkahest import DAE, ExprPool, rosenfeld_groebner + +pool = ExprPool() +t, x, dx = pool.symbol("t"), pool.symbol("x"), pool.symbol("dx/dt") + +dae = DAE.new([dx - x], [x], [dx], t) # x' = x +result = rosenfeld_groebner(dae, max_prolong_rounds=1) + +result.consistent # True — the unit ideal was never reached +result.truncated # True — stopped at the round budget, not at saturation +result.variables() # [t, x, dx/dt, ddx/dt/dt] — jets, in exponent-slot order + +[str(e) for e in result.final_basis().to_exprs()] +# ['(x + (-1 * ddx/dt/dt))', '(dx/dt + (-1 * ddx/dt/dt))'] +``` + +`final_basis()` returns a `GroebnerBasis` that knows its variable ordering, so `to_exprs()` reads the relations back as `Expr` (each meaning `= 0`). See [Solving](./solving.md) for the rest of the `GroebnerBasis` surface. + +`truncated=True` means prolongation stopped because `max_prolong_rounds` ran out, not because the chain saturated. A truncated basis is a *sound* set of consequences but need not be complete — "not in the basis" then does not mean "not a consequence". Nonlinear jets often do not saturate in finitely many algebraic steps, so this is the common case. + +`final_basis()` returns `None` when `consistent` is `False`. ## Sensitivity analysis -Sensitivity analysis computes how solutions depend on parameters: +Sensitivity analysis computes how solutions depend on parameters. `sensitivity_system(ode, params)` augments the state with `∂x/∂p`; `adjoint_system(ode, objective_grad)` takes the gradient of the objective with respect to the state, as a list parallel to `ode.state_vars()`: ```python -from alkahest import sensitivity_system, adjoint_system +from alkahest import ExprPool, ODE, adjoint_system, sensitivity_system pool = ExprPool() -p = pool.symbol("p") # parameter +t, y, k = pool.symbol("t"), pool.symbol("y"), pool.symbol("k") + +ode = ODE.new([y], [pool.integer(-1) * k * y], t) # y' = -k*y -# Forward sensitivity: generates ∂x/∂p equations alongside the ODE -sens = sensitivity_system(ode, [p]) +sens = sensitivity_system(ode, [k]) +sens.original_dim # 1 +sens.n_params # 1 +sens.extended_ode.state_vars() # [y, dS_y_k] -# Adjoint method: more efficient for many parameters, one output -adj = adjoint_system(ode, output_expr, [p]) +# Objective J = y(T)**2 → dJ/dy = 2y +adj = adjoint_system(ode, [pool.integer(2) * y]) +adj.state_vars() # [lambda_y] +adj.rhs() # [(k * lambda_y)] ``` ## Acausal modeling @@ -113,6 +183,7 @@ circuit.connect(cap.port("C1.n"), src.port("V1.n")) # Flatten into a DAE dae = circuit.flatten(t) +dae.n_equations # 10 ``` Built-in components (`resistor`, `capacitor`, `voltage_source`, and others registered via the component API) generate their constitutive equations automatically. `AcausalSystem.flatten` assembles them, plus the Kirchhoff-style connection equations, into a `DAE` that Pantelides can reduce. See `examples/acausal_and_laplace.py` for a runnable end-to-end example. @@ -122,24 +193,32 @@ Built-in components (`resistor`, `capacitor`, `voltage_source`, and others regis The Laplace transform lives in `alkahest.experimental` (the calculus/ODE/transform surface is not yet semver-frozen): ```python +from alkahest import ExprPool from alkahest.experimental import inverse_laplace_transform, laplace_transform +pool = ExprPool() s, t = pool.symbol("s"), pool.symbol("t") -F = laplace_transform(pool.integer(1), t, s) # 1/s -f = inverse_laplace_transform(F, s, t) # back to 1 (Heaviside(t)) +F = laplace_transform(pool.integer(1), t, s) # s^-1 +f = inverse_laplace_transform(F, s, t) # back to 1 ``` ## Hybrid systems -`HybridODE` adds event handling to an ODE: at a crossing event, the state is reset and integration resumes with a new ODE: +`HybridODE` adds event handling to an ODE: at a crossing event, the state is reset and integration resumes. `Event.new(name, condition, reset_map)` takes the reset map as a list of `(variable, new_value)` pairs, and both `HybridODE.new` and `add_event` return new objects: ```python -from alkahest import HybridODE, Event +from alkahest import Event, ExprPool, HybridODE, ODE -# Bouncing ball: velocity reverses at floor contact -bounce_event = Event( - condition=x, # fires when x = 0 - reset={v: pool.integer(-1) * v}, # reverse velocity -) -hybrid = HybridODE(ode=base_ode, events=[bounce_event]) +pool = ExprPool() +t, x, v = pool.symbol("t"), pool.symbol("x"), pool.symbol("v") + +# Bouncing ball: x' = v, v' = -1; velocity reverses at floor contact +base_ode = ODE.new([x, v], [v, pool.integer(-1)], t) +bounce = Event.new("bounce", x, [(v, pool.integer(-1) * v)]) + +hybrid = HybridODE.new(base_ode).add_event(bounce) +hybrid.n_events # 1 +hybrid.guards() # [x] — fires when x crosses 0 ``` + +`Event.rising()` and `Event.falling()` restrict an event to one crossing direction; by default it fires in both. diff --git a/docs/mdbook/src/representations.md b/docs/mdbook/src/representations.md index d4d4e749..101d9daa 100644 --- a/docs/mdbook/src/representations.md +++ b/docs/mdbook/src/representations.md @@ -42,7 +42,7 @@ from alkahest import UniPoly # x^3 - 2x + 1 p = UniPoly.from_symbolic(x**3 + pool.integer(-2) * x + pool.integer(1), x) -print(p.degree()) # 3 +print(p.degree) # 3 print(p.coefficients()) # [1, -2, 0, 1] (constant first) print(p.leading_coeff()) # 1 @@ -71,7 +71,7 @@ from alkahest import MultiPoly expr = x**2 * y + x * y**2 + pool.integer(-1) mp = MultiPoly.from_symbolic(expr, [x, y]) -print(mp.total_degree()) # 3 +print(mp.total_degree) # 3 print(mp.integer_content()) # 1 # Arithmetic diff --git a/docs/mdbook/src/solving.md b/docs/mdbook/src/solving.md index 8dc634f7..ae807e40 100644 --- a/docs/mdbook/src/solving.md +++ b/docs/mdbook/src/solving.md @@ -53,19 +53,53 @@ Pass `numeric=True` to return float values directly: `solve(eqs, vars, numeric=T A `GroebnerBasis` can be constructed directly for ideal-theoretic operations: ```python -from alkahest import GroebnerBasis, GbPoly +from alkahest import GroebnerBasis -# Compute a Gröbner basis under GrLex order +# Compute a Gröbner basis (lex by default) polys = [x**2 + y**2 - pool.integer(1), x - y] gb = GroebnerBasis.compute(polys, [x, y]) # Check ideal membership print(gb.contains(x - pool.rational(1, 2))) # False -# Reduce a polynomial modulo the ideal +# Reduce a polynomial modulo the ideal — the remainder is a GbPoly reduced = gb.reduce(x**3 + y**3) +print(reduced.to_expr()) # y ``` +### Reading a basis + +A `GroebnerBasis` is a sequence of `GbPoly`, and each `GbPoly` converts back to an `Expr`. This is how you read an elimination result — the generators of a `Lex` basis that are free of the eliminated variables *are* the eliminated relations: + +```python +len(gb) # 2 — number of generators +gb.order # "lex" +gb.variables() # [x, y] — what exponent slots 0, 1 refer to + +for g in gb: + print(g.to_expr(), "= 0") +# (y^2 + -1/2) = 0 +# (x + (y * -1)) = 0 + +gb.to_exprs() # the same list in one call +gb[0].terms() # [((0, 0), Fraction(-1, 2)), ((0, 2), 1)] +``` + +`terms()` gives `(exponent tuple, coefficient)` pairs with the coefficient as an exact Python `int` or `fractions.Fraction`; the exponent tuple is parallel to `variables()`. + +The conversion runs the other way with `expr_to_gbpoly`, which is what `reduce` and `contains` accept alongside plain `Expr`: + +```python +from alkahest import expr_to_gbpoly + +p = expr_to_gbpoly(x**2 + y**2 - pool.integer(1), [x, y]) +p.n_terms # 3 +gb.contains(p) # True +GroebnerBasis.compute_raw([p]) +``` + +A `GbPoly` stores exponent vectors, not names, so converting one back needs the variable list its slots refer to. Every basis Alkahest hands out carries that list — including the ones from `solve` (solve variables followed by the free parameters), `triangularize` and `rosenfeld_groebner` — so `to_expr()` normally takes no arguments. Naming too few variables raises `ValueError` rather than quietly misreading the exponent slots. + ### Monomial orders @@ -85,14 +119,18 @@ Because the Rust entry point falls back to CPU row reduction when no device is p ## Elimination ideals -`GroebnerBasis.eliminate` computes the elimination ideal by dropping generators involving specified variables: +`GroebnerBasis.eliminate` computes the elimination ideal `I ∩ k[remaining vars]` by dropping every generator whose support mentions one of the given variables. Under a `lex` basis with the eliminated variables ordered **first**, what is left is a Gröbner basis for that ideal: ```python -# Eliminate y to get a univariate ideal in x -x_ideal = gb.eliminate([y]) +# Implicitize the parametric curve (t, t**2): eliminate the parameter t. +gb = GroebnerBasis.compute([x - t, y - t**2], [t, x, y]) +gb.to_exprs() # [(t + (x * -1)), ((y * -1) + x^2)] + +implicit = gb.eliminate([t]) +implicit.to_exprs() # [((y * -1) + x^2)] — y = x**2 ``` -This is the algebraic geometry operation underlying implicitization of parametric curves and surfaces. +Note the variable order passed to `compute`: `t` comes first, so `lex` eliminates it. `eliminate` requires the basis to know its variables (`gb.variables()`), and rejects a variable it is not written over. ## Performance diff --git a/docs/mdbook/src/validated-bounds.md b/docs/mdbook/src/validated-bounds.md index a20d2250..42a730e0 100644 --- a/docs/mdbook/src/validated-bounds.md +++ b/docs/mdbook/src/validated-bounds.md @@ -88,6 +88,33 @@ enclosures lies strictly on the wrong side. Without it, `x > 0` on `[-1,1]` — plainly false — could only ever be `"undecided"`, since the range enclosure straddles zero by construction. +### Proving a root *exists* + +`verified_no_roots` returns `"false"` only with a proof in hand, and the proof +is the intermediate value theorem. The full-box enclosure succeeding is already +a continuity certificate — a Taylor model is only ever built where every +elementary step stayed strictly inside its domain — and a box is convex. So if +two points of the box can be found where `f` is *rigorously proven* to have +opposite signs, the segment between them stays in the box and `f` vanishes +somewhere on it. + +The two points do not have to be the box's own endpoints, and that is what makes +the test usable: the search subdivides the box, records the sign of any sub-box +whose enclosure has a determined one, and samples the centres of the rest. + +| Box | Roots inside | Endpoint signs | Verdict | +|---|---|---|---| +| `x²−2` on `[-2,0]` | 1 | + → − | `"false"` | +| `x²−2` on `[-2,2]` | 2 | + → + | `"false"` | +| `x²−2` on `[-10,10]` | 2 | + → + | `"false"` | +| `(x²−2)(x²+1)` on `[-2,2]` | 2 | + → + | `"false"` | +| `x−y` on `[-1,1]²` | a whole line | — | `"false"` | +| `(x−1)²` on `[0,2]` | 1 (double) | + → + | `"undecided"` | + +The last row is the honest limit. A double root never changes sign, so no +witness pair exists; `"undecided"` is the answer, and it is not upgraded to +`"false"` on the strength of an enclosure that merely touches zero. + ## Refusals | Code | Meaning | @@ -104,6 +131,62 @@ genuine interior pole. Only after the box has been bisected far enough for that explanation to be exhausted does it refuse — which is the right answer for a real interior singularity, where the range is not a bounded interval at all. +## Removable singularities in `verified_integral` + +`∫₀¹ ln(1+x)/x dx = π²/12` has nothing singular about it — only the *expression* +is singular at `x = 0`, and the integrand extends continuously to 1 there. A +Taylor model still refuses, because the reciprocal's enclosure contains zero. + +`verified_integral` recognises this shape. If the integrand splits as `N(x)/D(x)` +and there is a point `p` of the offending sub-interval at which `N` and `D` both +vanish, it enclosures that piece with **Cauchy's mean value theorem** instead: + +```text +N(p) = D(p) = 0, D' ≠ 0 on J ⟹ ∀ x ∈ J\{p} : N(x)/D(x) = N'(ξ)/D'(ξ) for some ξ ∈ J + ⟹ ∫_J N/D dx ∈ |J| · range(N'/D' on J) +``` + +so the piece is bounded by an enclosure of `N'/D'`, which is perfectly regular. +The number returned is the integral of the continuous extension. + +Three guards keep this from swallowing a genuine pole: + +- `N(p) = 0` and `D(p) = 0` are established **symbolically** (substitute the + exact rational `p`, simplify, require a literal zero). No numeric enclosure + can prove a value is exactly zero, so none is asked to. +- `D'` must be *certified non-vanishing* on the sub-interval. That is what fails + for `sin(x)/x²`, where the denominator has a double zero and the integral does + not converge. +- `N` and `D` must each have a successful enclosure over the whole + sub-interval, which certifies they are analytic — and hence that the symbolic + derivatives really are their derivatives. + +```python +ak.verified_integral(ak.log(pool.integer(1) + x) / x, x, 0.0, 1.0) # ≈ π²/12 +ak.verified_integral(ak.sin(x) / x, x, -1.0, 1.0) # ≈ 1.8921661 +ak.verified_integral(pool.integer(1) / x, x, -1.0, 1.0) # refuses: N(0) ≠ 0 +``` + +### What is still refused + +An **integrable but non-removable** singularity is refused, and the message says +so rather than implying the integral does not exist: + +| Integral | Value | Status | +|---|---|---| +| `∫₀¹ ln(1+x)/x dx` | `π²/12` | enclosed (removable) | +| `∫_{-1}^{1} sin(x)/x dx` | `2·Si(1)` | enclosed (removable) | +| `∫₀¹ −ln x dx` | 1 | refused — `log` enclosure reaches 0, not a `0/0` quotient | +| `∫₀¹ (ln x)² dx` | 2 | refused, same reason | +| `∫₀¹ dx/√(1−x²) dx` | `π/2` | refused — endpoint singularity, numerator does not vanish | +| `∫₀¹ xˣ dx` | 0.78343… | refused — `log` enclosure reaches 0 | +| `∫₀¹ ln(x)·ln(1−x) dx` | `2 − π²/6` | refused — singular at both ends | + +These need an integrable-tail bound or a singularity-removing substitution, +neither of which can be derived rigorously from the expression alone today. The +refusal is the honest answer; widening an enclosure to make them pass would +break the contract that makes the module worth using. + ## Relation to the rest of the stack This is the slow, certifying half of *falsify fast, certify slow*: @@ -123,8 +206,9 @@ antiderivative. Shipped: Taylor model arithmetic over a box (arithmetic, powers, division, and the elementary functions with rigorous remainders), range enclosure by -branch-and-bound, verified 1-D definite integrals, root absence and sign -predicates. +branch-and-bound, verified 1-D definite integrals including removable +singularities, root absence, root existence and sign predicates. Not shipped: multivariate verified quadrature (`verified_integral` is 1-D), -improper or singular integrals, and Taylor-model-based ODE enclosures. +improper integrals, integrable-but-not-removable singularities, and +Taylor-model-based ODE enclosures. diff --git a/docs/sphinx/api/solve.rst b/docs/sphinx/api/solve.rst index c9b9b46f..d7e25ffb 100644 --- a/docs/sphinx/api/solve.rst +++ b/docs/sphinx/api/solve.rst @@ -62,31 +62,84 @@ GroebnerBasis A Gröbner basis for a polynomial ideal. - .. classmethod:: compute(polys: list[Expr], vars: list[Expr], order: str = "GRevLex") -> GroebnerBasis + A basis is a sequence: ``len(gb)``, ``gb[i]`` and iteration all yield its + generators as :class:`GbPoly`. + + .. classmethod:: compute(polys: list[Expr], vars: list[Expr], order: str = "lex") -> GroebnerBasis Compute a Gröbner basis using the F4 algorithm. - :param order: Monomial order — ``"Lex"``, ``"GrLex"``, or ``"GRevLex"``. - Use ``"Lex"`` for elimination; ``"GRevLex"`` is generally fastest. + :param order: Monomial order — ``"lex"`` (default), ``"grlex"``, or + ``"grevlex"``. Use ``"lex"`` for elimination; ``"grevlex"`` is + generally fastest. For ``"lex"`` on a zero-dimensional ideal the + grevlex-then-FGLM strategy is used automatically. + + .. attribute:: order + + The monomial order the generators are reduced under, as a string. + + .. method:: variables() -> list[Expr] + + The variables naming exponent slots 0, 1, … of the generators. - Computes the basis using the F4 algorithm with product-criterion pruning. + .. method:: polynomials() -> list[GbPoly] - .. method:: reduce(expr: Expr) -> Expr + The generators. Equivalent to ``list(gb)``. - Reduce *expr* modulo the ideal. + .. method:: to_exprs(vars: list[Expr] | None = None) -> list[Expr] - .. method:: contains(expr: Expr) -> bool + The generators as expressions, each meaning ``g = 0``. This is the read + path for elimination results. + + .. method:: reduce(p: GbPoly | Expr) -> GbPoly + + Reduce *p* modulo the ideal and return the remainder. Call + :meth:`GbPoly.to_expr` on it to read it back; it is zero exactly when + :meth:`contains` is true. + + .. method:: contains(p: GbPoly | Expr) -> bool Test ideal membership. .. method:: eliminate(vars: list[Expr]) -> GroebnerBasis - Compute the elimination ideal by removing generators that involve - the specified variables. + The elimination ideal ``I ∩ k[remaining vars]``: drops every generator + whose support mentions one of *vars*. Under a ``"lex"`` basis with the + eliminated variables ordered **first**, what is left is a Gröbner basis + for that ideal. Useful for implicitization of parametric curves and surfaces. .. class:: GbPoly - A polynomial element of a Gröbner basis computation, with rational - coefficients represented as FLINT rationals. + A sparse multivariate polynomial over ℚ, as used by the Gröbner machinery. + It stores exponent *vectors*, so reading one back needs the variable list + its slots refer to — polynomials Alkahest hands out carry that list. + + .. attribute:: is_zero + .. attribute:: n_vars + .. attribute:: n_terms + + .. method:: variables() -> list[Expr] + + The variables naming this polynomial's exponent slots, in order. + + .. method:: terms() -> list[tuple[tuple[int, ...], int | Fraction]] + + ``(exponents, coefficient)`` pairs, in ascending exponent-vector order. + Coefficients are exact. + + .. method:: to_expr(vars: list[Expr] | None = None) -> Expr + + Convert back to an :class:`Expr`. Defaults to :meth:`variables`; raises + ``ValueError`` if *vars* names fewer variables than the polynomial uses. + +.. function:: expr_to_gbpoly(expr: Expr, vars: list[Expr]) -> GbPoly + + Convert a polynomial expression into the :class:`GbPoly` representation — + the inverse of :meth:`GbPoly.to_expr`. Exponent slot ``i`` refers to + ``vars[i]``, and the result remembers *vars*, so it can be passed straight + to :meth:`GroebnerBasis.reduce`, :meth:`GroebnerBasis.contains` or + ``GroebnerBasis.compute_raw``. + + Raises ``ValueError`` if *expr* is not polynomial in *vars*. diff --git a/examples/acausal_and_laplace.py b/examples/acausal_and_laplace.py index d8858723..a7c5403e 100644 --- a/examples/acausal_and_laplace.py +++ b/examples/acausal_and_laplace.py @@ -60,8 +60,8 @@ def part1_acausal_rc_circuit() -> None: dae = sys.flatten(t) - print(f"\nFlattened DAE: {dae.n_equations()} equations, " - f"{dae.n_variables()} variables") + print(f"\nFlattened DAE: {dae.n_equations} equations, " + f"{dae.n_variables} variables") print(dae) diff --git a/examples/library_tour.ipynb b/examples/library_tour.ipynb index 7865602d..10ff7b8d 100644 --- a/examples/library_tour.ipynb +++ b/examples/library_tour.ipynb @@ -231,7 +231,7 @@ "# ── UniPoly: dense univariate ─────────────────────────────────────────────\n", "p = ak.UniPoly.from_symbolic(x**3 + pool.integer(-2)*x + pool.integer(1), x)\n", "print(\"UniPoly p =\", p)\n", - "print(\"degree =\", p.degree())\n", + "print(\"degree =\", p.degree)\n", "print(\"coefficients =\", p.coefficients()) # constant-first\n", "\n", "# GCD\n", @@ -247,7 +247,7 @@ "print()\n", "mp = ak.MultiPoly.from_symbolic(x**2 * y + x * y**2, [x, y])\n", "print(\"MultiPoly =\", mp)\n", - "print(\"total_degree =\", mp.total_degree())\n", + "print(\"total_degree =\", mp.total_degree)\n", "\n", "# ── RationalFunction: automatic GCD normalization ─────────────────────────\n", "print()\n", diff --git a/examples/ode_modeling.py b/examples/ode_modeling.py index 4548f973..08e8ada2 100644 --- a/examples/ode_modeling.py +++ b/examples/ode_modeling.py @@ -39,7 +39,7 @@ ode_decay = ODE.new([y], [pool.integer(-1) * y], t) print(f"ODE : {ode_decay}") -print(f"order : {ode_decay.order()}") +print(f"order : {ode_decay.order}") print(f"autonomous : {ode_decay.is_autonomous()}") print(f"state_vars : {ode_decay.state_vars()}") print(f"rhs : {ode_decay.rhs()}") @@ -66,7 +66,7 @@ print(f"Lowered ODE:") print(f" state_vars : {ode_harmonic.state_vars()}") print(f" rhs : {ode_harmonic.rhs()}") -print(f" order : {ode_harmonic.order()}") +print(f" order : {ode_harmonic.order}") # --------------------------------------------------------------------------- # 3. Parametric ODE + forward sensitivity analysis diff --git a/examples/phase3_polynomials.md b/examples/phase3_polynomials.md index 4e92703a..5974ca8d 100644 --- a/examples/phase3_polynomials.md +++ b/examples/phase3_polynomials.md @@ -69,7 +69,7 @@ p = ak.UniPoly.from_symbolic( ) print(p) # x^2+2*x+1 -print(p.degree()) # 2 +print(p.degree) # 2 print(p.coefficients()) # [1, 2, 1] ``` @@ -95,7 +95,7 @@ x2m1 = ak.UniPoly.from_symbolic( ) g = x2m1.gcd(xm1) print(g) # x-1 -print(g.degree()) # 1 +print(g.degree) # 1 ``` GCD is computed by FLINT's `fmpz_poly_gcd` and is normalised to a positive @@ -128,7 +128,7 @@ a = ak.MultiPoly.from_symbolic( pool, ) print(a) # x1^2 + x0x1 + x0^2 -print(a.total_degree()) # 2 +print(a.total_degree) # 2 ``` ### Integer content diff --git a/examples/phase3_polynomials.py b/examples/phase3_polynomials.py index f08ffdb3..3eeee6c9 100644 --- a/examples/phase3_polynomials.py +++ b/examples/phase3_polynomials.py @@ -34,7 +34,7 @@ p = ak.UniPoly.from_symbolic(expr, x, pool) print(f"p = {p}") # x^2+2*x+1 -print(f" degree : {p.degree()}") # 2 +print(f" degree : {p.degree}") # 2 print(f" coeffs : {p.coefficients()}") # [1, 2, 1] (ascending degree) # q = x + 1 @@ -62,7 +62,7 @@ pool.add([pool.pow(x, pool.integer(2)), pool.integer(-1)]), x, pool ) g = x2m1.gcd(xm1) -print(f"\ngcd(x^2-1, x-1) = {g} (degree {g.degree()})") +print(f"\ngcd(x^2-1, x-1) = {g} (degree {g.degree})") # --------------------------------------------------------------------------- # MultiPoly demo @@ -84,7 +84,7 @@ ]) a = ak.MultiPoly.from_symbolic(expr_a, [x2, y2], pool2) print(f"a = x^2 + xy + y^2 = {a}") -print(f" total_degree : {a.total_degree()}") # 2 +print(f" total_degree : {a.total_degree}") # 2 print(f" integer_content : {a.integer_content()}") # 1 # 6x + 4 → content = 2 diff --git a/examples/polynomials.py b/examples/polynomials.py index 3d3f5919..c423106b 100644 --- a/examples/polynomials.py +++ b/examples/polynomials.py @@ -20,7 +20,7 @@ def main(): expr = x ** 3 + pool.integer(-2) * x + pool.integer(1) p = UniPoly.from_symbolic(expr, x) print(f"p = {p}") - print(f"degree = {p.degree()}") + print(f"degree = {p.degree}") print(f"coeffs = {p.coefficients()}") # GCD: gcd(x^2 - 1, x - 1) = x - 1 @@ -46,7 +46,7 @@ def main(): expr2 = x2y_expr + xy2_expr + pool.integer(-1) mp = MultiPoly.from_symbolic(expr2, [x, y]) print(f"mp = {mp}") - print(f"total_degree = {mp.total_degree()}") + print(f"total_degree = {mp.total_degree}") print(f"int_content = {mp.integer_content()}") # Bivariate addition diff --git a/pyproject.toml b/pyproject.toml index 93f8e59b..23f6b6ea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "alkahest" -version = "3.8.0" +version = "3.9.0" description = "A high-performance computer algebra system for Python" readme = "README.md" license = { text = "Apache-2.0" } diff --git a/python/alkahest/__init__.py b/python/alkahest/__init__.py index 03054e27..bd54fb78 100644 --- a/python/alkahest/__init__.py +++ b/python/alkahest/__init__.py @@ -314,6 +314,7 @@ RosenfeldGroebnerResult, dae_index_reduce, diophantine, + expr_to_gbpoly, primary_decomposition, radical, rosenfeld_groebner, @@ -1843,6 +1844,13 @@ def solve_side_conditions(*_args, **_kwargs): "See alkahest.solve.__doc__ for details." ) + def expr_to_gbpoly(*_args, **_kwargs): + """Expr → GbPoly conversion (groebner feature missing from this build).""" + raise ImportError( + "alkahest.expr_to_gbpoly is unavailable — groebner feature missing. " + "See alkahest.solve.__doc__ for details." + ) + class GroebnerBasis: """Gröbner basis type (groebner feature missing from this build). @@ -2253,6 +2261,8 @@ def wrapper(*args, **kwargs): "eval_expr", "evaluate", "exp", + # V1-16 — Expr ↔ GbPoly round trip + "expr_to_gbpoly", "factor_univariate_mod_p", "flatten_exprs", "floor", diff --git a/tests/test_accessor_convention.py b/tests/test_accessor_convention.py new file mode 100644 index 00000000..51e78f2f --- /dev/null +++ b/tests/test_accessor_convention.py @@ -0,0 +1,208 @@ +"""The property-vs-method convention for zero-argument accessors. + +The rule (see CONTRIBUTING.md § "Accessors: property or method?"): + + A zero-argument, O(1), non-allocating accessor that returns a scalar or a + flag is a **property**. Anything that returns a collection, allocates, or + does real work is a **method**. + +Two guards live here: + +``test_converted_accessors_are_properties`` + Pins the accessors converted in 3.9.0. Runs against the *installed* + extension, so it also catches a wheel built from stale sources. + +``test_no_zero_arg_scalar_accessor_is_a_method`` + Scans ``alkahest-py/src/lib.rs`` and fails on any *new* zero-argument + scalar/flag ``#[pymethods]`` entry that is not a ``#[getter]``. A runtime + version of this check is not possible: PyO3 does not expose the Rust + return type, and calling every zero-argument method to inspect its result + would need a live instance of every class. The source scan is exact and + catches the mistake where it is made. +""" + +from __future__ import annotations + +import inspect +import re +from pathlib import Path + +import alkahest.alkahest as _native +import pytest + +# -------------------------------------------------------------------------- +# Part 1 — the accessors converted in 3.9.0 stay properties +# -------------------------------------------------------------------------- + +CONVERTED_ACCESSORS = [ + ("UniPoly", "degree"), + ("UniPoly", "is_zero"), + ("MultiPoly", "is_zero"), + ("MultiPoly", "total_degree"), + ("RationalFunction", "is_zero"), + ("OdeTrajectory", "t_final"), + ("ODE", "order"), + ("DAE", "n_equations"), + ("DAE", "n_variables"), + ("HybridODE", "n_events"), + ("Component", "n_equations"), + ("Component", "n_ports"), + ("ArbBall", "is_exact"), + ("GbPoly", "is_zero"), + ("GbPoly", "n_vars"), + ("MultiPolyFp", "is_zero"), + ("MultiPolyFp", "total_degree"), +] + +# Accessors that were already properties before 3.9.0 and must stay that way — +# the convention is only useful if it holds in both directions. +PRE_EXISTING_PROPERTIES = [ + ("Enclosure", "lower"), + ("Enclosure", "upper"), + ("Enclosure", "width"), + ("Enclosure", "subdivisions"), + ("Matrix", "rows"), + ("Matrix", "cols"), + ("RegularChain", "n_vars"), + ("RosenfeldGroebnerResult", "consistent"), + ("ArbBall", "mid"), + ("ArbBall", "rad"), + ("EvaluationResult", "is_enclosure"), +] + + +def _descriptor(class_name: str, attr: str): + cls = getattr(_native, class_name, None) + if cls is None: + pytest.skip(f"{class_name} is not in this build (feature-gated)") + desc = inspect.getattr_static(cls, attr, None) + assert desc is not None, f"{class_name}.{attr} does not exist" + return desc + + +@pytest.mark.parametrize(("class_name", "attr"), CONVERTED_ACCESSORS) +def test_converted_accessors_are_properties(class_name: str, attr: str) -> None: + desc = _descriptor(class_name, attr) + assert not callable(desc), ( + f"{class_name}.{attr} is a method; the convention makes zero-argument " + f"scalar/flag accessors properties. Note that `x.{attr}` in a boolean " + f"or format context is a *silent* bug when it is a bound method." + ) + assert type(desc).__name__ == "getset_descriptor" + + +@pytest.mark.parametrize(("class_name", "attr"), PRE_EXISTING_PROPERTIES) +def test_pre_existing_properties_stay_properties(class_name: str, attr: str) -> None: + assert not callable(_descriptor(class_name, attr)) + + +def test_a_converted_accessor_actually_reads_as_a_value() -> None: + """End-to-end: the value is the scalar, not a bound method.""" + pool = _native.ExprPool() + x = pool.symbol("x") + p = _native.UniPoly.from_coefficients([-1, 0, 1], x) + assert p.degree == 2 + assert p.is_zero is False + # The failure mode this guards against: a bound method is always truthy. + assert not p.is_zero + + +# -------------------------------------------------------------------------- +# Part 2 — no *new* zero-argument scalar accessor may be a method +# -------------------------------------------------------------------------- + +_LIB_RS = Path(__file__).resolve().parent.parent / "alkahest-py" / "src" / "lib.rs" + +_SCALAR_RETURN = re.compile( + r"->\s*(?:PyResult<)?\s*" + r"(?:usize|u8|u16|u32|u64|i8|i16|i32|i64|f32|f64|bool" + r"|Option<\s*(?:usize|u32|u64|i64|f64|bool)\s*>)\s*>?\s*$" +) + +# Zero-argument scalar-returning methods that are *correctly* methods because +# they do real work, and so are exempt. Adding an entry here is a deliberate +# claim that the call is not O(1); say why. +REAL_WORK_EXEMPTIONS = { + # Gaussian elimination over the symbolic entries. + ("PyMatrix", "rank"), + # Walks every RHS expression looking for the time variable. + ("PyODE", "is_autonomous"), + # Re-runs the SOS identity check in exact arithmetic. + ("PyPositivityCertificate", "verify"), +} + + +def _pymethods_fns(source: str): + """Yield ``(class, fn_name, attrs, signature)`` for every ``#[pymethods]`` fn.""" + lines = source.split("\n") + n = len(lines) + i = 0 + while i < n: + if lines[i].strip() != "#[pymethods]": + i += 1 + continue + j = i + 1 + while j < n and not re.search(r"\bimpl\b", lines[j]): + j += 1 + m = re.search(r"impl\s+(?:\w+\s+for\s+)?([A-Za-z0-9_]+)", lines[j]) + cls = m.group(1) if m else "?" + depth = 0 + started = False + k = j + attrs: list[str] = [] + while k < n: + line = lines[k] + if not started: + depth += line.count("{") - line.count("}") + if "{" in line: + started = True + k += 1 + continue + stripped = line.strip() + fn = re.match(r"(?:pub\s+)?fn\s+([A-Za-z0-9_]+)\s*\(", stripped) + if fn and depth == 1: + sig = stripped + kk = k + while "{" not in sig and ";" not in sig and kk + 1 < n: + kk += 1 + sig += " " + lines[kk].strip() + yield cls, fn.group(1), list(attrs), re.sub(r"\s+", " ", sig.split("{")[0]) + attrs = [] + elif stripped.startswith("#["): + attrs.append(stripped) + elif stripped.startswith("//") or not stripped: + pass + elif depth == 1: + attrs = [] + depth += line.count("{") - line.count("}") + if depth <= 0: + break + k += 1 + i = k + + +@pytest.mark.skipif(not _LIB_RS.exists(), reason="running against an installed wheel") +def test_no_zero_arg_scalar_accessor_is_a_method() -> None: + offenders = [] + for cls, name, attrs, sig in _pymethods_fns(_LIB_RS.read_text()): + if name.startswith("__") or any("getter" in a for a in attrs): + continue + if any(x in a for a in attrs for x in ("staticmethod", "classmethod", "#[new]")): + continue + if (cls, name) in REAL_WORK_EXEMPTIONS: + continue + args = re.match(r".*?fn\s+\w+\s*\((.*)\)\s*(?:->.*)?$", sig) + if not args: + continue + rest = re.sub(r"^\s*(?:&\s*self|slf\s*:\s*[^,]+|self)\s*,?", "", args.group(1)).strip() + rest = re.sub(r"^\s*_?py\s*:\s*Python<[^>]*>\s*,?", "", rest).strip() + if rest.strip(","): + continue + if _SCALAR_RETURN.search(sig): + offenders.append(f"{cls}::{name} — {sig.strip()}") + + assert not offenders, ( + "These zero-argument scalar/flag accessors are exposed as methods. " + "Add `#[getter]` (see CONTRIBUTING.md § 'Accessors: property or " + "method?'), or add them to REAL_WORK_EXEMPTIONS with a reason:\n " + "\n ".join(offenders) + ) diff --git a/tests/test_api.py b/tests/test_api.py index 541fa401..f2491688 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -311,11 +311,11 @@ def test_from_symbolic_quadratic(self): def test_degree(self): p = self._make([0, 0, 1]) - assert p.degree() == 2 + assert p.degree == 2 def test_is_zero(self): p = self._make([0]) - assert p.is_zero() + assert p.is_zero def test_add(self): p = self._make([1, 2]) @@ -363,12 +363,12 @@ def test_from_symbolic_bivariate(self): x, y = self.x, self.y expr = x * y p = MultiPoly.from_symbolic(expr, [x, y]) - assert not p.is_zero() - assert p.total_degree() == 2 + assert not p.is_zero + assert p.total_degree == 2 def test_is_zero(self): p = MultiPoly.from_symbolic(self.pool.integer(0), [self.x]) - assert p.is_zero() + assert p.is_zero def test_repr_str(self): x, y = self.x, self.y diff --git a/tests/test_api_consistency_rough_edges.py b/tests/test_api_consistency_rough_edges.py index fe8ce4d2..41e6f259 100644 --- a/tests/test_api_consistency_rough_edges.py +++ b/tests/test_api_consistency_rough_edges.py @@ -15,7 +15,7 @@ def pool_x(): def test_unipoly_from_coefficients_accepts_python_ints(pool_x): _pool, x = pool_x poly = ak.UniPoly.from_coefficients([-1, 0, 1], x) - assert poly.degree() == 2 + assert poly.degree == 2 assert poly.coefficients() == [-1, 0, 1] @@ -36,7 +36,7 @@ def test_multipoly_from_symbolic_infers_vars(pool_x): pool, x = pool_x y = pool.symbol("y") mp = ak.MultiPoly.from_symbolic(x**2 * y + x) - assert not mp.is_zero() + assert not mp.is_zero def test_error_message_includes_stable_code(pool_x): diff --git a/tests/test_diffalg_v213.py b/tests/test_diffalg_v213.py index 7dc5e168..011fa917 100644 --- a/tests/test_diffalg_v213.py +++ b/tests/test_diffalg_v213.py @@ -44,4 +44,4 @@ def test_dae_index_reduce_prefers_pantelides(): assert out.used_pantelides is True assert out.used_rosenfeld_groebner is False reduced = alkahest.pantelides(dae) - assert out.dae().n_equations() == reduced.n_equations() + assert out.dae().n_equations == reduced.n_equations diff --git a/tests/test_docs_ode_dae.py b/tests/test_docs_ode_dae.py new file mode 100644 index 00000000..30839f5c --- /dev/null +++ b/tests/test_docs_ode_dae.py @@ -0,0 +1,245 @@ +"""Every code block in ``docs/mdbook/src/ode-dae.md``, executed. + +The page previously documented keyword constructors (``ODE(state=..., ...)``, +``DAE(equations=..., ...)``) and a one-argument ``lower_to_first_order`` that +never existed, plus ``pantelides(...).differentiated``, which the returned +object does not have. Following it failed on the first line. These tests are +the guard: each one is a doc snippet run verbatim, with the values the page +claims asserted. +""" + +import alkahest +import pytest +from alkahest import ( + DAE, + ODE, + AcausalSystem, + Event, + ExprPool, + HybridODE, + adjoint_system, + capacitor, + lower_to_first_order, + pantelides, + resistor, + sensitivity_system, + voltage_source, +) + + +def test_ode_constructor(): + pool = ExprPool() + t = pool.symbol("t") + x = pool.symbol("x") + v = pool.symbol("v") + + ode = ODE.new([x, v], [v, pool.integer(-1) * x], t) + + assert ode.order == 2 + assert [str(s) for s in ode.state_vars()] == ["x", "v"] + assert [str(r) for r in ode.rhs()] == ["v", "(x * -1)"] + assert ode.is_autonomous() is True + + # ODE(...) and ODE.new(...) agree. + assert str(ODE([x, v], [v, pool.integer(-1) * x], t)) == str(ode) + + +def test_ode_with_ic_returns_a_new_ode(): + pool = ExprPool() + t, x, v = pool.symbol("t"), pool.symbol("x"), pool.symbol("v") + ode = ODE.new([x, v], [v, pool.integer(-1) * x], t) + + ode_with_ic = ode.with_ic(x, pool.integer(1)).with_ic(v, pool.integer(0)) + + assert isinstance(ode_with_ic, ODE) + assert ode_with_ic is not ode + + +def test_lower_to_first_order_takes_four_arguments(): + pool = ExprPool() + t, x = pool.symbol("t"), pool.symbol("x") + + ode = lower_to_first_order(x, pool.integer(-4) * x, 2, t) + + assert [str(s) for s in ode.state_vars()] == ["x", "x_1"] + assert [str(r) for r in ode.rhs()] == ["x_1", "(x * -4)"] + + # The documented one-argument form does not exist. + with pytest.raises(TypeError): + lower_to_first_order(ode) + + +def pendulum_dae(): + """The index-3 Cartesian pendulum from the DAE section of the page.""" + pool = ExprPool() + t = pool.symbol("t") + x, y, u, w = (pool.symbol(n) for n in ("x", "y", "u", "w")) + lam = pool.symbol("lam") + dx, dy, du, dw = (pool.symbol(n) for n in ("dx/dt", "dy/dt", "du/dt", "dw/dt")) + one, two = pool.integer(1), pool.integer(2) + + dae = DAE.new( + [dx - u, dy - w, du + lam * x, dw + lam * y + one, x**two + y**two - one], + [x, y, u, w, lam], + [dx, dy, du, dw], + t, + ) + return pool, t, dae + + +def test_dae_new_and_read_back(): + _pool, t, dae = pendulum_dae() + + assert dae.n_equations == 5 + assert dae.n_variables == 5 + assert len(dae.equations()) == 5 + assert [str(d) for d in dae.derivatives()] == ["dx/dt", "dy/dt", "du/dt", "dw/dt"] + assert str(dae.time_var) == "t" + assert dae.index is None + + # There is no keyword constructor. + with pytest.raises(TypeError): + DAE(equations=dae.equations(), variables=dae.variables(), independent=t) + + +def test_pantelides_returns_a_reduced_dae(): + _, _, dae = pendulum_dae() + + reduced = pantelides(dae) + + assert isinstance(reduced, DAE) + assert reduced.index == 1 + assert reduced.n_equations == 6 + # The appended equation is the differentiated constraint 2x·x' + 2y·y' = 0. + assert str(reduced.equations()[-1]) == "((x * dx/dt * 2) + (y * dy/dt * 2))" + # Prolongation introduced second-order jets. + assert "ddx/dt/dt" in [str(d) for d in reduced.derivatives()] + + # The page used to claim a `.differentiated` attribute; it never existed. + assert not hasattr(reduced, "differentiated") + + +def test_pantelides_index_zero_when_already_matched(): + pool = ExprPool() + t, x, dx = pool.symbol("t"), pool.symbol("x"), pool.symbol("dx/dt") + reduced = pantelides(DAE.new([dx - x], [x], [dx], t)) + assert reduced.index == 0 + + +@pytest.mark.skipif( + not hasattr(alkahest, "rosenfeld_groebner"), + reason="native module built without groebner feature", +) +def test_rosenfeld_groebner_snippet(): + from alkahest import rosenfeld_groebner + + pool = ExprPool() + t, x, dx = pool.symbol("t"), pool.symbol("x"), pool.symbol("dx/dt") + + dae = DAE.new([dx - x], [x], [dx], t) + result = rosenfeld_groebner(dae, max_prolong_rounds=1) + + assert result.consistent is True + assert result.truncated is True + assert [str(v) for v in result.variables()] == ["t", "x", "dx/dt", "ddx/dt/dt"] + assert [str(e) for e in result.final_basis().to_exprs()] == [ + "(x + (-1 * ddx/dt/dt))", + "(dx/dt + (-1 * ddx/dt/dt))", + ] + + +@pytest.mark.skipif( + not hasattr(alkahest, "rosenfeld_groebner"), + reason="native module built without groebner feature", +) +def test_rosenfeld_groebner_inconsistent_has_no_final_basis(): + from alkahest import rosenfeld_groebner + + pool = ExprPool() + t, y, dy = pool.symbol("t"), pool.symbol("y"), pool.symbol("dy/dt") + dae = DAE.new([dy - y, dy - y - pool.integer(1)], [y], [dy], t) + + result = rosenfeld_groebner(dae, max_prolong_rounds=1) + + assert result.consistent is False + assert result.final_basis() is None + + +def test_sensitivity_and_adjoint_snippet(): + pool = ExprPool() + t, y, k = pool.symbol("t"), pool.symbol("y"), pool.symbol("k") + + ode = ODE.new([y], [pool.integer(-1) * k * y], t) + + sens = sensitivity_system(ode, [k]) + assert sens.original_dim == 1 + assert sens.n_params == 1 + assert [str(s) for s in sens.extended_ode.state_vars()] == ["y", "dS_y_k"] + + adj = adjoint_system(ode, [pool.integer(2) * y]) + assert [str(s) for s in adj.state_vars()] == ["lambda_y"] + assert [str(r) for r in adj.rhs()] == ["(k * lambda_y)"] + + +def test_acausal_snippet(): + pool = ExprPool() + t = pool.symbol("t") + + src = voltage_source("V1", pool.symbol("Vs"))["component"] + res = resistor("R1", pool.symbol("R"))["component"] + cap = capacitor("C1", pool.symbol("C"))["component"] + + circuit = AcausalSystem(pool) + circuit.add_component(src) + circuit.add_component(res) + circuit.add_component(cap) + + circuit.connect(src.port("V1.p"), res.port("R1.p")) + circuit.connect(res.port("R1.n"), cap.port("C1.p")) + circuit.connect(cap.port("C1.n"), src.port("V1.n")) + + dae = circuit.flatten(t) + assert dae.n_equations == 10 + + +def test_laplace_snippet(): + from alkahest.experimental import inverse_laplace_transform, laplace_transform + + pool = ExprPool() + s, t = pool.symbol("s"), pool.symbol("t") + + F = laplace_transform(pool.integer(1), t, s) + f = inverse_laplace_transform(F, s, t) + + assert str(F) == "s^-1" + assert str(f) == "(1)" + + +def test_hybrid_snippet(): + pool = ExprPool() + t, x, v = pool.symbol("t"), pool.symbol("x"), pool.symbol("v") + + base_ode = ODE.new([x, v], [v, pool.integer(-1)], t) + bounce = Event.new("bounce", x, [(v, pool.integer(-1) * v)]) + + hybrid = HybridODE.new(base_ode).add_event(bounce) + assert hybrid.n_events == 1 + assert [str(g) for g in hybrid.guards()] == ["x"] + + +def test_documented_apis_carry_docstrings(): + """The page was the only documentation because none of these had one.""" + for obj in ( + DAE.new, + DAE.equations, + DAE.index, + ODE.new, + pantelides, + lower_to_first_order, + ): + assert obj.__doc__, obj + + if hasattr(alkahest, "rosenfeld_groebner"): + assert alkahest.rosenfeld_groebner.__doc__ + assert alkahest.RosenfeldGroebnerResult.final_basis.__doc__ + assert alkahest.RosenfeldGroebnerResult.variables.__doc__ diff --git a/tests/test_gbpoly_read_path.py b/tests/test_gbpoly_read_path.py new file mode 100644 index 00000000..7f1ed18a --- /dev/null +++ b/tests/test_gbpoly_read_path.py @@ -0,0 +1,277 @@ +"""A computed basis must be readable — `GbPoly`/`GroebnerBasis` round trips. + +Before this, `GbPoly` exposed only `is_zero`/`n_vars` and `GroebnerBasis` only +its constructors plus `reduce`/`contains`, so nothing that returned a basis — +`rosenfeld_groebner`, `triangularize`, `primary_decomposition`, parametric +`solve` — could be inspected at all. `expr_to_gbpoly` was named in +`compute_raw`'s own docstring but was not exported, so `reduce()` could not +even be called: a caller had no way to build its argument. +""" + +import alkahest +import pytest +from alkahest import DAE, ExprPool + +pytestmark = pytest.mark.skipif( + not hasattr(alkahest, "GroebnerBasis"), + reason="native module built without groebner feature", +) + + +@pytest.fixture +def circle_line(): + """`x**2 + y**2 = 1` intersected with `x = y`, and its lex basis.""" + pool = ExprPool() + x, y = pool.symbol("x"), pool.symbol("y") + two, one = pool.integer(2), pool.integer(1) + polys = [x**two + y**two - one, x - y] + gb = alkahest.GroebnerBasis.compute(polys, [x, y]) + return pool, x, y, gb + + +# --------------------------------------------------------------------------- +# GroebnerBasis is a readable sequence +# --------------------------------------------------------------------------- + + +def test_basis_is_a_sequence_of_gbpoly(circle_line): + _, _, _, gb = circle_line + + assert len(gb) == 2 + assert len(list(gb)) == len(gb) + assert all(isinstance(g, alkahest.GbPoly) for g in gb) + assert isinstance(gb[0], alkahest.GbPoly) + assert gb.polynomials()[0].terms() == gb[0].terms() + + +def test_basis_indexing_bounds(circle_line): + _, _, _, gb = circle_line + + assert gb[-1].terms() == gb[len(gb) - 1].terms() + with pytest.raises(IndexError): + gb[len(gb)] + with pytest.raises(IndexError): + gb[-len(gb) - 1] + + +def test_basis_reports_its_order_and_variables(circle_line): + _, x, y, gb = circle_line + + assert gb.order == "lex" + assert [str(v) for v in gb.variables()] == ["x", "y"] + assert [str(v) for v in gb[0].variables()] == ["x", "y"] + + grevlex = alkahest.GroebnerBasis.compute([x - y], [x, y], "grevlex") + assert grevlex.order == "grevlex" + + +# --------------------------------------------------------------------------- +# GbPoly ↔ Expr round trip +# --------------------------------------------------------------------------- + + +def test_expr_to_gbpoly_round_trips(circle_line): + pool, x, y, _ = circle_line + expr = x ** pool.integer(2) + y ** pool.integer(2) - pool.integer(1) + + p = alkahest.expr_to_gbpoly(expr, [x, y]) + + assert p.n_vars == 2 + assert p.n_terms == 3 + # `to_expr` rebuilds a flat sum rather than the caller's exact tree, so the + # round trip is closed on the canonical side: Expr → GbPoly → Expr → GbPoly. + assert alkahest.expr_to_gbpoly(p.to_expr(), [x, y]).terms() == p.terms() + assert str(p.to_expr()) == "(x^2 + y^2 + -1)" + + +def test_gbpoly_terms_are_exact(circle_line): + pool, x, y, _ = circle_line + # x**2 - (3/2)*y + expr = x ** pool.integer(2) - pool.rational(3, 2) * y + + terms = dict(alkahest.expr_to_gbpoly(expr, [x, y]).terms()) + + from fractions import Fraction + + assert terms[(2, 0)] == 1 + assert terms[(0, 1)] == Fraction(-3, 2) + + +def test_basis_generators_convert_back_to_expr(circle_line): + _pool, _x, _y, gb = circle_line + + exprs = gb.to_exprs() + + assert len(exprs) == len(gb) + assert [str(e) for e in exprs] == [str(g.to_expr()) for g in gb] + # The lex basis eliminates x: one generator is univariate in y, and it is + # the relation 2y**2 = 1 that the intersection actually satisfies. + assert "(y^2 + -1/2)" in [str(e) for e in exprs] + + +def test_round_trip_expr_to_gbpoly_to_expr_is_in_the_ideal(circle_line): + """Read a generator back out, feed it in again, and it is still in `I`.""" + _pool, x, y, gb = circle_line + + for g in gb: + recovered = g.to_expr() + assert gb.contains(recovered) + assert alkahest.expr_to_gbpoly(recovered, [x, y]).n_terms == g.n_terms + + +def test_to_expr_accepts_an_explicit_variable_list(circle_line): + pool, _x, _y, gb = circle_line + a, b = pool.symbol("a"), pool.symbol("b") + + # Same exponent vectors, different names. + renamed = gb[0].to_expr([a, b]) + assert "y" not in str(renamed) + + +def test_to_expr_refuses_a_short_variable_list(circle_line): + _pool, x, _y, gb = circle_line + + with pytest.raises(ValueError, match=r"only 1 were named|more variables"): + gb.to_exprs([x]) + + +def test_zero_polynomial_converts_to_zero(circle_line): + _pool, x, y, gb = circle_line + + # Anything in the ideal reduces to zero. + remainder = gb.reduce(x - y) + assert remainder.is_zero + assert str(remainder.to_expr()) == "0" + + +# --------------------------------------------------------------------------- +# reduce() / contains() are usable on a basis the caller has in hand +# --------------------------------------------------------------------------- + + +def test_reduce_accepts_expr_and_gbpoly(circle_line): + pool, x, y, gb = circle_line + cube = x ** pool.integer(3) + + from_expr = gb.reduce(cube) + from_poly = gb.reduce(alkahest.expr_to_gbpoly(cube, [x, y])) + + assert from_expr.terms() == from_poly.terms() + # x**3 = y*y**2 = y/2 on the intersection. + assert str(from_expr.to_expr()) == "(y * 1/2)" + + +def test_compute_raw_keeps_the_variable_context(circle_line): + pool, x, y, _ = circle_line + p = alkahest.expr_to_gbpoly(x ** pool.integer(2) - y, [x, y]) + + gb = alkahest.GroebnerBasis.compute_raw([p]) + + assert [str(e) for e in gb.to_exprs()] == [str(p.to_expr())] + assert gb.contains(x ** pool.integer(2) - y) + + +def test_empty_vars_falls_back_to_the_stored_context(circle_line): + """An empty `vars` list means "use what you know", not "name nothing".""" + _, _, _, gb = circle_line + + assert [str(e) for e in gb.to_exprs([])] == [str(e) for e in gb.to_exprs()] + + +# --------------------------------------------------------------------------- +# eliminate() — documented in docs/mdbook/src/solving.md, previously unbound +# --------------------------------------------------------------------------- + + +def test_eliminate_implicitizes_a_parametric_curve(): + pool = ExprPool() + t, x, y = pool.symbol("t"), pool.symbol("x"), pool.symbol("y") + + # (t, t**2), with the parameter ordered first so lex eliminates it. + gb = alkahest.GroebnerBasis.compute([x - t, y - t ** pool.integer(2)], [t, x, y]) + implicit = gb.eliminate([t]) + + assert [str(e) for e in implicit.to_exprs()] == ["((y * -1) + x^2)"] + assert [str(v) for v in implicit.variables()] == ["t", "x", "y"] + + +def test_eliminate_rejects_an_unknown_variable(circle_line): + pool, _x, _y, gb = circle_line + + with pytest.raises(ValueError, match="not written over"): + gb.eliminate([pool.symbol("q")]) + + +# --------------------------------------------------------------------------- +# The three APIs issue #11 called write-only +# --------------------------------------------------------------------------- + + +def test_rosenfeld_groebner_basis_is_readable(): + pool = ExprPool() + t, x, dx = pool.symbol("t"), pool.symbol("x"), pool.symbol("dx/dt") + + result = alkahest.rosenfeld_groebner(DAE.new([dx - x], [x], [dx], t), max_prolong_rounds=1) + basis = result.final_basis() + + assert len(basis) > 0 + assert [str(v) for v in result.variables()] == [str(v) for v in basis.variables()] + exprs = basis.to_exprs() + assert len(exprs) == len(basis) + # x' = x prolongs to x = x'' and x' = x''. + assert [str(e) for e in exprs] == [ + "(x + (-1 * ddx/dt/dt))", + "(dx/dt + (-1 * ddx/dt/dt))", + ] + + +def test_rosenfeld_working_dae_variables_can_be_named(): + """The prolonged system reports extra variables; they now have names.""" + pool = ExprPool() + t, x, dx = pool.symbol("t"), pool.symbol("x"), pool.symbol("dx/dt") + + result = alkahest.rosenfeld_groebner(DAE.new([dx - x], [x], [dx], t), max_prolong_rounds=1) + working = result.working_dae() + + assert len(working.variables()) == working.n_variables + assert [str(v) for v in working.derivatives()] == ["dx/dt", "ddx/dt/dt"] + + +def test_triangularize_chain_is_readable(): + pool = ExprPool() + x, y = pool.symbol("x"), pool.symbol("y") + two, one = pool.integer(2), pool.integer(1) + + chains = alkahest.triangularize([x**two + y**two - one, x - y], [x, y]) + + assert chains + for chain in chains: + assert [str(v) for v in chain.variables()] == ["x", "y"] + assert [str(e) for e in chain.to_exprs()] == [str(p.to_expr()) for p in chain.polys()] + + +def test_primary_decomposition_components_are_readable(): + pool = ExprPool() + x, y = pool.symbol("x"), pool.symbol("y") + two = pool.integer(2) + + comps = alkahest.primary_decomposition([x**two, y], [x, y]) + + assert comps + for c in comps: + assert c.primary().to_exprs() + assert c.associated_prime().to_exprs() + + +def test_parametric_solve_basis_is_readable(): + """`solve` returns a basis for a positive-dimensional ideal; read it.""" + pool = ExprPool() + x, y, z = pool.symbol("x"), pool.symbol("y"), pool.symbol("z") + + result = alkahest.solve([x ** pool.integer(2) - y * z], [x, y]) + + assert isinstance(result, alkahest.GroebnerBasis) + # Free parameters are appended after the solve variables, and the basis + # names all of them — otherwise `to_exprs()` could not run at all. + assert [str(v) for v in result.variables()] == ["x", "y", "z"] + assert [str(e) for e in result.to_exprs()] == ["(x^2 + (y * z * -1))"] diff --git a/tests/test_modular.py b/tests/test_modular.py index 6957444d..ccf3b217 100644 --- a/tests/test_modular.py +++ b/tests/test_modular.py @@ -35,21 +35,21 @@ def test_basic(self): fp = modular.reduce_mod(poly, 5) assert isinstance(fp, MultiPolyFp) assert fp.modulus == 5 - assert not fp.is_zero() + assert not fp.is_zero def test_negative_coeff(self): # -3x mod 7 → 4x poly, _pool, _x, _y = make_poly(lambda p, x, y: p.integer(-3) * x) fp = modular.reduce_mod(poly, 7) assert fp.modulus == 7 - assert not fp.is_zero() + assert not fp.is_zero def test_vanishing_term(self): # 5x + 7 mod 5 → 2 (x term vanishes) poly, _pool, _x, _y = make_poly(lambda p, x, y: p.integer(5) * x + p.integer(7)) fp = modular.reduce_mod(poly, 5) # degree falls (x^1 term is gone) — the result is just the constant 2 - assert fp.total_degree() == 0 + assert fp.total_degree == 0 def test_zero_poly(self): pool = alkahest.ExprPool() @@ -58,7 +58,7 @@ def test_zero_poly(self): zero = pool.integer(0) poly = MultiPoly.from_symbolic(zero, [x, y]) fp = modular.reduce_mod(poly, 7) - assert fp.is_zero() + assert fp.is_zero def test_invalid_modulus_composite(self): poly, *_ = make_poly(lambda p, x, y: x) @@ -89,7 +89,7 @@ def test_total_degree(self): x2 = x**2 poly = MultiPoly.from_symbolic(x2 + x + pool.integer(1), [x, y]) fp = modular.reduce_mod(poly, 5) - assert fp.total_degree() == 2 + assert fp.total_degree == 2 # --------------------------------------------------------------------------- diff --git a/tests/test_numeric_ode.py b/tests/test_numeric_ode.py index 1488088d..54636d8f 100644 --- a/tests/test_numeric_ode.py +++ b/tests/test_numeric_ode.py @@ -67,7 +67,9 @@ def test_trajectory_attributes(self): traj = ex.ode_integrate_rk4(ode, [1.0], 0.0, 0.5, h=0.1) assert hasattr(traj, "t") assert hasattr(traj, "y") - assert callable(traj.t_final) + # `t_final` is a property (a scalar read); `y_final` is a method (it + # allocates a list). See CONTRIBUTING.md "Accessors: property or method?". + assert not callable(traj.t_final) assert callable(traj.y_final) def test_trajectory_len(self): @@ -85,7 +87,7 @@ def test_t_starts_at_t_start(self): def test_t_final_approx_t_end(self): ode = make_exp_ode() traj = ex.ode_integrate_rk4(ode, [1.0], 0.0, 1.0, h=0.1) - assert traj.t_final() == pytest.approx(1.0, abs=1e-10) + assert traj.t_final == pytest.approx(1.0, abs=1e-10) def test_y_initial_matches_ic(self): ode = make_harmonic_ode() diff --git a/tests/test_poly_factor.py b/tests/test_poly_factor.py index 078f38b0..9155d9d4 100644 --- a/tests/test_poly_factor.py +++ b/tests/test_poly_factor.py @@ -14,7 +14,7 @@ def test_unipoly_factor_quadratic(): assert len(fac.factor_list()) == 2 for base, exp in fac.factor_list(): assert exp == 1 - assert base.degree() == 1 + assert base.degree == 1 assert fac.verification == { "status": "exactly_verified", "evidence": "factor_product", diff --git a/tests/test_resultant.py b/tests/test_resultant.py index 19e4e97f..b04e46b3 100644 --- a/tests/test_resultant.py +++ b/tests/test_resultant.py @@ -147,7 +147,7 @@ def test_bivariate_sanity(self): # Convert the result to UniPoly in x and check coefficients. res_poly = UniPoly.from_symbolic(res_expr, x) coeffs = res_poly.coefficients() - assert res_poly.degree() == 2, f"expected degree 2, got {res_poly.degree()}" + assert res_poly.degree == 2, f"expected degree 2, got {res_poly.degree}" # Coefficients in ascending order: [-1, 0, 2] assert coeffs[0] == -1, f"constant term should be -1, got {coeffs[0]}" assert coeffs[2] == 2, f"leading coeff should be 2, got {coeffs[2]}" @@ -218,8 +218,8 @@ def test_prs_gcd_x2m1_xm1(self): seq = subresultant_prs(p_expr, q_expr, x) last = seq[-1] last_poly = UniPoly.from_symbolic(last, x) - assert last_poly.degree() == 1, ( - f"last element should be degree 1 (matching gcd), got {last_poly.degree()}" + assert last_poly.degree == 1, ( + f"last element should be degree 1 (matching gcd), got {last_poly.degree}" ) def test_prs_coprime_terminates_at_constant(self): @@ -232,8 +232,8 @@ def test_prs_coprime_terminates_at_constant(self): assert len(seq) >= 2 last = seq[-1] last_poly = UniPoly.from_symbolic(last, x) - assert last_poly.degree() == 0, ( - f"last element of coprime PRS should be degree 0, got {last_poly.degree()}" + assert last_poly.degree == 0, ( + f"last element of coprime PRS should be degree 0, got {last_poly.degree}" ) def test_prs_consistent_with_resultant(self): diff --git a/tests/test_smoke.py b/tests/test_smoke.py index b7372d42..7efd5fc5 100644 --- a/tests/test_smoke.py +++ b/tests/test_smoke.py @@ -228,7 +228,7 @@ def test_unipoly_degree(): p = pool() x = p.symbol("x") poly = UniPoly.from_symbolic(x**3 + p.integer(-2) * x + p.integer(1), x) - assert poly.degree() == 3 + assert poly.degree == 3 def test_unipoly_gcd(): @@ -245,7 +245,7 @@ def test_multipoly_total_degree(): x = p.symbol("x") y = p.symbol("y") mp = MultiPoly.from_symbolic(x**2 * y + x * y**2, [x, y]) - assert mp.total_degree() == 3 + assert mp.total_degree == 3 def test_rational_function_normalization(): diff --git a/tests/test_sparse_interp.py b/tests/test_sparse_interp.py index 828bb56f..3eeb5936 100644 --- a/tests/test_sparse_interp.py +++ b/tests/test_sparse_interp.py @@ -402,8 +402,8 @@ def test_univariate_coprime(self): g = _mp(x + pool.integer(1), [x]) h = gcd_sparse(f, g, term_bound=2, degree_bound=2) # Primitive 1 has total_degree 0 - assert h.total_degree() == 0, f"GCD of coprime polys should be constant, got {h}" - assert not h.is_zero(), "GCD should be 1 (non-zero constant)" + assert h.total_degree == 0, f"GCD of coprime polys should be constant, got {h}" + assert not h.is_zero, "GCD should be 1 (non-zero constant)" def test_univariate_identical_primitive_part(self): """gcd(2f, 2f) = primitive_part(f).""" @@ -424,8 +424,8 @@ def test_bivariate_common_factor(self): g = _mp((x + y) * (x + pool.integer(1)), [x, y]) h = gcd_sparse(f, g, term_bound=3, degree_bound=2) ref = f.gcd(g) - assert h.total_degree() == ref.total_degree(), ( - f"sparse GCD total_degree {h.total_degree()} != FLINT {ref.total_degree()}" + assert h.total_degree == ref.total_degree, ( + f"sparse GCD total_degree {h.total_degree} != FLINT {ref.total_degree}" ) assert str(h) == str(ref), f"bivariate sparse GCD {h} != FLINT {ref}" @@ -437,8 +437,8 @@ def test_bivariate_coprime(self): f = _mp(x, [x, y]) g = _mp(y, [x, y]) h = gcd_sparse(f, g, term_bound=2, degree_bound=2) - assert h.total_degree() == 0, f"gcd(x,y) should be constant, got {h}" - assert not h.is_zero(), "GCD should be 1 (non-zero)" + assert h.total_degree == 0, f"gcd(x,y) should be constant, got {h}" + assert not h.is_zero, "GCD should be 1 (non-zero)" def test_incompatible_vars_raises(self): """Different variable lists must raise SparseGcdError.""" diff --git a/tests/test_v02.py b/tests/test_v02.py index d8515cb2..23494f82 100644 --- a/tests/test_v02.py +++ b/tests/test_v02.py @@ -49,7 +49,7 @@ def test_sub(self): r1 = self.rf(self.x, self.one) r2 = self.rf(self.x, self.one) result = r1 - r2 - assert result.is_zero() + assert result.is_zero def test_div(self): # (x/1) / (x/1) = 1 @@ -65,7 +65,7 @@ def test_mul_cancels(self): ba = self.rf(two_expr, x_expr) result = ab * ba # Should reduce to 1 or at least is not zero - assert not result.is_zero() + assert not result.is_zero # --------------------------------------------------------------------------- diff --git a/tests/test_v03.py b/tests/test_v03.py index a4d940bb..4c628e65 100644 --- a/tests/test_v03.py +++ b/tests/test_v03.py @@ -244,7 +244,7 @@ def setup_method(self): def test_create_simple(self): # dx/dt = x ode = ODE.new([self.x], [self.x], self.t) - assert ode.order() == 1 + assert ode.order == 1 def test_autonomous(self): ode = ODE.new([self.x], [self.x], self.t) @@ -265,7 +265,7 @@ def test_with_ic(self): one = self.pool.integer(1) ode2 = ode.with_ic(self.x, one) # Should not raise; order unchanged - assert ode2.order() == 1 + assert ode2.order == 1 def test_state_vars(self): ode = ODE.new([self.x], [self.x], self.t) @@ -294,17 +294,17 @@ def setup_method(self): def test_first_order_passthrough(self): ode = lower_to_first_order(self.x, self.x, 1, self.t) - assert ode.order() == 1 + assert ode.order == 1 def test_second_order_harmonic_oscillator(self): # x'' = -x neg_x = self.pool.integer(-1) * self.x ode = lower_to_first_order(self.x, neg_x, 2, self.t) - assert ode.order() == 2 + assert ode.order == 2 def test_third_order_produces_3_states(self): ode = lower_to_first_order(self.x, self.pool.integer(0), 3, self.t) - assert ode.order() == 3 + assert ode.order == 3 # ───────────────────────────────────────────────────────────────────────────── @@ -323,8 +323,8 @@ def test_create_dae(self): neg_x = self.pool.integer(-1) * self.x eq = self.dx + neg_x # dx/dt - x = 0 dae = DAE.new([eq], [self.x], [self.dx], self.t) - assert dae.n_equations() == 1 - assert dae.n_variables() == 1 + assert dae.n_equations == 1 + assert dae.n_variables == 1 def test_pantelides_ode_no_steps(self): neg_x = self.pool.integer(-1) * self.x @@ -332,7 +332,7 @@ def test_pantelides_ode_no_steps(self): dae = DAE.new([eq], [self.x], [self.dx], self.t) reduced = pantelides(dae) # An ODE needs 0 differentiation steps - assert reduced.n_equations() >= 1 + assert reduced.n_equations >= 1 def test_pantelides_returns_dae(self): neg_x = self.pool.integer(-1) * self.x @@ -363,15 +363,15 @@ def test_acausal_system_flatten(self): sys = AcausalSystem(self.pool) dae = sys.flatten(self.t) # Empty system → empty DAE - assert dae.n_equations() == 0 + assert dae.n_equations == 0 def test_component_accessible_from_dict(self): R = self.pool.symbol("R") comp = resistor("R1", R)["component"] assert isinstance(comp, Component) assert comp.name == "R1" - assert comp.n_equations() == 1 - assert comp.n_ports() == 2 + assert comp.n_equations == 1 + assert comp.n_ports == 2 port = comp.port("R1.p") assert port is not None assert port.name == "R1.p" @@ -408,8 +408,8 @@ def test_acausal_rc_circuit_add_component_and_connect(self): dae = sys.flatten(self.t) # 1 (source) + 1 (resistor) + 2 (capacitor) + 2*3 (connections) = 10 - assert dae.n_equations() == 10 - assert dae.n_variables() > 0 + assert dae.n_equations == 10 + assert dae.n_variables > 0 # ───────────────────────────────────────────────────────────────────────────── @@ -433,14 +433,14 @@ def test_sensitivity_linear(self): assert sys.original_dim == 1 assert sys.n_params == 1 # Extended ODE has y + 1 sensitivity state - assert sys.extended_ode.order() == 2 + assert sys.extended_ode.order == 2 def test_sensitivity_constant(self): # dy/dt = p, param = p → dS/dt = 1 p = self.pool.symbol("p") ode = ODE.new([self.y], [p], self.t) sys = sensitivity_system(ode, [p]) - assert sys.extended_ode.order() == 2 + assert sys.extended_ode.order == 2 # The sensitivity RHS should be 1 s_rhs = sys.extended_ode.rhs()[1] assert str(s_rhs) == "1" @@ -451,7 +451,7 @@ def test_sensitivity_two_params(self): ode = ODE.new([self.y], [rhs], self.t) sys = sensitivity_system(ode, [self.a, b]) assert sys.n_params == 2 - assert sys.extended_ode.order() == 3 + assert sys.extended_ode.order == 3 def test_adjoint_system(self): # dy/dt = -y, objective ∂g/∂y = 1 @@ -459,7 +459,7 @@ def test_adjoint_system(self): ode = ODE.new([self.y], [neg_y], self.t) one = self.pool.integer(1) adj = adjoint_system(ode, [one]) - assert adj.order() == 1 + assert adj.order == 1 # dλ/dt = λ (adjoint of -y system) lam = adj.state_vars()[0] assert str(adj.rhs()[0]) == str(lam) @@ -487,7 +487,7 @@ def _bouncing_ball(self): def test_create_hybrid(self): h = self._bouncing_ball() - assert h.n_events() == 1 + assert h.n_events == 1 def test_guards(self): h = self._bouncing_ball() @@ -513,7 +513,7 @@ def test_event_direction(self): def test_no_events(self): ode = ODE.new([self.y], [self.v], self.t) h = HybridODE.new(ode) - assert h.n_events() == 0 + assert h.n_events == 0 assert h.guards() == [] def test_multiple_events(self): @@ -521,7 +521,7 @@ def test_multiple_events(self): ev1 = Event.new("ev1", self.y, []) ev2 = Event.new("ev2", self.v, []) h = HybridODE.new(ode).add_event(ev1).add_event(ev2) - assert h.n_events() == 2 + assert h.n_events == 2 guards = h.guards() assert str(guards[0]) == str(self.y) assert str(guards[1]) == str(self.v) diff --git a/tests/test_v04.py b/tests/test_v04.py index aff00451..7bfddefb 100644 --- a/tests/test_v04.py +++ b/tests/test_v04.py @@ -190,7 +190,7 @@ def test_construction(self): def test_default_exact(self): b = ArbBall(2.0) - assert b.is_exact() + assert b.is_exact def test_contains_midpoint(self): b = ArbBall(3.0, 0.5) diff --git a/tests/test_validated_bounds.py b/tests/test_validated_bounds.py index b8f1d5fb..3fbf9f30 100644 --- a/tests/test_validated_bounds.py +++ b/tests/test_validated_bounds.py @@ -110,6 +110,126 @@ def test_verified_integral_is_an_enclosure_not_an_estimate(): assert r.lower <= 2.0 <= r.upper +# --------------------------------------------------------------------------- +# Removable singularities: the integrand is singular, the integral is not +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + ("build", "a", "b", "exact", "label"), + [ + ( + lambda pool, x: ak.log(pool.integer(1) + x) / x, + 0.0, + 1.0, + math.pi**2 / 12, + "log(1+x)/x", + ), + (lambda pool, x: ak.sin(x) / x, -1.0, 1.0, 1.8921661407343664, "sin(x)/x"), + (lambda pool, x: ak.sin(x) / x, 0.0, 1.0, 0.9460830703671832, "sin(x)/x on [0,1]"), + ( + lambda pool, x: (ak.exp(x) - pool.integer(1)) / x, + 0.0, + 1.0, + 1.3179021514544038, + "(exp(x)-1)/x", + ), + ( + lambda pool, x: (pool.integer(1) - ak.cos(x)) / x, + 0.0, + 1.0, + 0.23981174200056, + "(1-cos x)/x = Cin(1)", + ), + ], +) +def test_removable_singularity_encloses_the_exact_value(build, a, b, exact, label): + """The enclosure has to *bracket* the truth — returning one is worthless + if it is wrong. These integrands are all undefined at a point of the + interval and extend continuously across it.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + r = ak.verified_integral(build(pool, x), x, a, b) + + assert r.lower <= exact <= r.upper, f"{label}: {exact} not in [{r.lower}, {r.upper}]" + assert r.width < 1e-6, f"{label}: enclosure width {r.width}" + + +def test_removable_branch_agrees_with_ordinary_quadrature_on_the_regular_part(): + """Splitting the interval so that only one half touches the singularity + must not move the answer: a biased removable branch would show up here.""" + pool = ak.ExprPool() + x = pool.symbol("x") + f = ak.log(pool.integer(1) + x) / x + + whole = ak.verified_integral(f, x, 0.0, 1.0) + left = ak.verified_integral(f, x, 0.0, 0.25) + right = ak.verified_integral(f, x, 0.25, 1.0) + + assert whole.lower <= left.upper + right.upper + assert whole.upper >= left.lower + right.lower + + +def test_a_genuine_pole_is_still_refused(): + """`1/x` on [-1,1]: the numerator does not vanish, so nothing is removable.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + with pytest.raises(ak.ValidatedError) as excinfo: + ak.verified_integral(pool.integer(1) / x, x, -1.0, 1.0) + + assert excinfo.value.code == "E-VALIDATED-003" + + +def test_a_double_pole_with_a_simple_numerator_zero_is_still_refused(): + """`sin(x)/x**2 ~ 1/x` does not converge; the numerator's zero is only + order 1 against the denominator's order 2, and `D' = 2x` vanishes.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + with pytest.raises(ak.ValidatedError): + ak.verified_integral(ak.sin(x) / (x * x), x, -1.0, 1.0) + + +def test_a_second_order_removable_singularity_is_refused_not_guessed(): + """`(1-cos x)/x**2` really is removable (it tends to 1/2), but the proof + needs a *second*-order argument: `D' = 2x` vanishes at 0, so Cauchy's mean + value theorem does not apply and the enclosure is declined rather than + stretched to fit.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + with pytest.raises(ak.ValidatedError): + ak.verified_integral((pool.integer(1) - ak.cos(x)) / (x * x), x, 0.0, 1.0) + + +@pytest.mark.parametrize( + ("build", "a", "b", "label"), + [ + (lambda pool, x: -ak.log(x), 0.0, 1.0, "-log x"), + (lambda pool, x: ak.log(x) * ak.log(x), 0.0, 1.0, "(log x)^2"), + (lambda pool, x: ak.exp(x * ak.log(x)), 0.0, 1.0, "x^x"), + ], +) +def test_integrable_but_not_removable_singularities_refuse_with_an_honest_message( + build, a, b, label +): + """These integrals all exist. What does not exist is a rigorous enclosure + of the *integrand*, and the error text must say which of the two it means + rather than implying the integral is undefined.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + with pytest.raises(ak.ValidatedError) as excinfo: + ak.verified_integral(build(pool, x), x, a, b) + + message = str(excinfo.value) + assert "integrand is singular" in message, f"{label}: {message}" + assert "integrable singularity" in message, f"{label}: {message}" + assert excinfo.value.remediation + + # --------------------------------------------------------------------------- # Three-valued predicates — the third value is never collapsed # --------------------------------------------------------------------------- @@ -124,6 +244,94 @@ def test_no_roots_verified_true(): assert ak.verified_no_roots(f, _box(x, 0.0, 1.0)) == "true" +@pytest.mark.parametrize( + ("lo", "hi", "roots"), + [ + (-2.0, 0.0, "1 root, endpoint signs + -> -"), + (0.0, 2.0, "1 root, endpoint signs - -> +"), + (1.3, 1.5, "1 root, endpoint signs - -> +"), + (-2.0, 2.0, "2 roots, endpoint signs + -> +"), + (-10.0, 10.0, "2 roots, endpoint signs + -> +"), + ], +) +def test_no_roots_false_whatever_the_root_count(lo, hi, roots): + """`x**2 - 2` has provable roots on every one of these boxes. An even + number of them used to defeat the test, because only the box's *own* + endpoints were checked for a sign change.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + assert ak.verified_no_roots(x * x - pool.integer(2), [(x, lo, hi)]) == "false", roots + + +def test_no_roots_false_when_the_roots_hide_behind_a_positive_factor(): + """`(x**2-2)(x**2+1)` has the same two roots; the second factor never + vanishes and never changes the sign pattern at the endpoints.""" + pool = ak.ExprPool() + x = pool.symbol("x") + f = (x * x - pool.integer(2)) * (x * x + pool.integer(1)) + + assert ak.verified_no_roots(f, [(x, -2.0, 2.0)]) == "false" + + +def test_no_roots_false_across_a_multivariate_box(): + """A box is convex, so two points of opposite proven sign certify a root + anywhere in it — `x - y` is +2 at (1,-1) and -2 at (-1,1).""" + pool = ak.ExprPool() + x, y = pool.symbol("x"), pool.symbol("y") + + assert ak.verified_no_roots(x - y, [(x, -1.0, 1.0), (y, -1.0, 1.0)]) == "false" + + +@pytest.mark.parametrize( + ("build", "lo", "hi"), + [ + (lambda pool, x: x * x + pool.integer(1), -10.0, 10.0), + (lambda pool, x: x * x + pool.integer(2), -2.0, 2.0), + (lambda pool, x: ak.exp(x), -5.0, 5.0), + ( + lambda pool, x: (x - pool.integer(5)) * (x - pool.integer(5)) + pool.integer(1), + 0.0, + 1.0, + ), + ], +) +def test_no_roots_true_cases_stay_true(build, lo, hi): + """The existence search must never be able to turn a proven `"true"` into + anything else — it only runs once the enclosure already contains zero.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + assert ak.verified_no_roots(build(pool, x), [(x, lo, hi)]) == "true" + + +@pytest.mark.parametrize( + ("build", "lo", "hi", "why"), + [ + ( + lambda pool, x: (x - pool.integer(1)) * (x - pool.integer(1)), + 0.0, + 2.0, + "double root at x=1: no sign change anywhere", + ), + ( + lambda pool, x: (x * x - pool.integer(1)) * (x * x - pool.integer(1)), + -2.0, + 2.0, + "two double roots at x=+-1", + ), + ], +) +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.""" + pool = ak.ExprPool() + x = pool.symbol("x") + + assert ak.verified_no_roots(build(pool, x), [(x, lo, hi)]) == "undecided", why + + def test_sign_positive_verified_true(): pool = ak.ExprPool() x = pool.symbol("x") diff --git a/tests/textbook_gate/test_tg_polynomials.py b/tests/textbook_gate/test_tg_polynomials.py index 25736e91..db548352 100644 --- a/tests/textbook_gate/test_tg_polynomials.py +++ b/tests/textbook_gate/test_tg_polynomials.py @@ -14,7 +14,7 @@ compared back to the original polynomial's coefficients. - ``MultiPoly``: has no coefficient accessor, so equality is checked via arithmetic: build a *reference* ``MultiPoly`` from a known expression, - subtract, and assert ``.is_zero()`` — this only depends on ``MultiPoly``'s + subtract, and assert ``.is_zero`` — this only depends on ``MultiPoly``'s own (canonical, content-addressed) arithmetic, not on any printed form. - ``real_roots``: each analytic root is checked against exactly one returned bracketing interval (``lo``/``hi``), with a small epsilon tolerance. @@ -94,7 +94,7 @@ def test_unipoly_factor_x4_minus_1(pool, x): factors = fz.factor_list() assert fz.unit == "1" assert len(factors) == 3 - degrees = sorted(f.degree() for f, _ in factors) + degrees = sorted(f.degree for f, _ in factors) assert degrees == [1, 1, 2] assert all(m == 1 for _, m in factors) _assert_factorization_reconstructs(p, [-1, 0, 0, 0, 1]) @@ -107,7 +107,7 @@ def test_unipoly_factor_x3_minus_x_three_linear_factors(pool, x): factors = fz.factor_list() assert fz.unit == "1" assert len(factors) == 3 - assert all(f.degree() == 1 and m == 1 for f, m in factors) + assert all(f.degree == 1 and m == 1 for f, m in factors) _assert_factorization_reconstructs(p, [0, -1, 0, 1]) @@ -120,7 +120,7 @@ def test_unipoly_factor_x2_plus_1_irreducible_over_z(pool, x): assert fz.unit == "1" assert len(factors) == 1 factor, multiplicity = factors[0] - assert factor.degree() == 2 + assert factor.degree == 2 assert multiplicity == 1 assert factor.coefficients() == [1, 0, 1] @@ -132,7 +132,7 @@ def test_unipoly_factor_cubic_rational_root_plus_irreducible_quadratic(pool, x): factors = fz.factor_list() assert fz.unit == "1" assert len(factors) == 2 - degrees = sorted(f.degree() for f, _ in factors) + degrees = sorted(f.degree for f, _ in factors) assert degrees == [1, 2] assert all(m == 1 for _, m in factors) _assert_factorization_reconstructs(p, [-2, 1, -2, 1]) @@ -145,7 +145,7 @@ def test_unipoly_factor_repeated_linear_factor(pool, x): factors = fz.factor_list() assert fz.unit == "1" assert len(factors) == 2 - assert all(f.degree() == 1 for f, _ in factors) + assert all(f.degree == 1 for f, _ in factors) multiplicities = sorted(m for _, m in factors) assert multiplicities == [1, 2] _assert_factorization_reconstructs(p, [3, -5, 1, 1]) @@ -159,15 +159,15 @@ def test_unipoly_factor_quartic_four_linear_factors(pool, x): factors = fz.factor_list() assert fz.unit == "1" assert len(factors) == 4 - assert all(f.degree() == 1 and m == 1 for f, m in factors) + assert all(f.degree == 1 and m == 1 for f, m in factors) _assert_factorization_reconstructs(p, [4, 0, -5, 0, 1]) def test_unipoly_degree_and_coefficients_roundtrip(pool, x): - """Sanity check: from_coefficients/.degree()/.coefficients() round-trip + """Sanity check: from_coefficients/.degree/.coefficients() round-trip exactly for a simple polynomial (2x^3 - 7 = 2x^3 + 0x^2 + 0x - 7).""" p = _unipoly(pool, x, [-7, 0, 0, 2]) - assert p.degree() == 3 + assert p.degree == 3 assert p.coefficients() == [-7, 0, 0, 2] @@ -179,7 +179,7 @@ def test_unipoly_gcd_shared_linear_factor(pool, x): a = _unipoly(pool, x, [-1, 0, 1]) # x^2 - 1 = (x-1)(x+1) b = _unipoly(pool, x, [2, -3, 1]) # x^2 - 3x + 2 = (x-1)(x-2) g = a.gcd(b) - assert g.degree() == 1 + assert g.degree == 1 # coefficients are only defined up to a unit; normalize so the constant # term's sign matches (x - 1) has coefficients [-1, 1] coeffs = g.coefficients() @@ -194,7 +194,7 @@ def test_unipoly_gcd_shared_quadratic_factor(pool, x): a = _unipoly(pool, x, _convolve(common, [-5, 1])) b = _unipoly(pool, x, _convolve(common, [7, 1])) g = a.gcd(b) - assert g.degree() == 2 + assert g.degree == 2 coeffs = g.coefficients() if coeffs[-1] < 0: coeffs = [-c for c in coeffs] @@ -206,7 +206,7 @@ def test_unipoly_gcd_coprime_polys_is_constant(pool, x): a = _unipoly(pool, x, [-1, 1]) b = _unipoly(pool, x, [-2, 1]) g = a.gcd(b) - assert g.degree() == 0 + assert g.degree == 0 # --- ak.real_roots --------------------------------------------------------- @@ -261,31 +261,31 @@ def test_real_roots_no_real_roots(pool, x): assert len(roots) == 0 -# --- MultiPoly.total_degree() --------------------------------------------- +# --- MultiPoly.total_degree --------------------------------------------- def test_multipoly_total_degree_binomial_cube(pool, x, y): """(x+y)^3 has total degree 3.""" mp = ak.MultiPoly.from_symbolic((x + y) ** 3, [x, y]) - assert mp.total_degree() == 3 + assert mp.total_degree == 3 def test_multipoly_total_degree_trinomial_square(pool, x, y, z): """(x+y+z)^2 has total degree 2.""" mp = ak.MultiPoly.from_symbolic((x + y + z) ** 2, [x, y, z]) - assert mp.total_degree() == 2 + assert mp.total_degree == 2 def test_multipoly_total_degree_mixed(pool, x, y): """x^2*y + y^3 -- both monomials have total degree 3.""" mp = ak.MultiPoly.from_symbolic(x**2 * y + y**3, [x, y]) - assert mp.total_degree() == 3 + assert mp.total_degree == 3 def test_multipoly_total_degree_linear(pool, x, y): """x + 2y + 1 has total degree 1.""" mp = ak.MultiPoly.from_symbolic(x + pool.integer(2) * y + pool.integer(1), [x, y]) - assert mp.total_degree() == 1 + assert mp.total_degree == 1 # --- MultiPoly.gcd() -------------------------------------------------------- @@ -294,7 +294,7 @@ def test_multipoly_total_degree_linear(pool, x, y): def _assert_multipoly_equal(a: ak.MultiPoly, b: ak.MultiPoly) -> None: """Compare two MultiPoly values via arithmetic (no coefficient accessor exists): subtract and check the result is the zero polynomial.""" - assert (a - b).is_zero() + assert (a - b).is_zero def test_multipoly_gcd_shared_linear_factor(pool, x, y): @@ -303,7 +303,7 @@ def test_multipoly_gcd_shared_linear_factor(pool, x, y): b = ak.MultiPoly.from_symbolic((x + y) * (x + pool.integer(2) * y), [x, y]) g = a.gcd(b) expected = ak.MultiPoly.from_symbolic(x + y, [x, y]) - assert g.total_degree() == 1 + assert g.total_degree == 1 _assert_multipoly_equal(g, expected) @@ -313,7 +313,7 @@ def test_multipoly_gcd_shared_trivariate_factor(pool, x, y, z): b = ak.MultiPoly.from_symbolic((x + y + z) * (y + z), [x, y, z]) g = a.gcd(b) expected = ak.MultiPoly.from_symbolic(x + y + z, [x, y, z]) - assert g.total_degree() == 1 + assert g.total_degree == 1 _assert_multipoly_equal(g, expected) @@ -322,7 +322,7 @@ def test_multipoly_gcd_coprime_is_constant(pool, x, y): a = ak.MultiPoly.from_symbolic(x, [x, y]) b = ak.MultiPoly.from_symbolic(y, [x, y]) g = a.gcd(b) - assert g.total_degree() == 0 + assert g.total_degree == 0 # --- MultiPoly.factor_z() --------------------------------------------------- @@ -334,7 +334,7 @@ def test_multipoly_factor_z_difference_of_squares(pool, x, y): fz = mp.factor_z() factors = fz.factor_list() assert len(factors) == 2 - assert all(f.total_degree() == 1 and m == 1 for f, m in factors) + assert all(f.total_degree == 1 and m == 1 for f, m in factors) # reconstruct the product and check it matches the original (up to unit) product = None