Polynomial long division in Python
Repeatedly subtracts multiples of the divisor to eliminate leading terms; quadratic in the degrees.
- File
polynomial_long_division.py, 55 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/polynomial-long-division/. The authors are listed in prov/0Q78BZ.jsonl.
"""Polynomial long division with remainder.
A polynomial is a list of coefficients, lowest degree first: [a0, a1, ..., an].
The zero polynomial is []. The coefficients can be int, Fraction, or any type
with +, -, *, == 0 and 1 / a for a unit a (for example an element of Z/pZ).
"""
from typing import Any, Sequence
def _strip(a: Sequence[Any]) -> list[Any]:
"""Return a copy of a without zero coefficients at the top."""
b = list(a)
while b and b[-1] == 0:
b.pop()
return b
def _unit_inverse(a: Any) -> Any:
"""Return the inverse of the unit a. For an int, only 1 and -1 are units."""
if a == 1 or a == -1:
return a # a * a = 1, and a keeps its type
if isinstance(a, int):
raise ValueError("the leading coefficient of the divisor is not a unit")
return 1 / a
def polynomial_long_division(f: Sequence[Any], g: Sequence[Any]) -> tuple[list[Any], list[Any]]:
"""Return (q, r) with f = q*g + r and deg r < deg g.
Precondition: g is not zero and its leading coefficient is a unit of the
coefficient ring (every nonzero element over a field; 1 or -1 over Z).
Cost: at most (2n + 1)(m - n + 1) ring operations and one inversion, for
m = deg f and n = deg g.
"""
f = _strip(f) # @line L1
g = _strip(g)
if not g:
raise ZeroDivisionError("division by the zero polynomial")
m, n = len(f) - 1, len(g) - 1
if m < n: # @line L2
return [], f
u = _unit_inverse(g[n]) # @line L3
r = list(f) # @line L4
q: list[Any] = [None] * (m - n + 1)
for k in range(m - n, -1, -1): # @line L5
c = r[n + k] * u # @line L6
q[k] = c # @line L7
for j in range(n): # @line L8
r[j + k] = r[j + k] - c * g[j] # @line L9
r[n + k] = c - c # @line L10: the top term cancels, since c * g[n] = r[n + k]
return _strip(q), _strip(r[:n]) # @line L11
test_polynomial_long_division.py 189 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/polynomial-long-division/. The authors are listed in prov/0Q78BZ.jsonl.
import random
from fractions import Fraction
import pytest
from polynomial_long_division import polynomial_long_division
class Fp:
"""An element of Z/pZ, for the tests."""
__slots__ = ("v", "p")
def __init__(self, v: int, p: int) -> None:
self.v = v % p
self.p = p
def _c(self, o: object) -> int:
return o.v if isinstance(o, Fp) else o # type: ignore[return-value]
def __add__(self, o: object) -> "Fp":
return Fp(self.v + self._c(o), self.p)
__radd__ = __add__
def __sub__(self, o: object) -> "Fp":
return Fp(self.v - self._c(o), self.p)
def __rsub__(self, o: object) -> "Fp":
return Fp(self._c(o) - self.v, self.p)
def __mul__(self, o: object) -> "Fp":
return Fp(self.v * self._c(o), self.p)
__rmul__ = __mul__
def __neg__(self) -> "Fp":
return Fp(-self.v, self.p)
def __truediv__(self, o: object) -> "Fp":
return Fp(self.v * pow(self._c(o), -1, self.p), self.p)
def __rtruediv__(self, o: object) -> "Fp":
return Fp(self._c(o) * pow(self.v, -1, self.p), self.p)
def __eq__(self, o: object) -> bool:
return (self.v - self._c(o)) % self.p == 0
def __hash__(self) -> int:
return hash(self.v)
def __repr__(self) -> str:
return f"Fp({self.v}, {self.p})"
def strip(a):
a = list(a)
while a and a[-1] == 0:
a.pop()
return a
def mul(a, b):
"""Schoolbook product, the reference."""
if not a or not b:
return []
c = [a[0] * 0] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
for j, y in enumerate(b):
c[i + j] = c[i + j] + x * y
return strip(c)
def add(a, b):
n = max(len(a), len(b))
return strip([(a[i] if i < len(a) else 0) + (b[i] if i < len(b) else 0) for i in range(n)])
def ref_divmod(f, g, inv_lc):
"""Independent reference: the quotient from the reversal and a power series inverse."""
f, g = strip(f), strip(g)
m, n = len(f) - 1, len(g) - 1
if m < n:
return [], f
k = m - n + 1
rf, rg = f[::-1], g[::-1]
h = [inv_lc]
for i in range(1, k):
s = 0
for j in range(1, min(i, n) + 1):
s = s + rg[j] * h[i - j]
h.append(-s * inv_lc)
rq = [0] * k
for i in range(k):
for j in range(i + 1):
if i - j < len(h) and j < len(rf):
rq[i] = rq[i] + rf[j] * h[i - j]
q = strip(rq[::-1])
r = add(f, [-c for c in mul(q, g)])
return q, r
def check(f, g, inv_lc):
q, r = polynomial_long_division(f, g)
assert q == strip(q) and r == strip(r)
assert len(r) < len(strip(g))
assert add(mul(q, g), r) == strip(f)
rq, rr = ref_divmod(f, g, inv_lc)
assert q == rq and r == rr
def test_example() -> None:
q, r = polynomial_long_division([5, -1, 0, 3, 2], [1, 1, 1])
assert q == [-3, 1, 2] and r == [8, 1]
def test_example_f7() -> None:
p = 7
f = [Fp(c, p) for c in [1, 0, 3, 0, 1, 4]]
g = [Fp(c, p) for c in [3, 0, 2]]
q, r = polynomial_long_division(f, g)
assert q == [6, 4, 4, 2] and r == [4, 2]
def test_edges() -> None:
assert polynomial_long_division([], [1, 1]) == ([], [])
assert polynomial_long_division([1, 1], [0, 0, 0, 1]) == ([], [1, 1])
assert polynomial_long_division([4, 5, 6], [1]) == ([4, 5, 6], [])
assert polynomial_long_division([4, 5, 6], [-1]) == ([-4, -5, -6], [])
assert polynomial_long_division([3, 2, 0, 0], [1, 1, 0]) == ([2], [1])
assert polynomial_long_division([7], [7 * Fraction(1)]) == ([1], [])
q, r = polynomial_long_division([Fraction(1), 0, 1], [Fraction(2), Fraction(3)])
assert add(mul(q, [2, 3]), r) == [1, 0, 1] and len(r) <= 1
# a divisor of degree 0 that is a zero after stripping
with pytest.raises(ZeroDivisionError):
polynomial_long_division([1, 2], [0, 0])
with pytest.raises(ZeroDivisionError):
polynomial_long_division([1, 2], [])
# over Z the leading coefficient must be a unit
with pytest.raises(ValueError):
polynomial_long_division([0, 0, 1], [0, 2])
def test_division_by_linear_gives_value() -> None:
f = [5, -1, 0, 3, 2]
for a in range(-5, 6):
_, r = polynomial_long_division(f, [-a, 1])
val = sum(c * a**i for i, c in enumerate(f))
assert r == strip([val])
def test_random_integers_unit_leading() -> None:
rng = random.Random(4061)
for _ in range(600):
f = [rng.randint(-20, 20) for _ in range(rng.randint(0, 12))]
g = [rng.randint(-20, 20) for _ in range(rng.randint(0, 7))] + [rng.choice([1, -1])]
g += [0] * rng.randint(0, 2) # zero top coefficients that must be stripped
check(f, g, strip(g)[-1])
def test_random_fractions() -> None:
rng = random.Random(4062)
for _ in range(400):
f = [Fraction(rng.randint(-9, 9), rng.randint(1, 4)) for _ in range(rng.randint(0, 10))]
g = [Fraction(rng.randint(-9, 9), rng.randint(1, 4)) for _ in range(rng.randint(0, 6))]
g.append(Fraction(rng.choice([-3, -2, 1, 2, 5]), rng.randint(1, 3)))
check(f, g, 1 / g[-1])
@pytest.mark.parametrize("p", [2, 3, 7, 998244353])
def test_random_mod_p(p: int) -> None:
rng = random.Random(4063 + p)
for _ in range(300):
f = [Fp(rng.randrange(p), p) for _ in range(rng.randint(0, 12))]
g = [Fp(rng.randrange(p), p) for _ in range(rng.randint(0, 6))]
g.append(Fp(rng.randrange(1, p), p))
check(f, g, 1 / g[-1])
def test_very_different_degrees() -> None:
rng = random.Random(4064)
f = [Fraction(rng.randint(-5, 5)) for _ in range(200)] + [Fraction(1)]
check(f, [Fraction(3), Fraction(-2)], Fraction(-1, 2))
g = [Fraction(rng.randint(-5, 5)) for _ in range(150)] + [Fraction(7)]
check(f, g, Fraction(1, 7))
check(g, f, Fraction(1))
What the next level needs
- the tests did not run in this build