"""Quadratic formula: exact check of the derivation, and the cancellation table.  Standard library only.
Each derivation step is an identity between left and right sides; we verify with exact fractions that
every step keeps the same solution set, and that the final formula solves a x^2 + b x + c = 0."""
from fractions import Fraction as F
from decimal import Decimal, getcontext
import math

def steps(a, b, c, x):
    """(left, right) of each frame 0..7 evaluated at x; the difference L-R must be (a x^2+b x+c)/a^k for fixed k."""
    a, b, c, x = map(F, (a, b, c, x))
    p = b / (2 * a)
    return [
        (a*x*x + b*x + c, 0),
        ((a*x*x + b*x + c) / a, 0),
        (x*x + b/a*x + c/a, 0),
        (x*x + b/a*x, -c/a),
        (x*x + b/a*x + p*p, -c/a + p*p),
        ((x + p)**2, -c/a + p*p),
        ((x + p)**2, -4*a*c/(4*a*a) + b*b/(4*a*a)),
        ((x + p)**2, (b*b - 4*a*c) / (4*a*a)),
    ]

def check(a, b, c):
    a, b, c = map(F, (a, b, c))
    for x in (F(0), F(1), F(-3, 7), F(5, 2)):
        f = a*x*x + b*x + c
        for i, (L, R) in enumerate(steps(a, b, c, x)):
            want = f if i == 0 else f / a
            assert L - R == want, (i, x)
    return True

def exact_roots(a, b, c):
    """Roots as Decimals (60 digits) through exact discriminant; None for complex."""
    getcontext().prec = 60
    a, b, c = map(F, (a, b, c)); D = b*b - 4*a*c
    if D < 0: return None
    s = Decimal(D.numerator).sqrt() / Decimal(D.denominator).sqrt()
    A, B = Decimal(a.numerator) / Decimal(a.denominator), Decimal(b.numerator) / Decimal(b.denominator)
    return (-B - s) / (2*A), (-B + s) / (2*A)

def naive(a, b, c):
    s = math.sqrt(b*b - 4*a*c); return (-b - s) / (2*a), (-b + s) / (2*a)

def stable(a, b, c):
    s = math.sqrt(b*b - 4*a*c); q = -0.5 * (b + math.copysign(s, b)); return q / a, c / q

if __name__ == "__main__":
    for abc in [(1, -5, 6), (2, 3, -7), (F(1, 3), F(-2, 5), F(1, 7))]:
        print(abc, "derivation checks:", check(*abc))
    print("\nroots of x^2 + 10^k x + 1: relative error of the small root")
    print("   k   naive        stable")
    for k in range(2, 15, 2):
        b = 10.0**k; ex = exact_roots(1, int(b), 1)[1]
        n = min(naive(1, b, 1), key=abs); s = min(stable(1, b, 1), key=abs)
        err = lambda v: abs((Decimal(v) - ex) / ex)
        print(f"{k:4d}   {float(err(n)):.2e}   {float(err(s)):.2e}")