Calculus API

Differentiate expr with respect to var.

exprExpr or DerivedResult

Expression to differentiate. If a DerivedResult is passed (e.g. the output of a previous diff() call) its .value is used automatically.

varExpr

The differentiation variable.

DerivedResult

Contains the derivative in .value and a step log in .steps.

>>> p = ExprPool(); x = p.symbol("x")
>>> d = diff(x**3, x)   # DerivedResult; d.value == 3*x^2
>>> diff(d, x).value    # second derivative — passes DerivedResult directly

Compute the derivative using forward-mode automatic differentiation (dual numbers). Produces the same result as diff through a different path; useful for cross-validation.

Symbolic partial derivatives of expr with respect to each variable in vars.

This is not the same as grad() from _transform, which differentiates a TracedFn produced by trace() and returns a GradTracedFn for numeric evaluation (often composed with jit()).

exprExpr

Expression to differentiate.

varslist[Expr]

Variables; one partial derivative per entry.

list[Expr]

[∂expr/∂vars[0], ∂expr/∂vars[1], …] (plain expressions, no step log).

DiffError

E-DIFF-001 when a function that depends on one of vars has no derivative rule (e.g. floor). A partial is never silently 0.

See also

• diff — single-variable derivative with DerivedResult steps.
• grad — JAX-style gradient of a traced Python function.
• jacobian — matrix of partials for a vector-valued function.

>>> p = ExprPool(); x, y = p.symbol("x"), p.symbol("y")
>>> symbolic_grad(x**2 + y**2, [x, y])  # [2*x, 2*y]

alkahest.euler_lagrange(L, coords, t) — the Euler–Lagrange equations.

coords are applied undefined functions [q1(t), q2(t), …] (built with pool.func("q1", [t])); L is written with them and their time derivatives (diff(q1, t)). Returns one expression per coordinate, to be read as = 0:

d/dt(∂L/∂q̇ᵢ) − ∂L/∂qᵢ

(higher derivatives in L contribute the Euler–Poisson terms). Pass the result through functions_to_symbols() to solve for the accelerations.

Raises DiffError: E-DIFF-005 when a coordinate is not an undefined function applied to arguments, E-DIFF-007 when it is not a function of t alone, E-DIFF-001 when L contains a function with no derivative rule.

alkahest.functions_to_symbols(exprs, funcs, t, domain=None) — replace q(t), q'(t), q''(t), … by plain symbols.

exprs is an Expr or a list of them; funcs the applied undefined functions [q1(t), …]. Every qᵢ(t) and every derivative of it that occurs becomes a fresh symbol named q, q_dot, q_ddot, q_d3, … (a trailing _ is added while the name is taken by a symbol already in exprs), in domain (default real).

Returns (new_exprs, mapping) where new_exprs has the shape of exprs and mapping is a dict {node: symbol}; subs(e, {s: n for n, s in mapping.items()}) undoes the rewrite.

Raises DiffError (E-DIFF-005 / E-DIFF-007) for an entry of funcs that is not an undefined function of t alone.

alkahest.replace_function(expr, pattern, body) — substitute a definition for an undefined function everywhere, derivatives included.

pattern is f(x1, …, xn) with distinct variables and body an expression in them. Every f(u…) becomes body at u… and every derivative of f becomes the matching derivative of body at u… (SymPy's subs(f, Lambda(x, body)).doit()). subs() with the key f(x…) does the same, but only at that exact argument.

Raises DiffError E-DIFF-005 when pattern is not an undefined function applied to distinct variables.

Indefinite or definite integral of expr with respect to var.

With two arguments this returns the indefinite integral (antiderivative). When both bounds a and b are supplied it returns the definite integral ∫_a^b expr d(var) via the fundamental theorem of calculus, F(b) − F(a) where F is the antiderivative.

exprExpr or DerivedResult

Integrand. DerivedResult is auto-coerced to .value.

varExpr

Integration variable.

a, bExpr, optional

Lower and upper bounds. Provide both for a definite integral.

DerivedResult

ValueError

If exactly one of a, b is provided.

∫_{-∞}^{∞} of a rational function takes the residue theorem rather than the fundamental theorem: 2πi·Σ_{Im α > 0} Res(P/Q, α), with both convergence conditions checked exactly (deg Q ≥ deg P + 2, and Q with no real root) and the value cross-checked against a rigorous enclosure of the same integral before it is returned — so integrate(1/(x**4+1), x, -oo, oo).value is π/√2. A divergent integral is reported as divergent and never given a finite value, including when a Cauchy principal value exists (x/(x**2+1)). The route covers denominators whose Hurwitz spectral factor is rational, and in radicals every denominator of degree ≤ 4 after even-normalisation; outside that it declines explicitly.

Everything else is the elementary FTC wrapper: it needs an elementary antiderivative whose value at each bound can be established — by substitution for a finite bound, by limit for an infinite one. An endpoint value that still mentions the integration variable (what an unevaluated RootSum limit or substitution produces) is refused rather than pushed into F(b) − F(a), where the two ends would cancel to a finite-looking but meaningless number.

A singularity of the integrand inside the interval is detected and refused (E-INT-001) rather than pushed through F(b) − F(a), which would produce a finite-looking but wrong number. Rational integrands are checked exactly (real root isolation of the reduced denominator); other integrands are checked by confirming numeric blow-up, which catches e.g. ∫_0^2 sec²x dx (pole at π/2; the FTC difference tan 2 is negative for a positive integrand). Finally, when both bounds are numeric the closed form must denote a finite real: results such as −log(−1) or 2·(0^{1/2})^{−1} are refused, not returned.

Where the antiderivative passes through a pole of tan/cot/sec/ csc (the tan(x/2) of a half-angle substitution), the endpoint value is the one-sided limit, and an interior pole splits the interval, so ∫_0^{2π} dx/(2 + cos x) is 2π/√3; when the limit cannot be established the integral is refused. An antiderivative carrying constants lifted from floating point is refused rather than presented as exact. (A + B·sin²x)^{∓1/2} integrates exactly to EllipticF/EllipticE (parameter m = k²), so ∫_0^{π/2} dx/√(1 − sin²x/4) is EllipticK(1/4).

Integrands with symbolic parameters or in unexpanded form — x*exp(-k*x), 1/(k**2+m**2), (s*z)**2*exp(-z**2/2), exp(-a*t)*exp(-b*(u-t)), x*x**(n-1) — are normalised (powers of sums expanded, exponentials merged, powers combined) and integrated termwise when the direct routes decline. A division by a parameter, and the 4ac - b**2 > 0 an atan of a quadratic with an unknown-sign discriminant needs, are listed in result.verification["side_conditions"]. integrate(f, x, 0, oo) of Σ c·x**p·exp(-α·x) is Σ c·Γ(p+1)/α**(p+1) when α > 0 is established (a Positive symbol, a numeric rate or assumptions); with the sign unknown it is refused (E-INT-001). So integrate(exp(-y/2)/sqrt(y), y, 0, oo) is √(2π).

For an indefinite integral, result.verification reports "exactly_verified" only when the kernel proves the symbolic residual diff(result, var) - expr is zero. A successful result with "unverified" may instead have passed the integrator's separate numeric soundness gate. Definite integral results have no antiderivative verification metadata.

>>> p = ExprPool(); x = p.symbol("x")
>>> integrate(x**2, x).value         # x^3/3
>>> integrate(x**2, x, 0, 1).value   # 1/3

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

No docstring.

Build a piecewise expression from a list of (condition, value) pairs and a default value.

Conditions must be Predicate expressions built with the pool's predicate constructors (pool.lt, pool.le, pool.gt, pool.ge, etc.). SymPy's (value, condition) order is also accepted: a pair whose second element is a predicate and whose first is not is read that way. A pair in which neither element is a predicate raises TypeError — it used to be accepted with the value in the condition slot, which diff and eval_expr then mis-read.

Compute the limit of expr as var → point.

exprExpr or DerivedResult

varExpr

pointExpr

dirstr or None

"+" for right limit, "-" for left limit.

DerivedResult

The limit value is .value. .certificate is Lean Filter.Tendsto source when the emitter recognises the pattern (currently a small x → +∞ fragment: exp(-x) → 0, xⁿ exp(-x) → 0, exp(x) → +∞, and a crude exp-ratio). Unrecognised shapes — including finite and one-sided limits — withhold rather than emitting sorry.

A limit that does not exist raises rather than returning a value. An infinite limit is returned as the canonical ∞ symbol (or −1·∞); an oscillating one — lim_{x→∞} sin x, lim_{x→0} sin(1/x) — raises LimitError (E-LIMIT-005), and a two-sided limit at an odd-order pole raises E-LIMIT-003 asking for a direction.

A limit whose value turns on the sign of a free parameter — lim_{t→∞} exp(-k·t) is 0 for k > 0, +∞ for k < 0 and 1 for k = 0 — raises E-LIMIT-006 naming the parameter, unless the sign is stated: either inside alkahest.context(assumptions=...) or by declaring the symbol with Domain.Positive.

Undefined functions. f(t) and Derivative(f(t), t) are opaque: a limit that contains one is refused (E-LIMIT-007) unless every such node is listed in continuous=[f(t), ...], the caller's statement that it is continuous at point (a finite point). The value is then the substitution f(t) -> f(point), and only for expressions that are jointly continuous there (sums, products, non-negative integer powers, exp, sin, cos, sinh, cosh, atan, tanh and negative powers of a factor that is a non-zero constant at the point). A 0/0 form such as (f(t) - f(a))/(t - a) is a derivative quotient; continuity does not give it, so it is refused rather than guessed.

Nested calls that expect an Expr still work: _coerce_expr() unwraps a DerivedResult automatically. Use .value when comparing against an Expr directly.

Truncated Taylor/Laurent series expansion.

exprExpr or DerivedResult

varExpr

Expansion variable.

pointExpr

Expansion point.

orderint

Truncation order: every term with exponent < order is present and the remainder is O(h**order), for a Laurent series (a pole at point) as for a Taylor one.

Series

The expansion is bounded: it honours an active Budget and, with none, an internal work ceiling. When every coefficient is rational (sin, exp, log(1 + x), rational functions and their compositions, where no irrational constant such as sin(1) appears) they come from truncated power-series arithmetic over the rationals, which is fast at any order. Otherwise coefficients are formed by repeated differentiation without re-simplifying, so an expression whose derivatives do not close — a nested radical such as sqrt(t**-2 + t**-1) — grows by a constant factor per coefficient, and a high order is unreachable rather than merely slow.

Running out of room raises SeriesError with code E-SERIES-003 (or BudgetExceededError when a budget stopped it). It never returns a shorter series: the O(h**order) term is a claim about the remainder, and attaching it to fewer coefficients than were asked for would be a false one that no caller could audit.

A removable singularity at point is expanded, not refused: the expression is put over a common denominator and the two power series are divided, so series(sin(x)/x, x, 0, 4) is 1 - x**2/6 + O(x**4) and series(1/sin(x), x, 0, 3) is x**-1 + x/6 + O(x**3). Substituting the expansion point into repeated derivatives would give 0/0 for the first of those, which is why the direct route cannot see it.

A singularity that is not removable and not a pole — a branch point (sqrt(x)), a logarithmic singularity (log(x)), an essential one (exp(1/x)) — raises SeriesError with code E-SERIES-004. There is no Laurent expansion to return, and a Series whose coefficients are 0**-1 would report success while evaluating to NaN.

The branch point is the one of those three that has an expansion this function cannot hold: sqrt(x) = x**(1/2) is a Puiseux series, and experimental.puiseux_series() returns it (together with the ramification index, and only after verifying it). Series carries no ramification and every consumer of the underlying expansion reads an integer valuation, so the fractional case lives in a sibling type rather than widening this one — the refusal here is deliberate, not a gap. The refusal says so: for series(sqrt(x), x, 0, 3), sqrt(x)*exp(x) or x**(1/3)*sin(x) its message reads "not a power series at the expansion point: it is a branch point …", names puiseux_series and ends with the verified Puiseux expansion itself (x^(1/2) + O(x^(3))):

>>> from alkahest.experimental import puiseux_series
>>> px = puiseux_series(sqrt(sin(x)), x, pool.integer(0), 5)  # doctest: +SKIP
>>> px.ramification, px.valuation                            # doctest: +SKIP
(2, Fraction(1, 2))

At an infinite point (pool.pos_infinity() or pool.neg_infinity()) the expansion is in powers of 1/var, by substituting var = ±1/t and expanding about t = 0+. order keeps its meaning in t: every term var**-e with e < order is present and the remainder is O(var**-order), as in SymPy's series(f, x, oo, n):

>>> series(1/x + 1/(x**2 + 1), x, pool.pos_infinity(), 4)  # doctest: +SKIP
x**-1 + x**-2 + O(x**-4)

A function with no Laurent expansion in 1/var there (exp(x), log(x), sqrt(x), sin(x)), or a point that mentions infinity without being ±oo, raises SeriesError with code E-SERIES-007 — it is never expanded by substituting infinity into the Taylor formula.

Use Series.truncated() to get the expansion as an Expr with the O(.) term dropped, which is what eval_expr() and simplify() accept.

No docstring.

expr

truncated()

No docstring.

Definite symbolic sum of expr for k from lo to hi (inclusive).

Uses Gosper / hypergeometric summation, which covers the standard first-course cases directly: Faulhaber sums (Σ k, Σ k², Σ k³, …), geometric series (Σ r^k), and telescoping products.

hi may also be ExprPool.pos_infinity (or parse("oo")) for an infinite sum. Gosper's algorithm doesn't apply there (its antidifference has no reason to vanish at ∞); instead these families are summed exactly, each knowing its own convergence condition:

• polynomial × geometric — Σ_{k=a}^∞ P(k)·q^k for any q free of k (r, exp(-ke*tau), 1/2), e.g. Σ k·rᵏ = r/(1-r)². Convergence needs |q| < 1. When the assumptions (symbol domains, alkahest.context(assumptions=...)) establish it, nothing is attached; when they refute it the sum is refused; otherwise the result carries the hypothesis as a side condition:

s = ak.sum_definite(r**k, k, pool.integer(0), pool.pos_infinity())
s.value                               # -1 * (r - 1)^-1
s.steps[0]["side_conditions"]         # ['(1 + (-1 * abs(r))) > 0']

• rational series whose poles are linear in k — partial fractions, telescoping between poles that differ by an integer, and the Hurwitz-zeta values left over: Σ 1/(k²+k) = 1, Σ 1/(4k²−1) = 1/2, Σ_{k≥1} 1/k² = π²/6 (however 1/k² is spelled), Σ 1/(2k+1)² = π²/8, Σ 1/(k(k+½)) = 4 − 4·log 2; digamma / trigamma at non-special points, with the positivity of a symbolic shift stated as a side condition. Irreducible quadratic denominators k² + p·k + s (p rational) are summed through Im ψ(x + iβ): Σ_{k≥1} 1/(k²+a²) = (π·a·coth(π·a) − 1)/(2a²) (printed with tanh; side condition a² > 0 unless known), Σ_{k≥1} k/(k⁴+4) = 3/8. When the non-elementary Re ψ(x + iβ) would be needed (Σ 1/((k+1)(k²+1))) the sum is refused (E-SUM-002), not guessed.

• rational × geometric with integer poles (the polylogarithm family): Σ rᵏ/k = −log(1 − r), Σ rᵏ/k^s = polylog(s, r) (dilog for s = 2), Σ rᵏ/(k(k+1)) elementary; |r| < 1 as for the geometric series, and r = −1 when the rational factor has no polynomial part (Σ (−1)ᵏ/k = −log 2).

• hypergeometric terms — t_{k+1}/t_k rational in k (powers of k-polynomials with rational roots, gamma/factorial/ binomial/pochhammer of integer-slope linear arguments, q^k, exp(c·k)): the sum from the first non-zero term t₀ is t₀·hyper(a, b, z), rewritten by expand_hyper() when a closed form is known — Σ qᵏ/k! = exp(q), Σ x^{2k}/(2k)! = cosh x, Σ (−1)ᵏ/(2k+1) = π/4, Σ binomial(2k, k)·xᵏ = (1 − 4x)^{−1/2} — and returned as the exact hyper node otherwise (Σ (k!)³/(3k)!). Convergence: always when p ≤ q or the series terminates; |z| < 1 (stated or decided) when p = q + 1, with z = 1 needing Σb − Σa > 0 (Gauss) and z = −1 needing Σb − Σa > −1; p > q + 1 diverges (Σ k!·xᵏ, E-SUM-003).

Divergent series (Σ 1/k, Σ 2ᵏ, Grandi's Σ (−1)ᵏ, Σ rᵏ under r > 1) raise SumError (E-SUM-003). Odd zeta values are named, not approximated: Σ 1/k³ = zeta(3), Σ 1/(2k+1)³ = 7/8·zeta(3). An infinite lower bound, −∞, or ∞ − 1 is refused with E-SUM-003.

Reversed bounds. When hi < lo the result follows the reversal convention Σ_{k=a}^{b} f = −Σ_{k=b+1}^{a−1} f, the same convention SymPy uses: sum_definite(k, k, 5, 1) is −9, not 0. Systems that treat a reversed range as empty (Mathematica) return 0 here, so do not rely on the two agreeing. hi == lo − 1 is the empty sum and is 0 under both conventions.

A pole in the range. If a term is undefined for some integer between the bounds — Σ_{k=1}^{10} 1/(k²−1) at k = 1 — the telescoped difference is not the sum, and SumError (E-SUM-003) is raised instead of returning the division by zero embedded in an otherwise plausible expression.

A geometric ratio that is a symbol. Σ_{k=0}^{n} rᵏ is summed to (r^{n+1} − 1)/(r − 1), which Gosper cannot reach: its certificate lives in ℚ(k) and a symbolic ratio is not in ℚ. r = 1 is a genuine second branch — the sum is n + 1 there and the formula is 0/0 — so it is not assumed away: the derivation step carries the side condition r − 1 ≠ 0, readable as:

s = ak.sum_definite(r**k, k, pool.integer(0), n)
s.steps[0]["side_conditions"]        # ['(r + -1) ≠ 0']

and the returned expression declines to evaluate at r = 1 rather than producing a number there. The same holds for Σ k^p·rᵏ (finite bounds) and for ratios such as exp(-ke*tau) — the multiple-dosing superposition Σ_{k=0}^{n-1} e^{-ke·k·τ}.

Indefinite symbolic sum of expr with respect to index k.

Same algorithm family as sum_definite() (Gosper / hypergeometric), which covers polynomial terms (Faulhaber sums, e.g. Σ k, Σ k²) and geometric terms (Σ 2^k) as special cases — both are hypergeometric in k, so no separate handler is needed.

Definite symbolic product of expr for k from lo to hi (inclusive).

Supports q ∈ ℚ(k) whose numerator and denominator split into ℤ-linear factors; the result is emitted as a ratio of gamma values times integer powers. Anything outside that class raises ProductError (E-PROD-001…E-PROD-003) rather than being approximated.

Empty range. hi < lo gives 1, the empty product, whatever the term is — including a term that is identically zero, since no factor is ever taken. (Note this differs from sum_definite(), which follows the reversal convention for hi < lo rather than treating it as empty.)

Indefinite symbolic product of expr with respect to index k.

Unevaluated discrete product — SymPy-shaped Product(term, (k, lo, hi)).

Call doit() for a closed form when term(k) reduces to ℚ(k) whose numerator and denominator factor into ℤ-linear polynomials in k.

Bounds are inclusive on both ends, matching sum_definite().

import alkahest as ah p = ah.ExprPool() k, n = p.symbol("k"), p.symbol("n") prod = ah.Product(k, (k, p.integer(1), n)) out = ah.simplify(prod.doit().value)

doit()

Solve a linear recurrence; equation may be DerivedResult.

equation is read as equation == 0 and may be written with either spelling of the shift — f(n+1) - f(n) - n**2 and f(n) - f(n-1) - (n-1)**2 are the same equation and give the same answer. Both are solved as stated: the right-hand side is re-indexed together with the sequence terms, so substituting the result back into the equation you wrote reproduces it.

Without initials the general solution is returned with fresh symbols C0, C1, …. For an order-2 equation with a repeated characteristic root the basis is {r**n, n*r**n}, so the family really is two-parameter and two independent initial conditions can be met.

No docstring.

alkahest.zeilberger(term, n, k, *, limits=None, max_order=4, max_degree=16, minimal=False) -> ZeilbergerCertificate

Zeilberger's algorithm (creative telescoping) for a proper hypergeometric term F(n, k). Returns a verified certificate: the recurrence Σ_i a_i(n)·F(n+i,k) = ΔG with G = R·F is re-checked as an exact identity in Q(n)(k) before it is returned.

The verified statement is that identity in k. A recurrence for the sum S(n) = Σ_{k=k_lo}^{k_hi} F(n,k) is a separate claim, and it is decided here rather than left to the caller: boundary is "vanishes" (the homogeneous recurrence holds), "nonzero" (the inhomogeneous one does, with b(n) in boundary_rhs) or "unknown" (nothing follows about the sum).

limits is the summation range as (k_lo, k_hi), each an Expr or an int. It defaults to (0, n) — Σ_{k=0}^{n}, the convention the classical identities and the OEIS formula field use — and the range actually used is echoed back on limits, so the assumption is on the record rather than silent. The verdict is about that range and changes with it: truncating Σ_{k=0}^{n} to Σ_{k=0}^{n-1} generally turns "vanishes" into "nonzero". A range this analysis cannot place — endpoints that are not integer-affine in n — is "unknown", never "vanishes".

max_order and max_degree are upper bounds, not starting points. The search visits the (order, degree) grid by iterative deepening, cheapest candidate first, and returns the first relation that passes exact verification — so raising either bound widens what can be found without slowing down an input that was already decided at a low order and degree.

The order it returns is not claimed to be minimal. Cheapest-first means a cheap order-2 probe can be reached before an expensive order-1 one, so an order-2 result does not rule out an order-1 relation; order_is_minimal is False to say so rather than leaving it to be assumed. Pass minimal=True to search order-ascending instead — every degree 0..=max_degree at order J is refused before order J+1 is tried — and the flag is set.

What that establishes is minimality at certificate degree <= max_degree, not minimality outright: a lower-order relation whose certificate needs a higher degree than the bound is never probed, and order–degree trade-offs are real in creative telescoping. So order_is_minimal is a statement about the grid that was swept, and max_degree is part of it — raise the bound if the minimality is what the caller is after. Order-ascending is the hopeless low-order sweep the default plan exists to avoid, and it costs accordingly; ask for it when minimality is the result, not as a habit.

Raises HolonomicError rather than guessing when term is outside the proper hypergeometric class (E-HOLO-001) or when the bounded search is exhausted (E-HOLO-002).

A verified Zeilberger certificate, returned by zeilberger().

Carries the recurrence coefficients a_0(n), …, a_J(n) and the rational certificate R(n, k) satisfying, as an exact identity in Q(n)(k),

Σ_i a_i(n)·F(n+i, k) = G(n, k+1) − G(n, k) with G = R·F.

That identity is re-checked exactly before this object is constructed — a returned certificate is a proof, not a numerical match.

It is an identity in k, and only that. A recurrence for the sum S(n) = Σ_{k=k_lo}^{k_hi} F(n,k) is a second statement, and it is decided separately over the range in limits — see boundary, which is "vanishes", "nonzero" or "unknown":

• "vanishes" — proved; Σ_i a_i(n)·S(n+i) = 0 holds for the sum.
• "nonzero" — proved; the recurrence is the inhomogeneous Σ_i a_i(n)·S(n+i) = b(n) with b(n) returned as boundary_rhs. Still a theorem, just not the homogeneous one.
• "unknown" — neither was established. Nothing follows about the sum; boundary_reason says what stopped the proof.

The distinction is not cosmetic. For Σ_{k=0}^{n} C(n,k)/(k+1) the certificate is perfectly valid and the homogeneous recurrence is false — the true relation is (n+2)·S(n+1) − (2n+2)·S(n) = 1. Reading a recurrence off a certificate without this verdict is how a valid certificate becomes a false theorem.

A verdict is also not a statement about every n. The range in limits need not be a range at every n — k = 3..n−3 runs backwards at n = 3, 4 — and the telescoping needs G = R·F to be finite at every integer k in it. boundary_valid_from and certificate_poles carry those two bounds, rather than leaving a bare b(n) with an implied "for every n" attached to it.

order

order_is_minimal

coeffs

n

certificate

boundary_term

limits

boundary

boundary_valid_from

certificate_poles

boundary_rhs

boundary_reason

implies_sum_recurrence

side_conditions

derivation

boundary_at(k_lo, k_hi)

to_lean()

z = alkahest.zeilberger(binomial(n, k)**2, n, k)
print(alkahest.to_lean(z))   # same as z.to_lean()

Check the WZ identity F(n+1,k) - F(n,k) = G(n,k+1) - G(n,k).

Returns True only when the identity is proved and False only when it is refuted. The test is exact: both sides are normalised in Q(n)(k) (for proper hypergeometric F, G -- gamma, factorial, binomial, pochhammer -- the identity is divided by F and the ratios of factorials reduced to rational functions), and a counterexample search at small integers can also refute.

Raises SumError (E-SUM-001) when neither a proof nor a refutation is available, instead of answering False.

Fit a P-recursive recurrence to terms, or refuse.

Searches for the smallest order (then smallest degree) such that

Σ_{i=0}^{order} p_i(n) · u(n+i) = 0

holds on every equation the terms provide, with deg p_i <= degree and u(start) = terms[0]. Everything is exact rational arithmetic: terms must be Python int or fractions.Fraction (arbitrary size), and a float is refused rather than rounded.

A candidate is fitted only when it is over-determined: a recurrence with U = (order+1)(degree+1) unknowns is tried only where the terms give at least U + min_surplus equations, min_surplus defaulting to U. Below that a nullspace vector exists no matter what the data is, so a fit would be interpolation wearing a recurrence's clothes.

Returns None only for a genuine negative — every (order, degree) in bounds was fitted with adequate surplus and none of them worked. If the terms were too few to sweep the bounds, the call raises HolonomicError (E-HOLO-005) saying how many terms the remaining candidates need, because a loop that reads "not holonomic" off a grid it never swept has closed a branch it never explored.

Check GuessedRecurrence.status on what comes back. A returned fit satisfies every equation the terms supplied; that it is the sequence's recurrence is a separate question, and only status == "confirmed" (equivalently confirmed is True) says the data is entitled to suggest it. The two undecided outcomes are returned rather than refused, because in both the relation genuinely holds and the caller can act on it: "underdetermined" (several independent relations — read basis) and "singular" (the leading coefficient vanishes inside the data at singular_indices, which is what a corrupted term looks like — recompute those terms).

>>> import alkahest as ak
>>> motzkin = [1, 1, 2, 4, 9, 21, 51, 127, 323, 835, 2188,
...            5798, 15511, 41835, 113634, 310572, 853467,
...            2356779, 6536382, 18199284, 50852019]
>>> guess = ak.guess_holonomic(motzkin)
>>> guess.order, guess.degree
(2, 1)
>>> guess.coeffs
((-3, -3), (-5, -2), (4, 1))
>>> guess.confirmed, guess.surplus_terms
(True, 14)

which is (n+4)·M(n+2) = (2n+5)·M(n+1) + (3n+3)·M(n), Motzkin's recurrence. Seven terms suffice to determine it and confirm nothing, so they are refused rather than fitted:

>>> try:
...     ak.guess_holonomic(motzkin[:7])
... except ak.HolonomicError as exc:
...     print(exc.code, str(exc).startswith("7 terms are not enough"))
E-HOLO-005 True

A sequence with no such relation comes back None — but only once every candidate in bounds has actually been tested:

>>> ak.guess_holonomic([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
...                     43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
...                     101, 103, 107, 109, 113, 127, 131, 137, 139, 149,
...                     151, 157, 163, 167, 173, 179, 181, 191, 193, 197,
...                     199, 211, 223, 227, 229, 233, 239, 241, 251, 257,
...                     263, 269, 271, 277, 281]) is None
True

One wrong term does not stop a fit — it is absorbed into roots of the leading coefficient, which is what status and singular_indices exist to say. The relation returned holds on every equation these terms provide and is not Motzkin's recurrence:

>>> spoiled = [1, 1]
>>> for i in range(1, 35):
...     spoiled.append(((2 * i + 3) * spoiled[-1] + 3 * i * spoiled[-2]) // (i + 3))
>>> spoiled[12] += 1
>>> fit = ak.guess_holonomic(spoiled)
>>> fit.status, fit.confirmed
('singular', None)
>>> fit.singular_indices
(10, 11, 12)

terms

the first terms of the sequence, exact int / Fraction.

max_order

largest recurrence order to try.

max_degree

largest coefficient-polynomial degree to try.

start

index n of terms[0]; the coefficient polynomials are polynomials in that n.

min_surplus

surplus equations demanded of a fit. Defaults to the candidate's own unknown count, i.e. the data must be twice what the ansatz needs. Setting it to 0 disables the requirement without disabling the reporting.

check_evidence

when False, every candidate is fitted regardless of surplus and the first fit is returned with GuessedRecurrence.status set honestly. This is the escape hatch, not the default — the same role check_precision=False plays on guess_relation().

HolonomicError

E-HOLO-005 when the terms cannot support the search, or when the only fit found had no surplus left to confirm it.

TypeError

when a term is not an exact rational.

ValueError

when the bounds are not positive.

A P-recursive recurrence fitted to a finite sequence, with its evidence.

The relation is

Σ_{i=0}^{order} p_i(n) · u(n+i) = 0

where u(start) = terms[0] and p_i is the polynomial whose coefficients are coeffs[i], lowest degree first. The coefficients are integers, content-free, with the leading one positive, so the same data always produces the same object.

A fitted recurrence is a conjecture. It is exact on the terms it was shown — that is checked, not assumed — and says nothing about the terms it was not. surplus_terms is how much of the data agreed without being asked to; hand the result to zeilberger() when the sequence has a hypergeometric summand, and to holds_for() when more terms turn up.

Read status rather than bool(guess.confirmed) when the difference between "the data does not support this" and "the data cannot decide" matters — confirmed is True / False / None.

order: int

degree: int

start: int

coeffs: tuple[tuple[int, ...], ...]

n_terms: int

n_equations: int

equations_used: int

surplus_terms: int

dimension: int

basis: tuple[tuple[tuple[int, ...], ...], ...]

singular_indices: tuple[int, ...]

min_surplus: int

untested_candidates: int

status: str

means: str

confirmed: bool | None

evidence() → dict

holds_for(terms: Sequence[Any]) → bool

to_exprs(pool: ExprPool, var: Expr) → list[Expr]

GUESS_STATUSES = ('confirmed', 'singular', 'underdetermined', 'unconfirmed')

No docstring.

GUESS_STATUS_MEANINGS = {'confirmed': 'the terms over-determined the fit, singled it out, and the operator is non-singular at every index they constrain — still a conjecture about t…

No docstring.

Check S(index(p)) ≡ expect(p) (mod p**k) at every prime in primes.

This is the loop a supercongruence investigation runs, with the residue coming from ModularRecurrence.value_mod() rather than from big-integer arithmetic — so the cost per prime is O(index(p)) machine-word multiplications instead of O(index(p)) operations on integers with Θ(index(p)) digits.

>>> import alkahest as ak
>>> apery = ak.ModularRecurrence(
...     [[1, 3, 3, 1], [-117, -231, -153, -34], [8, 12, 6, 1]], [1, 5]
... )
>>> sweep = ak.supercongruence_sweep(
...     apery, [5, 7, 11, 13, 17, 19, 23, 29, 31], k=3, expect=1
... )
>>> sweep.holds, sweep.n_tested, sweep.sharp
(True, 9, True)

sharp being True is the sweep earning its keep: some prime has v_p(A(p−1) − 1) exactly 3, so the mod-p⁴ version of the same statement is false and the p³ in Beukers' theorem is best possible.

recurrence

the sequence, as a ModularRecurrence.

primes

the primes to test. Checked for primality here, and a composite raises HolonomicError (E-HOLO-006) rather than being skipped, because a sweep that silently drops its inputs reports a range it did not cover. The check is up front so that the other E-HOLO-006 — p**(k + extra_precision) past the machine-word ceiling, a fact about one prime rather than about the call — can be recorded in CongruenceSweep.skipped() and the sweep continue.

k

the claimed exponent. a(p-1) ≡ 1 (mod p**4) is k=4.

index

p -> n, the index to evaluate at. Defaults to p - 1, which is the shape of nearly every Apéry-like supercongruence.

expect

p -> value, or a constant. Defaults to 0.

extra_precision

how many digits beyond k to compute, so that v_p(LHS − RHS) can be measured rather than merely bounded below by k. This is what makes CongruenceSweep.sharp and the valuation histogram possible; 0 disables both. One is enough — a residue that is 0 mod p**k and non-zero mod p**(k+1) has v_p exactly k — and it is the default because every extra digit costs modulus headroom: the evaluation runs at p**(k + extra + loss) and refuses past 2**62, so a generous extra_precision buys a finer histogram at the price of the larger primes in the range.

max_counterexamples

stop after this many failures. A conjecture that fails at the first ten primes does not need the eleventh.

The outcome of supercongruence_sweep().

Everything the sweep learned, including the parts that argue against the conjecture. The residues are kept so a caller can re-derive any of it.

holds: bool

n_tested: int

largest: int

claimed_exponent: int

sharp: bool

n_skipped: int

residues() → dict[int, int]

counterexamples() → list[tuple[int, int, int]]

valuations() → dict[int, int]

skipped() → list[tuple[int, str]]

A P-recursive recurrence prepared for evaluation modulo prime powers.

Holds Σ_{i=0}^{J} a_i(n)·S(n+i) = b(n) with integer polynomial coefficients, plus the J initial values S(start) … S(start+J-1). Nothing about p is fixed at construction, so one object serves an entire supercongruence sweep.

coeffs[i] is a_i written lowest-degree coefficient first — the same convention as GuessedRecurrence.coeffs, so a fitted recurrence is handed straight over.

Why this is Rust and not Python: it is exact modular arithmetic in a hot loop over machine words, and the p-adic precision accounting that makes a singular index safe has to happen inside that loop. Per CONTRIBUTING.md § Rust vs Python, points 2 and 5 of the Rust column.

The recurrence is a hypothesis about your sequence. This class checks that it is well formed and that every forward step is determined p-adically; it cannot check that your sequence satisfies it. Certify with zeilberger(), or fit and confirm with guess_holonomic().

>>> import alkahest as ak
>>> # Apéry A005259: (n+2)³A(n+2) = (34n³+153n²+231n+117)A(n+1) − (n+1)³A(n)
>>> apery = ak.ModularRecurrence(
...     [[1, 3, 3, 1], [-117, -231, -153, -34], [8, 12, 6, 1]],
...     [1, 5],
... )
>>> apery.value_mod(12, 13, 3)              # A(p−1) ≡ 1 (mod p³), p = 13
1
>>> apery.value_mod(10006, 10007, 3)        # …and at p = 10007, in ~5 ms
1

order

degree

start

is_homogeneous

coeffs()

rhs()

initial()

value_mod(n, p, k)

values_mod(indices, p, k)

evaluate(indices, p, k)