Annihilator

The content on this page was written by AI under human supervision.

Annihilator takes the opening terms of a sequence, supplied exactly as integers or fractions, and finds the lowest-order linear recurrence with polynomial coefficients that those terms satisfy. From the recurrence it writes down the linear differential operator that annihilates the sequence's generating function. For the series expansion of a period integral, that operator is its Picard–Fuchs operatorthe linear differential operator of lowest order, with polynomial coefficients, whose solutions include the period; its order and singular points encode the geometry behind the integral. When no recurrence within the requested range fits the data it says so, and that negative answer is often the useful one.

What it does

Many sequences met in practice are holonomicsatisfying a linear recurrence whose coefficients are polynomials in the index; equivalently, the generating function satisfies a linear differential equation with polynomial coefficients: Taylor coefficients of a Feynman integral on its maximal cut, moments of a random walk on a lattice, counting sequences such as the Catalan numbers. Such a sequence obeys

$$\sum_{j=0}^{r} c_j(n)\, a_{n+j} = 0, \qquad \deg c_j \le s,$$

with the order $r$, the degree bound $s$ and the polynomials $c_j$ unknown in advance. Annihilator finds them from the numbers alone: it needs no integrand and no formula for $a_n$, only enough exact terms. Fitting a recurrence to series data by linear algebra is the standard method of the computer-algebra packages gfun (Salvy and Zimmermann) and Guess (Kauers). This package adds speed on long integer sequences and a search order that returns the true minimal differential order.

For a trial pair $(r,s)$ there are $(r+1)(s+1)$ unknown coefficients and one linear equation per index $n$, and the search insists on at least five spare equations so that a nonzero solution signals a real recurrence. The systems are solved modulo a prime $p \lt 2^{31}$ by Gaussian elimination on 64-bit integer numpy arrays, ten to a hundred times faster than the same modular nullspace computed with sympy, according to the package notes. Every candidate is checked against all supplied terms. The degree $s$, which becomes the order of the differential operator, is the outer loop of the scan, so the first hit has minimal order; the cheaper scan by system size can return an inflated order carrying spurious apparent-singularity factors.

With $(r,s)$ fixed, one exact rational nullspace at that single pair recovers the coefficients as fractions, which are then verified on trailing terms kept out of the fit (15 by default). The result is (r, s, coeffs, n_verified), where coeffs maps (j, k) to the fraction multiplying $n^k$ in $c_j(n)$; rec_to_theta converts it to the operator $\sum_j x^{r-j} c_j(\theta - j)$ in the Euler derivative $\theta = x\,d/dx$.

The output is a recurrence verified on every supplied term, not a proof that the sequence is holonomic. A pair found at one prime can occasionally be an artifact of that prime, so rerun with a second prime (the p argument) before relying on it. In mode='modp' the call returns (r, s, None, 0) as soon as the pair is found. A return value of None means that no recurrence with $r \le$ rmax and $s \le$ smax is consistent with the data modulo $p$, which rules out every operator in that range. The modular kernel returns a single null vector, so it does not measure the dimension of the solution space (nullspace_basis_modp in the factor/ module returns a full basis and the exact nullity when that is needed). The prime must be below $2^{31}$, the sequence at most a few thousand terms long, and the terms exact integers or fractions, not floating-point numbers. When the integrand is available, the pf_rank probe in Dipstick computes the operator's order from the integrand polynomial for comparison.

Three optional modules work with operators after one is found, for factorizing large operators modulo a prime. ore/ multiplies operators, applies them to series, tests whether one operator divides another from the right, and computes the least common left multiplethe lowest-order operator that each of the given operators divides from the right; its solutions are spanned by the solutions of all of them of several operators, searching orders upward so that the first verified result is minimal. factor/ extracts an order-4 right factor of an operator in $\theta$ form from its logarithmic series solutions at the origin, verifies it by right division at several primes, and reconstructs its rational coefficients by the Chinese remainder theorem. eigenring/ separates right factors that right division alone cannot split, using the maps of the operator's solution space to itself; a candidate factor is accepted only after it annihilates an independently generated series. These modules take an operator as a list of coefficient lists, lowest derivative order first, and most of them need python-flint.

Examples

Run the self-tests. From the root of BootLoops' bootloops-dev repository (GitHub organization BootLoops-ai):

python3 tools/annihilator/nullspace_fast/test_nullspace_fast.py
python3 tools/annihilator/annihilator.py --selftest

The first script checks the modular nullspace on a small rank-deficient matrix, then recovers two known recurrences from 100 terms each: Fibonacci (order 2, constant coefficients) and Catalan, $(n+2)C_{n+1} = (4n+2)C_n$. It then runs the second command three ways: unmodified, which must exit with status 0 and print ALL PASS; with an unknown flag, which must be rejected; and on a temporary copy with one series term altered, which must exit with status 1 and report FAIL. It prints one [PASS] line per check, fifteen in all, then ALL PASS, in a few seconds. The second command builds the exact return-probability moments of the simple-cubic lattice walk in dimensions 2 through 5, whose operators have differential orders 2 through 5. It prints

[series_pf self-test] d-dim simple-cubic LGF, mode='auto' (ode_order-first)
  d=2: (r,s)=(1,2)  held-out=15  [0.03s]  PASS
  d=3: (r,s)=(2,3)  held-out=15  [0.07s]  PASS
  d=4: (r,s)=(2,4)  held-out=15  [0.13s]  PASS
  d=5: (r,s)=(3,5)  held-out=15  [0.26s]  PASS
ALL PASS

Each line gives the order $r$ and coefficient degree $s$ found, the number of trailing terms (excluded from the fit) on which the exact coefficients were verified, and the run time. A line ends in PASS when $(r,s)$ is the known minimal pair for that dimension, whose $s$ equals the dimension, and the exact coefficients were recovered. The script exits with status 0 after ALL PASS and status 1 after SELFTEST FAILED; the self-test is also what runs when no flag is given.

A recurrence modulo a prime. The Catalan check from the test file is a template for any integer sequence already reduced modulo $p$:

from nullspace_fast import nullspace_mod_fast, find_recurrence_fast
p = 2147483629  # < 2^31, prime
cat = [1]
for n in range(99):
    cat.append(cat[-1] * (4 * n + 2) * pow(n + 2, p - 2, p) % p)
rec = find_recurrence_fast(cat, p, rmax=3, dmax=3)
r, d, w = rec

Run it from inside nullspace_fast/, as the test does, or import the same two functions from annihilator instead. The call returns r = 1, d = 1 and a flat list w of $(r+1)(d+1) = 4$ residues, ordered so that $c_0(n) = w_0 + w_1 n$ and $c_1(n) = w_2 + w_3 n$. Dividing through by $w_3$ modulo $p$ gives $(-2, -4, 2, 1)$, that is $-(4n+2)C_n + (n+2)C_{n+1} = 0$. find_recurrence_fast returns None when nothing up to rmax, dmax fits.

An exact operator from a moment sequence. The calls the built-in self-test makes, here for the three-dimensional lattice, followed by a printout of the recurrence:

from annihilator import moments_sc, pf_from_series, pretty_rec, rec_to_theta
d = 3
a = moments_sc(d, 100)
res = pf_from_series(a, rmax=14, smax=8, nverify=15, mode='auto')
r, s, coeffs, nver = res
print(pretty_rec(r, s, coeffs))

moments_sc returns 101 exact fractions beginning $1, \tfrac16, \tfrac{5}{72}, \dots$; pf_from_series finds $(r,s) = (2,3)$ modulo $2^{31}-1$, reconstructs the rational coefficients and verifies them on 15 unused terms (nver = 15). pretty_rec prints

(4*n^3+12*n^2+11*n+3)*a[n+0] + (-40*n^3+-180*n^2+-272*n+-138)*a[n+1] + (36*n^3+216*n^2+432*n+288)*a[n+2] = 0

which is $(n+1)(2n+1)(2n+3)\,a_n - (40n^3+180n^2+272n+138)\,a_{n+1} + 36(n+2)^3\,a_{n+2} = 0$. rec_to_theta(r, s, coeffs) then gives the third-order operator in $\theta$. With mode='modp' the same call returns (2, 3, None, 0) almost instantly; repeating it with p=2147483629 and getting the same pair is the second-prime confirmation.

Routines

Core (numpy; sympy only for rec_to_theta)

Operator modules (python-flint unless noted; import explicitly, for example from annihilator.ore import lclm_modp)

Used on this site

Requirements and source

Python 3 with numpy; sympy is needed by rec_to_theta, ore/ore_ops.py and ore/known_factors.py. The ore/lclm_modp.py, ore/ore_rdiv_modp.py, factor/ and eigenring/ modules need python-flint; without it the core still imports and annihilator.ore_rightdiv_modp is None. Nothing is compiled. From the repository root the self-tests are

python3 tools/annihilator/annihilator.py --selftest
python3 tools/annihilator/nullspace_fast/test_nullspace_fast.py
python3 tools/annihilator/eigenring/s4_eigenring_pilot.py --selftest
python3 tools/annihilator/ore/ore_ops.py
python3 tools/annihilator/ore/ore_rdiv_modp.py
python3 tools/annihilator/ore/known_factors.py
python3 tools/annihilator/ore/lclm_modp.py

The first two are the commands of the first example; each of the others finishes in seconds except lclm_modp.py, which takes about half a minute. The code is in tools/annihilator/ in the bootloops-dev repository, released under the MIT license.

← back to the tools index