Laguerre's method in Python
A polynomial root method that uses p, p' and p'' at each step. Cubic convergence to simple roots; global and monotone when all roots are real; finds complex roots from real starts.
- File
laguerre_method.py, 116 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/laguerre-method/. The authors are listed in prov/K13QF2.jsonl.
"""Laguerre's method for one root of a polynomial, and all roots with deflation."""
import cmath
import math
from dataclasses import dataclass
from typing import List, Sequence
U = 2.0 ** -53 # unit roundoff of binary64
@dataclass
class LaguerreResult:
"""The returned point, why the loop ended, and the work done."""
root: complex
reason: str # "root", "residual", "converged", "failure" or "kmax reached"
iterations: int
evaluations: int # Horner sweeps, each giving p, p' and p''
@property
def converged(self) -> bool:
return self.reason in ("root", "residual", "converged")
def horner2(coeffs: Sequence[complex], x: complex):
"""Return p(x), p'(x), p''(x) and sum |a_j| |x|^j in one sweep.
coeffs lists a_n, ..., a_0 from the leading coefficient down.
"""
p = coeffs[0]
d = dd = 0
s = abs(coeffs[0])
ax = abs(x)
for c in coeffs[1:]:
dd = dd * x + d
d = d * x + p
p = p * x + c
s = s * ax + abs(c)
return p, d, 2 * dd, s
def _finite(v) -> bool:
return cmath.isfinite(complex(v))
def laguerre(coeffs: Sequence[complex], x0: complex, tol_abs: float = 0.0,
tol_rel: float = 2.0 ** -50, max_iter: int = 100) -> LaguerreResult:
"""Run Laguerre's method on the polynomial with coefficients a_n..a_0.
Precondition: degree n >= 1 and a_n != 0; tol_abs >= 0, tol_rel >= 0;
max_iter >= 1. Cost per pass: one Horner sweep of about 3n complex
multiply-adds and one complex square root.
"""
n = len(coeffs) - 1
if n < 1 or coeffs[0] == 0:
raise ValueError("need a polynomial of degree at least 1")
x = complex(x0) # @line L1
for k in range(1, max_iter + 1): # @line L2
p, d1, d2, s = horner2(coeffs, x) # @line L3
if p == 0: # @line L4
return LaguerreResult(x, "root", k, k)
if not (_finite(p) and _finite(d1) and _finite(d2)):
return LaguerreResult(x, "failure", k, k)
if abs(p) <= 8 * n * U * s: # @line L5
return LaguerreResult(x, "residual", k, k)
g = d1 / p # @line L6
h = g * g - d2 / p
r = cmath.sqrt((n - 1) * (n * h - g * g)) # @line L7
dp, dm = g + r, g - r # @line L8
den = dp if abs(dp) >= abs(dm) else dm
if not _finite(den):
return LaguerreResult(x, "failure", k, k)
if den == 0: # @line L9
# p' = p'' = 0 here: take a fixed-size step in a varying direction
a = (1 + abs(x)) * cmath.exp(1j * k)
else:
a = n / den # @line L10
x_new = x - a
if not _finite(x_new):
return LaguerreResult(x, "failure", k, k)
if abs(a) <= tol_abs + tol_rel * abs(x_new): # @line L11
return LaguerreResult(x_new, "converged", k, k)
x = x_new
return LaguerreResult(x, "kmax reached", max_iter, max_iter) # @line L12
def deflate(coeffs: Sequence[complex], z: complex) -> List[complex]:
"""Divide by (x - z) with synthetic division and drop the remainder."""
out = [coeffs[0]]
for c in coeffs[1:-1]:
out.append(out[-1] * z + c)
return out
def laguerre_all_roots(coeffs: Sequence[complex], polish: bool = True,
max_iter: int = 100) -> List[complex]:
"""All n roots: Laguerre from 0 on the deflated polynomials, then polish.
Each root found on a deflated polynomial is refined by a few Laguerre
steps on the original polynomial, which removes most of the error that
deflation introduces. Cost: O(n^2) per sweep of passes over all roots.
"""
coeffs = [complex(c) for c in coeffs]
work = list(coeffs)
roots = []
while len(work) > 1:
r = laguerre(work, 0.0, max_iter=max_iter)
z = r.root
if polish:
z = laguerre(coeffs, z, max_iter=10).root
roots.append(z)
work = deflate(work, z)
return roots
test_laguerre_method.py 141 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/laguerre-method/. The authors are listed in prov/K13QF2.jsonl.
import math
import random
from decimal import Decimal, getcontext
import pytest
from laguerre_method import U, deflate, horner2, laguerre, laguerre_all_roots
Q4 = [1, -10, 35, -50, 24] # (x - 1)(x - 2)(x - 3)(x - 4)
def poly_from_roots(zs):
c = [1 + 0j]
for z in zs:
c = [a - z * b for a, b in zip(c + [0], [0] + c)]
return c
def test_example_real_roots_from_zero():
r = laguerre(Q4, 0)
assert r.reason == "root" and r.iterations == 4 and r.evaluations == 4
assert r.root == 0.9999999999999998 and r.root.imag == 0
def test_observed_cubic_order_in_high_precision():
getcontext().prec = 200
c = [Decimal(a) for a in Q4]
n = 4
x = Decimal(0)
errs = []
for _ in range(5):
errs.append(abs(x - 1))
p = d = dd = Decimal(0)
p = c[0]
for a in c[1:]:
dd = dd * x + d
d = d * x + p
p = p * x + a
g = d / p
h = g * g - 2 * dd / p
r = ((n - 1) * (n * h - g * g)).sqrt()
den = g + r if abs(g + r) >= abs(g - r) else g - r
x = x - n / den
ratios = [errs[k + 1] / errs[k] ** 3 for k in range(2, 4)]
assert all(abs(q - Decimal("0.12037")) < Decimal("1e-5") for q in ratios)
assert errs[1] < errs[0] and x < 1 # monotone from the left
def test_real_rooted_moves_to_neighbouring_root():
# start between 2 and 3: the iterates stay in [2, 3]
for x0, want in [(2.4, 2), (2.6, 3), (10.0, 4), (-7.0, 1), (3.5, None)]:
r = laguerre(Q4, x0)
assert r.converged
z = r.root.real
assert abs(r.root.imag) == 0
if want is not None:
assert abs(z - want) < 1e-12
else:
assert abs(z - 3) < 1e-12 or abs(z - 4) < 1e-12
def test_complex_root_from_real_start():
r = laguerre([1, 0, 1], 0)
assert r.root == 1j and r.reason == "root" and r.evaluations == 2
def test_zero_denominator_fallback():
# at 0, p' = p'' = 0 for x^3 - 1
r = laguerre([1, 0, 0, -1], 0)
assert r.converged
assert abs(r.root ** 3 - 1) < 1e-14
def test_root_of_power_is_exact_in_one_step():
for n in range(1, 8):
c = [math.comb(n, j) * (-1) ** j for j in range(n + 1)] # (x - 1)^n
r = laguerre(c, -3.0)
assert r.root == 1 and r.reason == "root" and r.iterations == 2
def test_multiple_root_converges_linearly():
# (x - 1)^2 (x + 2): double root at 1, linear convergence
c = [1, 0, -3, 2]
r = laguerre(c, 3.0, max_iter=200)
assert r.converged and abs(r.root - 1) < 1e-7
def test_root_at_start_and_kmax_and_bad_input():
r = laguerre(Q4, 3)
assert r.reason == "root" and r.evaluations == 1
r = laguerre(Q4, 100.0, max_iter=1)
assert r.reason == "kmax reached" and r.evaluations == 1
with pytest.raises(ValueError):
laguerre([0, 1, 2], 0)
with pytest.raises(ValueError):
laguerre([5], 0)
def test_nan_and_overflow():
r = laguerre([1, float("nan"), 1], 0.5)
assert r.reason == "failure"
r = laguerre([1, 0, 1], 1e200)
assert r.reason in ("failure", "converged", "residual", "root")
assert r.evaluations <= 100
def test_horner2_and_deflate():
p, d1, d2, s = horner2(Q4, 5)
assert (p, d1, d2, s) == (24, 50, 70, 1 * 625 + 10 * 125 + 35 * 25 + 50 * 5 + 24)
assert deflate(Q4, 1) == [1, -9, 26, -24]
def test_all_roots_wilkinson_small_and_cubic():
rs = laguerre_all_roots([1, 0, 0, -1])
want = [1, complex(-0.5, math.sqrt(3) / 2), complex(-0.5, -math.sqrt(3) / 2)]
for w in want:
assert min(abs(z - w) for z in rs) < 1e-15
rs = laguerre_all_roots(poly_from_roots(range(1, 11)))
assert sorted(round(z.real) for z in rs) == list(range(1, 11))
assert max(min(abs(z - j) for z in rs) for j in range(1, 11)) < 1e-7
def test_random_polynomials():
rng = random.Random(20261108)
for _ in range(500):
n = rng.randint(1, 10)
zs = [complex(rng.uniform(-1, 1), rng.uniform(-1, 1)) for _ in range(n)]
c = poly_from_roots(zs)
rs = laguerre_all_roots(c)
assert len(rs) == n
rem = list(rs)
for z in zs:
j = min(range(len(rem)), key=lambda i: abs(rem[i] - z))
assert abs(rem[j] - z) < 1e-8
rem.pop(j)
for z in rs:
p, _, _, s = horner2(c, z)
assert abs(p) <= 1e-12 * s
What the next level needs
- the tests did not run in this build