Ferrari's method in Python

Closed form for quartic roots by reduction to a resolvent cubic. Shift away the cube term, solve a cubic for m, and split the quartic into two quadratics with the root m of largest modulus.

File: ferrari_method.py, 74 lines
Level: I0 not checked yet
Toolchain: Python 3.14, pytest 9
Test framework: pytest
Run the tests: python3 -m pytest

ferrari_method.pyRaw GitHub

# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/ferrari-method/. The authors are listed in prov/82X66E.jsonl.

"""Ferrari's method for the four roots of a quartic, in complex arithmetic."""

import cmath
from typing import Tuple


def _cardano(a: complex, b: complex, c: complex, d: complex) -> Tuple[complex, complex, complex]:
    """The three roots of a x^3 + b x^2 + c x + d by Cardano's formula (a != 0)."""
    omega = complex(-0.5, 3 ** 0.5 / 2)
    s, cc, dd = b / (3 * a), c / a, d / a
    p, q = cc - 3 * s * s, dd - s * cc + 2 * s ** 3
    rt = cmath.sqrt(q * q / 4 + p ** 3 / 27)
    w = -q / 2 + rt
    if abs(-q / 2 - rt) > abs(w):
        w = -q / 2 - rt
    if w == 0:
        return (-s, -s, -s)
    u = complex(w) ** (1 / 3)
    out = []
    for k in range(3):
        uk = u * omega ** k
        out.append(uk - p / (3 * uk) - s)
    return (out[0], out[1], out[2])


def _quad(beta: complex, gamma: complex) -> Tuple[complex, complex]:
    """The two roots of y^2 + beta y + gamma by the stable quadratic formula."""
    w = cmath.sqrt(beta * beta - 4 * gamma)  # @line L10
    t1, t2 = -(beta + w) / 2, -(beta - w) / 2
    t = t1 if abs(t1) >= abs(t2) else t2
    if t == 0:  # @line L11
        return (0j, 0j)
    return (t, gamma / t)


def ferrari(a: complex, b: complex, c: complex, d: complex, e: complex) -> Tuple[complex, complex, complex, complex]:
    """Return the four roots, with multiplicity, of a x^4 + b x^3 + c x^2 + d x + e.

    Precondition: a != 0. The coefficients may be real or complex. The root
    of the resolvent cubic of largest modulus is used. Cost: a fixed number of
    complex operations, five square roots and one cube root, O(1) space.
    """
    if a == 0:
        raise ValueError("a must be nonzero")
    s, C, D, E = b / (4 * a), c / a, d / a, e / a  # @line L1
    p = C - 6 * s * s  # @line L2
    q = D - 2 * C * s + 8 * s ** 3
    r = E - D * s + C * s * s - 3 * s ** 4
    ms = _cardano(8, 8 * p, 2 * p * p - 8 * r, -q * q)  # @line L3
    m = max(ms, key=abs)  # @line L4
    if m == 0:  # @line L5
        return (complex(-s), complex(-s), complex(-s), complex(-s))
    sigma = cmath.sqrt(2 * m)  # @line L6
    A = p / 2 + m + q / (2 * sigma)  # @line L7
    B = p / 2 + m - q / (2 * sigma)
    y1, y2 = _quad(-sigma, A)  # @line L8
    y3, y4 = _quad(sigma, B)
    return (y1 - s, y2 - s, y3 - s, y4 - s)  # @line L9


def polish(a: complex, b: complex, c: complex, d: complex, e: complex, x: complex, steps: int = 2) -> complex:
    """A few Newton steps on the quartic from x, by Horner's rule, to remove rounding errors."""
    for _ in range(steps):
        f = df = 0j
        for coef in (a, b, c, d, e):
            df = df * x + f
            f = f * x + coef
        if df == 0:
            break
        x = x - f / df
    return x

test_ferrari_method.py 109 lines

# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/ferrari-method/. The authors are listed in prov/82X66E.jsonl.

import random
import unittest
from fractions import Fraction

from ferrari_method import ferrari, polish


def coeffs(roots):
    """Coefficients of prod (x - r), leading first."""
    c = [1]
    for x in roots:
        c = [(c[i] if i < len(c) else 0) - x * (c[i - 1] if i >= 1 else 0) for i in range(len(c) + 1)]
    return c


def match_error(got, roots):
    """Largest distance from a true root to its matched computed root (greedy matching)."""
    left = list(got)
    worst = 0.0
    for x in roots:
        j = min(range(len(left)), key=lambda i: abs(left[i] - x))
        worst = max(worst, abs(left[j] - x))
        left.pop(j)
    return worst


class TestFerrari(unittest.TestCase):
    def test_ferrari_equation(self):
        got = ferrari(1, 0, 6, -60, 36)
        reals = sorted(z.real for z in got if abs(z.imag) < 1e-12)
        self.assertAlmostEqual(reals[0], 0.6443988642268, places=12)
        self.assertAlmostEqual(reals[1], 3.0998744240188, places=12)
        for z in got:
            self.assertLess(abs(z ** 4 + 6 * z * z - 60 * z + 36), 1e-11)

    def test_integer_roots(self):
        got = ferrari(1, -10, 35, -50, 24)
        self.assertLess(match_error(got, [1, 2, 3, 4]), 1e-12)

    def test_biquadratic(self):
        got = ferrari(1, 0, -5, 0, 4)
        self.assertLess(match_error(got, [1, -1, 2, -2]), 1e-12)

    def test_complex_roots(self):
        got = ferrari(1, 0, 0, 0, 4)  # roots +-1 +-i
        self.assertLess(match_error(got, [1 + 1j, 1 - 1j, -1 + 1j, -1 - 1j]), 1e-12)

    def test_quadruple_root(self):
        got = ferrari(*coeffs([2, 2, 2, 2]))
        self.assertLess(match_error(got, [2, 2, 2, 2]), 1e-12)
        got = ferrari(1, 0, 0, 0, 0)
        self.assertEqual(got, (0j, 0j, 0j, 0j))

    def test_multiple_roots(self):
        self.assertLess(match_error(ferrari(*coeffs([1, 1, 2, 2])), [1, 1, 2, 2]), 1e-8)
        self.assertLess(match_error(ferrari(*coeffs([2, 2, 2, -3])), [2, 2, 2, -3]), 1e-4)

    def test_scaling_and_complex_coefficients(self):
        roots = [1j, 2 - 1j, -0.5, 3 + 2j]
        c = [3 * x for x in coeffs(roots)]
        self.assertLess(match_error(ferrari(*c), roots), 1e-11)

    def test_largest_modulus_choice(self):
        roots = [1, 2, 3, 4 + 1e-7]
        self.assertLess(match_error(ferrari(*coeffs(roots)), roots), 1e-12)

    def test_polish_small_root(self):
        roots = [1e-4, 1, 2, 1000]
        c = coeffs(roots)
        got = ferrari(*c)
        raw = min(abs(z - 1e-4) for z in got) / 1e-4
        self.assertGreater(raw, 1e-8)  # the shift loses digits
        fixed = [polish(*c, z) for z in got]
        for x in roots:
            self.assertLess(min(abs(z - x) for z in fixed) / x, 1e-14)

    def test_exact_reference(self):
        # roots chosen as fractions; compare against exact values
        roots = [Fraction(1, 3), Fraction(-5, 7), Fraction(9, 4), Fraction(-2, 1)]
        c = [float(v) for v in coeffs(roots)]
        got = ferrari(*c)
        self.assertLess(match_error(got, [float(x) for x in roots]), 1e-12)

    def test_zero_leading(self):
        with self.assertRaises(ValueError):
            ferrari(0, 1, 2, 3, 4)

    def test_random(self):
        rng = random.Random(20261107)
        for _ in range(600):
            roots = []
            while len(roots) < 4:
                if rng.random() < 0.5 and len(roots) <= 2:
                    z = complex(rng.uniform(-3, 3), rng.uniform(0.1, 3))
                    roots += [z, z.conjugate()]
                else:
                    roots.append(complex(rng.uniform(-3, 3), 0))
            gaps = min(abs(roots[i] - roots[j]) for i in range(4) for j in range(i))
            got = ferrari(*coeffs(roots))
            pol = [polish(*coeffs(roots), z) for z in got]
            err = match_error(pol, roots)
            self.assertLess(err, 1e-9 / min(1.0, gaps), msg=str(roots))


if __name__ == "__main__":
    unittest.main()

What the next level needs

the tests did not run in this build

Read the full page about Ferrari's method