Aberth-Ehrlich method in Python
Simultaneous iteration for all roots of a polynomial; each approximation takes a Newton step on p divided by the factors of the others. Order 3 at simple roots.
- File
aberth_ehrlich_method.py, 98 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/aberth-ehrlich-method/. The authors are listed in prov/4VYZ43.jsonl.
"""Aberth-Ehrlich iteration: all roots of a polynomial at once, order 3 at simple roots."""
import cmath
import math
from dataclasses import dataclass, field
from typing import List, Optional, Sequence
U = 2.0 ** -53
@dataclass
class AberthResult:
roots: List[complex]
iterations: int
evaluations: int # evaluations of p; p' is evaluated alongside, the same number of times
converged: bool
history: List[List[complex]] = field(default_factory=list)
def horner2(coeffs: Sequence[complex], z: complex):
"""Return p(z), p'(z) and a running rounding bound 4 u mu for p(z) (Higham, Section 5.1)."""
p = 0j
dp = 0j
mu = 0.0
az = abs(z)
for a in coeffs:
dp = dp * z + p
p = p * z + a
mu = mu * az + abs(p)
return p, dp, 4 * U * mu
def start_points(coeffs: Sequence[complex]) -> List[complex]:
"""n points on a circle about c = -a_1/(n a_0), radius 1.5 (|p(c)|/|a_0|)^(1/n), angles offset by 0.4."""
n = len(coeffs) - 1
c = -coeffs[1] / (n * coeffs[0])
r = abs(horner2(coeffs, c)[0] / coeffs[0]) ** (1.0 / n)
r = 1.5 * r if r > 0 else 1.0
return [c + r * cmath.exp(1j * (2 * math.pi * k / n + 0.4)) for k in range(n)]
def aberth(coeffs: Sequence[complex], tol: float = 1e-14, max_iter: int = 500,
start: Optional[Sequence[complex]] = None, keep_history: bool = False) -> AberthResult:
"""All roots of p(z) = coeffs[0] z^n + ... + coeffs[n] by the Aberth correction.
The correction is p/(p' - p S) with S = sum over j != i of 1/(z_i - z_j),
the same as N/(1 - N S) with the Newton step N = p/p', but defined when
p'(z_i) = 0. A root is frozen when |p(z_i)| is at the rounding floor of
Horner's rule or its correction is at most tol |z_i|. Raises ValueError
on bad input and ArithmeticError on non-finite values, coinciding
approximations or a zero denominator.
"""
coeffs = [complex(a) for a in coeffs]
if len(coeffs) < 2 or coeffs[0] == 0:
raise ValueError("need degree >= 1 and a nonzero leading coefficient")
if any(not cmath.isfinite(a) for a in coeffs):
raise ArithmeticError("coefficients must be finite")
n = len(coeffs) - 1
z = start_points(coeffs) if start is None else [complex(s) for s in start] # @line L1
if len(z) != n or len(set(z)) != n:
raise ValueError("need n distinct starting points")
active = [True] * n
evals = 0
history = [list(z)] if keep_history else []
for it in range(1, max_iter + 1): # @line L2
new = list(z)
for i in range(n): # @line L3
if not active[i]:
continue
p, dp, floor = horner2(coeffs, z[i]) # @line L4
evals += 1
if abs(p) <= floor:
active[i] = False
continue
s = 0j # @line L5
for j in range(n):
if j != i:
if z[i] == z[j]:
raise ArithmeticError("two approximations coincide")
s += 1 / (z[i] - z[j])
den = dp - p * s # @line L6
if den == 0:
raise ArithmeticError("zero denominator in the Aberth correction")
w = p / den
new[i] = z[i] - w
if not cmath.isfinite(new[i]):
raise ArithmeticError("non-finite iterate")
if abs(w) <= tol * abs(new[i]): # @line L7
active[i] = False
z = new # @line L8
if keep_history:
history.append(list(z))
if not any(active): # @line L9
return AberthResult(z, it, evals, True, history)
return AberthResult(z, max_iter, evals, False, history) # @line L10
test_aberth_ehrlich_method.py 139 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/aberth-ehrlich-method/. The authors are listed in prov/4VYZ43.jsonl.
import cmath
import math
import random
import unittest
from fractions import Fraction
from aberth_ehrlich_method import aberth, horner2, start_points
CUBIC_ROOTS = [1.5213797068045676, complex(-0.7606898534022838, 0.8578736265951786),
complex(-0.7606898534022838, -0.8578736265951786)]
def poly(roots):
c = [1 + 0j]
for r in roots:
c = c + [0j]
for i in range(len(c) - 1, 0, -1):
c[i] -= r * c[i - 1]
return c
def max_err(z, roots):
return max(max(min(abs(a - b) for a in z) for b in roots), max(min(abs(a - b) for b in roots) for a in z))
class TestAberth(unittest.TestCase):
def test_integer_roots(self):
r = aberth(poly([1, 2, 3, 4, 5]))
self.assertTrue(r.converged)
self.assertEqual(r.iterations, 6)
self.assertEqual(r.evaluations, 30)
self.assertLess(max_err(r.roots, [1, 2, 3, 4, 5]), 1e-12)
def test_cubic_with_kerner_start_is_cubically_convergent(self):
r = aberth([1, 0, -1, -2], start=[(0.4 + 0.9j) ** k for k in range(3)], keep_history=True)
self.assertTrue(r.converged)
e = [max_err(h, CUBIC_ROOTS) for h in r.history]
c1 = e[3] / e[2] ** 3
c2 = e[4] / e[3] ** 3
self.assertLess(e[4], 1e-12)
self.assertLess(max(c1, c2) / min(c1, c2), 3)
def test_roots_of_unity_and_complex(self):
for n in (2, 5, 8, 15):
roots = [cmath.exp(2j * math.pi * k / n) for k in range(n)]
r = aberth(poly(roots))
self.assertTrue(r.converged)
self.assertLess(max_err(r.roots, roots), 1e-12)
roots = [1 + 2j, 1 - 2j, -3, 0.5j]
self.assertLess(max_err(aberth(poly(roots)).roots, roots), 1e-12)
def test_widely_spread_roots(self):
roots = [1e-3, 1, 1e3]
r = aberth(poly(roots))
for x in roots:
self.assertLess(min(abs(z - x) for z in r.roots), 1e-12 * x)
def test_multiple_root_ratio(self):
roots = [1, 1, 1, -2]
r = aberth(poly(roots), keep_history=True)
self.assertTrue(r.converged)
self.assertLess(max_err(r.roots, roots), 1e-4)
e = [max_err(h, roots) for h in r.history]
for k in range(5, 15):
self.assertAlmostEqual(e[k + 1] / e[k], 0.5, delta=0.05)
def test_exact_coefficients_from_fractions(self):
roots = [Fraction(1, 3), Fraction(-5, 7), Fraction(9, 4)]
c = [Fraction(1)]
for x in roots:
c = c + [Fraction(0)]
for i in range(len(c) - 1, 0, -1):
c[i] -= x * c[i - 1]
r = aberth([float(a) for a in c])
self.assertLess(max_err(r.roots, [float(x) for x in roots]), 1e-14)
def test_horner_derivative(self):
p, dp, _ = horner2([1, 0, -1, -2], 2.0)
self.assertEqual((p, dp), (4, 11))
def test_zero_derivative_at_iterate(self):
# p'(0) = 0 for z^2 - 1; the Newton step is undefined but p/(p' - pS) is not.
r = aberth([1, 0, -1], start=[0.0, 3.0])
self.assertTrue(r.converged)
self.assertLess(max_err(r.roots, [1, -1]), 1e-14)
def test_exact_root_in_start(self):
r = aberth([1, -3, 2], start=[1.0, 5.0])
self.assertEqual(r.roots[0], 1.0)
self.assertLess(abs(r.roots[1] - 2), 1e-14)
def test_degree_one(self):
r = aberth([2, -3])
self.assertAlmostEqual(r.roots[0], 1.5)
def test_bad_input(self):
with self.assertRaises(ValueError):
aberth([0, 1, 2])
with self.assertRaises(ValueError):
aberth([5])
with self.assertRaises(ValueError):
aberth([1, 0, -1], start=[1.0, 1.0])
with self.assertRaises(ArithmeticError):
aberth([1, math.nan, 2])
with self.assertRaises(ArithmeticError):
aberth([1, 2, math.inf])
def test_max_iter(self):
r = aberth(poly(list(range(1, 11))), max_iter=2)
self.assertFalse(r.converged)
self.assertEqual((r.iterations, r.evaluations), (2, 20))
def test_wilkinson_20_backward_floor(self):
c = poly(list(range(1, 21)))
r = aberth(c)
self.assertTrue(r.converged)
self.assertLess(max_err(r.roots, list(range(1, 21))), 0.5)
def test_start_points_distinct(self):
self.assertEqual(len(set(start_points(poly([1, 2, 3])))), 3)
def test_random_polynomials(self):
rng = random.Random(7)
for _ in range(500):
n = rng.randrange(2, 9)
roots = [complex(rng.gauss(0, 1), rng.gauss(0, 1)) for _ in range(n)]
r = aberth(poly(roots))
self.assertTrue(r.converged)
self.assertLessEqual(r.evaluations, n * r.iterations)
sep = min(abs(a - b) for i, a in enumerate(roots) for b in roots[i + 1:])
if sep > 0.1:
self.assertLess(max_err(r.roots, roots), 1e-8)
if __name__ == "__main__":
unittest.main()
What the next level needs
- the tests did not run in this build