Graeffe's method in Python
Repeatedly squares the roots so that they separate in magnitude and can be read off the coefficients. The moduli converge quadratically once the squared roots separate.
- File
graeffe_method.py, 111 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/graeffe-method/. The authors are listed in prov/4YAGDA.jsonl.
"""Graeffe's root squaring method: the moduli of all roots of a real polynomial.
Coefficients are real and listed highest degree first: a[0] x^n + ... + a[n].
One squaring step maps p to the polynomial q with q(x^2) = (-1)^n p(x) p(-x),
whose roots are the squares of the roots of p. After k steps the roots are
raised to the power N = 2^k, the moduli separate, and the modulus of a root
can be read off two coefficients.
The coefficients grow like r^N, far beyond the range of a float, so they are
kept as Decimal numbers with a huge exponent range. A coefficient index j is
"regular" when the step squares it: (-1)^j c_j / b_j^2 is within tol of 1. Between
two consecutive regular indices i < l lie l - i roots of one common modulus
|b_l / b_i|^(1 / (N (l - i))).
"""
import math
from collections import namedtuple
from decimal import Context, Decimal, MAX_EMAX, MIN_EMIN
GraeffeResult = namedtuple("GraeffeResult", "groups roots steps status")
# groups: list of (modulus, count), largest modulus first, zeros last.
# roots: the real roots found for groups of count 1 (a lone modulus must
# belong to a real root when the coefficients are real).
# status: "separated" when every index became regular, else "max-iter"; then
# the groups merge roots whose moduli were not yet separated.
CTX = Context(prec=40, Emax=MAX_EMAX, Emin=MIN_EMIN)
def square_step(b):
"""One Graeffe step: c_j = (-1)^j (b_j^2 + 2 sum_i (-1)^i b_(j-i) b_(j+i))."""
n = len(b) - 1
c = []
for j in range(n + 1):
s = CTX.multiply(b[j], b[j])
for i in range(1, min(j, n - j) + 1):
t = CTX.multiply(CTX.multiply(b[j - i], b[j + i]), Decimal(2))
s = CTX.add(s, t) if i % 2 == 0 else CTX.subtract(s, t)
c.append(s if j % 2 == 0 else CTX.minus(s))
return c
def _regular(b, c):
"""The deviations |(-1)^j c_j / b_j^2 - 1|, infinite where b_j or c_j is 0."""
dev = []
for j in range(len(b)):
if b[j] == 0 or c[j] == 0:
dev.append(math.inf) # cancelled coefficient: irregular
continue
r = float(CTX.divide(c[j], CTX.multiply(b[j], b[j])))
dev.append(abs((r if j % 2 == 0 else -r) - 1.0))
return dev
def _groups(b, reg, N):
"""Moduli of the groups between consecutive regular indices."""
out = []
for i, l in zip(reg, reg[1:]):
lg = CTX.subtract(CTX.ln(CTX.abs(b[l])), CTX.ln(CTX.abs(b[i])))
out.append((math.exp(float(lg) / (N * (l - i))), l - i))
return out
def _horner(a, x):
s = 0.0
for c in a:
s = s * x + c
return s
def graeffe(a, tol=1e-12, max_iter=30):
a = [float(c) for c in a]
if any(math.isnan(c) or math.isinf(c) for c in a):
raise ValueError("coefficients must be finite")
while a and a[0] == 0.0: # @line L1
a.pop(0)
if not a:
raise ValueError("the zero polynomial has no well defined roots")
zeros = 0
while len(a) > 1 and a[-1] == 0.0: # @line L1
a.pop()
zeros += 1
n = len(a) - 1
if n == 0:
groups = [(0.0, zeros)] if zeros else []
return GraeffeResult(groups, [0.0] * zeros, 0, "separated")
b = [Decimal(c) for c in a] # @line L2
k, status, dev = 0, "max-iter", [0.0] * (n + 1)
while k < max_iter: # @line L3
c = square_step(b) # @line L4
k += 1 # @line L5
dev = _regular(b, c) # @line L6
b = c # @line L7
if max(dev) <= tol: # @line L8
status = "separated"
break
reg = [j for j in range(n + 1) if dev[j] <= tol or j in (0, n)] # @line L9
N = 2 ** k
groups = _groups(b, reg, N) # @line L10
groups.sort(key=lambda g: -g[0])
roots = []
for rho, cnt in groups: # @line L11
if cnt == 1:
roots.append(rho if abs(_horner(a, rho)) <= abs(_horner(a, -rho)) else -rho)
if zeros:
groups.append((0.0, zeros))
roots.extend([0.0] * zeros)
return GraeffeResult(groups, roots, k, status) # @line L12
test_graeffe_method.py 129 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/graeffe-method/. The authors are listed in prov/4YAGDA.jsonl.
import math
import random
import unittest
from decimal import Decimal
from graeffe_method import graeffe, square_step
def poly(roots):
c = [1.0]
for r in roots:
c = [x - r * y for x, y in zip(c + [0.0], [0.0] + c)]
return c
class TestGraeffe(unittest.TestCase):
def test_square_step_squares_roots(self):
b = square_step([Decimal(x) for x in poly([1, 2, 3])])
self.assertEqual(b, [Decimal(1), Decimal(-14), Decimal(49), Decimal(-36)])
b = square_step(b)
self.assertEqual(b, [Decimal(1), Decimal(-98), Decimal(1393), Decimal(-1296)])
def test_example_converges(self):
r = graeffe(poly([1, 2, 3]))
self.assertEqual(r.status, "separated")
self.assertEqual(r.steps, 8)
for x, y in zip(r.roots, [3.0, 2.0, 1.0]):
self.assertAlmostEqual(x, y, delta=1e-14)
def test_signs_of_real_roots(self):
r = graeffe(poly([1, -2, 3, -4.5]))
for x, y in zip(sorted(r.roots), [-4.5, -2.0, 1.0, 3.0]):
self.assertAlmostEqual(x, y, delta=1e-14)
def test_quadratic_convergence_of_deviation(self):
# the deviation (-1)^j c_j / b_j^2 - 1 roughly squares at every step
b = [Decimal(x) for x in poly([1, 2, 3])]
devs = []
for _ in range(6):
c = square_step(b)
devs.append(abs(float(c[2] / (b[2] * b[2])) - 1)) # j = 2 is even
b = c
ratios = [math.log(devs[k + 1]) / math.log(devs[k]) for k in range(3, 5)]
for q in ratios:
self.assertGreater(q, 1.8)
self.assertLess(q, 3.0)
def test_complex_pair_gives_group(self):
r = graeffe([1, 0, 3, 6, 10]) # roots -1 +- i and 1 +- 2i
self.assertEqual(r.status, "max-iter")
self.assertEqual(len(r.groups), 2)
self.assertAlmostEqual(r.groups[0][0], math.sqrt(5), delta=1e-14)
self.assertAlmostEqual(r.groups[1][0], math.sqrt(2), delta=1e-14)
self.assertEqual([g[1] for g in r.groups], [2, 2])
self.assertEqual(r.roots, [])
def test_no_real_root(self):
r = graeffe([1, 0, 1])
self.assertEqual(r.groups, [(1.0, 2)])
self.assertEqual(r.roots, [])
def test_opposite_roots_and_multiple_roots(self):
self.assertEqual(graeffe([1, 0, -4]).groups, [(2.0, 2)])
for m in (2, 3, 4):
r = graeffe(poly([1.5] * m + [0.5]))
self.assertEqual([g[1] for g in r.groups], [m, 1])
self.assertAlmostEqual(r.groups[0][0], 1.5, delta=1e-14)
self.assertAlmostEqual(r.roots[0], 0.5, delta=1e-15)
def test_root_exactly_zero(self):
r = graeffe(poly([1, -2, 3, 0, 0]))
self.assertEqual(r.groups[-1], (0.0, 2))
for x, y in zip(sorted(r.roots), [-2.0, 0.0, 0.0, 1.0, 3.0]):
self.assertAlmostEqual(x, y, delta=1e-14)
def test_close_moduli(self):
r = graeffe(poly([1.0, 1.01]))
self.assertEqual(r.status, "separated")
self.assertAlmostEqual(r.roots[0], 1.01, delta=1e-12)
self.assertAlmostEqual(r.roots[1], 1.0, delta=1e-12)
def test_max_iter(self):
r = graeffe(poly([1, 2, 3]), max_iter=3)
self.assertEqual(r.status, "max-iter")
self.assertEqual(r.steps, 3)
self.assertEqual(r.groups[0][1], 3) # not yet separated: one group
def test_wilkinson_20(self):
r = graeffe(poly(range(1, 21)))
self.assertEqual(r.status, "separated")
self.assertEqual(r.steps, 11)
err = max(abs(x - k) / k for x, k in zip(sorted(r.roots), range(1, 21)))
self.assertLess(err, 1e-4) # the rounded coefficients move the roots this much
def test_degenerate_inputs(self):
self.assertEqual(graeffe([0, 0, 2, -4]).roots, [2.0])
self.assertEqual(graeffe([5]).groups, [])
with self.assertRaises(ValueError):
graeffe([0, 0])
with self.assertRaises(ValueError):
graeffe([1, float("nan"), 1])
with self.assertRaises(ValueError):
graeffe([1, float("inf")])
def test_random_distinct_moduli(self):
rng = random.Random(20240601)
worst, steps = 0.0, 0
for _ in range(500):
n = rng.randint(1, 10)
while True:
mods = sorted(rng.uniform(0.1, 10) for _ in range(n))
if all(mods[i + 1] / mods[i] > 1.05 for i in range(n - 1)):
break
roots = [m * rng.choice([-1, 1]) for m in mods]
r = graeffe(poly(roots))
self.assertEqual(r.status, "separated")
self.assertEqual(len(r.roots), n)
steps += r.steps
for x, y in zip(sorted(r.roots), sorted(roots)):
worst = max(worst, abs(x - y) / abs(y))
self.assertLess(worst, 1e-9)
self.assertLess(steps / 500, 12)
if __name__ == "__main__":
unittest.main()
What the next level needs
- the tests did not run in this build