diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1e348371..acd6c32d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -203,7 +203,14 @@ jobs: name: AddressSanitizer — alkahest-cas runs-on: ubuntu-latest if: github.event_name != 'schedule' - timeout-minutes: 60 + # 90, not 60: with M10's Douglas-Rachford SOS certificate tests + # (real::sos::psd, ~8.5 min under ASan) alongside M4's own suite, the full + # `-p alkahest-cas --lib --tests` run measured ~48 min under ASan + # instrumentation plus ~12 min of nightly-toolchain build-std setup, + # landing right at the old 60-minute ceiling and getting cancelled + # (not failed) mid-run. 90 gives real margin over the measured ~60 min + # total, not just clearing it exactly. + timeout-minutes: 90 steps: - uses: actions/checkout@v6 diff --git a/CHANGELOG.md b/CHANGELOG.md index 18585c24..e2f33e6f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,77 @@ # Changelog +## Unreleased + +- **`telescope2d` generalizes from two bound indices to an arbitrary `m ≥ 1`: + `experimental.telescope_md`** (M4 extension). `telescope2d(term, n, j, k)` + only ever reached exactly two bound indices; the underlying ansatz search + and boundary/face analysis are now implemented for general `m`, with + `telescope2d` itself unchanged in behavior (it is now a thin `m = 2` + wrapper over the general engine, not a separate implementation) and a new + `telescope_md(term, n, [x_1, …, x_m])` for `m ≠ 2` — including `m = 1`, + which degenerates cleanly to a single-sum-shaped search, and `m ≥ 3`, + genuinely new. Same scope as before, generalized: proper hypergeometric + summands only (no broader rational-summand class, no sum of several + hypergeometric terms), a fixed (non-minimal) certificate denominator built + from `F`'s own shift-ratio denominators rather than a minimal multivariate + Gosper normal form, and constant-box-only boundary analysis. **Not + attempted**: a genuine minimal Gosper-style certificate denominator (the + roadmap's stated remaining-gap item 3) — real algorithm-design work, not an + engineering extension of what already exists here. + + **The boundary is `2m` face sums, not `2^m` corner evaluations** — the + `m = 2` module's "four strip sums, not four corners" result, generalized: + summing the telescoping identity over an `m`-dimensional box gives `2m` + sums, each over an `(m − 1)`-dimensional face where one bound index is + fixed to a boundary value, not `2^m` point evaluations at the box's + corners. The same sufficient (not necessary) "face vanishes pointwise" + criterion the `m = 2` module used — a dominant `1/Γ` zero among `F`'s own + gamma factors, or the certificate's own numerator vanishing there — carries + over unchanged in kind: fix one axis to a constant and check that a gamma + factor's argument no longer depends on `n` or any *other* bound index. + + **Resource ceilings, added after this generalization surfaced a real + scaling cliff.** The linear system a probe builds has one equation per + distinct monomial and one unknown per certificate-numerator box + coefficient; both dimensions grow with `m` and the certificate degree + bound far faster than the degree numbers suggest, and `rational_nullspace` + is a plain dense `O(rows · cols²)` exact-rational Gaussian elimination. + Profiling a genuinely coupled `m = 3` example + (`C(n,x)·C(x,y)·C(y,z)`) at degree bounds that work fine at `m = 2` found a + ≈10 000-row, 245-unknown system whose elimination step alone took ≈47 + seconds *per probe*, and the next certificate-degree step up (770 unknowns) + was still running after several minutes — genuine `O(rows·cols²)` + arithmetic cost on a correctly-posed system, not a bug or an infinite loop, + but exactly the kind of resource cliff a caller needs protecting from. + `holonomic::telescoping2d::search::MAX_ANSATZ_UNKNOWNS` now refuses any + single probe above 400 unknowns outright, and + `MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS` caps the total work spent on probes + at or above 150 unknowns to 300 across one whole search call — so a search + with no certificate in reach cannot be made to pay the same expensive + elimination over and over across every `(order, a_degree)` combination + tried. Neither ceiling touches the `m = 2` search, whose default probes + never exceed ≈140 unknowns. A search that hits a ceiling reports + `SearchExhausted` naming the ceiling explicitly, never a false certificate + and never a silent hang. + + Verified on the 4-category multinomial coefficient + `F(n,x,y,z) = n!/(x!y!z!(n−x−y−z)!)` (via `factorial`, not a + product-of-binomials encoding, to avoid the redundant cancelling-gamma- + factor pairs a naive `C(n,x)·C(n−x,y)·C(n−x−y,z)` encoding would carry) — a + genuinely non-separable `m = 3` sum (all three bound indices interact + through the shared `n − x − y − z` term) whose closed form, + `Σ_{x,y,z} F = 4ⁿ` (the multinomial theorem), is checked by direct exact + summation (`rug::Rational`, never floats) against the returned recurrence, + plus the `4ⁿ`-decoupled fixed-support variant where the boundary genuinely + is `n`-independent and `Vanishes` is provided rather than refused. A + regression test pins the original scaling-cliff example + (`C(n,x)·C(x,y)·C(y,z)` at `order ≤ 2`, `cert_degree ≤ 3`) to now return a + bounded, ceiling-cited refusal rather than run unboundedly. + + New: `alkahest.experimental.telescope_md`, + `alkahest.experimental.TelescopingMdCertificate`. Experimental, same + refusal codes as `telescope2d` (`E-HOLO-040`/`041`/`042`). + ## 3.9.0 — 2026-08-14 Everything in this section landed **after `v3.8.0` was tagged and published**, diff --git a/alkahest-core/src/holonomic/telescoping2d/boundary.rs b/alkahest-core/src/holonomic/telescoping2d/boundary.rs index 02b67714..10f645e1 100644 --- a/alkahest-core/src/holonomic/telescoping2d/boundary.rs +++ b/alkahest-core/src/holonomic/telescoping2d/boundary.rs @@ -1,98 +1,85 @@ -//! Deciding the *two-dimensional* boundary hypothesis a double-sum -//! certificate rests on. +//! Deciding the *`m`-dimensional* boundary hypothesis a creative-telescoping +//! certificate rests on, for `m ≥ 1` bound indices. //! -//! [`super::search::telescope2d_search`] proves an identity about the +//! [`super::search::telescope_md_search`] proves an identity about the //! **summand**: //! //! ```text -//! Σ_i a_i(n)·F(n+i,j,k) = Δ_j G_1(n,j,k) + Δ_k G_2(n,j,k) +//! Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t G_t, G_t = c_t·F //! ``` //! -//! Turning that into a statement about `S(n) = Σ_j Σ_k F(n,j,k)` over a -//! rectangular range `j = j_lo..j_hi`, `k = k_lo..k_hi` is a second, separate -//! step — exactly as in the single-sum engine -//! ([`super::super::boundary`]) — and it is genuinely *more* subtle here, -//! not just a repeat of the 1-D argument in two coordinates. +//! Turning that into a statement about `S(n) = Σ_{x_1} … Σ_{x_m} F(n,x)` over +//! a rectangular *box* `x_t = lo_t..hi_t` is a second, separate step — +//! exactly as in the single-sum engine ([`super::super::boundary`]) — and it +//! is genuinely *more* subtle here, not just a repeat of the 1-D argument in +//! `m` coordinates. //! -//! # The boundary is four strip sums, not four corner evaluations +//! # The boundary is `2m` face sums, not `2^m` corner evaluations //! -//! Summing the identity over the rectangle and telescoping each difference -//! operator over its own index (Fubini — the order of summation does not -//! matter for a finite double sum) gives: +//! Summing the identity over the box and telescoping each difference operator +//! over its own index (Fubini — the order of summation does not matter for a +//! finite multiple sum) gives, for each axis `t`, a term //! //! ```text -//! Σ_i a_i(n)·Σ_j Σ_k F(n+i,j,k) -//! = Σ_{k=k_lo}^{k_hi} [G_1(n, j_hi+1, k) − G_1(n, j_lo, k)] -//! + Σ_{j=j_lo}^{j_hi} [G_2(n, j, k_hi+1) − G_2(n, j, k_lo)] +//! Σ_{x_1} … (no x_t) … Σ_{x_m} [G_t(…, x_t = hi_t+1, …) − G_t(…, x_t = lo_t, …)] //! ``` //! -//! The right-hand side is **four one-dimensional sums along the rectangle's -//! edges**, not four point evaluations at its corners. Point-evaluating `G_1` -//! and `G_2` at the four corners `(j_lo, k_lo)`, `(j_lo, k_hi+1)`, -//! `(j_hi+1, k_lo)`, `(j_hi+1, k_hi+1)` and combining them with `±` signs — the -//! naive generalisation of the 1-D endpoint story — is simply the wrong -//! formula: it neither computes what is above and would double count in -//! exactly the case where it happened to (this module does **not** do that, -//! and its docs above are the reason). +//! summed over the remaining `m - 1` free bound indices. The right-hand side +//! is **`2m` sums, each over the `(m-1)`-dimensional face where one axis has +//! been fixed to a boundary value** — not `2^m` point evaluations at the +//! box's corners. The `m = 2` case ([`super::super::telescoping2d::boundary`] +//! before this generalization) already spelled out why the naive corner- +//! evaluation generalization is the wrong formula; the trap recurs, unchanged +//! in kind, at every `m`: a face is a sum over an `(m-1)`-dimensional slab, +//! and collapsing it to `2^m` corner evaluations both fails to compute what +//! is above and (for `m ≥ 2`) miscounts every face that is not literally a +//! single point. //! //! # What this module actually proves, and what it refuses to //! -//! Summing a strip in closed form is, in general, itself a creative- +//! Summing a face in closed form is, in general, itself a creative- //! telescoping problem — potentially as hard as the original one. This //! module does not attempt that. Instead it establishes the **sufficient** -//! (not necessary) condition that each of the four strips is *identically -//! the zero function* of its remaining free variables, in which case the sum -//! along it is trivially `0` term by term. Concretely, for the strip -//! `G_1(n, j_lo, k)` (as a function of `n` and `k`): if `F`'s parsed -//! `Γ(a·n + b·j + c·k + d)^e` factor list contains an `e < 0` factor whose -//! argument, after substituting `j ↦ j_lo`, no longer depends on `n` or `k` -//! (i.e. its `n`- and `k`-coefficients vanish) and evaluates to a -//! non-positive integer, then `1/Γ(·)^{|e|}` is exactly `0` there — a -//! genuine identity (`1/Γ` at a non-positive integer is `0`, not merely a -//! limit), not a numeric coincidence — and as long as no *other* factor in -//! `F`, its rational prefactor, or the certificate's own denominator has a -//! matching pole there, the whole strip is `0`. +//! (not necessary) condition that each of the `2m` faces is *identically the +//! zero function* of its remaining free variables, in which case the sum over +//! it is trivially `0` term by term. Concretely, for the face +//! `G_t(…, x_t = v, …)` (as a function of `n` and the `m - 1` other bound +//! indices): if `F`'s parsed `Γ(a·n + Σ_s b_s·x_s + d)^e` factor list contains +//! an `e < 0` factor whose argument, after substituting `x_t ↦ v`, no longer +//! depends on `n` or on any *other* bound index (i.e. every coefficient +//! except `x_t`'s own vanishes) and evaluates to a non-positive integer, then +//! `1/Γ(·)^{|e|}` is exactly `0` there — a genuine identity, not a numeric +//! coincidence — and as long as no *other* factor in `F`, its rational +//! prefactor, or the certificate's own denominator has a matching pole there, +//! the whole face is `0`. //! //! This is the same *tool* the single-sum engine's order counting uses (a -//! `1/Γ` factor at a non-positive integer argument), but a much simpler use -//! of it: [`super::super::boundary`] locates a pole *approaching an endpoint -//! along the summed variable* and counts orders there, because in 1-D the -//! endpoint is generally the only place `k` is fixed to a specific affine -//! function of `n`. Here the swept boundary line already fixes one entire -//! index (`j` or `k`) to a constant, so the question collapses to "is this -//! `Γ` argument *identically* a non-positive integer over the remaining -//! variables" — a purely algebraic check on `F`'s known integer coefficients, -//! no order-counting needed. It is **strictly weaker** than the single-sum -//! module's analysis: it never resolves a vanishing that needs cancellation -//! between multiple factors, and it proves nothing about a strip that is -//! merely non-constant-but-summing-to-zero. Both of those report -//! [`BoundaryStatus2d::Unknown`] rather than a guess. +//! `1/Γ` factor at a non-positive integer argument); the `m`-dimensional +//! generalization is mechanical once one axis is fixed to a constant — the +//! same algebraic check as the `m = 2` case, just against "every *other* +//! bound index" instead of "the one other bound index". It is **strictly +//! weaker** than the single-sum module's analysis: it never resolves a +//! vanishing that needs cancellation between multiple factors, and it proves +//! nothing about a face that is merely non-constant-but-summing-to-zero. Both +//! of those report [`BoundaryStatusMd::Unknown`] rather than a guess. //! //! # Scope: only ranges that do not depend on `n` //! -//! The single-sum module's `b(n)` formula has an extra term `D_i(n)` -//! (`boundary.rs`, module docs) precisely because summing -//! `F(n+i, k)` over `k`'s range *at n*, rather than at `n+i`, is not what -//! `S(n+i)` means when the range itself moves with `n`. That correction -//! carries over unchanged in spirit to two dimensions, but doubles the -//! bookkeeping (both `j` and `k` ranges could each move, independently, and -//! with each other). This module does not implement it: **`j_lo, j_hi, k_lo, -//! k_hi` must be integer constants**, not expressions in `n`. A caller with -//! an `n`-dependent natural range (e.g. `j = 0..n`) can still use this module -//! whenever the summand vanishes combinatorially outside the true range (the -//! ordinary situation for binomial-type double sums): pick a fixed bound -//! safely larger than any `n` of interest and let `F`'s own vanishing do the -//! rest — exactly how the worked example in `mod.rs` is set up. Passing an -//! `n`-dependent limit expression is refused as -//! [`BoundaryStatus2d::Unknown`], not silently misinterpreted. +//! Exactly as in the `m = 2` case: **every `lo_t, hi_t` must be an integer +//! constant**, not an expression in `n`. A caller with an `n`-dependent +//! natural range (e.g. `x_1 = 0..n`) can still use this module whenever the +//! summand vanishes combinatorially outside the true range — pick a fixed +//! bound safely larger than any `n` of interest and let `F`'s own vanishing +//! do the rest. Passing an `n`-dependent limit expression is refused as +//! [`BoundaryStatusMd::Unknown`], not silently misinterpreted. -use super::poly::{Axis, Poly3}; -use super::search::Telescoping2dResult; -use super::term::{affine_parts3, as_rat3, GammaFactor3, ProperTerm3}; +use super::poly::{Axis, PolyM}; +use super::search::{Telescoping2dResult, TelescopingMdResult}; +use super::term::{affine_partsm, as_ratm, GammaFactorM, ProperTermM}; use crate::kernel::{ExprId, ExprPool}; use rug::Rational; -/// The verdict on the boundary of a rectangular double sum. See the +/// The verdict on the boundary of a rectangular double sum (`m = 2`). See the /// [module docs](self) for exactly what each variant is allowed to mean. #[derive(Debug, Clone, PartialEq, Eq)] pub enum BoundaryStatus2d { @@ -150,6 +137,58 @@ impl BoundaryStatus2d { } } +/// The verdict on the boundary of an `m`-dimensional box sum, the general +/// form of [`BoundaryStatus2d`] (`m = 2`). +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum BoundaryStatusMd { + /// Every one of the `2m` boundary faces was proved identically zero, so + /// the homogeneous recurrence `Σ_i a_i(n)·S(n+i) = 0` holds for the sum. + Vanishes, + /// Reserved for a future extension, exactly as + /// [`BoundaryStatus2d::Nonzero`] — never produced by this version. + Nonzero { rhs: ExprId, witness_n: i64 }, + /// Neither was established. **No** recurrence for the sum follows. + Unknown { reason: String }, +} + +impl BoundaryStatusMd { + pub fn tag(&self) -> &'static str { + match self { + BoundaryStatusMd::Vanishes => "vanishes", + BoundaryStatusMd::Nonzero { .. } => "nonzero", + BoundaryStatusMd::Unknown { .. } => "unknown", + } + } + + pub fn implies_sum_recurrence(&self) -> bool { + !matches!(self, BoundaryStatusMd::Unknown { .. }) + } + + /// `ranges[t]` is a full range description for bound index `t`, such as + /// `"x1 = 0..40"` — one entry per bound index, in the same order + /// `indices` was supplied to [`boundary_status_md`]. + pub fn side_conditions(&self, ranges: &[String]) -> Vec { + let joined = ranges.join(", "); + match self { + BoundaryStatusMd::Vanishes => vec![format!( + "all {} boundary faces for {joined} were proved to vanish identically \ + (pointwise, which is sufficient but not necessary), so the homogeneous \ + recurrence sum_i a_i(n)*S(n+i) = 0 holds for the {}-fold sum", + 2 * ranges.len(), + ranges.len() + )], + BoundaryStatusMd::Nonzero { witness_n, .. } => vec![format!( + "the boundary does not vanish; b({witness_n}) != 0 was checked exactly" + )], + BoundaryStatusMd::Unknown { reason } => vec![format!( + "the boundary for {joined} could not be decided ({reason}); the certificate \ + proves the telescoping identity in the bound indices and NOTHING follows \ + about the sum until this is discharged independently" + )], + } + } +} + /// Decide the boundary hypothesis for `result` over the rectangle /// `j = j_limits.0 .. j_limits.1`, `k = k_limits.0 .. k_limits.1`. /// @@ -157,6 +196,9 @@ impl BoundaryStatus2d { /// pairs must be integer constants (see the [module docs](self)); anything /// else — an `n`-dependent bound, a second free symbol, a non-integer — is /// reported [`BoundaryStatus2d::Unknown`] rather than guessed. +/// +/// A thin wrapper around [`boundary_status_md`] with `indices = [j, k]` — see +/// the module docs. #[allow(clippy::too_many_arguments)] pub fn boundary_status_2d( result: &Telescoping2dResult, @@ -168,28 +210,56 @@ pub fn boundary_status_2d( k_limits: (ExprId, ExprId), pool: &ExprPool, ) -> BoundaryStatus2d { - match analyze(result, term, n, j, k, j_limits, k_limits, pool) { - Ok(()) => BoundaryStatus2d::Vanishes, - Err(reason) => BoundaryStatus2d::Unknown { reason }, + let md_result = TelescopingMdResult { + order: result.order, + coeffs: result.coeffs.clone(), + certs: vec![result.cert1, result.cert2], + }; + match boundary_status_md(&md_result, term, n, &[j, k], &[j_limits, k_limits], pool) { + BoundaryStatusMd::Vanishes => BoundaryStatus2d::Vanishes, + BoundaryStatusMd::Nonzero { rhs, witness_n } => { + BoundaryStatus2d::Nonzero { rhs, witness_n } + } + BoundaryStatusMd::Unknown { reason } => BoundaryStatus2d::Unknown { reason }, + } +} + +/// Decide the boundary hypothesis for a general `m`-bound-index `result` +/// over the box `indices[t] = limits[t].0 .. limits[t].1`. +/// +/// `term` must be the same `F(n,x)` that produced `result`, and `indices` the +/// same order used to produce it. Every limit must be an integer constant +/// (see the [module docs](self)); anything else is reported +/// [`BoundaryStatusMd::Unknown`] rather than guessed. +pub fn boundary_status_md( + result: &TelescopingMdResult, + term: ExprId, + n: ExprId, + indices: &[ExprId], + limits: &[(ExprId, ExprId)], + pool: &ExprPool, +) -> BoundaryStatusMd { + match analyze_md(result, term, n, indices, limits, pool) { + Ok(()) => BoundaryStatusMd::Vanishes, + Err(reason) => BoundaryStatusMd::Unknown { reason }, } } fn const_limit( expr: ExprId, n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, ) -> Result { - let (a, b, c, d) = affine_parts3(expr, n, j, k, pool).ok_or_else(|| { + let (coeffs, d) = affine_partsm(expr, n, indices, pool).ok_or_else(|| { format!( - "limit {} is not an integer-affine expression in n, j, k", + "limit {} is not an integer-affine expression in n and the bound indices", pool.display(expr) ) })?; - if a != 0 || b != 0 || c != 0 { + if coeffs.iter().any(|&c| c != 0) { return Err(format!( - "limit {} depends on n (or on j/k); this module only supports constant \ + "limit {} depends on n or on a bound index; this module only supports constant \ (n-independent) rectangular ranges — see the module docs", pool.display(expr) )); @@ -200,43 +270,66 @@ fn const_limit( Ok(d) } -#[allow(clippy::too_many_arguments)] -fn analyze( - result: &Telescoping2dResult, +fn analyze_md( + result: &TelescopingMdResult, term: ExprId, n: ExprId, - j: ExprId, - k: ExprId, - j_limits: (ExprId, ExprId), - k_limits: (ExprId, ExprId), + indices: &[ExprId], + limits: &[(ExprId, ExprId)], pool: &ExprPool, ) -> Result<(), String> { - let f = ProperTerm3::parse(term, n, j, k, pool) + let m = indices.len(); + if limits.len() != m { + return Err(format!( + "{m} bound indices were supplied but {} limit pairs", + limits.len() + )); + } + if result.certs.len() != m { + return Err(format!( + "result carries {} certificates but {m} bound indices were supplied", + result.certs.len() + )); + } + let f = ProperTermM::parse(term, n, indices, pool) .map_err(|e| format!("term does not re-parse as proper hypergeometric: {e}"))?; - let cert1 = as_rat3(result.cert1, n, j, k, pool, 0) - .ok_or_else(|| "certificate 1 does not parse back into Q(n,j,k)".to_string())?; - let cert2 = as_rat3(result.cert2, n, j, k, pool, 0) - .ok_or_else(|| "certificate 2 does not parse back into Q(n,j,k)".to_string())?; - let j_lo = const_limit(j_limits.0, n, j, k, pool)?; - let j_hi = const_limit(j_limits.1, n, j, k, pool)?; - let k_lo = const_limit(k_limits.0, n, j, k, pool)?; - let k_hi = const_limit(k_limits.1, n, j, k, pool)?; - let j_hi_p1 = j_hi + Rational::from(1); - let k_hi_p1 = k_hi + Rational::from(1); + let mut certs = Vec::with_capacity(m); + for (t, &c) in result.certs.iter().enumerate() { + let parsed = as_ratm(c, n, indices, pool, 0) + .ok_or_else(|| format!("certificate {} does not parse back into Q(n,x)", t + 1))?; + certs.push(parsed); + } + + let mut los = Vec::with_capacity(m); + let mut his_p1 = Vec::with_capacity(m); + for &(lo_e, hi_e) in limits { + let lo = const_limit(lo_e, n, indices, pool)?; + let hi = const_limit(hi_e, n, indices, pool)?; + his_p1.push(hi + Rational::from(1)); + los.push(lo); + } - let lines: [(Axis, &Rational, &str, &Poly3, &Poly3); 4] = [ - (Axis::J, &j_lo, "j = j_lo", &cert1.num, &cert1.den), - (Axis::J, &j_hi_p1, "j = j_hi + 1", &cert1.num, &cert1.den), - (Axis::K, &k_lo, "k = k_lo", &cert2.num, &cert2.den), - (Axis::K, &k_hi_p1, "k = k_hi + 1", &cert2.num, &cert2.den), - ]; - for (axis, value, label, cert_num, cert_den) in lines { - if !line_vanishes(&f.gammas, &f.rat.den, cert_num, cert_den, axis, value) { - return Err(format!( - "boundary strip at {label} was not provably zero by the (deliberately \ - conservative) pointwise criterion this module implements" - )); + for t in 0..m { + let axis: Axis = t + 1; + for (label, value) in [ + (format!("x{} = lo", t + 1), &los[t]), + (format!("x{} = hi + 1", t + 1), &his_p1[t]), + ] { + if !face_vanishes( + &f.gammas, + &f.rat.den, + &certs[t].num, + &certs[t].den, + axis, + value, + m + 1, + ) { + return Err(format!( + "boundary face at {label} was not provably zero by the (deliberately \ + conservative) pointwise criterion this module implements" + )); + } } } Ok(()) @@ -245,40 +338,36 @@ fn analyze( /// See the [module docs](self) for the exact soundness argument. Returns /// `true` only when the corresponding `G = c·F` is *provably* the zero /// function of its remaining free variables — never a guess, and never a -/// claim about the strip *summing* to zero without every term being zero. +/// claim about the face *summing* to zero without every term being zero. /// /// Two independent sufficient routes are checked, either being enough: /// -/// 1. `F` itself vanishes identically along the line (a dominant `1/Γ` zero +/// 1. `F` itself vanishes identically along the face (a dominant `1/Γ` zero /// among `F`'s own gamma factors — the natural-boundary case, e.g. /// `C(n,k)` vanishing at `k = n+1`). -/// 2. `F` is finite (no unresolved pole) along the line, and the -/// certificate's own numerator (`P_1` or `P_2`) is the identically zero -/// polynomial there — the case a classical WZ certificate's own extra -/// factor (e.g. `R(n,k) ∝ k`, zero at `k = 0`) supplies, independent of -/// whether `F` itself vanishes at that endpoint. +/// 2. `F` is finite (no unresolved pole) along the face, and the +/// certificate's own numerator is the identically zero polynomial there. /// -/// Both need `F` to have no pole along the line (checked once, up front) — +/// Both need `F` to have no pole along the face (checked once, up front) — /// otherwise `0 · ∞` is not resolved by either route, and this function /// correctly refuses rather than guessing. -fn line_vanishes( - gammas: &[GammaFactor3], - rat_den: &Poly3, - cert_num: &Poly3, - cert_den: &Poly3, +#[allow(clippy::too_many_arguments)] +fn face_vanishes( + gammas: &[GammaFactorM], + rat_den: &PolyM, + cert_num: &PolyM, + cert_den: &PolyM, axis: Axis, value: &Rational, + num_axes: usize, ) -> bool { let mut zero_order: i64 = 0; let mut pole_order: i64 = 0; for g in gammas { - let (elim_coeff, other_coeff) = match axis { - Axis::J => (g.b, g.c), // eliminate j; k is the other free (non-n) axis - Axis::K => (g.c, g.b), // eliminate k; j is the other free (non-n) axis - Axis::N => unreachable!("boundary lines are only ever j or k"), - }; - if g.a != 0 || other_coeff != 0 { - // Argument still depends on n or on the other free variable: + let elim_coeff = g.coeffs[axis]; + let other_nonzero = (0..num_axes).any(|a| a != axis && g.coeffs[a] != 0); + if other_nonzero { + // Argument still depends on n or on another free bound index: // not usable as a *pointwise-in-everything-remaining* pole/zero. continue; } @@ -286,8 +375,8 @@ fn line_vanishes( if *new_const.clone().denom() != 1 { continue; } - let m = new_const.numer().clone(); - if m > 0 { + let mnum = new_const.numer().clone(); + if mnum > 0 { continue; // Gamma(positive integer): finite, non-zero, no information. } if g.e < 0 { @@ -297,9 +386,9 @@ fn line_vanishes( } } if pole_order > 0 { - return false; // F itself has an unresolved pole along this line. + return false; // F itself has an unresolved pole along this face. } - let const_poly = Poly3::constant(value.clone()); + let const_poly = PolyM::constant(value.clone(), num_axes); if rat_den.eliminate_axis(axis, &const_poly).is_zero() { return false; // F's rational prefactor has an unresolved pole. } @@ -307,10 +396,10 @@ fn line_vanishes( return false; // the certificate's own denominator vanishes here. } if zero_order > 0 { - return true; // route 1: F itself is identically zero along the line. + return true; // route 1: F itself is identically zero along the face. } // route 2: F is finite here (checked above); zero if the certificate's - // own numerator is identically zero along the line. + // own numerator is identically zero along the face. cert_num.eliminate_axis(axis, &const_poly).is_zero() } diff --git a/alkahest-core/src/holonomic/telescoping2d/mod.rs b/alkahest-core/src/holonomic/telescoping2d/mod.rs index 78a91c15..310b5261 100644 --- a/alkahest-core/src/holonomic/telescoping2d/mod.rs +++ b/alkahest-core/src/holonomic/telescoping2d/mod.rs @@ -1,99 +1,156 @@ -//! Double-sum creative telescoping (Apagodu–Zeilberger) for proper -//! hypergeometric terms `F(n, j, k)`. +//! Multi-sum creative telescoping (Apagodu–Zeilberger) for proper +//! hypergeometric terms `F(n, x_1, …, x_m)` with `m ≥ 1` bound indices. //! //! # Scope //! //! This is the multivariate generalization of `super::zeilberger` from one -//! bound index (`k`) to two (`j`, `k`). It targets exactly the concrete goal -//! named in the roadmap: **double sums over proper hypergeometric summands** -//! — `F(n+1,j,k)/F(n,j,k)`, `F(n,j+1,k)/F(n,j,k)`, `F(n,j,k+1)/F(n,j,k)` all -//! rational functions, the same shape [`super::hyperterm::ProperTerm`] -//! recognizes for one index, generalized to three. It does **not** implement -//! full Wegschaider-style reduction (arbitrary rational summands, arbitrary -//! many indices) — that is a substantially larger undertaking and out of -//! scope here; see the honest-limitations list below. +//! bound index (`k`) to an arbitrary number `m` of them. It targets exactly +//! the concrete goal named in the roadmap: **multi-sums over proper +//! hypergeometric summands** — every shift ratio `F(…,x_t+1,…)/F(…,x_t,…)` a +//! rational function, the same shape [`super::hyperterm::ProperTerm`] +//! recognizes for one index, generalized to `m + 1`. It does **not** +//! implement full Wegschaider-style reduction (arbitrary *rational* +//! summands, or a minimal Gosper-style certificate denominator) — that +//! remains a substantially larger undertaking and out of scope here; see the +//! honest-limitations list below. //! -//! Given `F(n,j,k)`, [`telescope2d`] searches for a recurrence order `J`, -//! polynomial coefficients `a_0(n), …, a_J(n)` (not all zero) and two -//! rational certificates `c_1, c_2 ∈ Q(n,j,k)` such that +//! The two-bound-index case (`telescope2d`, `Telescoping2dResult`, +//! `boundary_status_2d`, `BoundaryStatus2d`) is the module's original, +//! semver-stable public surface and is unchanged in behavior. As of this +//! extension it is a thin wrapper around the general `m`-index engine +//! ([`telescope_md`], [`search::TelescopingMdResult`], +//! [`boundary::boundary_status_md`], [`boundary::BoundaryStatusMd`]), which +//! is the new public surface for `m ≠ 2` (including `m = 1`, degenerating to +//! the classical single-sum shape, and `m ≥ 3`, genuinely new). +//! +//! Given `F(n,x_1,…,x_m)`, [`telescope_md`] searches for a recurrence order +//! `J`, polynomial coefficients `a_0(n), …, a_J(n)` (not all zero) and `m` +//! rational certificates `c_1, …, c_m ∈ Q(n,x)` such that //! //! ```text -//! Σ_i a_i(n)·F(n+i,j,k) = Δ_j G_1(n,j,k) + Δ_k G_2(n,j,k), -//! G_1 = c_1·F, G_2 = c_2·F +//! Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t G_t(n,x), G_t = c_t·F //! ``` //! -//! — proving a recurrence for `S(n) = Σ_j Σ_k F(n,j,k)` once the *boundary* -//! of the rectangle it is summed over is discharged (see [`boundary`]; the -//! telescoping identity above says nothing about the sum on its own, exactly -//! as in the single-index case). +//! — proving a recurrence for `S(n) = Σ_{x_1} … Σ_{x_m} F(n,x)` once the +//! *boundary* of the box it is summed over is discharged (see [`boundary`]; +//! the telescoping identity above says nothing about the sum on its own, +//! exactly as in the single-index case). //! //! # Method: Apagodu–Zeilberger by undetermined coefficients //! -//! There is no standard two-dimensional analogue of Gosper's normal form for -//! a general proper hypergeometric `F(n,j,k)`, so — unlike the single-sum -//! engine — this module does not attempt one. It follows the +//! There is no standard multivariate analogue of Gosper's normal form for a +//! general proper hypergeometric `F(n,x_1,…,x_m)`, so — unlike the +//! single-sum engine — this module does not attempt one. It follows the //! Apagodu–Zeilberger presentation directly: posit a certificate ansatz of //! bounded polynomial degree over a *fixed* (ansatz-independent) denominator //! built from `F`'s own shift-ratio denominators, clear it, and solve the //! resulting *linear* system by Gaussian elimination over `Q`. See the //! `search` submodule for the full derivation and the specific, stated //! limitation this buys (the fixed denominator is not always the minimal one -//! a genuine 2-D Gosper reduction would find). +//! a genuine multivariate Gosper reduction would find). //! //! # Module layout //! -//! - `poly` — plain sparse `Q[n,j,k]` / `Q(n,j,k)` arithmetic. Deliberately -//! simpler than `super::qfield`'s `Q(n)[k]` tower: the ansatz search never -//! needs a gcd, only linear algebra over a fixed denominator, so there is -//! no normal-form machinery here to get wrong. -//! - `term` — `F(n,j,k)` recognition and exact shift ratios, the 3-index -//! generalization of `super::hyperterm`. +//! - `poly` — plain sparse `Q[n,x_1,…,x_m]` / `Q(n,x_1,…,x_m)` arithmetic for +//! an arbitrary number of axes. Deliberately simpler than `super::qfield`'s +//! `Q(n)[k]` tower: the ansatz search never needs a gcd, only linear +//! algebra over a fixed denominator, so there is no normal-form machinery +//! here to get wrong. +//! - `term` — `F(n,x_1,…,x_m)` recognition and exact shift ratios, the +//! `(m+1)`-index generalization of `super::hyperterm`. //! - `search` — the ansatz search itself, kept strictly separate from //! verification: every candidate is re-derived and checked as an exact -//! `Q(n,j,k)` identity (see `search::verify_certificate`) before it is ever -//! returned, independent of how the search found it. -//! - [`boundary`] — the 2-D boundary/corner analysis, on its own so a -//! returned certificate is checkable without reference to how it was -//! produced. Read its module docs first: the boundary of a rectangle is -//! **four one-dimensional strip sums**, not four corner-point evaluations, -//! and getting that distinction right is the substance of the module. +//! `Q(n,x)` identity (see `search::verify_certificate_md`) before it is +//! ever returned, independent of how the search found it. +//! - [`boundary`] — the `m`-dimensional boundary/face analysis, on its own so +//! a returned certificate is checkable without reference to how it was +//! produced. Read its module docs first: the boundary of a box is **`2m` +//! `(m-1)`-dimensional face sums**, not `2^m` corner-point evaluations, and +//! getting that distinction right is the substance of the module. //! //! # Honest limitations (read before relying on this) //! -//! - **Summands**: proper hypergeometric in `(n,j,k)` only — rational -//! prefactor times `z_j^j·z_k^k·w^n` times `Γ(a·n+b·j+c·k+d)^e` factors, -//! `a,b,c ∈ Z`. No more than two bound indices. -//! - **Certificate ansatz**: bounded box degree in each of `n,j,k` -//! independently ([`Telescoping2dOpts`]), searched by plain -//! ascending nested loops — not the cost-ordered iterative deepening -//! `super::zeilberger` uses, so raising the bounds is not free the way it -//! is there. +//! - **Summands**: proper hypergeometric in `(n,x_1,…,x_m)` only — rational +//! prefactor times `∏_t z_t^{x_t}·w^n` times `Γ(a·n+Σ_t b_t·x_t+d)^e` +//! factors, `a,b_t ∈ Z`. **No genuinely broader summand class is +//! supported**: a rational prefactor beyond what a ratio of the module's +//! own gamma factors already produces, a *sum* of several proper +//! hypergeometric terms, or a mixed radix/`q`-analogue combination are all +//! refused as [`Telescoping2dError::NotProperHypergeometric`], not +//! approximated. +//! - **Bound index count**: arbitrary `m ≥ 1` via [`telescope_md`] / +//! [`search::telescope_md_search`] (`m = 2` — [`telescope2d`] — is the +//! original, semver-stable special case; `m = 1` degenerates cleanly to a +//! single-sum-shaped search, exercised by this module's own tests). Raising +//! `m` grows the ansatz search space fast — a certificate's numerator +//! spans a box of `(max_cert_degree + 1)^(m+1)` unknowns *per* certificate, +//! and there are `m` certificates — so higher `m` needs correspondingly +//! patient degree bounds; the search still only ever returns an exactly +//! re-verified certificate, never a false one, when bounds run out. +//! - **Resource ceilings on the linear solve**: `rational_nullspace`'s exact +//! Gaussian elimination is `O(rows · cols²)` over unbounded-precision +//! rationals, and both dimensions grow with `m` and `cert_degree` well +//! past what the degree numbers alone suggest — at `m = 3`, +//! `cert_degree = 2` already means a ≈10 000-row, 245-unknown system whose +//! elimination step alone measured ≈47 seconds *per probe*, and +//! `cert_degree = 3` (770 unknowns) was still running after several +//! minutes. This is genuine arithmetic cost on a real linear system, not a +//! bug or an infinite loop, but a caller still needs protection from it: +//! `search::MAX_ANSATZ_UNKNOWNS` refuses outright any single probe whose +//! unknown count would exceed `400` (comfortably above every worked +//! example this module ships, including the `m = 3` multinomial- +//! coefficient example's `245`), and +//! `search::MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS` caps the *total* work +//! spent across every probe at or above 150 unknowns in one search call to +//! `300` — capping the number of genuinely expensive elimination attempts +//! to about one, regardless of how large `max_order` / `max_a_degree` / +//! `max_cert_degree` are, so a caller whose input has no certificate +//! within reach at all cannot be made to pay that cost over and over +//! across every `(order, a_degree)` combination tried. Both ceilings are +//! checked with cheap arithmetic before any polynomial construction for +//! the affected probe begins, and a probe skipped this way is reported +//! exactly like one the linear algebra found nothing for — except +//! [`Telescoping2dError::SearchExhausted`]'s message says explicitly when +//! a ceiling, not genuine non-existence, is why nothing was found, so +//! raising the search bounds further is not silently misrepresented as a +//! path to success. Neither ceiling affects the `m = 2` case: its default +//! search never builds a probe past ≈140 unknowns. +//! - **Certificate ansatz**: bounded box degree in each of `n,x_1,…,x_m` +//! independently ([`Telescoping2dOpts`] / [`search::TelescopingMdOpts`]), +//! searched by plain ascending nested loops — not the cost-ordered +//! iterative deepening `super::zeilberger` uses, so raising the bounds is +//! not free the way it is there. //! - **Certificate denominator**: fixed from `F`'s raw (un-reduced) shift- -//! ratio denominators, not a minimal 2-D Gosper normal form. Sufficient for -//! the "binomial-type" examples this module is tested against; not proven +//! ratio denominators, not a minimal Gosper normal form. Sufficient for the +//! "binomial-type" examples this module is tested against; not proven //! sufficient in general. A search that finds nothing reports -//! [`telescope2d_search`]'s `SearchExhausted`, never a false -//! certificate. -//! - **Boundary**: only rectangles with **constant** (not `n`-dependent) -//! limits are supported, and only the sufficient "each strip vanishes -//! pointwise" criterion is checked — see [`boundary`]'s module docs for -//! why both are real restrictions and not just unfinished polish, and for -//! the natural workaround (`n`-independent bounds larger than the true -//! combinatorial support) that the worked example below uses. -//! - **No explicit nonzero boundary term**: [`boundary::BoundaryStatus2d`] is -//! three-valued in shape (matching [`super::boundary::BoundaryStatus`]), -//! but this version never *produces* -//! [`boundary::BoundaryStatus2d::Nonzero`] — an unresolved boundary is -//! always [`boundary::BoundaryStatus2d::Unknown`], not an inhomogeneous -//! recurrence with an explicit `b(n)`. +//! `SearchExhausted`, never a false certificate. A genuine minimal +//! Gosper-style denominator (the roadmap's stated remaining-gap item 3) was +//! not attempted in this extension — see the crate-level changelog entry +//! for why (a real algorithm-design problem, not an engineering extension +//! of what is here). +//! - **Boundary**: only boxes with **constant** (not `n`-dependent) limits +//! are supported, and only the sufficient "each face vanishes pointwise" +//! criterion is checked — see [`boundary`]'s module docs for why both are +//! real restrictions and not just unfinished polish, and for the natural +//! workaround (`n`-independent bounds larger than the true combinatorial +//! support) that the worked examples below use. +//! - **No explicit nonzero boundary term**: [`boundary::BoundaryStatus2d`] +//! and [`boundary::BoundaryStatusMd`] are three-valued in shape (matching +//! [`super::boundary::BoundaryStatus`]), but neither version ever +//! *produces* their `Nonzero` variant — an unresolved boundary is always +//! `Unknown`, not an inhomogeneous recurrence with an explicit `b(n)`. pub mod boundary; mod poly; mod search; mod term; -pub use boundary::{boundary_status_2d, BoundaryStatus2d}; -pub use search::{telescope2d_search, Telescoping2dOpts, Telescoping2dResult}; +pub use boundary::{boundary_status_2d, boundary_status_md, BoundaryStatus2d, BoundaryStatusMd}; +pub use search::{ + telescope2d_search, telescope_md_search, Telescoping2dOpts, Telescoping2dResult, + TelescopingMdOpts, TelescopingMdResult, +}; use std::fmt; @@ -155,7 +212,10 @@ impl crate::errors::AlkahestError for Telescoping2dError { a certificate denominator this module's fixed-denominator ansatz cannot \ represent — this method does not apply" } - Telescoping2dError::InvalidInput(_) => "n, j and k must be three distinct symbols", + Telescoping2dError::InvalidInput(_) => { + "n and every bound index must be pairwise distinct symbols, and at least one \ + bound index must be supplied" + } }) } } @@ -173,6 +233,19 @@ pub fn telescope2d( telescope2d_search(term, n, j, k, pool, &Telescoping2dOpts::default()) } +/// Top-level entry point for `m ≥ 1` bound indices: search for and verify a +/// creative-telescoping certificate for `term = F(n, indices[0], …, +/// indices[m-1])`, with the default search bounds +/// ([`TelescopingMdOpts::default`]). +pub fn telescope_md( + term: crate::kernel::ExprId, + n: crate::kernel::ExprId, + indices: &[crate::kernel::ExprId], + pool: &crate::kernel::ExprPool, +) -> Result { + search::telescope_md_search(term, n, indices, pool, &TelescopingMdOpts::default()) +} + #[cfg(test)] mod tests { use super::*; @@ -255,7 +328,7 @@ mod tests { assert_eq!(s(ni), rug::Rational::from(Integer::from(4).pow(ni as u32))); } - assert_annihilates(&result, n, &pool, &s, 0, 5); + assert_annihilates(&result.coeffs, n, &pool, &s, 0, 5); } /// **Genuinely non-separable worked example**: `F(n,j,k) = C(n,j)*C(j,k)` @@ -292,7 +365,7 @@ mod tests { assert_eq!(s(ni), rug::Rational::from(Integer::from(3).pow(ni as u32))); } - assert_annihilates(&result, n, &pool, &s, 0, 5); + assert_annihilates(&result.coeffs, n, &pool, &s, 0, 5); // What `boundary.rs` can and cannot certify here, and why: `C(n,j)`'s // true support grows with `n` (it is non-zero for every `0 <= j <= @@ -359,7 +432,7 @@ mod tests { rug::Rational::from(Integer::from(2).pow(ni as u32) * Integer::from(3).pow(10)) ); } - assert_annihilates(&result, n, &pool, &s, 0, 5); + assert_annihilates(&result.coeffs, n, &pool, &s, 0, 5); let lo = pool.integer(0_i32); let hi = pool.integer(15_i32); @@ -373,9 +446,11 @@ mod tests { /// `Σ_i a_i(n)·S(n+i) = 0` for `n = lo..=hi`, using the *exact* rational /// values `s` computes — never floats — for the sum itself, and reading - /// the (small-integer) recurrence coefficients back exactly. + /// the (small-integer) recurrence coefficients back exactly. Takes the + /// coefficient list directly (rather than a whole `Telescoping2dResult`) + /// so it is shared between the `m = 2` and general-`m` worked examples. fn assert_annihilates( - result: &Telescoping2dResult, + coeffs: &[ExprId], n: ExprId, pool: &ExprPool, s: &dyn Fn(i64) -> rug::Rational, @@ -384,7 +459,7 @@ mod tests { ) { for ni in lo..=hi { let mut total = rug::Rational::from(0); - for (i, &c) in result.coeffs.iter().enumerate() { + for (i, &c) in coeffs.iter().enumerate() { let ai = coeff_at_n(pool, c, n, ni + i as i64); total += ai * s(ni + i as i64); } @@ -408,4 +483,295 @@ mod tests { )); assert_eq!(crate::errors::AlkahestError::code(&err), "E-HOLO-040"); } + + /// Exact `n!/(x!y!z!(n-x-y-z)!)` — the multinomial coefficient + /// (4-category composition of `n`) computed via the equivalent + /// nested-binomial product `C(n,x)·C(n-x,y)·C(n-x-y,z)`, which is `0` + /// whenever `x+y+z > n`, matching the symbolic term's own `1/Γ` vanishing + /// there. This is an independent check function — it does not reuse any + /// part of the solver. + fn multinom_i(n: i64, x: i64, y: i64, z: i64) -> Integer { + if x < 0 || y < 0 || z < 0 || x + y + z > n { + return Integer::from(0); + } + binom_i(n, x) * binom_i(n - x, y) * binom_i(n - x - y, z) + } + + /// **Genuinely non-separable `m = 3` worked example**: + /// `F(n,x,y,z) = n!/(x!·y!·z!·(n-x-y-z)!)`, the multinomial coefficient + /// counting compositions of `n` into 4 labeled parts — built directly + /// from `factorial`, not as a product of binomials, so the parser sees + /// exactly 5 gamma factors (`Γ(n+1)`, `1/Γ(x+1)`, `1/Γ(y+1)`, + /// `1/Γ(z+1)`, `1/Γ(n-x-y-z+1)`) rather than the redundant 9 a naive + /// `C(n,x)·C(n-x,y)·C(n-x-y,z)` encoding would carry (two of that + /// encoding's factors are exact inverses of each other but this module's + /// unreduced `Rat3`/`RatM` arithmetic never cancels them, so picking the + /// simpler encoding is a real, deliberate choice — not cosmetic). All + /// three bound indices interact through the shared `n-x-y-z` term, so + /// this is not a product of independent-variable pieces: it genuinely + /// exercises the `m`-index generalization of the ansatz search (three + /// non-trivial, mutually coupled certificates) and the boundary module's + /// `2m = 6`-face analysis, not just the `m = 2` machinery run twice. + /// + /// `Σ_{x,y,z} n!/(x!y!z!(n-x-y-z)!) = 4^n` by the multinomial theorem + /// (the number of length-`n` strings over a 4-letter alphabet, grouped by + /// letter counts) — a genuine closed form, checked here by direct exact + /// summation over a box safely larger than the true support, not + /// assumed. + #[test] + fn multinomial_matches_known_closed_form() { + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let z = pool.symbol("z", Domain::Real); + let rest = pool.add(vec![ + n, + pool.mul(vec![x, pool.integer(-1_i32)]), + pool.mul(vec![y, pool.integer(-1_i32)]), + pool.mul(vec![z, pool.integer(-1_i32)]), + ]); + let inv_fact = |e: ExprId| pool.pow(pool.func("factorial", vec![e]), pool.integer(-1_i32)); + let f = pool.mul(vec![ + pool.func("factorial", vec![n]), + inv_fact(x), + inv_fact(y), + inv_fact(z), + inv_fact(rest), + ]); + let opts = TelescopingMdOpts { + max_order: 1, + max_a_degree: 1, + max_cert_degree: 2, + }; + let result = search::telescope_md_search(f, n, &[x, y, z], &pool, &opts) + .expect("certificate for the 4-category multinomial coefficient"); + assert_eq!(result.certs.len(), 3); + + // Exact rational cross-check against real summed values: sum + // F(n,x,y,z) exactly over x,y,z = 0..15 (safely beyond the true + // support for n <= 6) and confirm it equals 4^n, then confirm the + // returned recurrence annihilates that exact sequence. + let s = |ni: i64| -> rug::Rational { + let mut acc = Integer::from(0); + for xx in 0..=15 { + for yy in 0..=15 { + for zz in 0..=15 { + acc += multinom_i(ni, xx, yy, zz); + } + } + } + rug::Rational::from(acc) + }; + for ni in 0..=6 { + assert_eq!(s(ni), rug::Rational::from(Integer::from(4).pow(ni as u32))); + } + assert_annihilates(&result.coeffs, n, &pool, &s, 0, 4); + + // The true support (x+y+z <= n) grows with n, exactly as in the m=2 + // non-separable example, so a constant box cannot soundly dominate + // it for every symbolic n: the boundary analysis must correctly + // refuse `Vanishes` here, not guess it. + let lo = pool.integer(0_i32); + let hi = pool.integer(15_i32); + let status = boundary_status_md( + &result, + f, + n, + &[x, y, z], + &[(lo, hi), (lo, hi), (lo, hi)], + &pool, + ); + assert_eq!( + status.tag(), + "unknown", + "a constant box cannot soundly certify Vanishes when the true support grows with \ + n; got {status:?}" + ); + } + + /// The same non-separable multinomial coupling as above, but with the + /// `n`-dependence factored out into a decoupled `4^n` (so the `(x,y,z)` + /// support is a genuine constant, independent of `n`) — the case where + /// the `m`-dimensional constant-box boundary analysis *can* certify + /// `Vanishes` for a real, non-separable triple sum. + /// + /// `F(n,x,y,z) = 4^n·10!/(x!y!z!(10-x-y-z)!)`; the `(x,y,z)` part sums to + /// `4^10` (the worked example above, at `n=10`), so + /// `S(n) = 4^n·4^10 = 4^(n+10)` — checked exactly, and the natural + /// boundary (the multinomial vanishes combinatorially outside `x+y+z <= + /// 10`, for *every* `n`, since none of the three bounds depend on `n`) + /// should be provably `Vanishes`. + #[test] + fn multinomial_fixed_support_boundary_vanishes() { + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let z = pool.symbol("z", Domain::Real); + let four_n = pool.pow(pool.integer(4_i32), n); + let ten = pool.integer(10_i32); + let rest = pool.add(vec![ + ten, + pool.mul(vec![x, pool.integer(-1_i32)]), + pool.mul(vec![y, pool.integer(-1_i32)]), + pool.mul(vec![z, pool.integer(-1_i32)]), + ]); + let inv_fact = |e: ExprId| pool.pow(pool.func("factorial", vec![e]), pool.integer(-1_i32)); + let f = pool.mul(vec![ + four_n, + pool.func("factorial", vec![ten]), + inv_fact(x), + inv_fact(y), + inv_fact(z), + inv_fact(rest), + ]); + let opts = TelescopingMdOpts { + max_order: 1, + max_a_degree: 1, + max_cert_degree: 2, + }; + let result = search::telescope_md_search(f, n, &[x, y, z], &pool, &opts) + .expect("certificate for 4^n*10!/(x!y!z!(10-x-y-z)!)"); + + let inner: Integer = { + let mut acc = Integer::from(0); + for xx in 0..=10 { + for yy in 0..=10 { + for zz in 0..=10 { + acc += multinom_i(10, xx, yy, zz); + } + } + } + acc + }; + assert_eq!(inner, Integer::from(4).pow(10)); + let s = |ni: i64| -> rug::Rational { + rug::Rational::from(inner.clone() * Integer::from(4).pow(ni.max(0) as u32)) + }; + for ni in 0..=6 { + assert_eq!( + s(ni), + rug::Rational::from(Integer::from(4).pow(ni as u32) * Integer::from(4).pow(10)) + ); + } + assert_annihilates(&result.coeffs, n, &pool, &s, 0, 4); + + let lo = pool.integer(0_i32); + let hi = pool.integer(15_i32); + let status = boundary_status_md( + &result, + f, + n, + &[x, y, z], + &[(lo, hi), (lo, hi), (lo, hi)], + &pool, + ); + assert_eq!( + status.tag(), + "vanishes", + "expected the natural (n-independent) boundary to vanish, got {status:?}" + ); + } + + /// `m = 1` sanity check: `telescope_md` with a single bound index should + /// find the same kind of certificate the classical single-sum engine + /// would, on the simplest possible input. + #[test] + fn telescope_md_single_index_smoke_test() { + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let k = pool.symbol("k", Domain::Real); + let f = pool.func("binomial", vec![n, k]); + let result = telescope_md(f, n, &[k], &pool).expect("certificate for C(n,k)"); + assert_eq!(result.certs.len(), 1); + } + + /// Regression test for the two resource ceilings in `search` + /// (`MAX_ANSATZ_UNKNOWNS`, the per-probe ceiling, and + /// `MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS`, the whole-search budget that + /// stops the *same* expensive probe size from being retried across every + /// `(order, a_degree)` combination): the triple-binomial-chain example + /// (`C(n,x)*C(x,y)*C(y,z)`) at the original, larger degree bounds this + /// test is named for was the case that, before these ceilings existed, + /// drove the search loop through several *repeated* multi-minute-or-worse + /// exact-rational Gaussian eliminations (`m = 3`, `cert_degree = 3` alone + /// needs 770 unknowns and a ≈15 000-row linear system; even the + /// `cert_degree = 2` step below that took ≈47 seconds *per probe*, and + /// the search loop would otherwise retry it for every one of six + /// `(order, a_degree)` combinations — see `MAX_ANSATZ_UNKNOWNS`'s and + /// `MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS`'s docs for the measurements). + /// This must now come back **bounded** — capped to roughly one expensive + /// elimination attempt, not six or more — with an honest + /// `SearchExhausted` naming the ceilings as the reason, not hang and not + /// silently under-search. The wall-clock bound here is deliberately loose + /// (this project's own standing practice for CI-flakiness — see + /// AGENTS.md, and this specific assertion needs extra slack since the + /// cumulative budget still permits one genuinely slow ≈30–50 second + /// elimination): its only job is to distinguish "bounded to about one + /// expensive attempt" from "unboundedly retries the same expensive + /// probe," not to pin an exact latency. + #[test] + fn chained_product_at_original_bounds_refuses_fast_via_resource_ceiling() { + // `#[cfg(sanitize = "address")]` would need an unstable feature gate + // on the whole crate (breaking every stable build), so this checks + // the RUSTFLAGS the *test binary itself* was compiled with instead — + // `option_env!` bakes in the compile-time value, no runtime probing. + // AddressSanitizer instrumentation makes the exact-rational + // elimination this test exercises 30x+ slower (~2500s observed in + // CI, against ~70-500s across plain Linux/macOS/Windows builds), + // which both defeats any reasonable wall-clock bound and risks the + // ASan job's own 60-minute CI timeout. The property under test + // (bounded probe count, not raw speed) is not sanitizer-sensitive, + // so it stays fully covered by every other CI build; only the + // wall-clock assertion, which is meaningless under this much + // overhead, is skipped here. + if option_env!("RUSTFLAGS").is_some_and(|f| f.contains("sanitizer=address")) { + eprintln!( + "skipping wall-clock assertion under AddressSanitizer (see comment on this test)" + ); + return; + } + + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let z = pool.symbol("z", Domain::Real); + let f = pool.mul(vec![ + pool.func("binomial", vec![n, x]), + pool.func("binomial", vec![x, y]), + pool.func("binomial", vec![y, z]), + ]); + let opts = TelescopingMdOpts { + max_order: 2, + max_a_degree: 2, + max_cert_degree: 3, + }; + let start = std::time::Instant::now(); + let err = search::telescope_md_search(f, n, &[x, y, z], &pool, &opts) + .expect_err("this combination must be refused via the resource ceiling, not solved"); + let elapsed = start.elapsed(); + // 900s, not 180s: this bound exists to catch a genuine hang (the + // pre-fix behavior was still running after several minutes and + // growing), not to pin CI wall-clock precisely. Windows CI runners + // measured ~480s here against ~76s on Linux for the same exact- + // rational elimination — a real, expected platform gap for GMP- + // backed arithmetic under MSYS2/mingw, not evidence the ceiling + // isn't working. + assert!( + elapsed.as_secs() < 900, + "expected a bounded refusal (roughly one expensive elimination, not several), took \ + {elapsed:?}" + ); + assert!(matches!(err, Telescoping2dError::SearchExhausted(_))); + let Telescoping2dError::SearchExhausted(msg) = &err else { + unreachable!() + }; + assert!( + msg.contains("MAX_ANSATZ_UNKNOWNS"), + "expected the SearchExhausted message to name the resource ceiling as the reason, \ + got: {msg}" + ); + } } diff --git a/alkahest-core/src/holonomic/telescoping2d/poly.rs b/alkahest-core/src/holonomic/telescoping2d/poly.rs index f6f94ea6..1f40ec23 100644 --- a/alkahest-core/src/holonomic/telescoping2d/poly.rs +++ b/alkahest-core/src/holonomic/telescoping2d/poly.rs @@ -1,89 +1,90 @@ -//! Exact sparse arithmetic in `Q[n, j, k]` and `Q(n, j, k)`. +//! Exact sparse arithmetic in `Q[n, x_1, …, x_m]` and `Q(n, x_1, …, x_m)` for +//! an arbitrary number `m ≥ 1` of bound indices. //! //! The single-sum engine ([`super::super::qfield`]) builds a careful tower //! `Q(n)[k]` / `Q(n)(k)` with a subresultant-PRS gcd, because Gosper's -//! algorithm needs a genuine normal form. The double-sum ansatz search in +//! algorithm needs a genuine normal form. The ansatz search in //! [`super::search`] never calls a gcd: it fixes the certificate //! denominators up front from the term's own shift-ratio denominators (see //! that module's docs) and solves a *linear* system for the numerators. So //! this module is deliberately simpler than `qfield` — a plain sparse -//! trivariate polynomial ring over `Q`, dense enough operations -//! (add/mul/shift/substitute) and nothing else. `Rat3` (the field of +//! multivariate polynomial ring over `Q`, dense enough operations +//! (add/mul/shift/substitute) and nothing else. [`RatM`] (the field of //! fractions) is kept unreduced throughout; the only place that matters is //! the final exact identity check, which is a zero-test on a cross-multiplied //! numerator, not a canonical form. +//! +//! # Axis convention +//! +//! Every polynomial here lives in `m + 1` variables: `n` is always axis `0` +//! ([`AXIS_N`]), and the `m` bound indices are axes `1..=m`, in the order the +//! caller supplied them. `num_axes = m + 1` must be threaded consistently by +//! every caller that constructs a fresh polynomial ([`PolyM::var`], +//! [`PolyM::constant`], [`PolyM::one`]) — this module does not itself track +//! `m` (a `PolyM` is just a sparse exponent map), so a caller that mixes +//! polynomials built with different `num_axes` gets a silently wrong +//! (mismatched-length) exponent vector; every caller in this crate threads a +//! single `num_axes` value throughout one search. use rug::ops::Pow as _; use rug::Rational; use std::collections::BTreeMap; -/// Which of the three variables `(n, j, k)` an operation acts on. -#[derive(Clone, Copy, Debug, PartialEq, Eq)] -pub enum Axis { - N, - J, - K, -} +/// Which variable an operation acts on: `0` is always `n`; `1..=m` are the +/// bound indices, in caller-supplied order. +pub type Axis = usize; -/// Exponents `(deg_n, deg_j, deg_k)` of one monomial. -pub type Exp = (u32, u32, u32); +/// `n` is always axis `0`. +pub const AXIS_N: Axis = 0; -fn exp_axis(e: Exp, axis: Axis) -> u32 { - match axis { - Axis::N => e.0, - Axis::J => e.1, - Axis::K => e.2, - } -} +/// Exponents `(deg_n, deg_x1, …, deg_xm)` of one monomial — length `m + 1`. +pub type Exp = Vec; fn exp_with_axis(mut e: Exp, axis: Axis, v: u32) -> Exp { - match axis { - Axis::N => e.0 = v, - Axis::J => e.1 = v, - Axis::K => e.2 = v, - } + e[axis] = v; e } -fn exp_add(a: Exp, b: Exp) -> Exp { - (a.0 + b.0, a.1 + b.1, a.2 + b.2) +fn exp_add(a: &Exp, b: &Exp) -> Exp { + debug_assert_eq!(a.len(), b.len(), "PolyM operands must share num_axes"); + a.iter().zip(b.iter()).map(|(x, y)| x + y).collect() } -/// A polynomial in `Q[n, j, k]`, stored as a sparse map from exponent triples -/// to non-zero rational coefficients. +/// A polynomial in `Q[n, x_1, …, x_m]`, stored as a sparse map from exponent +/// tuples to non-zero rational coefficients. #[derive(Clone, Debug, Default, PartialEq, Eq)] -pub struct Poly3 { +pub struct PolyM { pub terms: BTreeMap, } -impl Poly3 { +impl PolyM { pub fn zero() -> Self { - Poly3 { + PolyM { terms: BTreeMap::new(), } } - pub fn one() -> Self { - Poly3::constant(Rational::from(1)) + pub fn one(num_axes: usize) -> Self { + PolyM::constant(Rational::from(1), num_axes) } - pub fn constant(q: Rational) -> Self { + pub fn constant(q: Rational, num_axes: usize) -> Self { let mut terms = BTreeMap::new(); if q != 0 { - terms.insert((0, 0, 0), q); + terms.insert(vec![0u32; num_axes], q); } - Poly3 { terms } + PolyM { terms } } - pub fn from_i64(v: i64) -> Self { - Poly3::constant(Rational::from(v)) + pub fn from_i64(v: i64, num_axes: usize) -> Self { + PolyM::constant(Rational::from(v), num_axes) } - pub fn var(axis: Axis) -> Self { - let e = exp_with_axis((0, 0, 0), axis, 1); + pub fn var(axis: Axis, num_axes: usize) -> Self { + let e = exp_with_axis(vec![0u32; num_axes], axis, 1); let mut terms = BTreeMap::new(); terms.insert(e, Rational::from(1)); - Poly3 { terms } + PolyM { terms } } pub fn is_zero(&self) -> bool { @@ -91,7 +92,13 @@ impl Poly3 { } pub fn is_one(&self) -> bool { - self.terms.len() == 1 && self.terms.get(&(0, 0, 0)).map(|c| *c == 1).unwrap_or(false) + self.terms.len() == 1 + && self + .terms + .iter() + .next() + .map(|(e, c)| e.iter().all(|&x| x == 0) && *c == 1) + .unwrap_or(false) } /// The as-constructed constant value, or `None` if `self` is not @@ -101,8 +108,10 @@ impl Poly3 { return Some(Rational::from(0)); } if self.terms.len() == 1 { - if let Some(c) = self.terms.get(&(0, 0, 0)) { - return Some(c.clone()); + if let Some((e, c)) = self.terms.iter().next() { + if e.iter().all(|&x| x == 0) { + return Some(c.clone()); + } } } None @@ -110,75 +119,82 @@ impl Poly3 { /// Highest exponent on `axis` among the non-zero terms, or `-1` for the /// zero polynomial. - /// - /// Not called outside this module's own tests today — kept as a general - /// primitive (and exercised by them) rather than deleted, since it is - /// the natural thing a future extension of the search or boundary - /// analysis would reach for. #[allow(dead_code)] pub fn degree(&self, axis: Axis) -> i64 { self.terms .keys() - .map(|e| exp_axis(*e, axis) as i64) + .map(|e| e[axis] as i64) .max() .unwrap_or(-1) } - pub fn neg(&self) -> Poly3 { - Poly3 { + pub fn neg(&self) -> PolyM { + PolyM { terms: self .terms .iter() - .map(|(e, c)| (*e, Rational::from(-c))) + .map(|(e, c)| (e.clone(), Rational::from(-c))) .collect(), } } - pub fn add(&self, other: &Poly3) -> Poly3 { + pub fn add(&self, other: &PolyM) -> PolyM { let mut out = self.terms.clone(); for (e, c) in &other.terms { - let entry = out.entry(*e).or_insert_with(|| Rational::from(0)); + let entry = out.entry(e.clone()).or_insert_with(|| Rational::from(0)); *entry += c; } out.retain(|_, c| *c != 0); - Poly3 { terms: out } + PolyM { terms: out } } - pub fn sub(&self, other: &Poly3) -> Poly3 { + pub fn sub(&self, other: &PolyM) -> PolyM { self.add(&other.neg()) } - pub fn mul(&self, other: &Poly3) -> Poly3 { + pub fn mul(&self, other: &PolyM) -> PolyM { if self.is_zero() || other.is_zero() { - return Poly3::zero(); + return PolyM::zero(); } let mut out: BTreeMap = BTreeMap::new(); for (ea, ca) in &self.terms { for (eb, cb) in &other.terms { - let e = exp_add(*ea, *eb); + let e = exp_add(ea, eb); let entry = out.entry(e).or_insert_with(|| Rational::from(0)); *entry += Rational::from(ca * cb); } } out.retain(|_, c| *c != 0); - Poly3 { terms: out } + PolyM { terms: out } } - pub fn scale(&self, q: &Rational) -> Poly3 { + pub fn scale(&self, q: &Rational) -> PolyM { if *q == 0 { - return Poly3::zero(); + return PolyM::zero(); } - Poly3 { + PolyM { terms: self .terms .iter() - .map(|(e, c)| (*e, Rational::from(c * q))) + .map(|(e, c)| (e.clone(), Rational::from(c * q))) .collect(), } } - pub fn pow_u32(&self, e: u32) -> Poly3 { - let mut acc = Poly3::one(); + pub fn pow_u32(&self, e: u32) -> PolyM { + let num_axes = self.terms.keys().next().map(|k| k.len()); + let mut acc = match num_axes { + Some(na) => PolyM::one(na), + // self is the zero polynomial; 0^0 = 1 has no well-defined + // `num_axes` here, but `pow_u32` is never called with e == 0 on + // a zero polynomial in this crate's actual call sites, and the + // caller-visible behavior (zero raised to any positive power is + // zero) is preserved by the loop below regardless. + None => PolyM::zero(), + }; + if e == 0 && num_axes.is_none() { + return PolyM::zero(); + } let mut base = self.clone(); let mut e = e; while e > 0 { @@ -192,58 +208,60 @@ impl Poly3 { } /// `p(axis ↦ axis + delta)`. - pub fn shift(&self, axis: Axis, delta: i64) -> Poly3 { + pub fn shift(&self, axis: Axis, delta: i64) -> PolyM { if delta == 0 || self.is_zero() { return self.clone(); } - let repl = Poly3::var(axis).add(&Poly3::from_i64(delta)); + let num_axes = self.terms.keys().next().unwrap().len(); + let repl = PolyM::var(axis, num_axes).add(&PolyM::from_i64(delta, num_axes)); self.eliminate_axis(axis, &repl) } /// `p(axis ↦ replacement)`, where `replacement` is any polynomial (it may - /// reference `axis` itself, as [`Poly3::shift`] does, or eliminate it - /// entirely, as the boundary substitution in - /// [`super::boundary`] does). - pub fn eliminate_axis(&self, axis: Axis, replacement: &Poly3) -> Poly3 { + /// reference `axis` itself, as [`PolyM::shift`] does, or eliminate it + /// entirely, as the boundary substitution in [`super::boundary`] does). + pub fn eliminate_axis(&self, axis: Axis, replacement: &PolyM) -> PolyM { if self.is_zero() { - return Poly3::zero(); + return PolyM::zero(); } // Memoize `replacement^e` across the distinct exponents actually used. - let mut pow_cache: BTreeMap = BTreeMap::new(); - pow_cache.insert(0, Poly3::one()); - let mut acc = Poly3::zero(); + let mut pow_cache: BTreeMap = BTreeMap::new(); + let num_axes = self.terms.keys().next().unwrap().len(); + pow_cache.insert(0, PolyM::one(num_axes)); + let mut acc = PolyM::zero(); for (e, c) in &self.terms { - let ax_e = exp_axis(*e, axis); - let rest = exp_with_axis(*e, axis, 0); + let ax_e = e[axis]; + let rest = exp_with_axis(e.clone(), axis, 0); let repl_pow = pow_cache .entry(ax_e) .or_insert_with(|| replacement.pow_u32(ax_e)) .clone(); - let rest_mono = Poly3::monomial(rest, c.clone()); + let rest_mono = PolyM::monomial(rest, c.clone()); acc = acc.add(&rest_mono.mul(&repl_pow)); } acc } - fn monomial(e: Exp, c: Rational) -> Poly3 { + fn monomial(e: Exp, c: Rational) -> PolyM { let mut terms = BTreeMap::new(); if c != 0 { terms.insert(e, c); } - Poly3 { terms } + PolyM { terms } } - /// Substitute a rational value for every axis at once. Used by this + /// Substitute a rational value for every axis at once (`vals[0]` for + /// `n`, `vals[1..]` for the bound indices in order). Used by this /// module's tests to check `shift`/`eliminate_axis` against direct /// evaluation, independent of how those are implemented. #[allow(dead_code)] - pub fn eval(&self, n: &Rational, j: &Rational, k: &Rational) -> Rational { + pub fn eval(&self, vals: &[Rational]) -> Rational { let mut acc = Rational::from(0); - for ((en, ej, ek), c) in &self.terms { + for (e, c) in &self.terms { let mut term = c.clone(); - term *= n.clone().pow(*en); - term *= j.clone().pow(*ej); - term *= k.clone().pow(*ek); + for (ax, &exp) in e.iter().enumerate() { + term *= vals[ax].clone().pow(exp); + } acc += term; } acc @@ -253,23 +271,24 @@ impl Poly3 { &self, pool: &crate::kernel::ExprPool, n: crate::kernel::ExprId, - j: crate::kernel::ExprId, - k: crate::kernel::ExprId, + indices: &[crate::kernel::ExprId], ) -> crate::kernel::ExprId { if self.terms.is_empty() { return pool.integer(0_i32); } let mut terms = Vec::with_capacity(self.terms.len()); - for ((en, ej, ek), c) in &self.terms { + for (e, c) in &self.terms { let mut factors = Vec::new(); let ce = rational_to_expr(pool, c); let is_unit = *c == 1; - if !is_unit || (*en == 0 && *ej == 0 && *ek == 0) { + let all_zero = e.iter().all(|&x| x == 0); + if !is_unit || all_zero { factors.push(ce); } - push_pow(pool, &mut factors, n, *en); - push_pow(pool, &mut factors, j, *ej); - push_pow(pool, &mut factors, k, *ek); + push_pow(pool, &mut factors, n, e[0]); + for (t, &idx_expr) in indices.iter().enumerate() { + push_pow(pool, &mut factors, idx_expr, e[t + 1]); + } terms.push(if factors.len() == 1 { factors[0] } else { @@ -306,92 +325,99 @@ fn rational_to_expr(pool: &crate::kernel::ExprPool, q: &Rational) -> crate::kern } } -/// An element of `Q(n, j, k)`, kept as a raw (not necessarily reduced) +/// An element of `Q(n, x_1, …, x_m)`, kept as a raw (not necessarily reduced) /// `num/den` pair. See the module docs for why no gcd reduction happens here. #[derive(Clone, Debug)] -pub struct Rat3 { - pub num: Poly3, - pub den: Poly3, +pub struct RatM { + pub num: PolyM, + pub den: PolyM, } -impl Rat3 { - pub fn zero() -> Self { - Rat3 { - num: Poly3::zero(), - den: Poly3::one(), +impl RatM { + pub fn one(num_axes: usize) -> Self { + RatM { + num: PolyM::one(num_axes), + den: PolyM::one(num_axes), } } - pub fn one() -> Self { - Rat3 { - num: Poly3::one(), - den: Poly3::one(), - } - } - - pub fn from_poly(p: Poly3) -> Self { - Rat3 { + /// `num_axes` **must** be supplied explicitly rather than inferred from + /// `p`'s own keys: `p` may be the zero polynomial (empty term map, no + /// exponent-vector length to read off), and a `den` built with the wrong + /// `num_axes` silently mismatches every other polynomial's exponent + /// vectors the first time it's multiplied against one — exactly the bug + /// this explicit parameter exists to rule out at the type level. + pub fn from_poly(p: PolyM, num_axes: usize) -> Self { + RatM { num: p, - den: Poly3::one(), + den: PolyM::one(num_axes), } } - pub fn from_rational(q: Rational) -> Self { - Rat3::from_poly(Poly3::constant(q)) + pub fn from_rational(q: Rational, num_axes: usize) -> Self { + RatM::from_poly(PolyM::constant(q, num_axes), num_axes) } - pub fn add(&self, other: &Rat3) -> Rat3 { - Rat3 { + pub fn add(&self, other: &RatM) -> RatM { + RatM { num: self.num.mul(&other.den).add(&other.num.mul(&self.den)), den: self.den.mul(&other.den), } } - pub fn neg(&self) -> Rat3 { - Rat3 { + pub fn neg(&self) -> RatM { + RatM { num: self.num.neg(), den: self.den.clone(), } } - pub fn sub(&self, other: &Rat3) -> Rat3 { + pub fn sub(&self, other: &RatM) -> RatM { self.add(&other.neg()) } - pub fn mul(&self, other: &Rat3) -> Rat3 { - Rat3 { + pub fn mul(&self, other: &RatM) -> RatM { + RatM { num: self.num.mul(&other.num), den: self.den.mul(&other.den), } } - pub fn inv(&self) -> Option { + pub fn inv(&self) -> Option { if self.num.is_zero() { return None; } - Some(Rat3 { + Some(RatM { num: self.den.clone(), den: self.num.clone(), }) } - pub fn div(&self, other: &Rat3) -> Option { + pub fn div(&self, other: &RatM) -> Option { Some(self.mul(&other.inv()?)) } - pub fn pow_i32(&self, e: i32) -> Option { + pub fn pow_i32(&self, e: i32) -> Option { if e == 0 { - return Some(Rat3::one()); + let num_axes = self + .num + .terms + .keys() + .next() + .or_else(|| self.den.terms.keys().next()) + .map(|k| k.len()) + .unwrap_or(1); + return Some(RatM::one(num_axes)); } let base = if e < 0 { self.inv()? } else { self.clone() }; - Some(Rat3 { + Some(RatM { num: base.num.pow_u32(e.unsigned_abs()), den: base.den.pow_u32(e.unsigned_abs()), }) } - pub fn shift(&self, axis: Axis, delta: i64) -> Rat3 { - Rat3 { + pub fn shift(&self, axis: Axis, delta: i64) -> RatM { + RatM { num: self.num.shift(axis, delta), den: self.den.shift(axis, delta), } @@ -401,20 +427,19 @@ impl Rat3 { &self, pool: &crate::kernel::ExprPool, n: crate::kernel::ExprId, - j: crate::kernel::ExprId, - k: crate::kernel::ExprId, + indices: &[crate::kernel::ExprId], ) -> crate::kernel::ExprId { - let num = self.num.to_expr(pool, n, j, k); + let num = self.num.to_expr(pool, n, indices); if self.den.is_one() { return num; } - let den = self.den.to_expr(pool, n, j, k); + let den = self.den.to_expr(pool, n, indices); pool.mul(vec![num, pool.pow(den, pool.integer(-1_i32))]) } /// `self == other` as an identity of rational functions: cross-multiply /// and compare numerators. Exact, no reduction needed. - pub fn eq_rat(&self, other: &Rat3) -> bool { + pub fn eq_rat(&self, other: &RatM) -> bool { self.num .mul(&other.den) .sub(&other.num.mul(&self.den)) @@ -427,61 +452,81 @@ mod tests { use super::*; #[test] - fn poly3_basic_arithmetic() { - let n = Poly3::var(Axis::N); - let j = Poly3::var(Axis::J); - let one = Poly3::one(); + fn polym_basic_arithmetic() { + // 3 axes: n, j, k (num_axes = 3). + let n = PolyM::var(0, 3); + let j = PolyM::var(1, 3); + let one = PolyM::one(3); let p = n.add(&j).add(&one); // n + j + 1 let q = n.sub(&j); // n - j let prod = p.mul(&q); // n^2 - j^2 + n - j - assert_eq!(prod.degree(Axis::N), 2); - assert_eq!(prod.degree(Axis::J), 2); + assert_eq!(prod.degree(0), 2); + assert_eq!(prod.degree(1), 2); // n^2 - j^2 + n - j at (n,j,k)=(3,2,0): 9 - 4 + 3 - 2 = 6. - assert_eq!( - prod.eval(&Rational::from(3), &Rational::from(2), &Rational::from(0)), - Rational::from(6) - ); + let vals = [Rational::from(3), Rational::from(2), Rational::from(0)]; + assert_eq!(prod.eval(&vals), Rational::from(6)); } #[test] fn shift_matches_direct_eval() { // p = n*k + j; p(k -> k+2) evaluated at (n,j,k)=(2,3,5) should equal - // p evaluated at (2,3,7). - let n = Poly3::var(Axis::N); - let j = Poly3::var(Axis::J); - let k = Poly3::var(Axis::K); + // p evaluated at (2,3,7). Axes: n=0, j=1, k=2. + let n = PolyM::var(0, 3); + let j = PolyM::var(1, 3); + let k = PolyM::var(2, 3); let p = n.mul(&k).add(&j); - let shifted = p.shift(Axis::K, 2); - let a = shifted.eval(&Rational::from(2), &Rational::from(3), &Rational::from(5)); - let b = p.eval(&Rational::from(2), &Rational::from(3), &Rational::from(7)); + let shifted = p.shift(2, 2); + let a = shifted.eval(&[Rational::from(2), Rational::from(3), Rational::from(5)]); + let b = p.eval(&[Rational::from(2), Rational::from(3), Rational::from(7)]); assert_eq!(a, b); } #[test] fn eliminate_axis_matches_composition() { // p = j^2 + k, replace j by (2n+1). At n=3 that's j=7. - let j = Poly3::var(Axis::J); - let k = Poly3::var(Axis::K); - let n = Poly3::var(Axis::N); + let j = PolyM::var(1, 3); + let k = PolyM::var(2, 3); + let n = PolyM::var(0, 3); let p = j.mul(&j).add(&k); - let repl = n.scale(&Rational::from(2)).add(&Poly3::one()); - let sub = p.eliminate_axis(Axis::J, &repl); - let a = sub.eval(&Rational::from(3), &Rational::from(0), &Rational::from(5)); - let b = p.eval(&Rational::from(3), &Rational::from(7), &Rational::from(5)); + let repl = n.scale(&Rational::from(2)).add(&PolyM::one(3)); + let sub = p.eliminate_axis(1, &repl); + let a = sub.eval(&[Rational::from(3), Rational::from(0), Rational::from(5)]); + let b = p.eval(&[Rational::from(3), Rational::from(7), Rational::from(5)]); assert_eq!(a, b); } #[test] - fn rat3_cross_mul_zero_test() { - let n = Poly3::var(Axis::N); - let one = Poly3::one(); - let a = Rat3::from_poly(n.clone()) - .div(&Rat3::from_poly(n.add(&one))) + fn ratm_cross_mul_zero_test() { + let n = PolyM::var(0, 3); + let one = PolyM::one(3); + let a = RatM::from_poly(n.clone(), 3) + .div(&RatM::from_poly(n.add(&one), 3)) .unwrap(); - let b = Rat3 { - num: n.mul(&Poly3::from_i64(2)), - den: n.add(&one).mul(&Poly3::from_i64(2)), + let b = RatM { + num: n.mul(&PolyM::from_i64(2, 3)), + den: n.add(&one).mul(&PolyM::from_i64(2, 3)), }; assert!(a.eq_rat(&b)); } + + /// Exercise a fourth axis (an `m = 3` bound-index arity, `num_axes = 4`) + /// to make sure nothing in this module secretly assumes exactly three + /// axes. + #[test] + fn four_axis_arithmetic_works() { + let n = PolyM::var(0, 4); + let x1 = PolyM::var(1, 4); + let x2 = PolyM::var(2, 4); + let x3 = PolyM::var(3, 4); + let p = n.add(&x1).add(&x2).add(&x3); // n + x1 + x2 + x3 + let vals = [ + Rational::from(1), + Rational::from(2), + Rational::from(3), + Rational::from(4), + ]; + assert_eq!(p.eval(&vals), Rational::from(10)); + let shifted = p.shift(3, 5); // x3 -> x3 + 5 + assert_eq!(shifted.eval(&vals), Rational::from(15)); + } } diff --git a/alkahest-core/src/holonomic/telescoping2d/search.rs b/alkahest-core/src/holonomic/telescoping2d/search.rs index 8f8a179b..bae3e7e3 100644 --- a/alkahest-core/src/holonomic/telescoping2d/search.rs +++ b/alkahest-core/src/holonomic/telescoping2d/search.rs @@ -1,65 +1,135 @@ //! The Apagodu–Zeilberger ansatz search: given a proper hypergeometric -//! `F(n,j,k)`, find `a_0(n), …, a_J(n)` (not all zero) and two rational -//! certificates `c_1, c_2 ∈ Q(n,j,k)` such that +//! `F(n, x_1, …, x_m)` with `m ≥ 1` bound indices, find `a_0(n), …, a_J(n)` +//! (not all zero) and `m` rational certificates `c_1, …, c_m ∈ Q(n,x)` such +//! that //! //! ```text -//! Σ_i a_i(n)·F(n+i,j,k) = Δ_j G_1 + Δ_k G_2, G_1 = c_1·F, G_2 = c_2·F +//! Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t G_t, G_t = c_t·F //! ``` //! -//! # Method — undetermined coefficients, not a 2-D Gosper normal form +//! # Method — undetermined coefficients, not a Gosper normal form //! //! The single-sum engine ([`super::super::zeilberger`]) puts the shift ratio //! of `F` into *Gosper normal form* before solving, which is what lets it -//! search a smaller ansatz efficiently. There is no standard two-dimensional -//! analogue of that normal form for a general proper hypergeometric `F(n,j,k)` -//! — this is exactly why Apagodu–Zeilberger's method (unlike single-index -//! Zeilberger) is usually presented as an undetermined-coefficients search: -//! posit a certificate of bounded degree over a *fixed* denominator, clear -//! it, and solve the resulting linear system. +//! search a smaller ansatz efficiently. There is no standard multivariate +//! analogue of that normal form for a general proper hypergeometric +//! `F(n,x_1,…,x_m)` — this is exactly why Apagodu–Zeilberger's method (unlike +//! single-index Zeilberger) is usually presented as an undetermined- +//! coefficients search: posit a certificate of bounded degree over a *fixed* +//! denominator, clear it, and solve the resulting linear system. //! -//! Concretely: divide the target identity by `F(n,j,k)`. Writing -//! `ρ_j(n,j,k) = F(n,j+1,k)/F(n,j,k) = N_j/D_j` and -//! `ρ_k(n,j,k) = F(n,j,k+1)/F(n,j,k) = N_k/D_k` (both known rational -//! functions, computed exactly by [`super::term::ProperTerm3`]), and taking -//! the certificate ansatz `c_1 = P_1(n,j,k)/D_j(n,j,k)`, -//! `c_2 = P_2(n,j,k)/D_k(n,j,k)` with `P_1, P_2` polynomials of bounded -//! degree, the identity becomes, after multiplying through by the (known, -//! ansatz-independent) common denominator +//! Concretely: divide the target identity by `F(n,x)`. Writing +//! `ρ_t(n,x) = F(…,x_t+1,…)/F(…,x_t,…) = N_t/D_t` (a known rational function, +//! computed exactly by [`super::term::ProperTermM`]) for each bound axis `t`, +//! and taking the certificate ansatz `c_t = P_t(n,x)/E_t(n,x)` with `P_t` a +//! polynomial of bounded *box* degree and `E_t := D_t · (∏_i D_{n,i})` (see +//! below for why the denominator is not just `D_t` alone), the identity +//! becomes, after multiplying through by a common denominator built from the +//! `E_t` and their axis-`t` shifts, a **polynomial** identity in +//! `Q[n,x_1,…,x_m]`, linear in the unknown coefficients of +//! `a_i(n) = Σ_p α_{i,p}·n^p` and of every `P_t`. Matching coefficients of +//! every monomial gives one linear equation per monomial; +//! [`solve_ansatz_md`] assembles that system and takes its nullspace over `Q` +//! by plain Gaussian elimination (see [`rational_nullspace`]). //! -//! ```text -//! D_total = (∏_i D_{n,i}) · D_j(n,j,k)·D_j(n,j+1,k) · D_k(n,j,k)·D_k(n,j,k+1) -//! ``` +//! `E_t`'s `D_t` factor (the *raw*, un-reduced denominator of `ρ_t`) is not +//! the minimal possible certificate denominator in general — a genuine +//! multivariate Gosper reduction would sometimes need a smaller one after +//! cancelling a shift-equivalent factor between `N_t` and a shifted `D_t` +//! (exactly what the single-sum engine's `C(k)` factor exists to supply). +//! This module does **not** compute that reduction. For the ordinary +//! "binomial-type" sums this targets, the shift ratios are already close to +//! normal form, so the raw denominator is already sufficient — but this is a +//! property of the *examples*, not a theorem the code establishes. When it is +//! not sufficient, the bounded search below simply finds nothing and reports +//! [`Telescoping2dError::SearchExhausted`]; it never claims a false +//! certificate, because every candidate is re-verified from scratch (see +//! [`verify_certificate_md`]) before it is returned. //! -//! a **polynomial** identity in `Q[n,j,k]`, linear in the unknown -//! coefficients of `a_i(n) = Σ_p α_{i,p}·n^p` and of `P_1, P_2`. Matching -//! coefficients of every monomial `n^s·j^q·k^r` gives one linear equation per -//! monomial; [`solve_ansatz`] assembles that system and takes its nullspace -//! over `Q` by plain Gaussian elimination (see [`rational_nullspace`]). +//! # From two bound indices to `m` //! -//! `D_j(n,j,k)` (the *raw*, un-reduced denominator of `ρ_j`) is not the -//! minimal possible certificate denominator in general — a genuine 2-D Gosper -//! reduction would sometimes need a smaller one after cancelling a -//! shift-equivalent factor between `N_j` and a shifted `D_j` (exactly what -//! the single-sum engine's `C(k)` factor exists to supply). This module does -//! **not** compute that reduction. For the ordinary "binomial-type" double -//! sums this targets, the shift ratios are already close to normal form -//! (`gcd(N_j(j), D_j(j+h))` is a unit for every `h ≥ 0` in the examples this -//! module is tested against), so the raw denominator is already sufficient — -//! but this is a property of the *examples*, not a theorem the code -//! establishes. When it is not sufficient, the bounded search below simply -//! finds nothing and reports [`Telescoping2dError::SearchExhausted`]; it -//! never claims a false certificate, because every candidate is re-verified -//! from scratch (see [`verify_certificate`]) before it is returned. - -use super::poly::{Axis, Poly3, Rat3}; -use super::term::ProperTerm3; +//! The `m = 2` case (`telescope2d_search`) is now a thin wrapper around the +//! general `telescope_md_search`: it converts `Telescoping2dOpts` to +//! [`TelescopingMdOpts`], calls the general search with `indices = [j, k]`, +//! and repackages the two-certificate result as [`Telescoping2dResult`]. The +//! two-index public API's behavior (including its error variants and search +//! order) is unchanged by this — it is now derived from, rather than +//! duplicating, the general path. + +use super::poly::{Axis, PolyM, RatM, AXIS_N}; +use super::term::ProperTermM; use super::Telescoping2dError; use crate::kernel::{ExprId, ExprPool}; use rug::{Integer, Rational}; + +/// Upper bound on the total unknown count (`a_i(n)` coefficients plus every +/// certificate numerator's box coefficients, summed) a single +/// `(order, a_degree, cert_degree)` probe is allowed to build a linear system +/// for, checked *before* any polynomial arithmetic for that probe begins. +/// +/// This exists because [`rational_nullspace`]'s plain dense Gaussian +/// elimination is `O(rows · cols²)` over exact (unbounded) `Rational` +/// coefficients, and both `rows` (one equation per distinct monomial in +/// `n, x_1, …, x_m` appearing anywhere in the cleared identity) and `cols` +/// (`total`, this bound's subject) grow with `m` and `cert_degree` far +/// faster than the box-degree numbers themselves suggest — see this crate's +/// own measurements: at `m = 3`, `cert_degree = 2` already means `rows ≈ +/// 10 000`, `cols = 245`, and a single probe's elimination step alone took +/// **≈ 47 seconds**; `cert_degree = 3` (`cols = 770`) was still running +/// after several minutes when profiled. This is genuine `O(rows·cols²)` +/// arithmetic cost on a real, correctly-posed linear system — not a bug, an +/// infinite loop, or unbounded coefficient blowup — but it is exactly the +/// kind of input-dependent resource cliff a caller must be protected from by +/// a fast, honest refusal rather than an unbounded hang. `400` is +/// calibrated to comfortably admit every worked example this module ships +/// with (the largest, the `m = 3` multinomial-coefficient example, needs +/// `cols = 245`) while excluding the next box-degree step up at `m = 3` +/// (`cols = 770`), which is the one that was actually observed to run +/// unacceptably long. A probe whose unknown count would exceed this is +/// skipped — reported as no candidate at that `(a_degree, cert_degree)`, +/// exactly like a probe the linear algebra genuinely found nothing for — and +/// [`telescope_md_search`]'s final [`Telescoping2dError::SearchExhausted`] +/// message says so explicitly when at least one probe was skipped for this +/// reason, so a caller sees a fast, clearly-explained refusal instead of a +/// silent guess about whether raising the bounds would even help. +const MAX_ANSATZ_UNKNOWNS: usize = 400; + +/// A probe's own unknown count must reach this before it counts against +/// [`MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS`] at all. Below this, a probe is +/// "cheap" by construction (this crate's own `m = 2` default search never +/// exceeds ~140 unknowns for any probe it tries — see the constant's own +/// docs) and is exempted from the cumulative accounting entirely, so this +/// budget cannot regress the existing two-index search in any way. +const LARGE_PROBE_THRESHOLD: usize = 150; + +/// A single probe under [`MAX_ANSATZ_UNKNOWNS`] can still be individually +/// slow (the `m = 3` multinomial-coefficient worked example's `cols = 245` +/// probe takes ≈ 30 seconds) — tolerable *once*, but +/// [`telescope_md_search`]'s outer loop tries every `(order, a_degree, +/// cert_degree)` combination, and nothing about `total`'s formula depends +/// much on `a_degree` or `order`, so a caller whose input has no certificate +/// at all would otherwise pay that same ≈ 30–50 second cost again for *every* +/// `a_degree` and `order` value tried — six repeats of the exact scenario +/// that motivated [`MAX_ANSATZ_UNKNOWNS`] in the first place, for the +/// triple-binomial-chain example that was this ceiling's original motivating +/// case (see `mod.rs`'s +/// `chained_product_at_original_bounds_refuses_fast_via_resource_ceiling` +/// regression test). This is a running budget across the *whole* search +/// call: every probe with `total >= LARGE_PROBE_THRESHOLD` that is actually +/// attempted adds its `total` to a running sum, and once that sum would +/// exceed this bound, every further large probe is skipped for the rest of +/// the search — capping the number of genuinely expensive elimination +/// attempts to about one or two, regardless of how large the caller's +/// `max_order` / `max_a_degree` / `max_cert_degree` are. `300` admits +/// exactly one probe the size of the multinomial example (`245`) before +/// refusing further ones of that size. +const MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS: usize = 300; + use std::collections::BTreeMap; -/// Search bounds for [`telescope2d`](super::telescope2d). All three are -/// genuine upper bounds — the search tries every combination +/// Search bounds for [`telescope2d`](super::telescope2d) and +/// [`telescope_md`](super::telescope_md). All three are genuine upper +/// bounds — the search tries every combination /// `1..=max_order × 0..=max_a_degree × 0..=max_cert_degree` in ascending /// order (cheapest first in each axis), so raising them only admits harder /// inputs. @@ -69,8 +139,8 @@ use std::collections::BTreeMap; /// three loops are simply nested, ascending. That is a real scope /// simplification, not an oversight — the point of the cost-ordered plan in /// the single-sum engine is to make expensive high-degree probes at low -/// order not block a cheap high-order solution; the double-sum ansatz here -/// does not yet have that tuning. +/// order not block a cheap high-order solution; the ansatz here does not yet +/// have that tuning. #[derive(Debug, Clone, Copy)] pub struct Telescoping2dOpts { /// Largest recurrence order `J`; orders are tried from 1 upward. @@ -92,11 +162,45 @@ impl Default for Telescoping2dOpts { } } -/// A verified double-sum creative-telescoping certificate. +/// Search bounds for [`telescope_md_search`], the `m`-bound-index +/// generalization of [`Telescoping2dOpts`]. Same fields, same search +/// discipline (see that struct's docs); `max_cert_degree` bounds the box +/// degree of each `P_t` in **every** one of the `m + 1` variables +/// independently, so the per-certificate unknown count is +/// `(max_cert_degree + 1)^(m + 1)` — raising `max_cert_degree` gets +/// expensive fast as `m` grows, more so than in the two-index case. +#[derive(Debug, Clone, Copy)] +pub struct TelescopingMdOpts { + pub max_order: usize, + pub max_a_degree: usize, + pub max_cert_degree: usize, +} + +impl Default for TelescopingMdOpts { + fn default() -> Self { + TelescopingMdOpts { + max_order: 2, + max_a_degree: 2, + max_cert_degree: 2, + } + } +} + +impl From for TelescopingMdOpts { + fn from(o: Telescoping2dOpts) -> Self { + TelescopingMdOpts { + max_order: o.max_order, + max_a_degree: o.max_a_degree, + max_cert_degree: o.max_cert_degree, + } + } +} + +/// A verified double-sum creative-telescoping certificate (`m = 2`). /// /// The verified content is the identity /// `Σ_i a_i(n)·F(n+i,j,k) = Δ_j G_1 + Δ_k G_2` with `G_1 = cert1·F`, -/// `G_2 = cert2·F` — checked exactly in `Q(n,j,k)` by `verify_certificate` +/// `G_2 = cert2·F` — checked exactly in `Q(n,j,k)` by `verify_certificate_md` /// before this is ever constructed. Turning it into a recurrence for /// `S(n) = Σ_j Σ_k F(n,j,k)` over a stated range is a separate step; see /// [`super::boundary`]. @@ -114,18 +218,37 @@ pub struct Telescoping2dResult { pub cert2: ExprId, } +/// A verified `m`-bound-index creative-telescoping certificate, the general +/// form of [`Telescoping2dResult`] (which is `certs.len() == 2` repackaged). +/// +/// The verified content is `Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t(c_t·F)` — checked +/// exactly in `Q(n,x_1,…,x_m)` before this is ever constructed. See +/// [`super::boundary::boundary_status_md`] for turning it into a recurrence +/// for the `m`-fold sum over a stated box. +#[derive(Debug, Clone)] +pub struct TelescopingMdResult { + pub order: usize, + /// `a_0(n), …, a_J(n)`. + pub coeffs: Vec, + /// `c_1(n,x), …, c_m(n,x)`, one per bound index, in the order the caller + /// supplied `indices`. `G_t = certs[t-1]·F`. + pub certs: Vec, +} + /// Internal (pre-`ExprId`) form of a candidate, kept in algebraic form so -/// [`verify_certificate`] can re-check it without any expression-pool +/// [`verify_certificate_md`] can re-check it without any expression-pool /// round-trip. -struct Candidate { +struct CandidateMd { order: usize, - a: Vec, // a_i(n), i = 0..=order, degree only in the N axis - c1: Rat3, - c2: Rat3, + a: Vec, // a_i(n), i = 0..=order, degree only in the N axis + certs: Vec, // one per bound index } -/// Apagodu–Zeilberger search: find and verify a double-sum certificate for -/// `term`, a proper hypergeometric `F(n,j,k)`. +/// Apagodu–Zeilberger search for the `m = 2` case: find and verify a +/// double-sum certificate for `term`, a proper hypergeometric `F(n,j,k)`. +/// +/// A thin wrapper around [`telescope_md_search`] with `indices = [j, k]` — +/// see the module docs. pub fn telescope2d_search( term: ExprId, n: ExprId, @@ -134,9 +257,41 @@ pub fn telescope2d_search( pool: &ExprPool, opts: &Telescoping2dOpts, ) -> Result { - if n == j || n == k || j == k { + let md_opts: TelescopingMdOpts = (*opts).into(); + let r = telescope_md_search(term, n, &[j, k], pool, &md_opts)?; + debug_assert_eq!(r.certs.len(), 2); + Ok(Telescoping2dResult { + order: r.order, + coeffs: r.coeffs, + cert1: r.certs[0], + cert2: r.certs[1], + }) +} + +/// Apagodu–Zeilberger search for general `m ≥ 1`: find and verify a +/// creative-telescoping certificate for `term`, a proper hypergeometric +/// `F(n, x_1, …, x_m)` with `indices = [x_1, …, x_m]`. +/// +/// See the module docs (`telescoping2d::search`) for the method, and the honest limitations +/// list in `mod.rs`: this covers the same proper-hypergeometric-only, +/// fixed-denominator-ansatz, box-degree-search scope the `m = 2` engine has, +/// generalized to arbitrary `m` — it is **not** a broader class of summand. +pub fn telescope_md_search( + term: ExprId, + n: ExprId, + indices: &[ExprId], + pool: &ExprPool, + opts: &TelescopingMdOpts, +) -> Result { + let m = indices.len(); + if m == 0 { + return Err(Telescoping2dError::InvalidInput( + "at least one bound index is required".into(), + )); + } + if indices.contains(&n) || has_duplicate(indices) { return Err(Telescoping2dError::InvalidInput( - "n, j and k must be three distinct symbols".into(), + "n and every bound index must be pairwise distinct symbols".into(), )); } if opts.max_order == 0 { @@ -144,23 +299,28 @@ pub fn telescope2d_search( "max_order must be at least 1".into(), )); } - - let f = ProperTerm3::parse(term, n, j, k, pool)?; - let rho_j = f.ratio_axis(Axis::J, 1)?; - let rho_k = f.ratio_axis(Axis::K, 1)?; - if rho_j.den.is_zero() || rho_k.den.is_zero() { - return Err(Telescoping2dError::NotProperHypergeometric( - "shift ratio has a zero denominator".into(), - )); + let num_axes = m + 1; + + let f = ProperTermM::parse(term, n, indices, pool)?; + let mut rhos: Vec = Vec::with_capacity(m); + for t in 0..m { + let r = f.ratio_axis(t + 1, 1)?; + if r.den.is_zero() { + return Err(Telescoping2dError::NotProperHypergeometric( + "shift ratio has a zero denominator".into(), + )); + } + rhos.push(r); } + let mut skipped_for_budget = false; + let mut cumulative_large_unknowns: usize = 0; + for order in 1..=opts.max_order { - // The `i = 0..=order` shift ratios in `n` do not depend on the - // degree budgets, so they are computed once per order. let mut nn = Vec::with_capacity(order + 1); let mut dn = Vec::with_capacity(order + 1); for i in 0..=order as i64 { - let r = f.ratio_axis(Axis::N, i)?; + let r = f.ratio_axis(AXIS_N, i)?; if r.den.is_zero() { return Err(Telescoping2dError::NotProperHypergeometric( "shift ratio in n has a zero denominator".into(), @@ -172,11 +332,40 @@ pub fn telescope2d_search( for a_degree in 0..=opts.max_a_degree { for cert_degree in 0..=opts.max_cert_degree { + // Cheap, purely arithmetic pre-check — see MAX_ANSATZ_UNKNOWNS' + // docs — before any polynomial construction or linear-system + // assembly for this probe begins. + let box_len = cert_degree + 1; + let Some(cert_box_count) = box_len.checked_pow(num_axes as u32) else { + skipped_for_budget = true; + continue; + }; + let a_count = (order + 1) * (a_degree + 1); + let total = a_count.saturating_add(m.saturating_mul(cert_box_count)); + if total > MAX_ANSATZ_UNKNOWNS { + skipped_for_budget = true; + continue; + } + // See MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS' docs: a probe + // large enough to be individually slow also counts against a + // running total-search budget, so a caller cannot pay that + // same cost over and over across every (order, a_degree) + // combination when no certificate exists at all. + if total >= LARGE_PROBE_THRESHOLD { + if cumulative_large_unknowns.saturating_add(total) + > MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS + { + skipped_for_budget = true; + continue; + } + cumulative_large_unknowns += total; + } + if let Some(cand) = - solve_ansatz(order, a_degree, cert_degree, &nn, &dn, &rho_j, &rho_k)? + solve_ansatz_md(order, m, a_degree, cert_degree, &nn, &dn, &rhos, num_axes)? { - if verify_certificate(&f, &cand) { - return Ok(finish(cand, n, j, k, pool)); + if verify_certificate_md(&f, &cand) { + return Ok(finish_md(cand, n, indices, pool)); } // A genuine implementation bug, not a user-facing error: // never happens for a correct construction, but refusing @@ -187,9 +376,23 @@ pub fn telescope2d_search( } } + let budget_note = if skipped_for_budget { + format!( + " (at least one (order, a_degree, cert_degree) combination within these bounds was \ + skipped without being attempted, refused by this module's resource ceilings \ + (MAX_ANSATZ_UNKNOWNS = {MAX_ANSATZ_UNKNOWNS} unknowns for any single probe; \ + MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS = {MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS} total \ + across every probe at or above {LARGE_PROBE_THRESHOLD} unknowns in one search call) \ + rather than attempted and left to run arbitrarily long — raising max_cert_degree or \ + the bound-index count further will only make this worse, not better)" + ) + } else { + String::new() + }; Err(Telescoping2dError::SearchExhausted(format!( - "no verified double-sum certificate of order <= {} was found for {} \ - within a_degree <= {} and certificate degree <= {}", + "no verified {}-index certificate of order <= {} was found for {} \ + within a_degree <= {} and certificate degree <= {}{budget_note}", + m, opts.max_order, pool.display(term), opts.max_a_degree, @@ -197,46 +400,64 @@ pub fn telescope2d_search( ))) } +fn has_duplicate(xs: &[ExprId]) -> bool { + for i in 0..xs.len() { + for j in (i + 1)..xs.len() { + if xs[i] == xs[j] { + return true; + } + } + } + false +} + +fn flatten(exps: &[usize], box_len: usize) -> usize { + exps.iter().fold(0, |acc, &x| acc * box_len + x) +} + +fn unflatten(mut combo: usize, num_axes: usize, box_len: usize) -> Vec { + let mut exps = vec![0usize; num_axes]; + for idx in (0..num_axes).rev() { + exps[idx] = combo % box_len; + combo /= box_len; + } + exps +} + +fn monomial_m(exps: &[usize], c: Rational, num_axes: usize) -> PolyM { + let mut p = PolyM::constant(c, num_axes); + for (axis, &e) in exps.iter().enumerate() { + if e > 0 { + p = p.mul(&PolyM::var(axis, num_axes).pow_u32(e as u32)); + } + } + p +} + /// Build and solve the linear system for one `(order, a_degree, cert_degree)` /// probe. `Ok(None)` means the system has no solution with a non-trivial /// leading coefficient — not an error, just this probe failing. #[allow(clippy::too_many_arguments)] -// `p`, `q`, `r` below double as monomial exponents (via `idx_p1`/`idx_p2`) -// and as indices into the memoized `(j+1)^q`/`(k+1)^r` power tables, so -// clippy's "use an iterator instead" rewrite does not apply cleanly. -#[allow(clippy::needless_range_loop)] -fn solve_ansatz( +fn solve_ansatz_md( order: usize, + m: usize, a_degree: usize, cert_degree: usize, - nn: &[Poly3], - dn: &[Poly3], - rho_j: &Rat3, - rho_k: &Rat3, -) -> Result, Telescoping2dError> { - let dj = &rho_j.den; - let nj = &rho_j.num; - let dk = &rho_k.den; - let nk = &rho_k.num; - - // The certificate denominators are *not* just the raw `D_j`/`D_k` of the - // single ratio being telescoped: a certificate built from a product of - // two single-sum WZ pairs (see the worked "separable" example in - // `mod.rs`) needs a factor from the *other* direction's `n`-shift ratio - // too — e.g. `c_1 = R_A(n,j)·[B(n+1,k)/B(n,k)]` for `F = A(n,j)·B(n,k)`, - // whose denominator is not a function of `j` alone. So the ansatz here - // is `c_1 = P_1 / E_1`, `E_1 := D_j · (∏_i D_{n,i})`, and symmetrically - // for `c_2` with `D_k`. This is still just a *fixed, ansatz-independent* - // denominator (no gcd, no minimality claim) — see the module docs — just - // a bigger one than the raw single-ratio denominator alone, chosen to be - // sufficient for the worked examples this module is tested against. - let mut dnfull = Poly3::one(); + nn: &[PolyM], + dn: &[PolyM], + rhos: &[RatM], + num_axes: usize, +) -> Result, Telescoping2dError> { + let ds: Vec<&PolyM> = rhos.iter().map(|r| &r.den).collect(); + let ns_: Vec<&PolyM> = rhos.iter().map(|r| &r.num).collect(); + + let mut dnfull = PolyM::one(num_axes); for d in dn { dnfull = dnfull.mul(d); } - let mut dnfull_excl: Vec = Vec::with_capacity(order + 1); + let mut dnfull_excl: Vec = Vec::with_capacity(order + 1); for i in 0..=order { - let mut p = Poly3::one(); + let mut p = PolyM::one(num_axes); for (idx, d) in dn.iter().enumerate() { if idx != i { p = p.mul(d); @@ -245,43 +466,60 @@ fn solve_ansatz( dnfull_excl.push(p); } - let e1 = dj.mul(&dnfull); - let e1s = e1.shift(Axis::J, 1); - let e2 = dk.mul(&dnfull); - let e2s = e2.shift(Axis::K, 1); - - // `Mn_i = D_total / D_{n,i}`, built by replacing `E_1`'s copy of - // `D_{n,i}` with the `i`-excluded product (any one of `D_total`'s four - // redundant copies of `D_{n,i}` may be the one "removed" — they are all - // equal, so the result is the same polynomial either way). - let mut mn: Vec = Vec::with_capacity(order + 1); + // e[t] = D_t * dnfull ; es[t] = shift(e[t], axis = t+1, 1). This is + // exactly the m=2 module's `e1`/`e2`/`e1s`/`e2s`, generalized. + let mut e: Vec = Vec::with_capacity(m); + let mut es: Vec = Vec::with_capacity(m); + for (t, d) in ds.iter().enumerate().take(m) { + let et = d.mul(&dnfull); + let ets = et.shift(t + 1, 1); + e.push(et); + es.push(ets); + } + // block[t] = e[t] * es[t] * D_t — the m=2 module's `e1*e1s*dj`/`e2*e2s*dk`. + let block: Vec = (0..m).map(|t| e[t].mul(&es[t]).mul(ds[t])).collect(); + + // mn[i]: the multiplier for a_i(n)'s column, built as `D_total / D_{n,i}` + // by replacing exactly one of the two `dnfull` copies inside `block[0]` + // with `dnfull_excl[i]` (any one of `D_total`'s `2m` redundant copies of + // `D_{n,i}` may be the one "removed" — always using axis 0's block for + // this, mirroring the m=2 module's asymmetric-but-valid choice of always + // using `e1`/`e1s`/`dj`). + let mut mn: Vec = Vec::with_capacity(order + 1); for excl in &dnfull_excl { - let e1_excl_i = dj.mul(excl); - mn.push(e1_excl_i.mul(&e1s).mul(dj).mul(&e2).mul(&e2s).mul(dk)); + let mut p = ds[0].mul(excl).mul(&es[0]).mul(ds[0]); + for blk in block.iter().skip(1) { + p = p.mul(blk); + } + mn.push(p); + } + // mt[t]: the multiplier for P_t's columns, `D_total / block[t]` — the + // product of every *other* axis's block. + let mut mt: Vec = Vec::with_capacity(m); + for t in 0..m { + let mut p = PolyM::one(num_axes); + for (tp, blk) in block.iter().enumerate() { + if tp != t { + p = p.mul(blk); + } + } + mt.push(p); } - let mj = e2.mul(&e2s).mul(dk); - let mk = e1.mul(&e1s).mul(dj); - // Unknown layout: a_{i,p} for i=0..=order, p=0..=a_degree; then P1 - // coefficients over the box n<=cert_degree, j<=cert_degree, k<=cert_degree; - // then P2 coefficients over the same box. let na = a_degree + 1; let box_len = cert_degree + 1; - let cert_count = box_len * box_len * box_len; + let cert_box_count = box_len.pow(num_axes as u32); let a_count = (order + 1) * na; - let total = a_count + 2 * cert_count; + let total = a_count + m * cert_box_count; let idx_a = |i: usize, p: usize| i * na + p; - let idx_p1 = |p: usize, q: usize, r: usize| a_count + (p * box_len + q) * box_len + r; - let idx_p2 = - |p: usize, q: usize, r: usize| a_count + cert_count + (p * box_len + q) * box_len + r; + let idx_cert = |t: usize, exps: &[usize]| a_count + t * cert_box_count + flatten(exps, box_len); - // Monomial -> row (dense-in-practice sparse map of coefficient vectors). - let mut rows: BTreeMap<(u32, u32, u32), Vec> = BTreeMap::new(); - let mut add_contribution = |mono: &BTreeMap<(u32, u32, u32), Rational>, col: usize| { + let mut rows: BTreeMap, Vec> = BTreeMap::new(); + let mut add_contribution = |mono: &BTreeMap, Rational>, col: usize| { for (exp, c) in mono { let row = rows - .entry(*exp) + .entry(exp.clone()) .or_insert_with(|| vec![Rational::from(0); total]); row[col] += c.clone(); } @@ -289,71 +527,82 @@ fn solve_ansatz( for i in 0..=order { for p in 0..=a_degree { - let basis = Poly3::var(Axis::N) + let basis = PolyM::var(AXIS_N, num_axes) .pow_u32(p as u32) .mul(&nn[i]) .mul(&mn[i]); add_contribution(&basis.terms, idx_a(i, p)); } } - // Memoized (j+1)^q and (k+1)^r. - let mut jp1_pow: Vec = Vec::with_capacity(box_len); - let mut kp1_pow: Vec = Vec::with_capacity(box_len); - { - let jp1 = Poly3::var(Axis::J).add(&Poly3::one()); - let kp1 = Poly3::var(Axis::K).add(&Poly3::one()); - let mut acc_j = Poly3::one(); - let mut acc_k = Poly3::one(); + + // Memoized plain (`x_t^q`) and shifted (`(x_t+1)^q`) powers, per bound + // axis, for q = 0..=cert_degree; and plain `n^p` powers. + let mut plain_pow: Vec> = Vec::with_capacity(m); + let mut shift_pow: Vec> = Vec::with_capacity(m); + for t in 0..m { + let axis = t + 1; + let xv = PolyM::var(axis, num_axes); + let xv1 = xv.add(&PolyM::one(num_axes)); + let mut pp = Vec::with_capacity(box_len); + let mut sp = Vec::with_capacity(box_len); + let mut accp = PolyM::one(num_axes); + let mut accs = PolyM::one(num_axes); for _ in 0..box_len { - jp1_pow.push(acc_j.clone()); - kp1_pow.push(acc_k.clone()); - acc_j = acc_j.mul(&jp1); - acc_k = acc_k.mul(&kp1); + pp.push(accp.clone()); + sp.push(accs.clone()); + accp = accp.mul(&xv); + accs = accs.mul(&xv1); } + plain_pow.push(pp); + shift_pow.push(sp); } - for p in 0..=cert_degree { - let np = Poly3::var(Axis::N).pow_u32(p as u32); - for q in 0..=cert_degree { - for r in 0..=cert_degree { - let kr = Poly3::var(Axis::K).pow_u32(r as u32); - // P1(n,j+1,k) term: n^p (j+1)^q k^r ; P1(n,j,k) term: n^p j^q k^r. - // P1(n,j+1,k)*Nj*E1 - P1(n,j,k)*E1s*Dj, times Mj — see the - // derivation in the module docs: this is - // `[c1(j+1)*rho_j - c1(j)] * (E1*E1s*Dj) * Mj`. - let mono_shift = np.mul(&jp1_pow[q]).mul(&kr); - let mono_plain = np.mul(&Poly3::var(Axis::J).pow_u32(q as u32)).mul(&kr); - let basis1 = mono_shift - .mul(nj) - .mul(&e1) - .sub(&mono_plain.mul(&e1s).mul(dj)) - .mul(&mj) - .neg(); - add_contribution(&basis1.terms, idx_p1(p, q, r)); - - let jq = Poly3::var(Axis::J).pow_u32(q as u32); - let mono_shift_k = np.mul(&jq).mul(&kp1_pow[r]); - let mono_plain_k = np.mul(&jq).mul(&kr); - let basis2 = mono_shift_k - .mul(nk) - .mul(&e2) - .sub(&mono_plain_k.mul(&e2s).mul(dk)) - .mul(&mk) - .neg(); - add_contribution(&basis2.terms, idx_p2(p, q, r)); + let n_pow: Vec = { + let nv = PolyM::var(AXIS_N, num_axes); + let mut acc = PolyM::one(num_axes); + let mut v = Vec::with_capacity(box_len); + for _ in 0..box_len { + v.push(acc.clone()); + acc = acc.mul(&nv); + } + v + }; + + let total_combos = box_len.pow(num_axes as u32); + for combo in 0..total_combos { + let exps = unflatten(combo, num_axes, box_len); + let p = exps[0]; + let np = n_pow[p].clone(); + for t in 0..m { + let mut mono_shift = np.clone(); + let mut mono_plain = np.clone(); + for tp in 0..m { + let q = exps[tp + 1]; + mono_plain = mono_plain.mul(&plain_pow[tp][q]); + mono_shift = mono_shift.mul(if tp == t { + &shift_pow[tp][q] + } else { + &plain_pow[tp][q] + }); } + // [c_t(x_t+1)*rho_t - c_t] * (D_total/block[t]) — see the module + // docs derivation, generalizing the m=2 module's `basis1`/`basis2`. + let basis = mono_shift + .mul(ns_[t]) + .mul(&e[t]) + .sub(&mono_plain.mul(&es[t]).mul(ds[t])) + .mul(&mt[t]) + .neg(); + add_contribution(&basis.terms, idx_cert(t, &exps)); } } let matrix: Vec> = rows.into_values().collect(); - let basis = rational_nullspace(matrix, total); - if basis.is_empty() { + let basis_vecs = rational_nullspace(matrix, total); + if basis_vecs.is_empty() { return Ok(None); } - for vec in &basis { - // Reject a basis vector whose top a_order(n) is identically zero: - // that is really a lower-order (or no) relation, and returning it - // under this `order` would misreport the recurrence's order. + for vec in &basis_vecs { let top_nonzero = (0..na).any(|p| vec[idx_a(order, p)] != 0); if !top_nonzero { continue; @@ -363,60 +612,45 @@ fn solve_ansatz( let mut a = Vec::with_capacity(order + 1); for i in 0..=order { - let mut p = Poly3::zero(); + let mut poly = PolyM::zero(); for pw in 0..=a_degree { let c = scaled[idx_a(i, pw)].clone(); if c != 0 { - p = p.add(&Poly3::var(Axis::N).pow_u32(pw as u32).scale(&c)); + poly = poly.add(&PolyM::var(AXIS_N, num_axes).pow_u32(pw as u32).scale(&c)); } } - a.push(p); + a.push(poly); } - let mut p1 = Poly3::zero(); - let mut p2 = Poly3::zero(); - for p in 0..=cert_degree { - for q in 0..=cert_degree { - for r in 0..=cert_degree { - let c1c = scaled[idx_p1(p, q, r)].clone(); - if c1c != 0 { - p1 = p1.add(&monomial3(p as u32, q as u32, r as u32, c1c)); - } - let c2c = scaled[idx_p2(p, q, r)].clone(); - if c2c != 0 { - p2 = p2.add(&monomial3(p as u32, q as u32, r as u32, c2c)); - } + + let mut certs_num: Vec = vec![PolyM::zero(); m]; + for combo in 0..total_combos { + let exps = unflatten(combo, num_axes, box_len); + for (t, cn) in certs_num.iter_mut().enumerate() { + let c = scaled[idx_cert(t, &exps)].clone(); + if c != 0 { + *cn = cn.add(&monomial_m(&exps, c, num_axes)); } } } - let c1 = Rat3 { - num: p1, - den: e1.clone(), - }; - let c2 = Rat3 { - num: p2, - den: e2.clone(), - }; - return Ok(Some(Candidate { order, a, c1, c2 })); + let certs: Vec = (0..m) + .map(|t| RatM { + num: certs_num[t].clone(), + den: e[t].clone(), + }) + .collect(); + return Ok(Some(CandidateMd { order, a, certs })); } Ok(None) } -fn monomial3(en: u32, ej: u32, ek: u32, c: Rational) -> Poly3 { - Poly3::var(Axis::N) - .pow_u32(en) - .mul(&Poly3::var(Axis::J).pow_u32(ej)) - .mul(&Poly3::var(Axis::K).pow_u32(ek)) - .scale(&c) -} - /// Scale a flat family of rationals to integer, content-1 (up to sign) — /// the same principle as [`super::super::qfield::make_primitive`], applied -/// to one combined vector spanning `a_i(n)` *and* both certificate -/// numerators, because a homogeneous linear relation stays a solution under +/// to one combined vector spanning `a_i(n)` *and* every certificate +/// numerator, because a homogeneous linear relation stays a solution under /// any single overall rescaling. Kept local rather than reusing /// `make_primitive` because that helper is specialized to a family of -/// single-variable `RatUniPoly`s; here the family spans three different -/// polynomial rings (one univariate, two trivariate) sharing only their +/// single-variable `RatUniPoly`s; here the family spans `1 + m` different +/// polynomial rings (one univariate, `m` multivariate) sharing only their /// scalar coefficients. fn primitive_scale_rationals(v: &mut [Rational]) { let mut den_lcm = Integer::from(1); @@ -448,10 +682,6 @@ fn primitive_scale_rationals(v: &mut [Rational]) { /// Nullspace basis of `mat` (rows = equations, `ncols` unknowns) over `Q`, by /// plain Gaussian elimination to row-echelon form. -// The pivot column `col` bounds the range read/written below; a plain slice -// iterator would need to skip a variable, non-slice-aligned prefix on both -// `mat[row]` and every other row simultaneously, so the index form is -// clearer here than clippy's suggested rewrite. #[allow(clippy::needless_range_loop)] fn rational_nullspace(mut mat: Vec>, ncols: usize) -> Vec> { let nrows = mat.len(); @@ -511,39 +741,43 @@ fn rational_nullspace(mut mat: Vec>, ncols: usize) -> Vec bool { - let mut lhs = Rat3::zero(); +fn verify_certificate_md(f: &ProperTermM, cand: &CandidateMd) -> bool { + let num_axes = f.m + 1; + let mut lhs = RatM::from_rational(Rational::from(0), num_axes); for (i, ai) in cand.a.iter().enumerate() { - let Ok(ratio) = f.ratio_axis(Axis::N, i as i64) else { + let Ok(ratio) = f.ratio_axis(AXIS_N, i as i64) else { return false; }; - lhs = lhs.add(&Rat3::from_poly(ai.clone()).mul(&ratio)); + lhs = lhs.add(&RatM::from_poly(ai.clone(), num_axes).mul(&ratio)); + } + let mut rhs = RatM::from_rational(Rational::from(0), num_axes); + for (t, cert_t) in cand.certs.iter().enumerate() { + let axis: Axis = t + 1; + let Ok(rho_t) = f.ratio_axis(axis, 1) else { + return false; + }; + let g_delta = cert_t.shift(axis, 1).mul(&rho_t).sub(cert_t); + rhs = rhs.add(&g_delta); } - let Ok(rho_j) = f.ratio_axis(Axis::J, 1) else { - return false; - }; - let Ok(rho_k) = f.ratio_axis(Axis::K, 1) else { - return false; - }; - let g1_delta = cand.c1.shift(Axis::J, 1).mul(&rho_j).sub(&cand.c1); - let g2_delta = cand.c2.shift(Axis::K, 1).mul(&rho_k).sub(&cand.c2); - let rhs = g1_delta.add(&g2_delta); lhs.eq_rat(&rhs) } -fn finish( - cand: Candidate, +fn finish_md( + cand: CandidateMd, n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, -) -> Telescoping2dResult { - let coeffs = cand.a.iter().map(|p| p.to_expr(pool, n, j, k)).collect(); - Telescoping2dResult { +) -> TelescopingMdResult { + let coeffs = cand.a.iter().map(|p| p.to_expr(pool, n, indices)).collect(); + let certs = cand + .certs + .iter() + .map(|c| c.to_expr(pool, n, indices)) + .collect(); + TelescopingMdResult { order: cand.order, coeffs, - cert1: cand.c1.to_expr(pool, n, j, k), - cert2: cand.c2.to_expr(pool, n, j, k), + certs, } } @@ -590,4 +824,30 @@ mod tests { telescope2d_search(n, n, n, n, &pool, &opts).expect_err("n, j, k must be distinct"); assert!(matches!(err, Telescoping2dError::InvalidInput(_))); } + + #[test] + fn md_refuses_zero_indices() { + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let opts = TelescopingMdOpts::default(); + let err = telescope_md_search(n, n, &[], &pool, &opts) + .expect_err("at least one bound index required"); + assert!(matches!(err, Telescoping2dError::InvalidInput(_))); + } + + /// `m = 1` should degenerate cleanly to the same shape as the classical + /// single-sum engine's easiest examples: `F(n,k) = C(n,k)`, order-1 + /// certificate, `Σ_k C(n,k) = 2^n`. + #[test] + fn md_single_index_matches_classical_binomial_sum() { + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let k = pool.symbol("k", Domain::Real); + let f = binom(&pool, n, k); + let opts = TelescopingMdOpts::default(); + let result = telescope_md_search(f, n, &[k], &pool, &opts) + .expect("a certificate should be found for C(n,k)"); + assert_eq!(result.certs.len(), 1); + assert!(result.order >= 1); + } } diff --git a/alkahest-core/src/holonomic/telescoping2d/term.rs b/alkahest-core/src/holonomic/telescoping2d/term.rs index c2d1436f..d8bf2eae 100644 --- a/alkahest-core/src/holonomic/telescoping2d/term.rs +++ b/alkahest-core/src/holonomic/telescoping2d/term.rs @@ -1,95 +1,111 @@ -//! Recognising *proper hypergeometric terms* `F(n, j, k)` in three affine -//! indices, and computing their exact shift ratios. +//! Recognising *proper hypergeometric terms* `F(n, x_1, …, x_m)` in `m + 1` +//! affine indices, for an arbitrary number `m ≥ 1` of bound indices, and +//! computing their exact shift ratios. //! -//! This is [`super::super::hyperterm`] generalized from two indices to three: +//! This is [`super::super::hyperterm`] generalized from one bound index to +//! `m`: //! //! ```text -//! F(n, j, k) = R(n,j,k) · z_j^j · z_k^k · w^n · ∏_i Γ(a_i·n + b_i·j + c_i·k + d_i)^(e_i) +//! F(n, x_1, …, x_m) = R(n,x) · ∏_t z_t^{x_t} · w^n · ∏_i Γ(a_i·n + Σ_t b_{i,t}·x_t + d_i)^(e_i) //! ``` //! -//! with `R ∈ Q(n,j,k)`, `z_j, z_k, w ∈ Q \ {0}`, `a_i, b_i, c_i ∈ Z`, -//! `d_i ∈ Q`, `e_i ∈ Z`. This is exactly the class for which the three shift -//! quotients `F(n+1,j,k)/F(n,j,k)`, `F(n,j+1,k)/F(n,j,k)` and -//! `F(n,j,k+1)/F(n,j,k)` — and, more generally, any `F(·+i,·,·)/F(·,·,·)` in a -//! single index — are rational functions in `Q(n,j,k)`, which is what the -//! double-sum ansatz search in [`super::search`] consumes. +//! with `R ∈ Q(n,x)`, `z_t, w ∈ Q \ {0}`, `a_i, b_{i,t} ∈ Z`, `d_i ∈ Q`, +//! `e_i ∈ Z`. This is exactly the class for which every shift quotient +//! `F(·+1 in one axis, ·) / F(·)` — and, more generally, any +//! `F(n+i,x)/F(n,x)` — is a rational function in `Q(n,x_1,…,x_m)`, which is +//! what the ansatz search in [`super::search`] consumes. //! //! The parser deliberately reuses the supported-function-head set of //! [`super::super::hyperterm`] (`gamma`, `factorial`, `binomial`, //! `pochhammer`) and its strictness: anything outside the class above is //! refused rather than approximated. -use super::poly::{Axis, Poly3, Rat3}; +use super::poly::{Axis, PolyM, RatM, AXIS_N}; use super::Telescoping2dError; use crate::kernel::{ExprData, ExprId, ExprPool}; use rug::Rational; -/// One `Γ(a·n + b·j + c·k + d)^e` factor. +/// One `Γ(a·n + Σ_t b_t·x_t + d)^e` factor. `coeffs[0]` is the coefficient of +/// `n`; `coeffs[1..]` are the coefficients of the `m` bound indices, in order +/// — length `m + 1`, matching [`super::poly`]'s axis convention. #[derive(Clone, Debug, PartialEq, Eq)] -pub struct GammaFactor3 { - pub a: i64, - pub b: i64, - pub c: i64, +pub struct GammaFactorM { + pub coeffs: Vec, pub d: Rational, pub e: i32, } -fn gamma_arg_poly3(g: &GammaFactor3) -> Poly3 { - Poly3::var(Axis::N) - .scale(&Rational::from(g.a)) - .add(&Poly3::var(Axis::J).scale(&Rational::from(g.b))) - .add(&Poly3::var(Axis::K).scale(&Rational::from(g.c))) - .add(&Poly3::constant(g.d.clone())) +fn gamma_arg_polym(g: &GammaFactorM, num_axes: usize) -> PolyM { + let mut p = PolyM::constant(g.d.clone(), num_axes); + for (axis, &c) in g.coeffs.iter().enumerate() { + if c != 0 { + p = p.add(&PolyM::var(axis, num_axes).scale(&Rational::from(c))); + } + } + p } -/// A parsed proper hypergeometric term `F(n, j, k)`. +/// A parsed proper hypergeometric term `F(n, x_1, …, x_m)`. #[derive(Clone, Debug)] -pub struct ProperTerm3 { - pub rat: Rat3, - pub zj: Rational, - pub zk: Rational, +pub struct ProperTermM { + /// Number of bound indices. + pub m: usize, + pub rat: RatM, pub w: Rational, - pub gammas: Vec, + /// `z_t` for `t = 0..m`, one per bound index (in caller order). + pub z: Vec, + pub gammas: Vec, } const MAX_POW: i32 = 32; const MAX_PARSE_DEPTH: usize = 64; const MAX_FOLD_BITS: u32 = 4096; -impl ProperTerm3 { - fn one() -> Self { - ProperTerm3 { - rat: Rat3::one(), - zj: Rational::from(1), - zk: Rational::from(1), +impl ProperTermM { + fn num_axes(&self) -> usize { + self.m + 1 + } + + fn one(m: usize) -> Self { + ProperTermM { + m, + rat: RatM::one(m + 1), w: Rational::from(1), + z: vec![Rational::from(1); m], gammas: Vec::new(), } } - fn from_rat3(r: Rat3) -> Self { - ProperTerm3 { + fn from_ratm(r: RatM, m: usize) -> Self { + ProperTermM { + m, rat: r, - zj: Rational::from(1), - zk: Rational::from(1), w: Rational::from(1), + z: vec![Rational::from(1); m], gammas: Vec::new(), } } - fn mul(&self, other: &ProperTerm3) -> ProperTerm3 { + fn mul(&self, other: &ProperTermM) -> ProperTermM { + debug_assert_eq!(self.m, other.m); let mut gammas = self.gammas.clone(); gammas.extend(other.gammas.iter().cloned()); - ProperTerm3 { + let z = self + .z + .iter() + .zip(other.z.iter()) + .map(|(a, b)| a.clone() * b.clone()) + .collect(); + ProperTermM { + m: self.m, rat: self.rat.mul(&other.rat), - zj: self.zj.clone() * other.zj.clone(), - zk: self.zk.clone() * other.zk.clone(), w: self.w.clone() * other.w.clone(), + z, gammas, } } - fn pow(&self, e: i32) -> Option { + fn pow(&self, e: i32) -> Option { if e.unsigned_abs() > MAX_POW as u32 { return None; } @@ -97,51 +113,51 @@ impl ProperTerm3 { .gammas .iter() .map(|g| { - Some(GammaFactor3 { - a: g.a, - b: g.b, - c: g.c, + Some(GammaFactorM { + coeffs: g.coeffs.clone(), d: g.d.clone(), e: g.e.checked_mul(e)?, }) }) .collect::>>()?; - Some(ProperTerm3 { + let z = self + .z + .iter() + .map(|zt| rat_pow(zt, e)) + .collect::>>()?; + Some(ProperTermM { + m: self.m, rat: self.rat.pow_i32(e)?, - zj: rat_pow(&self.zj, e)?, - zk: rat_pow(&self.zk, e)?, w: rat_pow(&self.w, e)?, + z, gammas, }) } - /// `F(·+i·axis, ·, ·) / F(·,·,·)` — the shift quotient in a single index, - /// as an exact element of `Q(n,j,k)`. - pub fn ratio_axis(&self, axis: Axis, i: i64) -> Result { + /// `F(·+i·axis, ·) / F(·)` — the shift quotient in a single axis (`0` for + /// `n`, `1..=m` for a bound index), as an exact element of + /// `Q(n,x_1,…,x_m)`. + pub fn ratio_axis(&self, axis: Axis, i: i64) -> Result { if i == 0 { - return Ok(Rat3::one()); + return Ok(RatM::one(self.num_axes())); } let shifted = self.rat.shift(axis, i); let mut acc = shifted.div(&self.rat).ok_or_else(|| { Telescoping2dError::NotProperHypergeometric("term vanishes identically".into()) })?; - let base = match axis { - Axis::N => &self.w, - Axis::J => &self.zj, - Axis::K => &self.zk, + let base = if axis == AXIS_N { + &self.w + } else { + &self.z[axis - 1] }; - acc = acc.mul(&Rat3::from_rational(rat_pow_i64(base, i)?)); + acc = acc.mul(&RatM::from_rational(rat_pow_i64(base, i)?, self.num_axes())); for g in &self.gammas { - let arg = gamma_arg_poly3(g); - let coeff = match axis { - Axis::N => g.a, - Axis::J => g.b, - Axis::K => g.c, - }; + let arg = gamma_arg_polym(g, self.num_axes()); + let coeff = g.coeffs[axis]; let shift = coeff.checked_mul(i).ok_or_else(|| { Telescoping2dError::SearchExhausted("gamma shift overflow".into()) })?; - let step = gamma_shift_ratio3(&arg, shift)?; + let step = gamma_shift_ratiom(&arg, shift, self.num_axes())?; acc = acc.mul(&step.pow_i32(g.e).ok_or_else(|| { Telescoping2dError::NotProperHypergeometric( "gamma factor is identically zero".into(), @@ -154,35 +170,39 @@ impl ProperTerm3 { pub fn parse( expr: ExprId, n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, - ) -> Result { - parse_rec(expr, n, j, k, pool, 0) + ) -> Result { + if indices.is_empty() { + return Err(Telescoping2dError::InvalidInput( + "at least one bound index is required".into(), + )); + } + parse_rec(expr, n, indices, pool, 0) } } /// `Γ(x + s) / Γ(x)` for integer `s`, as an exact rational function. -fn gamma_shift_ratio3(x: &Poly3, s: i64) -> Result { +fn gamma_shift_ratiom(x: &PolyM, s: i64, num_axes: usize) -> Result { if s == 0 { - return Ok(Rat3::one()); + return Ok(RatM::one(num_axes)); } if s.unsigned_abs() > 512 { return Err(Telescoping2dError::SearchExhausted(format!( "gamma argument shift {s} exceeds the supported limit of 512" ))); } - let mut prod = Poly3::one(); + let mut prod = PolyM::one(num_axes); if s > 0 { for t in 0..s { - prod = prod.mul(&x.add(&Poly3::from_i64(t))); + prod = prod.mul(&x.add(&PolyM::from_i64(t, num_axes))); } - Ok(Rat3::from_poly(prod)) + Ok(RatM::from_poly(prod, num_axes)) } else { for t in 1..=(-s) { - prod = prod.mul(&x.add(&Poly3::from_i64(-t))); + prod = prod.mul(&x.add(&PolyM::from_i64(-t, num_axes))); } - Rat3::from_poly(prod) + RatM::from_poly(prod, num_axes) .inv() .ok_or_else(|| Telescoping2dError::NotProperHypergeometric("gamma pole".into())) } @@ -220,30 +240,30 @@ fn rat_pow_i64(q: &Rational, e: i64) -> Result { fn parse_rec( expr: ExprId, n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, depth: usize, -) -> Result { +) -> Result { + let m = indices.len(); if depth > MAX_PARSE_DEPTH { return Err(Telescoping2dError::NotProperHypergeometric( "expression nests deeper than the parser supports".into(), )); } - if let Some(r) = as_rat3(expr, n, j, k, pool, 0) { - return Ok(ProperTerm3::from_rat3(r)); + if let Some(r) = as_ratm(expr, n, indices, pool, 0) { + return Ok(ProperTermM::from_ratm(r, m)); } match pool.get(expr) { ExprData::Mul(args) => { - let mut acc = ProperTerm3::one(); + let mut acc = ProperTermM::one(m); for a in args { - acc = acc.mul(&parse_rec(a, n, j, k, pool, depth + 1)?); + acc = acc.mul(&parse_rec(a, n, indices, pool, depth + 1)?); } Ok(acc) } ExprData::Pow { base, exp } => { if let Some(e) = as_i32(exp, pool) { - let b = parse_rec(base, n, j, k, pool, depth + 1)?; + let b = parse_rec(base, n, indices, pool, depth + 1)?; return b.pow(e).ok_or_else(|| { Telescoping2dError::NotProperHypergeometric(format!( "exponent {e} is outside the supported range (|e| <= {MAX_POW})" @@ -261,13 +281,12 @@ fn parse_rec( "0 raised to a symbolic power".into(), )); } - let (alpha, beta, delta, gamma_c) = - affine_parts3(exp, n, j, k, pool).ok_or_else(|| { - Telescoping2dError::NotProperHypergeometric(format!( - "exponent {} is not integer-affine in the three indices", - pool.display(exp) - )) - })?; + let (coeffs, gamma_c) = affine_partsm(exp, n, indices, pool).ok_or_else(|| { + Telescoping2dError::NotProperHypergeometric(format!( + "exponent {} is not integer-affine in the bound indices", + pool.display(exp) + )) + })?; if *gamma_c.clone().denom() != 1 { return Err(Telescoping2dError::NotProperHypergeometric( "constant part of an exponential exponent must be an integer".into(), @@ -277,15 +296,19 @@ fn parse_rec( .numer() .to_i64() .ok_or_else(|| Telescoping2dError::SearchExhausted("exponent too large".into()))?; - Ok(ProperTerm3 { - rat: Rat3::from_rational(rat_pow_i64(&c, gi)?), - zj: rat_pow_i64(&c, beta)?, - zk: rat_pow_i64(&c, delta)?, - w: rat_pow_i64(&c, alpha)?, + let z = coeffs[1..] + .iter() + .map(|&bt| rat_pow_i64(&c, bt)) + .collect::, _>>()?; + Ok(ProperTermM { + m, + rat: RatM::from_rational(rat_pow_i64(&c, gi)?, m + 1), + w: rat_pow_i64(&c, coeffs[0])?, + z, gammas: Vec::new(), }) } - ExprData::Func { name, args } => parse_func(&name, &args, n, j, k, pool), + ExprData::Func { name, args } => parse_func(&name, &args, n, indices, pool), other => Err(Telescoping2dError::NotProperHypergeometric(format!( "unsupported node {other:?} in {}", pool.display(expr) @@ -297,35 +320,30 @@ fn parse_func( name: &str, args: &[ExprId], n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, -) -> Result { - let gamma_of = |arg: ExprId, e: i32| -> Result { - let (a, b, c, d) = affine_parts3(arg, n, j, k, pool).ok_or_else(|| { +) -> Result { + let m = indices.len(); + let gamma_of = |arg: ExprId, e: i32| -> Result { + let (coeffs, d) = affine_partsm(arg, n, indices, pool).ok_or_else(|| { Telescoping2dError::NotProperHypergeometric(format!( - "gamma argument {} is not integer-affine in the three indices", + "gamma argument {} is not integer-affine in the bound indices", pool.display(arg) )) })?; - Ok(GammaFactor3 { a, b, c, d, e }) + Ok(GammaFactorM { coeffs, d, e }) }; let one_plus = |arg: ExprId| -> ExprId { pool.add(vec![arg, pool.integer(1_i32)]) }; + let plain = |gammas: Vec| ProperTermM { + m, + rat: RatM::one(m + 1), + w: Rational::from(1), + z: vec![Rational::from(1); m], + gammas, + }; match (name, args.len()) { - ("gamma", 1) => Ok(ProperTerm3 { - rat: Rat3::one(), - zj: Rational::from(1), - zk: Rational::from(1), - w: Rational::from(1), - gammas: vec![gamma_of(args[0], 1)?], - }), - ("factorial", 1) => Ok(ProperTerm3 { - rat: Rat3::one(), - zj: Rational::from(1), - zk: Rational::from(1), - w: Rational::from(1), - gammas: vec![gamma_of(one_plus(args[0]), 1)?], - }), + ("gamma", 1) => Ok(plain(vec![gamma_of(args[0], 1)?])), + ("factorial", 1) => Ok(plain(vec![gamma_of(one_plus(args[0]), 1)?])), ("binomial", 2) => { let top = one_plus(args[0]); let bot = one_plus(args[1]); @@ -334,23 +352,15 @@ fn parse_func( pool.mul(vec![args[1], pool.integer(-1_i32)]), pool.integer(1_i32), ]); - Ok(ProperTerm3 { - rat: Rat3::one(), - zj: Rational::from(1), - zk: Rational::from(1), - w: Rational::from(1), - gammas: vec![gamma_of(top, 1)?, gamma_of(bot, -1)?, gamma_of(rest, -1)?], - }) + Ok(plain(vec![ + gamma_of(top, 1)?, + gamma_of(bot, -1)?, + gamma_of(rest, -1)?, + ])) } ("pochhammer", 2) => { let sum = pool.add(vec![args[0], args[1]]); - Ok(ProperTerm3 { - rat: Rat3::one(), - zj: Rational::from(1), - zk: Rational::from(1), - w: Rational::from(1), - gammas: vec![gamma_of(sum, 1)?, gamma_of(args[0], -1)?], - }) + Ok(plain(vec![gamma_of(sum, 1)?, gamma_of(args[0], -1)?])) } _ => Err(Telescoping2dError::NotProperHypergeometric(format!( "function `{name}/{}` is not part of the proper hypergeometric class \ @@ -360,42 +370,41 @@ fn parse_func( } } -/// Evaluate an expression inside the field `Q(n, j, k)`, or `None` if it -/// leaves it. -pub fn as_rat3( +/// Evaluate an expression inside the field `Q(n, x_1, …, x_m)`, or `None` if +/// it leaves it. +pub fn as_ratm( expr: ExprId, n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, depth: usize, -) -> Option { +) -> Option { if depth > MAX_PARSE_DEPTH { return None; } + let num_axes = indices.len() + 1; if expr == n { - return Some(Rat3::from_poly(Poly3::var(Axis::N))); - } - if expr == j { - return Some(Rat3::from_poly(Poly3::var(Axis::J))); + return Some(RatM::from_poly(PolyM::var(AXIS_N, num_axes), num_axes)); } - if expr == k { - return Some(Rat3::from_poly(Poly3::var(Axis::K))); + for (t, &idx_expr) in indices.iter().enumerate() { + if expr == idx_expr { + return Some(RatM::from_poly(PolyM::var(t + 1, num_axes), num_axes)); + } } match pool.get(expr) { - ExprData::Integer(i) => Some(Rat3::from_rational(Rational::from(i.0.clone()))), - ExprData::Rational(r) => Some(Rat3::from_rational(r.0.clone())), + ExprData::Integer(i) => Some(RatM::from_rational(Rational::from(i.0.clone()), num_axes)), + ExprData::Rational(r) => Some(RatM::from_rational(r.0.clone(), num_axes)), ExprData::Add(args) => { - let mut acc = Rat3::zero(); + let mut acc = RatM::from_rational(Rational::from(0), num_axes); for a in args { - acc = acc.add(&as_rat3(a, n, j, k, pool, depth + 1)?); + acc = acc.add(&as_ratm(a, n, indices, pool, depth + 1)?); } Some(acc) } ExprData::Mul(args) => { - let mut acc = Rat3::one(); + let mut acc = RatM::one(num_axes); for a in args { - acc = acc.mul(&as_rat3(a, n, j, k, pool, depth + 1)?); + acc = acc.mul(&as_ratm(a, n, indices, pool, depth + 1)?); } Some(acc) } @@ -404,7 +413,7 @@ pub fn as_rat3( if e.unsigned_abs() > MAX_POW as u32 { return None; } - as_rat3(base, n, j, k, pool, depth + 1)?.pow_i32(e) + as_ratm(base, n, indices, pool, depth + 1)?.pow_i32(e) } _ => None, } @@ -468,43 +477,49 @@ fn as_i32(expr: ExprId, pool: &ExprPool) -> Option { } } -/// Decompose an expression as `a·n + b·j + c·k + d` with `a, b, c ∈ Z` and -/// `d ∈ Q`. -pub fn affine_parts3( +/// Decompose an expression as `a·n + Σ_t b_t·x_t + d` with `a, b_t ∈ Z` and +/// `d ∈ Q`. Returns `(coeffs, d)` with `coeffs[0] = a`, `coeffs[1..] = b_t`. +pub fn affine_partsm( expr: ExprId, n: ExprId, - j: ExprId, - k: ExprId, + indices: &[ExprId], pool: &ExprPool, -) -> Option<(i64, i64, i64, Rational)> { - let r = as_rat3(expr, n, j, k, pool, 0)?; +) -> Option<(Vec, Rational)> { + let r = as_ratm(expr, n, indices, pool, 0)?; let den_c = r.den.as_constant()?; if den_c == 0 { return None; } let inv = Rational::from(1) / den_c; let num = r.num.scale(&inv); - // Total degree at most 1: every stored exponent triple sums to <= 1. - if num.terms.keys().any(|(en, ej, ek)| en + ej + ek > 1) { + let num_axes = indices.len() + 1; + // Total degree at most 1: every stored exponent tuple sums to <= 1. + if num.terms.keys().any(|e| e.iter().sum::() > 1) { return None; } - let coeff_of = |axis_e: (u32, u32, u32)| -> Rational { + let coeff_of = |axis: usize| -> Rational { + let mut e = vec![0u32; num_axes]; + e[axis] = 1; num.terms - .get(&axis_e) + .get(&e) .cloned() .unwrap_or_else(|| Rational::from(0)) }; - let a = coeff_of((1, 0, 0)); - let b = coeff_of((0, 1, 0)); - let c = coeff_of((0, 0, 1)); - let d = coeff_of((0, 0, 0)); - if *a.clone().denom() != 1 || *b.clone().denom() != 1 || *c.clone().denom() != 1 { - return None; - } - let a_i = a.numer().to_i64()?; - let b_i = b.numer().to_i64()?; - let c_i = c.numer().to_i64()?; - Some((a_i, b_i, c_i, d)) + let mut coeffs = Vec::with_capacity(num_axes); + for axis in 0..num_axes { + let c = coeff_of(axis); + if *c.clone().denom() != 1 { + return None; + } + coeffs.push(c.numer().to_i64()?); + } + let zero_e = vec![0u32; num_axes]; + let d = num + .terms + .get(&zero_e) + .cloned() + .unwrap_or_else(|| Rational::from(0)); + Some((coeffs, d)) } #[cfg(test)] @@ -532,15 +547,15 @@ mod tests { let pool = ExprPool::new(); let (n, j, k) = njk(&pool); let f = pool.mul(vec![binom(&pool, n, j), binom(&pool, j, k)]); - let term = ProperTerm3::parse(f, n, j, k, &pool).expect("proper hypergeometric"); + let term = ProperTermM::parse(f, n, &[j, k], &pool).expect("proper hypergeometric"); - // ratio_j = F(n,j+1,k)/F(n,j,k) should equal + // ratio_j = F(n,j+1,k)/F(n,j,k) (axis 1 is j) should equal // (n-j)/(j+1) * (j+1)!/(j+1-k)! / (j!/(j-k)!) as a rational function; // check it numerically at a sample point instead of re-deriving the // closed form, since that *is* what the exact machinery computes. - let rj = term.ratio_axis(Axis::J, 1).expect("ratio_j"); - let (nn, jj, kk) = (Rational::from(6), Rational::from(3), Rational::from(2)); - let got = rj.num.eval(&nn, &jj, &kk) / rj.den.eval(&nn, &jj, &kk); + let rj = term.ratio_axis(1, 1).expect("ratio_j"); + let vals = [Rational::from(6), Rational::from(3), Rational::from(2)]; + let got = rj.num.eval(&vals) / rj.den.eval(&vals); // C(6,4)*C(4,2) / (C(6,3)*C(3,2)) = 15*6 / (20*3) = 90/60 = 3/2 assert_eq!(got, Rational::from((3, 2))); } @@ -550,10 +565,39 @@ mod tests { let pool = ExprPool::new(); let (n, j, k) = njk(&pool); let bad = pool.func("sin", vec![pool.mul(vec![n, j, k])]); - let err = ProperTerm3::parse(bad, n, j, k, &pool).expect_err("not hypergeometric"); + let err = ProperTermM::parse(bad, n, &[j, k], &pool).expect_err("not hypergeometric"); assert!(matches!( err, Telescoping2dError::NotProperHypergeometric(_) )); } + + /// `m = 3` bound indices: `F(n,x,y,z) = C(n,x)*C(x,y)*C(y,z)`, a chain of + /// three coupled binomial transforms — exercises a `num_axes = 4` parse. + #[test] + fn three_index_chain_parses_and_has_exact_ratios() { + let pool = ExprPool::new(); + let n = pool.symbol("n", Domain::Real); + let x = pool.symbol("x", Domain::Real); + let y = pool.symbol("y", Domain::Real); + let z = pool.symbol("z", Domain::Real); + let f = pool.mul(vec![ + binom(&pool, n, x), + binom(&pool, x, y), + binom(&pool, y, z), + ]); + let term = ProperTermM::parse(f, n, &[x, y, z], &pool).expect("proper hypergeometric"); + assert_eq!(term.m, 3); + let rz = term.ratio_axis(3, 1).expect("ratio_z"); + // F(n,x,y,z+1)/F(n,x,y,z) = C(y,z+1)/C(y,z) at a sample point. + let vals = [ + Rational::from(8), + Rational::from(5), + Rational::from(4), + Rational::from(1), + ]; + let got = rz.num.eval(&vals) / rz.den.eval(&vals); + // C(4,2)/C(4,1) = 6/4 = 3/2 + assert_eq!(got, Rational::from((3, 2))); + } } diff --git a/alkahest-py/src/lib.rs b/alkahest-py/src/lib.rs index f45aef7a..49e65f17 100644 --- a/alkahest-py/src/lib.rs +++ b/alkahest-py/src/lib.rs @@ -190,9 +190,11 @@ use alkahest_core::holonomic::qzeil::{ }; // M4 — double-sum (Apagodu–Zeilberger) creative telescoping use alkahest_core::holonomic::telescoping2d::{ - boundary_status_2d as core_boundary_status_2d, telescope2d_search as core_telescope2d_search, + boundary_status_2d as core_boundary_status_2d, boundary_status_md as core_boundary_status_md, + telescope2d_search as core_telescope2d_search, telescope_md_search as core_telescope_md_search, Telescoping2dError as CoreTelescoping2dError, Telescoping2dOpts as CoreTelescoping2dOpts, - Telescoping2dResult as CoreTelescoping2dResult, + Telescoping2dResult as CoreTelescoping2dResult, TelescopingMdOpts as CoreTelescopingMdOpts, + TelescopingMdResult as CoreTelescopingMdResult, }; // M6 — modular / p-adic evaluation of holonomic sequences use alkahest_core::holonomic::modular::{ @@ -5984,6 +5986,225 @@ fn py_telescope2d( }) } +/// A **verified** `m`-bound-index creative-telescoping certificate, returned +/// by :func:`alkahest.experimental.telescope_md`. The general form of +/// :class:`Telescoping2dCertificate` (`m = 2` is that class's special case, +/// unchanged and still returned by :func:`~alkahest.experimental.telescope2d` +/// itself). +/// +/// Carries the recurrence coefficients ``a_0(n), …, a_J(n)`` and `m` rational +/// certificates ``c_1, …, c_m`` satisfying, as an exact identity in +/// ``Q(n, x_1, …, x_m)``, ``Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t(c_t·F)``. That +/// identity is re-checked exactly before this object is constructed. It says +/// nothing on its own about the `m`-fold sum over a stated box — call +/// :meth:`boundary_status` to decide that, separately. +/// +/// See the Rust module docs (``alkahest_cas::holonomic::telescoping2d``) for +/// the complete, honestly-stated scope: proper hypergeometric summands only, +/// no general Wegschaider reduction, a fixed (non-minimal) certificate +/// denominator, constant-box-only boundary analysis, and the resource +/// ceilings (`MAX_ANSATZ_UNKNOWNS`, `MAX_CUMULATIVE_LARGE_PROBE_UNKNOWNS`) +/// that keep a search with no certificate in reach a fast, honest refusal +/// rather than an unbounded computation. +#[pyclass(name = "TelescopingMdCertificate")] +struct PyTelescopingMdCertificate { + order: usize, + coeff_ids: Vec, + cert_ids: Vec, + pool: Py, + term_id: ExprId, + n_id: ExprId, + index_ids: Vec, + result: CoreTelescopingMdResult, +} + +#[pymethods] +impl PyTelescopingMdCertificate { + /// Recurrence order ``J``; ``len(coeffs()) == order + 1``. + #[getter] + fn order(&self) -> usize { + self.order + } + + /// ``[a_0(n), …, a_J(n)]`` — polynomial coefficients of the recurrence. + fn coeffs(&self, py: Python<'_>) -> Vec { + self.coeff_ids + .iter() + .map(|&id| PyExpr { + id, + pool: self.pool.clone_ref(py), + }) + .collect() + } + + /// ``[c_1(n,x), …, c_m(n,x)]``, one per bound index in the order they + /// were supplied to :func:`~alkahest.experimental.telescope_md`. + fn certs(&self, py: Python<'_>) -> Vec { + self.cert_ids + .iter() + .map(|&id| PyExpr { + id, + pool: self.pool.clone_ref(py), + }) + .collect() + } + + /// Decide the boundary hypothesis for the `m`-fold sum + /// ``S(n) = Σ_{x_1} … Σ_{x_m} F(n,x)`` over the box + /// ``x_t = limits[t][0] .. limits[t][1]``. + /// + /// ``limits`` must have exactly one ``(lo, hi)`` pair per bound index, + /// in the same order they were supplied to + /// :func:`~alkahest.experimental.telescope_md`; each bound is an + /// **integer constant** (``Expr`` or plain ``int``), never an expression + /// in ``n`` — see the Rust module docs for + /// ``telescoping2d::boundary`` for why, and for the standard workaround + /// when the natural range is `n`-dependent. + /// + /// Returns a ``dict`` with keys ``status`` (``"vanishes"``, ``"nonzero"`` + /// or ``"unknown"`` — though this version never produces ``"nonzero"``, + /// see the class docs), ``implies_sum_recurrence`` and + /// ``side_conditions``. + #[pyo3(signature = (limits))] + fn boundary_status( + &self, + py: Python<'_>, + limits: Vec<(Bound<'_, PyAny>, Bound<'_, PyAny>)>, + ) -> PyResult> { + if limits.len() != self.index_ids.len() { + return Err(PyValueError::new_err(format!( + "expected {} (lo, hi) limit pairs, one per bound index, got {}", + self.index_ids.len(), + limits.len() + ))); + } + let mut coerced: Vec<(ExprId, ExprId)> = Vec::with_capacity(limits.len()); + for (i, (lo, hi)) in limits.iter().enumerate() { + let lo_id = coerce_limit(py, &self.pool, lo, &format!("limits[{i}][0]"))?; + let hi_id = coerce_limit(py, &self.pool, hi, &format!("limits[{i}][1]"))?; + coerced.push((lo_id, hi_id)); + } + let (status, ranges) = { + let pool = self.pool.borrow(py); + let status = core_boundary_status_md( + &self.result, + self.term_id, + self.n_id, + &self.index_ids, + &coerced, + &pool.inner, + ); + let ranges: Vec = self + .index_ids + .iter() + .zip(coerced.iter()) + .map(|(&idx, &(lo, hi))| { + format!( + "{} = {}..{}", + pool.inner.display(idx), + pool.inner.display(lo), + pool.inner.display(hi) + ) + }) + .collect(); + (status, ranges) + }; + let out = PyDict::new_bound(py); + out.set_item("status", status.tag())?; + out.set_item("implies_sum_recurrence", status.implies_sum_recurrence())?; + out.set_item("side_conditions", status.side_conditions(&ranges))?; + Ok(out.unbind()) + } + + fn __repr__(&self, py: Python<'_>) -> String { + let pool = self.pool.borrow(py); + let coeffs: Vec = self + .coeff_ids + .iter() + .map(|&id| pool.inner.display(id).to_string()) + .collect(); + let certs: Vec = self + .cert_ids + .iter() + .map(|&id| pool.inner.display(id).to_string()) + .collect(); + format!( + "TelescopingMdCertificate(order={}, coeffs=[{}], certs=[{}])", + self.order, + coeffs.join(", "), + certs.join(", ") + ) + } +} + +/// `alkahest.experimental.telescope_md(term, n, indices, *, max_order=2, max_a_degree=2, max_cert_degree=2) -> TelescopingMdCertificate` +/// +/// Creative telescoping (Apagodu–Zeilberger) for a proper hypergeometric term +/// ``F(n, x_1, …, x_m)`` with an arbitrary number `m ≥ 1` of bound indices — +/// the general form of :func:`telescope2d` (`m = 2`), which remains the +/// semver-stable special case with its own dedicated function. Returns a +/// **verified** certificate: the identity +/// ``Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t(c_t·F)`` is re-checked exactly in +/// ``Q(n,x_1,…,x_m)`` before it is returned. +/// +/// ``indices`` is the list of bound-index symbols ``[x_1, …, x_m]``, in the +/// order the returned certificate's ``certs()`` and any +/// :meth:`~alkahest.experimental.TelescopingMdCertificate.boundary_status` +/// call use. The supported summand class, the certificate ansatz's degree +/// budget, the fixed (non-minimal) certificate denominator and the +/// constant-box-only boundary analysis are exactly :func:`telescope2d`'s, +/// generalized from two indices to `m`; see the Rust module docs +/// (``alkahest_cas::holonomic::telescoping2d``) for the complete, honestly- +/// stated scope, including the resource ceilings that keep a search with no +/// certificate in reach a fast refusal rather than an unbounded computation +/// as `m` or `max_cert_degree` grow — raising `m` or the certificate degree +/// bound grows the search space (and the risk of hitting those ceilings) +/// much faster than in the two-index case. +/// +/// Raises :exc:`alkahest.HolonomicError` (``E-HOLO-040`` outside the +/// supported class, ``E-HOLO-041`` when the bounded search is exhausted — +/// including when a resource ceiling, not genuine non-existence, is the +/// reason — ``E-HOLO-042`` for a malformed call, e.g. an empty ``indices`` +/// or a repeated symbol) rather than guessing. +#[allow(clippy::too_many_arguments)] +#[pyfunction] +#[pyo3( + name = "telescope_md", + signature = (term, n, indices, *, max_order = 2, max_a_degree = 2, max_cert_degree = 2) +)] +fn py_telescope_md( + py: Python<'_>, + term: PyRef, + n: PyRef, + indices: Vec>, + max_order: usize, + max_a_degree: usize, + max_cert_degree: usize, +) -> PyResult { + let pool_py = term.pool.clone_ref(py); + let index_ids: Vec = indices.iter().map(|e| e.id).collect(); + let opts = CoreTelescopingMdOpts { + max_order, + max_a_degree, + max_cert_degree, + }; + let result = { + let pool = pool_py.borrow(py); + core_telescope_md_search(term.id, n.id, &index_ids, &pool.inner, &opts) + .map_err(telescoping2d_error_to_py)? + }; + Ok(PyTelescopingMdCertificate { + order: result.order, + coeff_ids: result.coeffs.clone(), + cert_ids: result.certs.clone(), + pool: pool_py, + term_id: term.id, + n_id: n.id, + index_ids, + result, + }) +} + // --------------------------------------------------------------------------- // P1 item 10 — asymptotic expansion at scale // --------------------------------------------------------------------------- @@ -15035,6 +15256,9 @@ fn alkahest(m: &Bound<'_, PyModule>) -> PyResult<()> { // M4 — double-sum (Apagodu–Zeilberger) creative telescoping m.add_class::()?; m.add_function(wrap_pyfunction!(py_telescope2d, m)?)?; + // M4 extension — arbitrary-m-bound-index generalization + m.add_class::()?; + m.add_function(wrap_pyfunction!(py_telescope_md, m)?)?; // M6 — modular / p-adic evaluation of holonomic sequences m.add_class::()?; m.add_class::()?; diff --git a/alkahest-skill/alkahest.md b/alkahest-skill/alkahest.md index add9af5f..7df084ab 100644 --- a/alkahest-skill/alkahest.md +++ b/alkahest-skill/alkahest.md @@ -1178,7 +1178,7 @@ All errors inherit `AlkahestError` and carry `.code`, `.remediation`, `.span`. | `EigenError` | `E-EIGEN-*` | *Subclass of `MatrixError`.* Eigen/Jordan; defective matrix (`E-EIGEN-005`) | | `CadError` | `E-CAD-*` | **`decide` refused** — outside the fragment, or an untestable irrational boundary point | | `SosError` | `E-SOS-*` | No positivity certificate of this shape/degree (`E-SOS-002` — **a refusal: record `unknown`, not "not SOS"**); proved negative with a witness point (`E-SOS-003` — the only SOS verdict) | -| `HolonomicError` | `E-HOLO-*` | `zeilberger` outside the proper-hypergeometric class; `q_zeilberger` outside the `q`-hypergeometric one (`E-HOLO-020`) or with a non-rational shift quotient (`E-HOLO-024` — **permanent, not a bounds problem**); `telescope2d` outside the proper-hypergeometric-in-`(n,j,k)` class (`E-HOLO-040`), search exhausted (`E-HOLO-041`), or a malformed call (`E-HOLO-042` — `n`, `j`, `k` not distinct); `guess_holonomic` given too few terms to confirm a fit (`E-HOLO-005` — **a refusal: record `unknown`, not "no recurrence"**); `ModularRecurrence` / `binomial_mod` given an unsupported prime-power modulus (`E-HOLO-006`), a step with no `p`-adic integer answer (`E-HOLO-007` — **permanent**) or a working precision past `2**62` (`E-HOLO-008` — **resource: record `unknown`**) | +| `HolonomicError` | `E-HOLO-*` | `zeilberger` outside the proper-hypergeometric class; `q_zeilberger` outside the `q`-hypergeometric one (`E-HOLO-020`) or with a non-rational shift quotient (`E-HOLO-024` — **permanent, not a bounds problem**); `telescope2d`/`telescope_md` outside the proper-hypergeometric-in-the-bound-indices class (`E-HOLO-040`), search exhausted — including a resource ceiling refusal, see item 31 — (`E-HOLO-041`), or a malformed call (`E-HOLO-042` — indices not pairwise distinct, or empty); `guess_holonomic` given too few terms to confirm a fit (`E-HOLO-005` — **a refusal: record `unknown`, not "no recurrence"**); `ModularRecurrence` / `binomial_mod` given an unsupported prime-power modulus (`E-HOLO-006`), a step with no `p`-adic integer answer (`E-HOLO-007` — **permanent**) or a working precision past `2**62` (`E-HOLO-008` — **resource: record `unknown`**) | | `ValidatedError` | `E-VALIDATED-*` | Rigorous-bounds request unsupported / singular / malformed | | `OdeError` | `E-ODE-*` | ODE construction failed | | `DaeError` | `E-DAE-*` | DAE index reduction failed | @@ -1438,4 +1438,4 @@ reg.coverage_report_markdown() # same, rendered as a Markdown table 29. **`cert.specialize_at_root_of_unity(d, n)` is the decision that carries a `q_zeilberger` verdict to `q = ζ_d`, and it is three-valued** (since 3.9). A proved `Q(q)` recurrence does not by itself license setting `q` to a primitive `d`-th root of unity — a coefficient or a sum value can have a pole there, and specialising anyway is the `q`-analogue of the A279013 failure (item 22): a certificate that re-checks perfectly while the specialised claim is false. The hypotheses (no pole in any `a_i(qⁿ)` or `S(n+i)` at `ζ_d`) are decided **exactly**, by polynomial divisibility by `Φ_d(q)` over `Q` in the cyclotomic field `Q(ζ_d) = Q[q]/(Φ_d(q))` — never numerically — and `cyclotomic_polynomial(pool, d)` exposes `Φ_d(q)` itself so a caller can redo the check by hand. `status` is `"specializes"` (proved, and re-checked as an exact identity in `Q(ζ_d)` before being returned), `"obstructed"` (a pole was **exhibited** — `sum_value`/`coefficient` raise, but `sum_valuation(i)` is still available since the negative valuation *is* the obstruction — and this is not a claim the specialised identity is false, only that this route is blocked), or `"unknown"` (the generic boundary verdict was already `"unknown"`, so there is nothing to specialise). Three things a `"specializes"` verdict does **not** by itself mean, each with its own accessor: `is_vacuous` (every coefficient died — always true at `d = 1`, the `q → 1` limit — so the recurrence is `0 = 0`, still true, but empty), `leading_coefficient_survives` (`False` means the specialised recurrence no longer determines the last value from the earlier ones), and `support_shrinks` (`q`-Lucas killing terms — `[2;1]_q = 1 + q` is non-zero in `Q(q)` and zero at `ζ_2` — reported via `effective_support`, which can shrink but never grow). `sum_valuation(i)` is the `q`-supercongruence content itself: the exact integer `v` with `Φ_d(q)^v ∥ S(n+i)`, so `v ≥ r` is precisely `Φ_d(q)^r | S(n)`. 30. **`sos_decompose` tries the full PSD Gram cone and a Reznick multiplier search before refusing, and now certifies Motzkin and Robinson's form too** (since 3.9). Past diagonal dominance (`E-SOS-002` from DSOS alone) it searches the general PSD Gram cone, and past that — when `p` itself is not SOS — tries `(x_1²+…+x_n²)^N·p` for `N = 1..4` and searches *that* cone; a witness for `p < 0` still refuses separately with `E-SOS-003`, unaffected. Every certificate this returns is exact end to end: the numeric search only ever proposes a Gram matrix, which is rounded to nearby rationals and re-expanded to check it equals the target exactly before anything is returned — a `Some`/returned certificate is always sound regardless of what the float search converged to. Budget exhaustion is still `E-SOS-002`, undecided, never "not SOS" — say so, don't paraphrase it as a disproof. **The textbook PSD-not-SOS examples whose multiplier certificates are *singular* Gram matrices sitting exactly on the boundary of the PSD cone** — Motzkin's polynomial and Robinson's form — used to be out of reach for the original annealed alternating-projection search (a diagnosed convergence limitation at tangential PSD-cone intersections, not a soundness bug); the search now also tries Douglas–Rachford splitting with over-relaxation and a facial-reduction step, and with them both examples are found and exactly re-verified. **What's still open:** the homogeneous 3-variable form of Motzkin at multiplier power `N = 2` (`(x²+y²+z²)²·Motzkin_hom`, not `N = 1` — this specific homogeneous ternary form is not classically expected to be SOS at `N = 1` at all, unlike the affine 2-variable case) is still not reached, and this has now been attempted to closure rather than left as a budget skip: a new symmetry-reduction fallback (`real::sos::psd::symmetry_reduced_search`) exploits the target's own signed-permutation symmetry (order 16 here) to shrink its 165 free parameters to 26, and an exact algebraic zero-vector restriction (no numerics — Motzkin's known zero at `(1,1,1)` forces a specific null vector on any witnessing Gram matrix) shrinks that again to 16 — but deep Douglas–Rachford on that 16-parameter family still leaves the minimum eigenvalue at roughly `−1.4·10⁻⁸` after 6,000,000 iterations, a genuinely slow (not budget-limited) convergence. So a boundary-only certificate is still not guaranteed to be found in general — `E-SOS-002` still means "not found within this search", never "not SOS". Raise `basis_degree`, or fall back to `alkahest.decide`, exactly as for any other `E-SOS-002`. -31. **Double sums need `experimental.telescope2d`, not `zeilberger`** (since 3.9). `zeilberger`/`q_zeilberger` reach a sum over *one* index; `telescope2d(term, n, j, k)` is the Apagodu–Zeilberger generalization to a proper hypergeometric `F(n,j,k)` with **two** bound indices `j`, `k`, returning `a_0(n), …, a_J(n)` and *two* certificates `cert1`, `cert2` with `Σ_i a_i(n)·F(n+i,j,k) = Δ_j(cert1·F) + Δ_k(cert2·F)`, re-checked exactly in `Q(n,j,k)`. Three real, stated scope limits, not unfinished polish: (1) the certificate ansatz uses a *fixed* denominator built from `F`'s own shift-ratio denominators rather than a minimal 2-D Gosper normal form, so a search that finds nothing raises `E-HOLO-041` and does not prove no certificate exists; (2) `cert.boundary_status(j_lo, j_hi, k_lo, k_hi)` only accepts **constant** (not `n`-dependent) rectangles — for a natural range like `j = 0..n`, pick a fixed bound safely larger than any `n` you check and let `F`'s own combinatorial vanishing do the rest, exactly as the module's own worked example does; (3) the boundary of a rectangle is **four one-dimensional strip sums along its edges, not four corner-point evaluations** — a naive corner-evaluation formula is simply wrong — and this version only proves the sufficient (not necessary) condition that each strip vanishes identically, so `boundary_status` can return `"unknown"` for a boundary that is genuinely `0` but not by that pointwise route; it never guesses `"vanishes"`. There is no inhomogeneous `"nonzero"` verdict yet — an unresolved strip is always `"unknown"`. `E-HOLO-040` is the class refusal (not proper hypergeometric in `n, j, k`), `E-HOLO-042` a malformed call (`n`, `j`, `k` not distinct). +31. **Multi-sums need `experimental.telescope2d` (two bound indices) or `experimental.telescope_md` (any number `m >= 1`), not `zeilberger`** (since 3.9; `telescope_md` since 3.10). `zeilberger`/`q_zeilberger` reach a sum over *one* index; `telescope2d(term, n, j, k)` is the Apagodu–Zeilberger generalization to a proper hypergeometric `F(n,j,k)` with **two** bound indices `j`, `k`, returning `a_0(n), …, a_J(n)` and *two* certificates `cert1`, `cert2` with `Σ_i a_i(n)·F(n+i,j,k) = Δ_j(cert1·F) + Δ_k(cert2·F)`, re-checked exactly in `Q(n,j,k)`. `telescope_md(term, n, [x_1, ..., x_m])` is the same engine generalized to arbitrary `m` — `m = 1` degenerates to a single-sum-shaped search, `m = 2` behaves identically to `telescope2d` (which is now a thin wrapper over it), `m >= 3` is genuinely new — returning `cert.certs()` (a list of `m` certificates, a method not a property since it's a collection) instead of `cert1`/`cert2`. Four real, stated scope limits, not unfinished polish: (1) the certificate ansatz uses a *fixed* denominator built from `F`'s own shift-ratio denominators rather than a minimal Gosper normal form, so a search that finds nothing raises `E-HOLO-041` and does not prove no certificate exists; (2) `cert.boundary_status(j_lo, j_hi, k_lo, k_hi)` / `cert.boundary_status([(lo_1, hi_1), ..., (lo_m, hi_m)])` only accept **constant** (not `n`-dependent) boxes — for a natural range like `j = 0..n`, pick a fixed bound safely larger than any `n` you check and let `F`'s own combinatorial vanishing do the rest, exactly as the module's own worked examples do; (3) the boundary of a box is **`2m` `(m-1)`-dimensional face sums, not `2^m` corner-point evaluations** — a naive corner-evaluation formula is simply wrong — and this version only proves the sufficient (not necessary) condition that each face vanishes identically, so `boundary_status` can return `"unknown"` for a boundary that is genuinely `0` but not by that pointwise route; it never guesses `"vanishes"`. There is no inhomogeneous `"nonzero"` verdict yet — an unresolved face is always `"unknown"`; (4) `telescope_md`'s underlying exact linear solve is `O(rows · cols²)` and both grow fast with `m` and the certificate degree bound (measured: `m = 3` at certificate degree 2 already means a ~10,000-row, 245-unknown system taking ~47s to solve *per probe*), so two resource ceilings apply — a single probe above 400 unknowns is refused outright, and total work across every probe at or above 150 unknowns in one search call is capped to 300 — meaning `E-HOLO-041` can also mean "refused by a resource ceiling, not searched and found nothing," which the error message states explicitly; raising `m` or `max_cert_degree` further will not help once a ceiling is the reason. `E-HOLO-040` is the class refusal (not proper hypergeometric in the bound indices), `E-HOLO-042` a malformed call (indices not pairwise distinct, or `indices` empty for `telescope_md`). diff --git a/docs/mdbook/src/telescoping.md b/docs/mdbook/src/telescoping.md index a87fe2af..e4ff4ea8 100644 --- a/docs/mdbook/src/telescoping.md +++ b/docs/mdbook/src/telescoping.md @@ -556,6 +556,88 @@ Two restrictions are real, not unfinished polish: is always `"unknown"`, never guessed as zero and never resolved to an explicit `b(n)`. +## Multi-sum telescoping for `m ≥ 1` bound indices (`alkahest.experimental.telescope_md`) + +`telescope2d` reaches exactly two bound indices. `telescope_md` is the same +engine generalized to an arbitrary number `m ≥ 1` — `m = 1` degenerates +cleanly to a single-sum-shaped search, `m = 2` behaves identically to +`telescope2d` (which is now a thin wrapper over the general engine, not a +separate implementation), and `m ≥ 3` is genuinely new capability. Given a +proper hypergeometric `F(n, x_1, …, x_m)`, it finds `a_0(n), …, a_J(n)` and +`m` rational certificates `c_1, …, c_m` such that + +```text +Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t(c_t·F) +``` + +checked as an exact identity in `Q(n,x_1,…,x_m)` before it is ever returned — +same discipline, same fixed (non-minimal) denominator ansatz, generalized +from two axes to `m + 1`. + +```python +import alkahest as ak +from alkahest.experimental import telescope_md + +pool = ak.ExprPool() +n, x, y, z = pool.symbol("n"), pool.symbol("x"), pool.symbol("y"), pool.symbol("z") + +def factorial(e): + return ak.gamma(e + pool.integer(1)) + +# F(n,x,y,z) = n! / (x! y! z! (n-x-y-z)!) — the 4-category multinomial +# coefficient, genuinely non-separable: all three bound indices interact +# through the shared (n-x-y-z) term. +rest = n - x - y - z +f = factorial(n) / (factorial(x) * factorial(y) * factorial(z) * factorial(rest)) +cert = telescope_md(f, n, [x, y, z]) +cert.order # 1 +cert.coeffs() # [-4, 1] -> S(n+1) = 4*S(n) +cert.certs() # [c_1, c_2, c_3] +``` + +`Σ_{x,y,z} F = 4ⁿ` by the multinomial theorem (the number of length-`n` +strings over a 4-letter alphabet, grouped by letter counts) — exactly what +the order-1 relation says, checked in the test suite by direct exact +summation against the actual sum, not just against the telescoping identity. + +### The boundary is `2m` face sums, not `2^m` corner evaluations + +`cert.boundary_status([(lo_1, hi_1), …, (lo_m, hi_m)])` generalizes +`telescope2d`'s "four strip sums, not four corners" result: telescoping an +`m`-dimensional box gives `2m` sums, each over an `(m − 1)`-dimensional +**face** where one bound index is fixed to a boundary value — not `2^m` +point evaluations at the box's corners. The same sufficient (not necessary) +pointwise-vanishing criterion applies, generalized mechanically: fix one axis +to a constant and check that a gamma factor's argument no longer depends on +`n` or on any *other* bound index. + +### A real scaling cliff, and the resource ceilings that bound it + +Raising `m` or the certificate degree bound grows the ansatz search space far +faster than the numbers suggest: a certificate numerator spans a box of +`(max_cert_degree + 1)^(m+1)` unknowns, and there are `m` certificates. Worse, +the underlying exact linear solve (`rational_nullspace`) is a plain dense +`O(rows · cols²)` Gaussian elimination over unbounded-precision rationals, +and both dimensions grow with `m` and the degree bound well past what the box +size alone implies. Measured directly: at `m = 3`, certificate degree 2 means +a ≈10,000-row, 245-unknown system whose elimination step alone took ≈47 +seconds *per probe*; certificate degree 3 (770 unknowns) was still running +after several minutes. This is genuine arithmetic cost on a real, +correctly-posed linear system — not a bug, an infinite loop, or unbounded +coefficient blowup — but a caller still needs protection from it, since the +search tries every `(order, a_degree, cert_degree)` combination within the +stated bounds and would otherwise pay that same cost repeatedly. Two ceilings +apply: a single probe above 400 total unknowns is refused outright, and the +total work spent on probes at or above 150 unknowns is capped to 300 across +one whole search call — capping the number of genuinely expensive +elimination attempts to about one, regardless of how large `max_order` / +`max_a_degree` / `max_cert_degree` are. Neither ceiling touches the `m = 2` +search, whose default probes never exceed ≈140 unknowns. A search that hits +a ceiling still reports `SearchExhausted`, exactly like one that genuinely +found nothing — except the message says explicitly when a ceiling, not +genuine non-existence, is the reason, so raising the bounds further is not +silently misrepresented as a path to success. + ## Method The implementation is the standard Gosper-style reduction (Petkovšek–Wilf– @@ -586,9 +668,14 @@ explicit minimal-order certification, `guess_holonomic` — recurrence guessing from finite data — the `q`-analogue `q_zeilberger` over `Q(q)(qⁿ)(q^k)` with its own two-valued boundary verdict, `specialize_at_root_of_unity` — the step from a `Q(q)` identity to `q = ζ_d`, decided exactly in the cyclotomic -field `Q(ζ_d)` with its own three-valued verdict — and `telescope2d`, the +field `Q(ζ_d)` with its own three-valued verdict — `telescope2d`, the Apagodu–Zeilberger generalization to **two** bound indices, with its own -2-D boundary analysis (four strip sums, not four corner evaluations). +2-D boundary analysis (four strip sums, not four corner evaluations) — and +`telescope_md`, the further generalization to an arbitrary number `m ≥ 1` of +bound indices, with the `m`-dimensional boundary analysis (`2m` face sums, +not `2^m` corner evaluations) and the two resource ceilings that keep a +search with no certificate in reach a fast, honest refusal rather than an +unbounded computation as `m` or the certificate degree bound grow. Not shipped on the `q` side: multivariate (`q`-)telescoping and an inhomogeneous boundary arm. A `q`-sum whose support cannot be bounded is @@ -599,12 +686,16 @@ the wider congruence statements (e.g. uniform-in-`n` supercongruences, or `p`-adic statements not phrased as `Φ_d`-adic valuations) that literature contains. -Not shipped on the double-sum side: more than two bound indices, arbitrary -rational (not proper hypergeometric) summands — i.e. no general Wegschaider -reduction — a minimal 2-D Gosper certificate denominator (the ansatz uses a -fixed, larger-than-necessary one built from `F`'s own shift ratios), an -`n`-dependent rectangle in the boundary analysis, and an inhomogeneous -(`"nonzero"`) boundary verdict. +Not shipped on the multi-sum side: arbitrary rational (not proper +hypergeometric) summands, or a sum of several proper hypergeometric terms — +i.e. no general Wegschaider reduction — a minimal multivariate Gosper +certificate denominator (the ansatz uses a fixed, larger-than-necessary one +built from `F`'s own shift ratios), an `n`-dependent box in the boundary +analysis, and an inhomogeneous (`"nonzero"`) boundary verdict. The bound-index +count itself is no longer capped at two (`telescope_md` reaches any `m ≥ 1`), +but the search's two resource ceilings mean a caller cannot simply raise `m` +or the certificate degree bound without limit and expect the search to keep +running longer — past a point it refuses fast instead. Not yet shipped, and tracked as follow-up work: Ore-operator closure properties for D-finite functions (sums and products of holonomic objects) and the diff --git a/python/alkahest/experimental/__init__.py b/python/alkahest/experimental/__init__.py index ca29d1ca..3add0440 100644 --- a/python/alkahest/experimental/__init__.py +++ b/python/alkahest/experimental/__init__.py @@ -71,10 +71,32 @@ for the derivation). Only **constant** (not ``n``-dependent) rectangles are supported, and only a sufficient "each strip vanishes pointwise" criterion is checked — real, stated restrictions, not unfinished polish +- :func:`telescope_md` / :class:`TelescopingMdCertificate` (since 3.10) — the + general form of :func:`telescope2d` for an arbitrary number ``m >= 1`` of + bound indices, not just two: ``telescope_md(term, n, [x_1, ..., x_m])`` + finds ``m`` rational certificates ``c_1, ..., c_m`` with + ``Σ_i a_i(n)·F(n+i,x) = Σ_t Δ_t(c_t·F)``, and + :meth:`~alkahest.experimental.TelescopingMdCertificate.boundary_status` + decides the box-sum boundary (``2m`` face sums, the ``m``-dimensional + generalization of the four-strip-sum result, **not** ``2**m`` corner + evaluations). Same proper-hypergeometric-only, fixed-denominator, + constant-box-only scope as ``telescope2d`` — no genuinely broader summand + class. Raising ``m`` or ``max_cert_degree`` grows the ansatz search space + fast (a certificate numerator is a box of + ``(max_cert_degree + 1)**(m + 1)`` unknowns, and there are ``m`` of them), + so this module also enforces two resource ceilings on the underlying exact + linear solve — a single search probe is refused outright above 400 + unknowns, and the total work spent on probes at or above 150 unknowns in + one search call is capped to 300 — so a search with no certificate in + reach at all comes back as a fast, honest ``SearchExhausted`` (naming the + ceiling when that, not genuine non-existence, is why) rather than running + unboundedly long. See the Rust module docs + (``alkahest_cas::holonomic::telescoping2d``) for the complete, honestly- + stated scope and the exact ceiling values - This is a genuinely scoped-down engine: proper hypergeometric summands - only, no more than two bound indices, no general Wegschaider reduction, a - bounded-degree ansatz search rather than a minimal 2-D Gosper normal form. - See the Rust module docs for the complete, honest limitations list + only, no general Wegschaider reduction, a bounded-degree ansatz search + rather than a minimal Gosper normal form. See the Rust module docs for the + complete, honest limitations list Novelty filtering (:mod:`alkahest.experimental.novelty`): - :class:`RecurrenceClaim` — a recurrence in a normal form two presentations @@ -140,6 +162,7 @@ QRootOfUnitySpecialization, QZeilbergerCertificate, Telescoping2dCertificate, + TelescopingMdCertificate, asymptotic_expand, # P1 item 10 — asymptotic expansion at scale coefficient_asymptotics, @@ -159,6 +182,7 @@ q_zeilberger, series_solve, telescope2d, + telescope_md, z_transform, ) @@ -223,6 +247,7 @@ "RecurrenceClaim", # M4 — double-sum (Apagodu-Zeilberger) creative telescoping "Telescoping2dCertificate", + "TelescopingMdCertificate", "arg", "asymptotic_expand", # M5 — recurrence -> asymptotics @@ -266,6 +291,7 @@ "solve", # M4 — double-sum (Apagodu-Zeilberger) creative telescoping "telescope2d", + "telescope_md", "to_jax", "to_lean", "to_stablehlo",