Gosper's algorithm in Python
Decides whether a hypergeometric term has a hypergeometric antidifference and finds it when it does.
- File
gosper_algorithm.py, 208 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/gosper-algorithm/. The authors are listed in prov/CAVQHD.jsonl.
"""Gosper's algorithm for indefinite hypergeometric summation over Q.
A hypergeometric term t(n) is given by its ratio r(n) = t(n+1)/t(n) as a pair
(num, den) of polynomials. Polynomials are lists of Fractions, lowest degree
first, with no trailing zero; [] is the zero polynomial.
gosper(num, den) returns a rational function y = (p, q) such that
z(n) = y(n) t(n) satisfies z(n+1) - z(n) = t(n), or None if t has no
hypergeometric antidifference.
Cost model: operations in Q. With d the larger degree of num and den, H the
bound on integer root differences used for the Gosper form, and D the degree
bound for the polynomial x, the cost is O(H d^2) for the form (gcds of degree
at most d) plus O((d + D)^3) for the linear system. H and D are integers read
off the input, so the bound is pseudo-polynomial in the bit size.
"""
from fractions import Fraction
from typing import List, Optional, Tuple
Poly = List[Fraction]
RatFun = Tuple[Poly, Poly]
def norm(a: Poly) -> Poly:
while a and a[-1] == 0:
a.pop()
return a
def P(*cs: object) -> Poly:
"""Polynomial from coefficients, lowest degree first."""
return norm([Fraction(c) for c in cs]) # type: ignore[arg-type]
def deg(a: Poly) -> int:
return len(a) - 1
def add(a: Poly, b: Poly) -> Poly:
n = max(len(a), len(b))
return norm([(a[i] if i < len(a) else Fraction(0)) + (b[i] if i < len(b) else Fraction(0)) for i in range(n)])
def scale(a: Poly, c: Fraction) -> Poly:
return norm([c * x for x in a])
def sub(a: Poly, b: Poly) -> Poly:
return add(a, scale(b, Fraction(-1)))
def mul(a: Poly, b: Poly) -> Poly:
if not a or not b:
return []
c = [Fraction(0)] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
if x:
for j, y in enumerate(b):
c[i + j] += x * y
return norm(c)
def divmod_poly(a: Poly, b: Poly) -> Tuple[Poly, Poly]:
if not b:
raise ZeroDivisionError("division by the zero polynomial")
r = list(a)
q = [Fraction(0)] * max(len(a) - len(b) + 1, 0)
inv = 1 / b[-1]
for k in range(len(a) - len(b), -1, -1):
t = r[k + len(b) - 1] * inv
q[k] = t
if t:
for j, y in enumerate(b):
r[k + j] -= t * y
return norm(q), norm(r[: len(b) - 1])
def exquo(a: Poly, b: Poly) -> Poly:
q, r = divmod_poly(a, b)
if r:
raise ArithmeticError("division is not exact")
return q
def monic(a: Poly) -> Poly:
return [x / a[-1] for x in a] if a else []
def gcd(a: Poly, b: Poly) -> Poly:
while b:
a, b = b, divmod_poly(a, b)[1]
return monic(a)
def shift(a: Poly, s: int) -> Poly:
"""a(n + s), by Horner's rule with the linear polynomial n + s."""
out: Poly = []
lin = P(s, 1)
for c in reversed(a):
out = add(mul(out, lin), [c] if c else [])
return out
def root_bound(a: Poly) -> int:
"""An integer bound on the absolute values of the complex roots of a (Fujiwara):
|z| <= 2 max_i |a_(d-i) / a_d|^(1/i), with the term for i = d halved."""
d = deg(a)
best = 0
for i in range(1, d + 1):
c = abs(a[d - i] / a[d])
if i == d:
c = c / 2
if c == 0:
continue
r = max(int(float(c) ** (1.0 / i)) - 1, 0)
while Fraction(r) ** i < c: # smallest integer r with r^i >= c
r += 1
best = max(best, r)
return 2 * best
def solve_linear(cols: List[Poly], rhs: Poly) -> Optional[List[Fraction]]:
"""A solution v of sum_j v_j cols[j] = rhs, or None. Gaussian elimination over Q."""
rows = max([len(c) for c in cols] + [len(rhs)])
m = len(cols)
A = [[(cols[j][i] if i < len(cols[j]) else Fraction(0)) for j in range(m)]
+ [rhs[i] if i < len(rhs) else Fraction(0)] for i in range(rows)]
piv: List[int] = []
r = 0
for c in range(m):
p = next((i for i in range(r, rows) if A[i][c] != 0), None)
if p is None:
continue
A[r], A[p] = A[p], A[r]
inv = 1 / A[r][c]
A[r] = [v * inv for v in A[r]]
for i in range(rows):
if i != r and A[i][c] != 0:
f = A[i][c]
A[i] = [u - f * v for u, v in zip(A[i], A[r])]
piv.append(c)
r += 1
if any(A[i][m] != 0 for i in range(r, rows)):
return None
v = [Fraction(0)] * m
for i, c in enumerate(piv):
v[c] = A[i][m]
return v
def degree_bound(a: Poly, B: Poly, c: Poly) -> int:
"""Largest possible degree of a polynomial x with a(n) x(n+1) - B(n) x(n) = c(n)."""
d1, d2 = deg(a), deg(B)
if d1 != d2:
return deg(c) - max(d1, d2)
if a[-1] != B[-1]:
return deg(c) - d1
e = deg(c) - d1 + 1
al = a[-2] if d1 >= 1 else Fraction(0)
be = B[-2] if d2 >= 1 else Fraction(0)
d0 = (be - al) / a[-1]
if d0.denominator == 1 and d0 >= 0:
e = max(e, int(d0))
return e
def gosper(num: Poly, den: Poly) -> Optional[RatFun]:
"""Antidifference of the term t with t(n+1)/t(n) = num/den, as y = z/t.
Precondition: num and den nonzero polynomials over Q. Returns (p, q) with
(p/q)(n+1) * num(n)/den(n) - (p/q)(n) = 1, or None if no rational y exists.
Cost: O(H d^2 + (d + D)^3) operations in Q, see the module docstring.
"""
if not num or not den:
raise ValueError("the ratio must be a nonzero rational function")
a = scale(monic(num), num[-1] / den[-1]) # @line L1
b = monic(den)
c = P(1)
H = root_bound(num) + root_bound(den) # every integer root difference is at most H
for h in range(H + 1): # @line L2
while True:
g = gcd(a, shift(b, h)) # @line L3
if deg(g) <= 0:
break
a = exquo(a, g) # @line L4
b = exquo(b, shift(g, -h))
for i in range(1, h + 1):
c = mul(c, shift(g, -i))
B = shift(b, -1) # @line L5
D = degree_bound(a, B, c)
if D < 0: # @line L6
return None
cols = []
for j in range(D + 1): # @line L7
e = P(*([0] * j + [1]))
cols.append(sub(mul(a, shift(e, 1)), mul(B, e)))
v = solve_linear(cols, c)
if v is None: # @line L8
return None
x = P(*v)
p, q = mul(B, x), c # @line L9
g = gcd(p, q)
p, q = exquo(p, g), exquo(q, g)
return scale(p, 1 / q[-1]), monic(q)
test_gosper_algorithm.py 218 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/gosper-algorithm/. The authors are listed in prov/CAVQHD.jsonl.
import random
from fractions import Fraction
from typing import List, Optional, Tuple
import pytest
from gosper_algorithm import P, Poly, RatFun, gosper
# Independent polynomial helpers for the oracles (dense lists, lowest degree first).
def trim(a: List[Fraction]) -> List[Fraction]:
a = list(a)
while a and a[-1] == 0:
a.pop()
return a
def pmul(a: List[Fraction], b: List[Fraction]) -> List[Fraction]:
if not a or not b:
return []
c = [Fraction(0)] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
for j, y in enumerate(b):
c[i + j] += x * y
return trim(c)
def padd(a: List[Fraction], b: List[Fraction], s: int = 1) -> List[Fraction]:
n = max(len(a), len(b))
return trim([(a[i] if i < len(a) else 0) + s * (b[i] if i < len(b) else 0) for i in range(n)]) # type: ignore[misc]
def pev(a: List[Fraction], x: Fraction) -> Fraction:
s = Fraction(0)
for c in reversed(a):
s = s * x + c
return s
def pshift(a: List[Fraction], s: int) -> List[Fraction]:
# binomial expansion of sum a_i (n + s)^i
out = [Fraction(0)] * len(a)
for i, c in enumerate(a):
binom = 1
for j in range(i + 1):
out[j] += c * binom * Fraction(s) ** (i - j)
binom = binom * (i - j) // (j + 1)
return trim(out)
def from_roots(z: Fraction, roots: List[Fraction]) -> List[Fraction]:
p = [Fraction(z)]
for u in roots:
p = pmul(p, [Fraction(-u), Fraction(1)])
return p
def property_holds(num: Poly, den: Poly, y: RatFun, rng: random.Random) -> bool:
"""r(n) y(n+1) - y(n) = 1, as polynomials and at random points."""
p, q = y
p1, q1 = pshift(p, 1), pshift(q, 1)
# num p1 q - den q1 p = den q1 q
lhs = padd(pmul(pmul(num, p1), q), pmul(pmul(den, q1), p), -1)
rhs = pmul(pmul(den, q1), q)
if lhs != rhs:
return False
for _ in range(3):
x = Fraction(rng.randint(-40, 40), rng.randint(1, 5))
if pev(den, x) != 0 and pev(q, x) != 0 and pev(q1, x) != 0:
if pev(num, x) / pev(den, x) * pev(p1, x) / pev(q1, x) - pev(p, x) / pev(q, x) != 1:
return False
return True
def solve_any(M: List[List[Fraction]], v: List[Fraction]) -> Optional[List[Fraction]]:
"""Some solution of M u = v by row reduction, or None."""
rows, cols = len(M), len(M[0]) if M else 0
A = [M[i][:] + [v[i]] for i in range(rows)]
r, piv = 0, []
for c in range(cols):
p = next((i for i in range(r, rows) if A[i][c] != 0), None)
if p is None:
continue
A[r], A[p] = A[p], A[r]
for i in range(rows):
if i != r and A[i][c] != 0:
f = A[i][c] / A[r][c]
A[i] = [x - f * y for x, y in zip(A[i], A[r])]
piv.append((r, c))
r += 1
if any(A[i][cols] != 0 for i in range(r, rows)):
return None
u = [Fraction(0)] * cols
for i, c in piv:
u[c] = A[i][cols] / A[i][c]
return u
def ansatz_summable(num_roots: List[Fraction], den_roots: List[Fraction], num: Poly, den: Poly) -> bool:
"""Oracle: search y = p/q with q = prod_{i=1..h} num(n - i), h the largest integer
difference of a denominator root and a numerator root, and a generous degree for p."""
h = max([int(b - a) for a in num_roots for b in den_roots if (b - a).denominator == 1 and b >= a] + [0])
q = [Fraction(1)]
for i in range(1, h + 1):
q = pmul(q, pshift(num, -i))
q1 = pshift(q, 1)
d = max(len(num), len(den)) - 1
m = (len(q) - 1) + d + h + 3
# num(n) p(n+1) q(n) - den(n) p(n) q(n+1) = den(n) q(n) q(n+1)
cols = []
for j in range(m + 1):
e = [Fraction(0)] * j + [Fraction(1)]
cols.append(padd(pmul(pmul(num, pshift(e, 1)), q), pmul(pmul(den, e), q1), -1))
rhs = pmul(pmul(den, q), q1)
n_rows = max([len(c) for c in cols] + [len(rhs)])
M = [[(c[i] if i < len(c) else Fraction(0)) for c in cols] for i in range(n_rows)]
return solve_any(M, [rhs[i] if i < len(rhs) else Fraction(0) for i in range(n_rows)]) is not None
def partial_sums_agree(num: Poly, den: Poly, y: RatFun, N: int) -> bool:
"""Brute force: sum_{k<n} t(k) = z(n) - z(0) with t(0) = 1, t(k+1) = r(k) t(k)."""
p, q = y
t, s = Fraction(1), Fraction(0)
ts = []
for k in range(N + 1):
ts.append(t)
if pev(den, Fraction(k)) == 0:
return True # the term is not defined further; nothing more to compare
t = t * pev(num, Fraction(k)) / pev(den, Fraction(k))
if any(pev(q, Fraction(k)) == 0 for k in range(N + 1)):
return True
z = [pev(p, Fraction(k)) / pev(q, Fraction(k)) * ts[k] for k in range(N + 1)]
for n in range(N + 1):
if s != z[n] - z[0]:
return False
s += ts[n]
return True
def test_worked_examples() -> None:
rng = random.Random(1)
# t = binom(2k, k) / 4^k: r = (2k + 1) / (2k + 2), y = 2k
y = gosper(P(1, 2), P(2, 2))
assert y == (P(0, 2), P(1))
# t = k * k!: r = (k + 1)^2 / k, y = 1/k, so z = k!
y = gosper(P(1, 2, 1), P(0, 1))
assert y == (P(1), P(0, 1)) and property_holds(P(1, 2, 1), P(0, 1), y, rng)
# t = 1/(k(k+1)): r = k/(k+2), y = -(k+1), so z = -1/k
assert gosper(P(0, 1), P(2, 1)) == (P(-1, -1), P(1))
# t = 1: y = n
assert gosper(P(1), P(1)) == (P(0, 1), P(1))
# k! and the harmonic term 1/(k+1) have no hypergeometric antidifference
assert gosper(P(1, 1), P(1)) is None
assert gosper(P(1, 1), P(2, 1)) is None
# t = 2^k: y = 1
assert gosper(P(2), P(1)) == (P(1), P(1))
def test_edge_cases() -> None:
rng = random.Random(2)
with pytest.raises(ValueError):
gosper([], P(1))
with pytest.raises(ValueError):
gosper(P(1), [])
# alternating sign: t = (-1)^k, y = -1/2
assert gosper(P(-1), P(1)) == (P(Fraction(-1, 2)), P(1))
# polynomial term t = k + 1: r = (k + 2)/(k + 1), z = k(k + 1)/2
y = gosper(P(2, 1), P(1, 1))
assert y is not None and property_holds(P(2, 1), P(1, 1), y, rng)
# many reduction steps: t = Delta of (k+6)!/k!-type terms, shift 6
num, den = P(1, 1), P(7, 1) # r = (k+1)/(k+7): t = 6!/((k+1)...(k+6)), summable
y = gosper(num, den)
assert y is not None and property_holds(num, den, y, rng)
assert partial_sums_agree(num, den, y, 30)
def rand_roots(rng: random.Random, k: int) -> List[Fraction]:
return [Fraction(rng.choice([-3, -2, -1, 0, 1, 2, 3, 4, 5]), rng.choice([1, 1, 1, 2])) for _ in range(k)]
def test_random_against_ansatz() -> None:
rng = random.Random(4242)
found = 0
for _ in range(270):
nr, dr = rand_roots(rng, rng.randint(0, 2)), rand_roots(rng, rng.randint(0, 2))
z = Fraction(rng.choice([1, 1, 1, 1, -1, 2]))
num, den = from_roots(z, nr), from_roots(Fraction(1), dr)
y = gosper(num, den)
assert (y is not None) == ansatz_summable(nr, dr, num, den)
if y is not None:
found += 1
assert property_holds(num, den, y, rng)
assert partial_sums_agree(num, den, y, 12)
assert found >= 20
def test_random_constructed_summable() -> None:
"""t = z(n+1) - z(n) for a random hypergeometric z, so t is Gosper-summable."""
rng = random.Random(77)
for _ in range(250):
while True:
u, v = rand_roots(rng, rng.randint(0, 2)), rand_roots(rng, rng.randint(0, 2))
zc = Fraction(rng.choice([1, 2, 3, -1, Fraction(1, 2)]))
rn, rd = from_roots(zc, u), from_roots(Fraction(1), v)
if rn != rd:
break
# rho = rn/rd, t(n) = z(n)(rho(n) - 1), r(n) = rho(n) (rho(n+1) - 1) / (rho(n) - 1)
num = pmul(rn, padd(pshift(rn, 1), pshift(rd, 1), -1))
den = pmul(padd(rn, rd, -1), pshift(rd, 1))
# an unreduced input: multiply both by the same factor
num, den = pmul(num, rd), pmul(den, rd)
y = gosper(num, den)
assert y is not None
assert property_holds(num, den, y, rng)
What the next level needs
- the tests did not run in this build