The Quadratic Formula Is Algebraically Right—and Numerically Fragile

Most of us learn the quadratic formula as a finished piece of algebra:

x = (-b ± sqrt(b² - 4ac)) / (2a)

For exact arithmetic, that formula is complete. For floating-point arithmetic, it is only the beginning.

The problem is not that JavaScript implements the formula incorrectly. The problem is that the two algebraically equivalent branches can have very different numerical behavior. When b and sqrt(b² - 4ac) are nearly equal, one numerator subtracts two large, nearly equal numbers. Most of their significant bits cancel, and the small remainder inherits the rounding error of the large operands.

That failure mode is called catastrophic cancellation. It can turn a well-conditioned large root and a representable small root into one plausible answer and one surprisingly inaccurate answer.

This article develops a stable real-root branch, explains the limits of the technique, and builds a validation strategy around residuals, Vieta's identities, root classification, and scale-sensitive cases.

A small equation that exposes a large error

Consider:

x² + 100000000x + 1 = 0

Here a = 1, b = 100000000, and c = 1. Vieta's identities tell us in advance that the roots have:

x1 + x2 = -b/a = -100000000
x1 * x2 =  c/a = 1

So one root is close to -100000000, while the other is close to -0.00000001.

A direct implementation looks harmless:

function naiveRealRoots(a, b, c) {
  const discriminant = b * b - 4 * a * c;
  const squareRoot = Math.sqrt(discriminant);

  return [
    (-b - squareRoot) / (2 * a),
    (-b + squareRoot) / (2 * a),
  ];
}

In binary64 arithmetic, that code produces approximately:

-100000000
-7.450580596923828e-9

The second value should be close to -1e-8. Its relative error is about 25%. The damage happens in this numerator:

-100000000 + sqrt(9999999999999996)

The two terms agree in almost all of their leading digits. Subtracting them discards those shared digits and leaves a tiny result whose useful precision has already been lost.

This is an important distinction: the subtraction did not make the absolute rounding error enormous. It made the error enormous relative to the small result we wanted.

Choose the non-cancelling numerator, then recover the other root

For two distinct real roots, define:

q = -1/2 * (b + sign(b) * sqrt(D))

where:

D = b² - 4ac

The sign choice makes the terms inside the parentheses reinforce each other instead of cancel. One root is then:

x1 = q/a

Vieta's product identity gives the other:

x1 * x2 = c/a
x2 = c/q

The second formula is not an approximation layered on top of the quadratic formula. It is an algebraically equivalent recovery step. Its value is numerical: it avoids evaluating the dangerous numerator.

In JavaScript, the real-root branch is compact:

function stableDistinctRealRoots(a, b, c, discriminant) {
  const squareRoot = Math.sqrt(discriminant);
  const signedSquareRoot = b >= 0 ? squareRoot : -squareRoot;
  const q = -0.5 * (b + signedSquareRoot);

  const first = q / a;
  const second = c / q;

  return first <= second ? [first, second] : [second, first];
}

For the cancellation example, this returns:

-100000000
-1e-8

The product of those represented values is exactly 1 in this run, matching c/a. The sum differs from -b/a only at the final representable digits.

A complete solver still needs three branches

The q method is for D > 0. A complete quadratic solver must distinguish three cases:

1. D > 0: two distinct real roots; use the stable q branch.
2. D === 0: one repeated real root, -b / (2a).
3. D < 0: a complex-conjugate pair.

Here is a self-contained implementation. It deliberately rejects non-finite coefficients and a = 0 rather than silently treating a linear equation as quadratic.

function assertFiniteNumber(value, name) {
  if (!Number.isFinite(value)) {
    throw new TypeError(`${name} must be a finite number`);
  }
}

function solveQuadratic(a, b, c) {
  assertFiniteNumber(a, "a");
  assertFiniteNumber(b, "b");
  assertFiniteNumber(c, "c");

  if (a === 0) {
    throw new RangeError("a must be non-zero for a quadratic equation");
  }

  const discriminant = b * b - 4 * a * c;
  if (!Number.isFinite(discriminant)) {
    throw new RangeError("the discriminant is outside the supported range");
  }

  if (discriminant > 0) {
    const squareRoot = Math.sqrt(discriminant);
    const q = -0.5 * (b + (b >= 0 ? squareRoot : -squareRoot));
    const roots = [q / a, c / q].sort((left, right) => left - right);

    if (!roots.every(Number.isFinite)) {
      throw new RangeError("a real root is outside the supported range");
    }

    return {
      discriminant,
      type: "two_real",
      roots: roots.map((real) => ({ real, imaginary: 0 })),
    };
  }

  if (discriminant === 0) {
    const root = -b / (2 * a);
    if (!Number.isFinite(root)) {
      throw new RangeError("the repeated root is outside the supported range");
    }

    return {
      discriminant,
      type: "repeated_real",
      roots: [
        { real: root, imaginary: 0 },
        { real: root, imaginary: 0 },
      ],
    };
  }

  const real = -b / (2 * a);
  const imaginary = Math.sqrt(-discriminant) / (2 * Math.abs(a));

  if (![real, imaginary].every(Number.isFinite)) {
    throw new RangeError("a complex root is outside the supported range");
  }

  return {
    discriminant,
    type: "complex_conjugate",
    roots: [
      { real, imaginary },
      { real, imaginary: -imaginary },
    ],
  };
}

Using Math.abs(a) for the displayed imaginary magnitude is safe because the result is the unordered pair real ± imaginary * i. Changing the sign of a swaps which conjugate is associated with the plus branch; it does not change the pair.

Test residuals, but scale them

Substituting a reported root into the polynomial is a useful test:

p(x) = ax² + bx + c

An absolute residual alone can be misleading. In the cancellation example, evaluating p(-100000000) gives an absolute residual of 1. That sounds bad until we notice that the two large terms being added have magnitude around 10^16.

A scale-aware residual is more informative:

function scaledResidual(a, b, c, x) {
  const residual = Math.abs(a * x * x + b * x + c);
  const scale = Math.abs(a * x * x) + Math.abs(b * x) + Math.abs(c);
  return scale === 0 ? residual : residual / scale;
}

For a = 1, b = 1e8, and c = 1, the observed values are:

Root calculation	Reported small root	Scaled residual
Direct (-b + sqrt(D)) / (2a)	-7.450580596923828e-9	about 1.46e-1
Stable c / q	-1e-8	about 5.55e-17

The stable large root has a scaled residual of about 5e-17 as well. The absolute residual did not disappear, but relative to the arithmetic scale it is near machine precision.

Residuals should be paired with structural checks:

function approximatelyEqual(actual, expected, absolute = 1e-12, relative = 1e-12) {
  return Math.abs(actual - expected) <=
    Math.max(absolute, relative * Math.max(Math.abs(actual), Math.abs(expected)));
}

function verifyRealPair(a, b, c, roots) {
  const [x1, x2] = roots;
  return {
    sumMatches: approximatelyEqual(x1 + x2, -b / a),
    productMatches: approximatelyEqual(x1 * x2, c / a),
    residuals: roots.map((x) => scaledResidual(a, b, c, x)),
  };
}

Vieta's sum and product can catch a broken pair even when one root happens to have a small residual. Conversely, residuals catch cases where a sum or product comparison looks acceptable only because the scale is huge.

A compact regression matrix

A solver should not be validated only on convenient factorable quadratics. The following matrix covers distinct real, repeated real, complex, cancellation-prone, scaled, and invalid inputs.

import assert from "node:assert/strict";

{
  const result = solveQuadratic(1, -3, 2);
  assert.equal(result.type, "two_real");
  assert.deepEqual(result.roots, [
    { real: 1, imaginary: 0 },
    { real: 2, imaginary: 0 },
  ]);
}

{
  const result = solveQuadratic(1, 2, 1);
  assert.equal(result.type, "repeated_real");
  assert.equal(result.discriminant, 0);
  assert.equal(result.roots[0].real, -1);
  assert.equal(result.roots[1].real, -1);
}

{
  const result = solveQuadratic(1, 2, 5);
  assert.equal(result.type, "complex_conjugate");
  assert.deepEqual(result.roots, [
    { real: -1, imaginary: 2 },
    { real: -1, imaginary: -2 },
  ]);
}

{
  const result = solveQuadratic(1, 1e8, 1);
  const realRoots = result.roots.map((root) => root.real);
  assert.ok(approximatelyEqual(realRoots[0], -1e8, 1e-7, 1e-15));
  assert.ok(approximatelyEqual(realRoots[1], -1e-8, 1e-20, 1e-12));
  assert.ok(scaledResidual(1, 1e8, 1, realRoots[1]) < 1e-15);
}

for (const [a, b, c, expectedRoot] of [
  [1e50, -2, 1e-50, 1e-50],
  [1e100, -2, 1e-100, 1e-100],
  [1e-100, -2e-100, 1e-100, 1],
]) {
  const result = solveQuadratic(a, b, c);
  assert.equal(result.discriminant, 0);
  assert.equal(result.type, "repeated_real");
  assert.ok(approximatelyEqual(result.roots[0].real, expectedRoot, 0, 1e-12));
}

{
  const positive = solveQuadratic(1e50, -2, (1 - 1e-12) * 1e-50);
  const negative = solveQuadratic(1e50, -2, (1 + 1e-12) * 1e-50);
  assert.equal(positive.type, "two_real");
  assert.ok(positive.discriminant > 0);
  assert.equal(negative.type, "complex_conjugate");
  assert.ok(negative.discriminant < 0);
}

assert.throws(() => solveQuadratic(0, 2, 1), /non-zero/);

The scaled repeated-root cases matter because a test suite made only of coefficients near 1 can hide accidental normalization bugs. In a bounded binary64 implementation, evaluating the discriminant in the original coefficient units can preserve an exact zero that an unnecessary rescaling step turns into a tiny signed value.

The positive and negative perturbations test the other side of that decision. A blanket epsilon rule such as "treat every sufficiently small discriminant as zero" can silently convert two real roots or a complex pair into a repeated root. There is no universal epsilon independent of coefficient scale and input uncertainty.

Parsing is part of numerical correctness

The solver can be stable while the input layer is not. A practical coefficient form should make its accepted grammar explicit.

The implementation that prompted this audit accepts:

• finite decimal values such as -3, .125, and 2.;
• scientific notation such as 1.2e2;
• one simple fraction such as 3/4 or +1.2e2/-3.

It rejects blank input, NaN, Infinity, division by zero, mixed numbers, chained fractions, and algebraic expressions in coefficient fields. It also rejects a = 0 with a message explaining that the equation is then linear.

Keeping parsing separate from root calculation has two benefits. First, the numerical core receives a small, documented input type. Second, the user interface does not need to duplicate sqrt, discriminant, or root logic. The display layer can render the core's discriminant, rootType, root1, and root2 values instead of becoming a second solver that might drift from the tested one.

What the stable branch does not solve

The q method is a targeted improvement, not arbitrary-precision arithmetic.

• If the input coefficients have already been rounded, a solver cannot recover information that was never represented.
• Coefficients outside a declared range can make b*b or 4*a*c overflow, or make a root overflow even when the discriminant is finite.
• A polynomial near a repeated root is intrinsically ill-conditioned. Tiny coefficient changes can legitimately change the root classification.
• Floating-point evaluation of the polynomial can itself cancel, which is why a scaled residual and Vieta checks are more useful together than either is alone.
• Sorting real roots is a presentation choice. It should happen after calculation and should not be confused with the algebraic plus and minus branch labels.

In the audited contract, each nonzero coefficient is bounded in magnitude between 1e-100 and 1e100, while zero remains valid for b and c. Those limits keep b² and 4ac finite and within the normal binary64 range. They are part of the algorithm's validity, not merely form-field decoration.

The engineering lesson

Algebra tells us which expressions are equivalent over exact numbers. Numerical analysis asks whether those expressions are equally useful on a finite machine.

For a quadratic with two real roots, the implementation pattern is:

1. validate finite, bounded coefficients and require a != 0;
2. compute and classify the discriminant;
3. choose the numerator that adds like-signed magnitudes;
4. compute one root as q/a and recover the other as c/q;
5. validate with scale-aware residuals and Vieta's identities;
6. test repeated, complex, cancellation-prone, perturbed, and cross-scale cases.

The implementation discussed here can also be exercised through this quadratic formula calculator, but the code and checks above are deliberately complete enough to reproduce without depending on a hosted interface.

The larger lesson applies far beyond quadratics: when a formula subtracts nearly equal quantities, algebraic elegance is not evidence of numerical stability. Look for an equivalent form that keeps the significant bits you actually need—and make the failure case a permanent regression test.

Comments

[Luis Cruz]

It's fascinating how you've addressed catastrophic cancellation in floating-point arithmetic—an issue that often goes unnoticed until it leads to significant errors. Your approach to use Vieta’s identities not only provides stability but also highlights the importance of choosing numerically favorable algorithms. It could be beneficial to extend this discussion to optimization techniques for cases involving larger coefficients or tighter tolerances. If you're looking for extra hands on further refining this solver or exploring other numerical methods, I’d be happy to discuss potential collaboration.

[Ekong Ikpe]

Don't blame the formula; fix the data pipeline and adapt to the environment. if abstracted properly this works at scale for other problems including AI use.
you did a great job in this post ✌️