Bairstow's method in Python
Newton's method on the remainder of p modulo a trial quadratic x^2 - ux - v. It finds real quadratic factors, so complex conjugate roots come out in real arithmetic, quadratically near a factor.
- File
bairstow_method.py, 149 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/bairstow-method/. The authors are listed in prov/J2HSDJ.jsonl.
"""Bairstow's method: Newton's method on the remainder of p modulo x^2 - u x - v.
Coefficients are real and listed highest degree first: a[0] x^n + ... + a[n].
"""
import cmath
import math
import random
from dataclasses import dataclass, field
from typing import List, Sequence, Tuple
@dataclass
class FactorResult:
u: float
v: float
quotient: List[float] # b_0, ..., b_(n-2): p = (x^2 - u x - v) q + remainder
iterations: int
evaluations: int # passes of the double synthetic division
status: str # "converged", "exact", "singular", "max-iter"
# "converged": the step test held, or the remainder is below its rounding floor
@dataclass
class RootsResult:
roots: List[complex]
factors: List[Tuple[float, float]] = field(default_factory=list)
evaluations: int = 0
restarts: int = 0
status: str = "converged" # or "failed"
def divide(a: Sequence[float], u: float, v: float) -> List[float]:
"""b_k = a_k + u b_(k-1) + v b_(k-2): p = (x^2 - u x - v)(b_0 x^(n-2) + ... + b_(n-2)) + b_(n-1)(x - u) + b_n."""
b = []
b1 = b2 = 0.0
for ak in a:
bk = ak + u * b1 + v * b2
b.append(bk)
b1, b2 = bk, b1
return b
def remainder_floor(a: Sequence[float], u: float, v: float) -> Tuple[float, float]:
"""Running bounds gamma * mu_k for the rounding in b_(n-1) and b_n, with
mu_k = |a_k| + |u| mu_(k-1) + |v| mu_(k-2) and gamma = 4 (n + 1) 2^-53."""
m1 = m2 = 0.0
mus = []
for ak in a:
mk = abs(ak) + abs(u) * m1 + abs(v) * m2
mus.append(mk)
m1, m2 = mk, m1
gamma = 4 * len(a) * 2.0 ** -53
return gamma * mus[-2], gamma * mus[-1]
def quadratic_roots(u: float, v: float) -> List[complex]:
"""Roots of x^2 - u x - v, computed without cancellation."""
disc = u * u + 4 * v
if disc >= 0:
x1 = (u + math.copysign(math.sqrt(disc), u)) / 2
x2 = -v / x1 if x1 != 0 else 0.0
return [complex(x1), complex(x2)]
im = math.sqrt(-disc) / 2
return [complex(u / 2, im), complex(u / 2, -im)]
def find_factor(a: Sequence[float], u: float, v: float, tol: float = 1e-14,
max_iter: int = 100) -> FactorResult:
"""Newton's method on (b_(n-1), b_n) as functions of (u, v). Needs deg p >= 3."""
n = len(a) - 1
if n < 3:
raise ValueError("find_factor needs degree at least 3")
evals = 0
for k in range(1, max_iter + 1): # @line L1
b = divide(a, u, v) # @line L2
evals += 1
if b[n - 1] == 0 and b[n] == 0: # @line L3
return FactorResult(u, v, b[:n - 1], k - 1, evals, "exact")
f1, f0 = remainder_floor(a, u, v)
if abs(b[n - 1]) <= f1 and abs(b[n]) <= f0:
return FactorResult(u, v, b[:n - 1], k - 1, evals, "converged")
c = divide(b[:n], u, v) # @line L4
det = c[n - 2] * c[n - 2] - c[n - 1] * c[n - 3] # @line L5
if det == 0 or not math.isfinite(det):
return FactorResult(u, v, b[:n - 1], k - 1, evals, "singular")
du = (b[n] * c[n - 3] - b[n - 1] * c[n - 2]) / det # @line L6
dv = (b[n - 1] * c[n - 1] - b[n] * c[n - 2]) / det
u, v = u + du, v + dv # @line L7
if not (math.isfinite(u) and math.isfinite(v)):
return FactorResult(u, v, b[:n - 1], k, evals, "singular")
if abs(du) + abs(dv) <= tol * (1 + abs(u) + abs(v)): # @line L8
b = divide(a, u, v)
evals += 1
return FactorResult(u, v, b[:n - 1], k, evals, "converged")
b = divide(a, u, v)
return FactorResult(u, v, b[:n - 1], max_iter, evals + 1, "max-iter") # @line L9
def cauchy_radius(a: Sequence[float]) -> float:
return 1 + max(abs(x / a[0]) for x in a[1:])
def bairstow_roots(coeffs: Sequence[float], tol: float = 1e-14, max_iter: int = 100,
max_restarts: int = 30, seed: int = 1) -> RootsResult:
"""All roots of a real polynomial by Bairstow's method with deflation.
The first start is Lin's (u, v) = (-a_(n-1)/a_(n-2), -a_n/a_(n-2)); a
failed factor search restarts from pseudo-random starts inside the
Cauchy radius, at most max_restarts times per factor.
"""
a = [float(x) for x in coeffs]
if not all(math.isfinite(x) for x in a):
raise ValueError("coefficients must be finite")
while a and a[0] == 0:
a.pop(0)
if len(a) == 0:
raise ValueError("the zero polynomial has every number as a root")
rng = random.Random(seed)
out = RootsResult([])
while len(a) - 1 >= 3: # @line L10
n = len(a) - 1
if a[n - 2] != 0:
u, v = -a[n - 1] / a[n - 2], -a[n] / a[n - 2]
else:
u, v = 0.0, 0.0
r = cauchy_radius(a)
for attempt in range(max_restarts + 1):
res = find_factor(a, u, v, tol, max_iter)
out.evaluations += res.evaluations
if res.status in ("converged", "exact"):
break
out.restarts += 1
z = cmath.rect(r * rng.random(), math.pi * rng.random())
u, v = 2 * z.real, -abs(z) ** 2
else:
out.status = "failed"
return out
out.factors.append((res.u, res.v))
out.roots.extend(quadratic_roots(res.u, res.v)) # @line L11
a = res.quotient # @line L12
if len(a) == 3: # @line L13
out.roots.extend(quadratic_roots(-a[1] / a[0], -a[2] / a[0]))
elif len(a) == 2:
out.roots.append(complex(-a[1] / a[0]))
return out
test_bairstow_method.py 142 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/bairstow-method/. The authors are listed in prov/J2HSDJ.jsonl.
import math
import random
import unittest
from bairstow_method import bairstow_roots, divide, find_factor, quadratic_roots
def from_roots(rs):
p = [1 + 0j]
for r in rs:
p = [x - r * y for x, y in zip(p + [0], [0] + p)]
return [x.real for x in p]
def match_error(found, expected):
left = list(found)
err = 0.0
for r in expected:
j = min(range(len(left)), key=lambda i: abs(left[i] - r))
err = max(err, abs(left[j] - r))
left.pop(j)
return err
QUARTIC = [1, 0, 3, 6, 10] # (x^2 + 2x + 2)(x^2 - 2x + 5)
class TestBairstow(unittest.TestCase):
def test_division_identity(self):
a = [2.0, -3.0, 5.0, 7.0, -11.0, 13.0]
u, v = 0.75, -1.5
b = divide(a, u, v)
n = len(a) - 1
for x in (-2.0, -0.5, 0.0, 1.0, 3.0):
p = sum(c * x ** (n - i) for i, c in enumerate(a))
q = sum(c * x ** (n - 2 - i) for i, c in enumerate(b[:n - 1]))
self.assertAlmostEqual(p, (x * x - u * x - v) * q + b[n - 1] * (x - u) + b[n], places=9)
def test_example_quartic(self):
r = find_factor(QUARTIC, -2.0, -10 / 3)
self.assertIn(r.status, ("converged", "exact"))
self.assertLess(abs(r.u + 2) + abs(r.v + 2), 1e-15)
self.assertLessEqual(r.evaluations, 8)
roots = bairstow_roots(QUARTIC)
self.assertLess(match_error(roots.roots, [-1 + 1j, -1 - 1j, 1 + 2j, 1 - 2j]), 1e-14)
def test_quadratic_convergence(self):
u, v = -2.0, -10 / 3
errs = []
for _ in range(5):
errs.append(abs(u + 2) + abs(v + 2))
r = find_factor(QUARTIC, u, v, tol=0, max_iter=1)
u, v = r.u, r.v
# e_{k+1} / e_k^2 stays bounded and the observed order approaches 2
for e0, e1, e2 in zip(errs, errs[1:], errs[2:]):
if e2 > 1e-13:
self.assertGreater(math.log(e2 / e1) / math.log(e1 / e0), 1.8)
self.assertLess(errs[4], 1e-10)
def test_double_root_factor_still_quadratic(self):
# factor (x - 1)^2 is coprime to the quotient, so the Jacobian is nonsingular
p = from_roots([1, 1, complex(-0.5, 0.75 ** 0.5), complex(-0.5, -0.75 ** 0.5), -3])
u, v = 1.9, -0.9
errs = []
for _ in range(3):
r = find_factor(p, u, v, tol=0, max_iter=1)
u, v = r.u, r.v
errs.append(abs(u - 2) + abs(v + 1))
self.assertLess(errs[1], 10 * errs[0] ** 2)
self.assertLess(errs[2], 1e-7)
def test_exact_factor_at_start(self):
r = find_factor(QUARTIC, -2.0, -2.0)
self.assertEqual((r.status, r.iterations, r.evaluations), ("exact", 0, 1))
self.assertEqual(r.quotient, [1.0, -2.0, 5.0])
def test_root_at_zero(self):
res = bairstow_roots([1, -6, 11, -6, 0]) # x (x - 1)(x - 2)(x - 3)
self.assertLess(match_error(res.roots, [0, 1, 2, 3]), 1e-12)
def test_singular_jacobian(self):
# u = v = 0 on x^4 + 1: b = a, c = (1, 0, 0, 0) gives det = 0
r = find_factor([1, 0, 0, 0, 1], 0.0, 0.0)
self.assertEqual(r.status, "singular")
res = bairstow_roots([1, 0, 0, 0, 1])
self.assertEqual(res.status, "converged")
self.assertGreaterEqual(res.restarts, 1)
s = 0.5 ** 0.5
self.assertLess(match_error(res.roots, [complex(s, s), complex(s, -s), complex(-s, s), complex(-s, -s)]), 1e-13)
def test_max_iter_and_failure(self):
r = find_factor([1, 0, 3, 6, 10], 40.0, 30.0, max_iter=2)
self.assertEqual((r.status, r.iterations, r.evaluations), ("max-iter", 2, 3))
res = bairstow_roots([1, 0, 3, 6, 10], max_iter=1, max_restarts=2)
self.assertEqual(res.status, "failed")
self.assertEqual(res.restarts, 3)
def test_low_degree_and_bad_input(self):
self.assertEqual(bairstow_roots([2, -4]).roots, [2 + 0j])
self.assertLess(match_error(bairstow_roots([0, 1, 0, 1]).roots, [1j, -1j]), 1e-15)
self.assertEqual(sorted(z.real for z in quadratic_roots(1e8, -1.0)), [1e-8, 1e8])
with self.assertRaises(ValueError):
bairstow_roots([1, float("nan"), 2])
with self.assertRaises(ValueError):
bairstow_roots([1, math.inf, 2])
with self.assertRaises(ValueError):
bairstow_roots([0, 0])
with self.assertRaises(ValueError):
find_factor([1, 2, 3], 0.0, 0.0)
def test_wilkinson_10(self):
res = bairstow_roots(from_roots(range(1, 11)))
self.assertEqual(res.status, "converged")
self.assertLess(match_error(res.roots, range(1, 11)), 1e-6)
def test_random_polynomials(self):
rng = random.Random(20261108)
total = 0
for _ in range(600):
n = rng.randint(3, 8)
rs = []
while len(rs) < n:
if n - len(rs) >= 2 and rng.random() < 0.5:
z = complex(rng.uniform(-3, 3), rng.uniform(0.3, 3))
cand = [z, z.conjugate()]
else:
cand = [complex(rng.uniform(-3, 3), 0)]
if all(abs(c - r) > 0.3 for c in cand for r in rs):
rs += cand
res = bairstow_roots(from_roots(rs))
self.assertEqual(res.status, "converged")
self.assertEqual(len(res.roots), n)
self.assertLess(match_error(res.roots, rs), 1e-9)
total += res.evaluations
self.assertLess(total / 600, 60)
if __name__ == "__main__":
unittest.main()
What the next level needs
- the tests did not run in this build