# 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

```python
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](https://mathssolvergpt.com/en) 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.