A Numerically Stable Quadratic Solver in Python (with Steps)

Most "quadratic formula" snippets online stop at two lines of code. That is fine until you hit a negative discriminant, a = 0, or floating-point noise that turns a double root into two slightly different numbers. Here is a small, careful solver in Python that handles all of those cases and prints each step, which is handy when you are building study tools or checking homework.

The math in one paragraph

For ax² + bx + c = 0 with a ≠ 0, the discriminant D = b² − 4ac decides everything. D > 0 gives two real roots, D = 0 gives one repeated root, and D < 0 gives a pair of complex conjugate roots. The roots are x = (−b ± √D) / 2a.

A step-by-step solver

import cmath
import math

def solve_quadratic(a: float, b: float, c: float, eps: float = 1e-12):
    steps = []
    if abs(a) < eps:
        if abs(b) < eps:
            return steps + ["a = 0 and b = 0: not an equation in x"], []
        x = -c / b
        return steps + [f"a = 0, so it is linear: x = -c/b = {x:g}"], [x]

    D = b * b - 4 * a * c
    steps.append(f"Discriminant D = b^2 - 4ac = {b:g}^2 - 4({a:g})({c:g}) = {D:g}")

    if abs(D) < eps:
        x = -b / (2 * a)
        steps.append(f"D = 0, one repeated root: x = -b/2a = {x:g}")
        return steps, [x]

    if D > 0:
        # Numerically stable form: avoid subtracting nearly equal numbers
        q = -0.5 * (b + math.copysign(math.sqrt(D), b))
        x1, x2 = q / a, c / q
        steps.append(f"D > 0, two real roots: x1 = {x1:g}, x2 = {x2:g}")
        return steps, sorted([x1, x2])

    r = cmath.sqrt(D)
    x1, x2 = (-b + r) / (2 * a), (-b - r) / (2 * a)
    steps.append(f"D < 0, complex roots: {x1:.4g} and {x2:.4g}")
    return steps, [x1, x2]

for coeffs in [(2, -7, 3), (1, -6, 9), (1, 2, 5), (0, 4, -8), (1, 1e8, 1)]:
    steps, roots = solve_quadratic(*coeffs)
    print(coeffs, "->", roots)
    for s in steps:
        print("   ", s)

Why the "stable" branch matters

Try a = 1, b = 10⁸, c = 1. The textbook formula computes (−b + √D) by subtracting two numbers that agree in almost every digit, and you lose most of your precision. The small root comes out as 0 or as garbage. Computing q = −½(b + sign(b)·√D) first and then using x₁ = q/a and x₂ = c/q (Vieta's formula, x₁·x₂ = c/a) keeps full precision for both roots.

Testing it

Check each root by plugging it back in: a·x² + b·x + c should be close to 0. For 2x² − 7x + 3 you should get x = 0.5 and x = 3. For x² − 6x + 9 you should get the single root 3. For x² + 2x + 5 you should get −1 ± 2i.

Checking your work

When you're studying by hand, it helps to compare each line of your working against a worked solution. I use the free step-by-step MathGPT AI math solver for that, alongside small scripts like this one.

Takeaways

• Handle a = 0 separately, because it's a linear equation.
• Use a tolerance for D = 0 instead of an exact comparison.
• Use the stable form for real roots when |b| is large.
• Use cmath for negative discriminants instead of raising an error.