Risch-Norman heuristic integration in Python
Heuristic symbolic integrator that is fast in practice but not complete.
- File
risch_norman_heuristic.py, 193 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-norman-heuristic/. The authors are listed in prov/XPCGZH.jsonl.
"""The Risch-Norman (parallel Risch) heuristic over Q.
The integrand lives in Q(x, t_1, ..., t_m), where each t_i is a monomial
(for example exp(a(x)), log(b(x)) or tan(x)) given by its derivative
t_i' = N_i / d with N_i in Q[x, t] and one common denominator d in Q[x].
The scaled derivation is D = d d/dx + sum_i N_i d/dt_i, so D = d * (d/dx).
Polynomials are dicts {exponent tuple (x, t_1, ..., t_m): Fraction}.
The integrand is P / Q with Q given as a list of factors (q, e), the q
pairwise coprime and squarefree. The method guesses
integral = U / B + sum_j c_j log(v_j),
with B from the factors of Q, v_j the factors of Q and any extra candidates,
U a polynomial of bounded total degree with unknown coefficients and c_j
unknown constants, and solves one linear system over Q. It returns
(U, B, [(c_j, v_j)]) or None. None does not prove that no elementary integral
exists: the guess may be too small.
Cost model: operations in Q. With N unknowns, S the largest number of terms
of a polynomial in the cleared equation and E <= the number of its
monomials, building the system takes O(N S^2) and Gaussian elimination
O(E N^2); both are proved bounds for this procedure. Whether it succeeds is
a matter of the guess, and it is observed to succeed on most integrals met
in practice.
"""
from __future__ import annotations
from fractions import Fraction
from itertools import product as iproduct
Poly = dict
def clean(p: Poly) -> Poly:
return {e: c for e, c in p.items() if c != 0}
def padd(*ps: Poly) -> Poly:
r: Poly = {}
for p in ps:
for e, c in p.items():
r[e] = r.get(e, Fraction(0)) + c
return clean(r)
def pscale(p: Poly, c) -> Poly:
return clean({e: v * c for e, v in p.items()})
def pmul(*ps: Poly) -> Poly:
if not ps:
raise ValueError("pmul needs at least one factor")
r = ps[0]
for p in ps[1:]:
out: Poly = {}
for a, x in r.items():
for b, y in p.items():
k = tuple(i + j for i, j in zip(a, b))
out[k] = out.get(k, Fraction(0)) + x * y
r = clean(out)
return r
def ppow(p: Poly, k: int, nv: int) -> Poly:
r = {tuple([0] * nv): Fraction(1)}
for _ in range(k):
r = pmul(r, p)
return r
def pdiff(p: Poly, i: int) -> Poly:
out: Poly = {}
for e, c in p.items():
if e[i]:
k = e[:i] + (e[i] - 1,) + e[i + 1:]
out[k] = out.get(k, Fraction(0)) + c * e[i]
return clean(out)
def tdeg(p: Poly) -> int:
return max((sum(e) for e in p), default=0)
def derive(p: Poly, d: Poly, N: list) -> Poly:
"""The scaled derivation D(p) = d dp/dx + sum N_i dp/dt_i."""
terms = [pmul(d, pdiff(p, 0))] if pdiff(p, 0) else []
for i, Ni in enumerate(N):
dp = pdiff(p, i + 1)
if dp:
terms.append(pmul(Ni, dp))
return padd(*terms) if terms else {}
def divides(a: Poly, b: Poly) -> bool:
"""True if a divides b in Q[x, t] (division by a single polynomial)."""
if not b:
return True
key = lambda e: e # lex order on (x, t_1, ..., t_m)
la = max(a, key=key)
r = dict(b)
while r:
lb = max(r, key=key)
if any(i < j for i, j in zip(lb, la)):
return False # the leading monomial of a does not divide it
q = tuple(i - j for i, j in zip(lb, la))
c = r[lb] / a[la]
r = padd(r, pscale(pmul(a, {q: Fraction(1)}), -c))
return True
def ansatz_denominator(Qf: list, d: Poly, N: list, nv: int) -> Poly:
"""B = prod q^(e-1) over normal q, times prod q^e over special q."""
B = {tuple([0] * nv): Fraction(1)}
for q, e in Qf:
special = divides(q, derive(q, d, N))
B = pmul(B, ppow(q, e if special else e - 1, nv))
return B
def _solve(cols: list, rhs: Poly):
"""Solve sum_k z_k cols[k] = rhs coefficientwise; None if inconsistent."""
mons = sorted(set(rhs).union(*[set(c) for c in cols]))
rows = [[c.get(m, Fraction(0)) for c in cols] + [rhs.get(m, Fraction(0))]
for m in mons]
nc = len(cols)
piv_cols = []
r = 0
for c in range(nc):
p = next((i for i in range(r, len(rows)) if rows[i][c]), None)
if p is None:
continue
rows[r], rows[p] = rows[p], rows[r]
inv = 1 / rows[r][c]
rows[r] = [v * inv for v in rows[r]]
for i in range(len(rows)):
if i != r and rows[i][c]:
f = rows[i][c]
rows[i] = [a - f * b for a, b in zip(rows[i], rows[r])]
piv_cols.append(c)
r += 1
if any(row[nc] for row in rows[r:]):
return None
z = [Fraction(0)] * nc # free unknowns are set to 0
for i, c in enumerate(piv_cols):
z[c] = rows[i][nc]
return z
def risch_norman(P: Poly, Qf: list, d: Poly, N: list,
extra_logs: tuple = (), deg_bound=None):
"""Try to integrate P / prod(q^e) in the tower given by d and N."""
nv = len(N) + 1
one = {tuple([0] * nv): Fraction(1)}
P = clean({e: Fraction(c) for e, c in P.items()})
if not P:
return {}, one, []
Q = one
for q, e in Qf:
Q = pmul(Q, ppow(q, e, nv))
B = ansatz_denominator(Qf, d, N, nv) # @line L1
V = [q for q, _ in Qf if tdeg(q) > 0] + list(extra_logs) # @line L2
if deg_bound is None: # @line L3
deg_bound = max(tdeg(P) - tdeg(Q), 0) + tdeg(B) + 1
mons = [e for e in iproduct(range(deg_bound + 1), repeat=nv)
if sum(e) <= deg_bound] # @line L4
W = one
for v in V:
W = pmul(W, v)
DB = derive(B, d, N)
QW = pmul(Q, W)
cols = [] # @line L5
for m in mons:
mono = {m: Fraction(1)}
cols.append(pmul(padd(pmul(derive(mono, d, N), B),
pscale(pmul(mono, DB), -1)) or {m: Fraction(0)}, QW))
B2Q = pmul(B, B, Q)
for j, v in enumerate(V):
rest = one
for k, w in enumerate(V):
if k != j:
rest = pmul(rest, w)
cols.append(pmul(derive(v, d, N) or {tuple([0] * nv): Fraction(0)}, B2Q, rest))
rhs = pmul(d, P, B, B, W) # @line L6
z = _solve([clean(c) for c in cols], rhs) # @line L7
if z is None:
return None
U = clean({m: z[k] for k, m in enumerate(mons)})
logs = [(z[len(mons) + j], v) for j, v in enumerate(V) if z[len(mons) + j]]
return U, B, logs # @line L8
test_risch_norman_heuristic.py 241 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/risch-norman-heuristic/. The authors are listed in prov/XPCGZH.jsonl.
"""Tests for the Risch-Norman heuristic, checked by differentiation.
Integrands are built as derivatives of known answers U0 / B0 + sum c log v
in several towers, so an elementary integral is known to exist. Every
answer is checked by differentiating it with code written here, and the
log coefficients of the normal candidates are compared with the known ones
(they are unique). Known non-elementary integrands must give None.
"""
import random
from fractions import Fraction
import pytest
from risch_norman_heuristic import ansatz_denominator, risch_norman
def F(d):
return {k: Fraction(v) for k, v in d.items() if v}
def add(*ps):
r = {}
for p in ps:
for e, c in p.items():
r[e] = r.get(e, Fraction(0)) + c
return {e: c for e, c in r.items() if c}
def mul(*ps):
r = ps[0]
for p in ps[1:]:
out = {}
for a, x in r.items():
for b, y in p.items():
k = tuple(i + j for i, j in zip(a, b))
out[k] = out.get(k, Fraction(0)) + x * y
r = {e: c for e, c in out.items() if c}
return r
def scale(p, c):
return {e: v * c for e, v in p.items() if v * c}
def deg(p):
return max((sum(e) for e in p), default=0)
def dx(p, tower):
"""The true derivative d/dx of p, as a pair (numerator, denominator d)."""
d, N = tower
out = {}
for e, c in p.items():
if e[0]:
k = (e[0] - 1,) + e[1:]
for m, v in d.items():
kk = tuple(i + j for i, j in zip(k, m))
out[kk] = out.get(kk, Fraction(0)) + c * e[0] * v
for i, Ni in enumerate(N):
if e[i + 1]:
k = e[:i + 1] + (e[i + 1] - 1,) + e[i + 2:]
for m, v in Ni.items():
kk = tuple(a + b for a, b in zip(k, m))
out[kk] = out.get(kk, Fraction(0)) + c * e[i + 1] * v
return {e: c for e, c in out.items() if c}
def check_derivative(P, Qf, tower, ans):
"""d/dx (U / B + sum c log v) == P / prod q^e, cleared of denominators."""
U, B, logs = ans
d = tower[0]
nv = len(tower[1]) + 1
one = {tuple([0] * nv): Fraction(1)}
Q = one
for q, e in Qf:
for _ in range(e):
Q = mul(Q, q)
W = one
for _, v in logs:
W = mul(W, v)
# d (U/B)' = (DU B - U DB) / B^2 and d (log v)' = Dv / v
lhs = mul(add(mul(dx(U, tower), B), scale(mul(U, dx(B, tower)), -1)) or {}, W) \
if U else {}
for j, (c, v) in enumerate(logs):
rest = one
for k, (_, w) in enumerate(logs):
if k != j:
rest = mul(rest, w)
lhs = add(lhs, scale(mul(dx(v, tower) or {tuple([0] * nv): 0}, B, B, rest), c))
lhs = mul(lhs, Q) if lhs else {}
rhs = mul(d, P, B, B, W)
assert add(lhs, scale(rhs, -1)) == {}
def X(nv, i=0, k=1):
return {tuple(k if j == i else 0 for j in range(nv)): Fraction(1)}
def C(nv, c):
return {tuple([0] * nv): Fraction(c)} if c else {}
# towers (d, [N_1, ...]): t_i' = N_i / d
RAT = (F({(0,): 1}), [])
EXP = (F({(0, 0): 1}), [F({(0, 1): 1})]) # t = exp(x)
EXP2 = (F({(0, 0): 1}), [F({(1, 1): 2})]) # t = exp(x^2)
LOG = (F({(1, 0): 1}), [F({(0, 0): 1})]) # t = log(x)
TAN = (F({(0, 0): 1}), [F({(0, 0): 1, (0, 2): 1})]) # t = tan(x)
LOGEXP = (F({(1, 0, 0): 1}), [F({(0, 0, 0): 1}), F({(1, 0, 1): 1})]) # log x, exp x
def test_examples():
# x exp(x^2) -> exp(x^2) / 2
U, B, logs = risch_norman(F({(1, 1): 1}), [], *EXP2)
assert U == F({(0, 1): Fraction(1, 2)}) and B == F({(0, 0): 1}) and logs == []
# log x -> x log x - x
U, B, logs = risch_norman(F({(0, 1): 1}), [], *LOG)
assert U == F({(1, 1): 1, (1, 0): -1}) and logs == []
# tan x -> log(1 + tan^2 x) / 2, with 1 + t^2 as a candidate
v = F({(0, 0): 1, (0, 2): 1})
U, B, logs = risch_norman(F({(0, 1): 1}), [], *TAN, extra_logs=(v,))
assert U == {} and logs == [(Fraction(1, 2), v)]
# 1 / x -> log x; 1 / (x^2 - 1) -> (log(x - 1) - log(x + 1)) / 2
U, B, logs = risch_norman(F({(0,): 1}), [(F({(1,): 1}), 1)], *RAT)
assert U == {} and logs == [(1, F({(1,): 1}))]
fa, fb = F({(1,): 1, (0,): -1}), F({(1,): 1, (0,): 1})
U, B, logs = risch_norman(F({(0,): 1}), [(fa, 1), (fb, 1)], *RAT)
assert U == {} and logs == [(Fraction(1, 2), fa), (Fraction(-1, 2), fb)]
# 1 / (x log x) -> log(log x)
t = F({(0, 1): 1})
U, B, logs = risch_norman(F({(0, 0): 1}), [(F({(1, 0): 1}), 1), (t, 1)], *LOG)
assert U == {} and logs == [(1, t)]
def test_failures():
# exp(x^2), exp(x) / x and 1 / (x^2 - 2) have no answer of the guessed form
for b in range(1, 7):
assert risch_norman(F({(0, 1): 1}), [], *EXP2, deg_bound=b) is None
assert risch_norman(F({(0, 1): 1}), [(F({(1, 0): 1}), 1)], *EXP, deg_bound=b) is None
assert risch_norman(F({(0,): 1}), [(F({(2,): 1, (0,): -2}), 1)], *RAT,
deg_bound=b) is None
def test_edges():
# the zero integrand, a constant, and a polynomial in x only
assert risch_norman({}, [], *EXP) == ({}, F({(0, 0): 1}), [])
U, B, logs = risch_norman(F({(0, 0): 3}), [], *EXP)
assert U == F({(1, 0): 3}) and logs == []
U, B, logs = risch_norman(F({(4,): 5}), [], *RAT)
assert U == F({(5,): 1})
# exp(-x) = 1 / t: t is special, so it stays in the denominator
U, B, logs = risch_norman(F({(0, 0): 1}), [(F({(0, 1): 1}), 1)], *EXP)
assert B == F({(0, 1): 1}) and U == F({(0, 0): -1}) and logs == []
# a deep pole (x - 1)^5 needs many terms in U
q = F({(1,): 1, (0,): -1})
U, B, logs = risch_norman(F({(0,): 1}), [(q, 6)], *RAT)
check_derivative(F({(0,): 1}), [(q, 6)], RAT, (U, B, logs))
# a degree bound that is too small gives None, a larger one succeeds
assert risch_norman(F({(1, 1): 1}), [], *EXP, deg_bound=1) is None
U, B, logs = risch_norman(F({(1, 1): 1}), [], *EXP, deg_bound=2)
assert U == F({(1, 1): 1, (0, 1): -1})
def pools(name):
"""Irreducible factors that are normal in each tower, coprime to d."""
if name == "RAT":
return [F({(1,): 1, (0,): a}) for a in (-3, -1, 2, 5)] + [F({(2,): 1, (0,): 1})]
if name == "LOGEXP":
return [F({(1, 0, 0): 1, (0, 0, 0): a}) for a in (-2, 1, 3)] + \
[F({(0, 1, 0): 1, (0, 0, 0): 1}), F({(0, 0, 1): 1, (1, 0, 0): 1})]
out = [F({(1, 0): 1, (0, 0): a}) for a in (-2, 1, 3)]
if name in ("EXP", "EXP2"):
out += [F({(0, 1): 1, (0, 0): 1}), F({(0, 1): 1, (1, 0): -1})]
else:
out += [F({(0, 1): 1, (0, 0): -2}), F({(0, 1): 1, (1, 0): 1})]
if name == "LOG":
out = out[1:] # x - 2 stays, x itself is d
return out
TOWERS = {"RAT": RAT, "EXP": EXP, "EXP2": EXP2, "LOG": LOG, "TAN": TAN, "LOGEXP": LOGEXP}
def build(rng, name):
d, N = TOWERS[name]
nv = len(N) + 1
one = C(nv, 1)
mons = [e for e in __import__("itertools").product(range(3), repeat=nv) if sum(e) <= 2]
U0 = {}
while deg(U0) == 0:
U0 = F({m: rng.randint(-3, 3) for m in rng.sample(mons, 3)})
pool = pools(name)
rng.shuffle(pool)
nb, nl = rng.randint(0, 2), rng.randint(0, 2)
bf = [(q, rng.randint(1, 2)) for q in pool[:nb]]
if name in ("EXP", "EXP2") and rng.random() < 0.5:
bf.append((X(nv, 1), rng.randint(1, 2))) # the special factor t
lf = [(Fraction(rng.choice([-2, -1, 1, 3])), v) for v in pool[nb:nb + nl]]
B0 = one
for q, e in bf:
for _ in range(e):
B0 = mul(B0, q)
W = one
for _, v in lf:
W = mul(W, v)
num = mul(add(mul(dx(U0, (d, N)), B0), scale(mul(U0, dx(B0, (d, N))), -1)) or {}, W)
for j, (c, v) in enumerate(lf):
rest = one
for k, (_, w) in enumerate(lf):
if k != j:
rest = mul(rest, w)
num = add(num, scale(mul(dx(v, (d, N)), B0, B0, rest), c))
# P / Q with Q = d B0^2 W, the factors listed with their multiplicities
Qf = [(q, 2 * e) for q, e in bf] + [(v, 1) for _, v in lf]
if deg(d) > 0:
Qf.append((d, 1))
return num, Qf, (d, N), U0, B0, lf, bf
@pytest.mark.parametrize("name", list(TOWERS))
def test_random_constructed(name):
rng = random.Random(sum(map(ord, name)))
for _ in range({"RAT": 30, "TAN": 25, "LOGEXP": 8}.get(name, 15)):
P, Qf, tower, U0, B0, lf, bf = build(rng, name)
if not P:
continue
nv = len(tower[1]) + 1
B = ansatz_denominator(Qf, tower[0], tower[1], nv)
bound = deg(U0) + deg(B) - deg(B0)
ans = risch_norman(P, Qf, *tower, deg_bound=bound)
assert ans is not None
check_derivative(P, Qf, tower, ans)
got = {tuple(sorted(v.items())): c for c, v in ans[2]}
for c, v in lf:
assert got.get(tuple(sorted(v.items()))) == c
for q, _ in bf:
if q != X(nv, 1):
assert tuple(sorted(q.items())) not in got
What the next level needs
- the tests did not run in this build