Risch algorithm in Python
Decision procedure that finds an elementary antiderivative or proves none exists for elementary functions.
- File
risch_algorithm.py, 178 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/risch-algorithm/. The authors are listed in prov/GY4RMM.jsonl.
"""The Risch algorithm for one exponential monomial over Q(x).
Input: b in Q[x] of degree at least 1, so that t = exp(b) is a monomial over
Q(x) with t' = b' t, and an integrand
f = sum over k of f_k t^k, f_k in Q(x), k in Z
(a Laurent polynomial in t). Output: None if f has no elementary integral,
otherwise (ys, R, (A, D)) with
integral of f = sum over k != 0 of y_k t^k + R + integral of A / D,
where y_k, R are in Q(x) and D is squarefree, so that the last integral is a
sum of logarithms (found by the Rothstein-Trager method, not done here).
Polynomials are lists of Fractions, low degree first; a rational function is
a pair (num, den). Cost model: operations in Q. The work is polynomial in the
degrees of b, of the numerators and of the denominators (proved); the size of
the rational numbers is not bounded here.
"""
from __future__ import annotations
from fractions import Fraction as Q
def norm(a):
a = [Q(c) for c in a]
while a and a[-1] == 0:
a.pop()
return a
def deg(a):
return len(a) - 1 # -1 for the zero polynomial
def add(a, b):
n = max(len(a), len(b))
return norm([(a[i] if i < len(a) else 0) + (b[i] if i < len(b) else 0) for i in range(n)])
def scale(a, c):
return norm([c * x for x in a])
def sub(a, b):
return add(a, scale(b, -1))
def mul(a, b):
if not a or not b:
return []
out = [Q(0)] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
if x:
for j, y in enumerate(b):
out[i + j] += x * y
return norm(out)
def divmod_p(a, b):
a, q = list(a), [Q(0)] * max(len(a) - len(b) + 1, 0)
while len(a) >= len(b) and a:
c, s = a[-1] / b[-1], len(a) - len(b)
q[s] = c
for i, y in enumerate(b):
a[i + s] -= c * y
a = norm(a)
return norm(q), a
def exquo(a, b):
q, r = divmod_p(a, b)
assert not r
return q
def gcd(a, b):
a, b = norm(a), norm(b)
while b:
a, b = b, divmod_p(a, b)[1]
return scale(a, 1 / a[-1]) if a else []
def deriv(a):
return norm([i * a[i] for i in range(1, len(a))])
def ext_euclid(a, b, c):
"""(s, t) with s a + t b = c and deg s < deg b; gcd(a, b) must divide c."""
r0, r1, s0, s1 = norm(a), norm(b), [Q(1)], []
while r1: # invariant: s_i a = r_i mod b
q, r = divmod_p(r0, r1)
r0, r1, s0, s1 = r1, r, s1, sub(s0, mul(q, s1))
g = r0
q, r = divmod_p(c, g)
assert not r
s = divmod_p(mul(s0, q), b)[1] if deg(b) > 0 else []
t = exquo(sub(c, mul(s, a)), b)
return s, t
def rat_add(x, y):
n = add(mul(x[0], y[1]), mul(y[0], x[1]))
d = mul(x[1], y[1])
g = gcd(n, d) if n else d
return (exquo(n, g), exquo(d, g)) if n else ([], [Q(1)])
def rat_deriv(x):
n, d = x
return rat_add((mul(deriv(n), d), mul(d, d)), (scale(mul(n, deriv(d)), -1), mul(d, d)))
def hermite_reduce(A, D):
"""(R, (A2, D2)) with A / D = R' + A2 / D2 and D2 squarefree (Mack's linear version)."""
P, A = divmod_p(A, D) # polynomial part, integrated directly
g = (norm([0] + [c / (i + 1) for i, c in enumerate(P)]), [Q(1)])
Dm = gcd(D, deriv(D))
Ds = exquo(D, Dm)
while deg(Dm) > 0:
Dm2 = gcd(Dm, deriv(Dm))
Dms = exquo(Dm, Dm2)
B, C = ext_euclid(scale(exquo(mul(Ds, deriv(Dm)), Dm), -1), Dms, A)
A = sub(C, exquo(mul(deriv(B), Ds), Dms))
g = rat_add(g, (B, Dm))
Dm = Dm2
return g, (A, Ds)
def rde_poly_coeff(fk, gk):
"""Rational y with y' + fk y = gk, fk in Q[x] of degree m >= 0; None if none.
The solution is unique because y' + fk y = 0 has no nonzero rational solution.
"""
a, d = gk
if not a:
return ([], [Q(1)])
g0 = gcd(a, d)
a, d = exquo(a, g0), exquo(d, g0)
h = gcd(d, deriv(d)) # @line L4
c, r = divmod_p(mul(a, mul(h, h)), d)
if r: # a simple pole in gk
return None
# h q' + (fk h - h') q = c, an operator that raises the degree by m + deg h
lead = sub(mul(fk, h), deriv(h))
shift, lc = deg(lead), lead[-1]
q = []
while c: # @line L5
j = deg(c) - shift
if j < 0:
return None
coef = c[-1] / lc
mono = [Q(0)] * j + [coef]
q = add(q, mono)
c = sub(c, add(mul(h, deriv(mono)), mul(lead, mono)))
gq = gcd(q, h) if q else h
return (exquo(q, gq), exquo(h, gq)) if q else ([], [Q(1)])
def risch_exp(b, f):
"""Integrate f = {k: (num, den)} in Q(x)[t, 1/t] with t = exp(b)."""
b = norm(b)
assert deg(b) >= 1
db = deriv(b)
ys = {}
for k in sorted(f): # @line L1
if k == 0 or not norm(f[k][0]):
continue
y = rde_poly_coeff(scale(db, k), (norm(f[k][0]), norm(f[k][1]))) # @line L2
if y is None:
return None # @line L3
ys[k] = y
A, D = f.get(0, ([], [Q(1)])) # @line L6
A, D = norm(A), norm(D)
if not A:
return ys, ([], [Q(1)]), ([], [Q(1)])
R, rest = hermite_reduce(A, D) # @line L7
return ys, R, rest # @line L8
test_risch_algorithm.py 166 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/risch-algorithm/. The authors are listed in prov/GY4RMM.jsonl.
"""Tests for the Risch algorithm with one exponential monomial t = exp(b).
Every answer is checked by differentiation in Q(x)[t, 1/t]. Integrable cases
are built from known antiderivatives, and the answer must match them (the
rational parts y_k are unique). A FAIL answer is compared with a direct
linear algebra search for a rational solution of y' + k b' y = f_k.
"""
import random
from fractions import Fraction as Q
from risch_algorithm import (add, deg, deriv, exquo, gcd, hermite_reduce,
mul, norm, rat_add, rat_deriv, risch_exp, scale, sub)
ONE = [Q(1)]
def rat_eq(x, y):
return not norm(sub(mul(x[0], y[1]), mul(y[0], x[1])))
def rat_mul(x, y):
return (mul(x[0], y[0]), mul(x[1], y[1]))
def check(b, f, out):
ys, R, rest = out
db = deriv(b)
for k in set(f) | set(ys):
want = f.get(k, ([], ONE))
if k == 0:
got = rat_add(rat_deriv(R), rest)
else:
y = ys.get(k, ([], ONE))
got = rat_add(rat_deriv(y), rat_mul((scale(db, k), ONE), y))
assert rat_eq(got, want), k
D = rest[1]
assert deg(gcd(D, deriv(D))) <= 0 # squarefree
def rand_poly(rng, d, lo=-3, hi=3):
return norm([rng.randint(lo, hi) for _ in range(d + 1)])
def rand_den(rng, maxpow=3):
d = ONE
for _ in range(rng.randint(0, 2)):
p = norm([rng.randint(-3, 3), rng.choice([1, 2])] + ([1] if rng.random() < 0.3 else []))
for _ in range(rng.randint(1, maxpow)):
d = mul(d, p)
return d
def brute_rde(fk, g):
"""Search y = q / d with deg q <= deg a + deg d + 2 by Gaussian elimination."""
a, d = g
N = deg(a) + deg(d) + 3
cols = []
for j in range(N): # image of q = x^j in q' d - q d' + fk q d
q = [Q(0)] * j + [Q(1)]
cols.append(add(sub(mul(deriv(q), d), mul(q, deriv(d))), mul(mul(fk, q), d)))
rhs = mul(a, d)
rows = max([len(c) for c in cols] + [len(rhs)])
M = [[(c[i] if i < len(c) else Q(0)) for c in cols] + [rhs[i] if i < len(rhs) else Q(0)]
for i in range(rows)]
r = 0
for c in range(N):
piv = next((i for i in range(r, rows) if M[i][c] != 0), None)
if piv is None:
continue
M[r], M[piv] = M[piv], M[r]
for i in range(rows):
if i != r and M[i][c] != 0:
t = M[i][c] / M[r][c]
M[i] = [x - t * y for x, y in zip(M[i], M[r])]
r += 1
return all(any(x != 0 for x in row[:N]) or row[N] == 0 for row in M)
def test_known():
# e^(x^2) and e^x / x and x^2 e^(-x^2) have no elementary integral
assert risch_exp([0, 0, 1], {1: ([1], ONE)}) is None
assert risch_exp([0, 1], {1: ([1], [0, 1])}) is None
assert risch_exp([0, 0, -1], {1: ([0, 0, 1], ONE)}) is None
# x e^x = ((x - 1) e^x)'
ys, R, rest = risch_exp([0, 1], {1: ([0, 1], ONE)})
assert rat_eq(ys[1], ([-1, 1], ONE))
# (2x^2 + 1) e^(x^2) = (x e^(x^2))'
ys, _, _ = risch_exp([0, 0, 1], {1: ([1, 0, 2], ONE)})
assert rat_eq(ys[1], ([0, 1], ONE))
# (x - 1) / x^2 e^x = (e^x / x)'
ys, _, _ = risch_exp([0, 1], {1: ([-1, 1], [0, 0, 1])})
assert rat_eq(ys[1], (ONE, [0, 1]))
def test_edges():
assert risch_exp([0, 1], {}) == ({}, ([], ONE), ([], ONE))
out = risch_exp([0, 1], {0: ([5], ONE)}) # a constant
check([Q(0), Q(1)], {0: ([5], ONE)}, out)
out = risch_exp([0, 1], {2: ([], ONE), 0: ([], [1, 1])}) # zero coefficients
assert out[0] == {}
f = {-3: ([7], ONE)} # negative power: 7 e^(-3x) = (-7/3 e^(-3x))'
ys, _, _ = risch_exp([0, 1], f)
assert rat_eq(ys[-3], ([Q(-7, 3)], ONE))
f = {0: ([1], mul(mul([1, 0, 1], [1, 0, 1]), mul([1, 0, 1], [1, 0, 1])))} # 1/(x^2+1)^4
out = risch_exp([0, 1], f)
check([Q(0), Q(1)], f, out)
assert rat_eq((out[2][1], ONE), ([1, 0, 1], ONE))
def test_random_integrable():
rng = random.Random(1)
for _ in range(60):
b = norm(rand_poly(rng, rng.randint(1, 3)))
if deg(b) < 1:
b = [Q(0), Q(1)]
db = deriv(b)
f, want = {}, {}
for k in rng.sample([-2, -1, 1, 2, 3], rng.randint(1, 3)):
y = (rand_poly(rng, rng.randint(0, 3)), rand_den(rng))
g = gcd(*y) if y[0] else y[1]
if not y[0]:
continue
y = (exquo(y[0], g), exquo(y[1], g))
want[k] = y
f[k] = rat_add(rat_deriv(y), rat_mul((scale(db, k), ONE), y))
if rng.random() < 0.7:
f[0] = (rand_poly(rng, rng.randint(0, 5)), rand_den(rng, 4))
out = risch_exp(b, f)
assert out is not None
check(b, f, out)
for k, y in want.items():
assert rat_eq(out[0][k], y)
def test_random_against_brute_force():
rng = random.Random(2)
fails = 0
for _ in range(80):
b = norm(rand_poly(rng, rng.randint(1, 2)))
if deg(b) < 1:
b = [Q(0), Q(1)]
k = rng.choice([-1, 1, 2])
g = (rand_poly(rng, rng.randint(0, 3)), rand_den(rng, 2))
if not g[0]:
continue
out = risch_exp(b, {k: g})
assert (out is not None) == brute_rde(scale(deriv(b), k), g)
if out is None:
fails += 1
else:
check(b, {k: g}, out)
assert fails > 20 # random integrands are usually not integrable
def test_hermite():
rng = random.Random(3)
for _ in range(40):
A, D = rand_poly(rng, rng.randint(0, 8)), rand_den(rng, 5)
if not A:
continue
R, rest = hermite_reduce(A, D)
assert rat_eq(rat_add(rat_deriv(R), rest), (A, D))
assert deg(gcd(rest[1], deriv(rest[1]))) <= 0
What the next level needs
- the tests did not run in this build