<!-- canonical: https://efficientnewlanguage.org/eml-p/projects/P114-quadratic-solver/ | updated: 2026-10-16 -->

# P114 Quadratic solver

ax^2 + bx + c = 0 with exact rational, irrational and complex roots, square roots and divisions worked digit by digit on whole numbers; and the small root of x^2 + bx + 1 in floats, by the textbook formula and a stable one, against its exact digits.

EML-P project `projects/quadratic-solver` in the EML language repo: 2 module(s), entry `main.eml`, terminal UI. There, `eml project run projects/quadratic-solver` runs it and `eml project verify projects/quadratic-solver` replays every session under CPython (two hash seeds) and in the interpreter; the site build replays every session in the interpreter again and publishes a session only if its screen matches.

Built on verified corpus cases: quadratic-solver (https://efficientnewlanguage.org/cases/034-quadratic-solver/), catastrophic-cancellation (https://efficientnewlanguage.org/cases/234-catastrophic-cancellation/).

## Sessions

### bad-input - interpreter: byte-equal to the golden

Input:

```text
0
x
1

1
abc
1
1.2345
1
1

1
1
2

1
1
2
1,5
1
9999999999
2

2
2
2
1000000000001
2
x
3
```

Screen:

```text
== Quadratic solver ==
ax^2 + bx + c = 0, solved exactly - and the same algebra in floats losing its small root.

1) solve  2) floats and cancellation  3) quit
choice> 0
Pick a number from 1 to 3.

1) solve  2) floats and cancellation  3) quit
choice> x
Pick a number from 1 to 3.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 
Cancelled.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> abc
Type a number with at most 9 digits before the point and 3 after, such as -2 or 0.125.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1.2345
Type a number with at most 9 digits before the point and 3 after, such as -2 or 0.125.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> 
Cancelled.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> 2
c> 
Cancelled.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> 2
c> 1,5
Type a number with at most 9 digits before the point and 3 after, such as -2 or 0.125.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 9999999999
Type a number with at most 9 digits before the point and 3 after, such as -2 or 0.125.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> 
Cancelled.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> 2
Type a whole number from 3 to 1000000000000.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> 1000000000001
Type a whole number from 3 to 1000000000000.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> x
Type a whole number from 3 to 1000000000000.

1) solve  2) floats and cancellation  3) quit
choice> 3
Bye.
```

### basic - interpreter: byte-equal to the golden

Input:

```text
1
1
-3
2
1
2
-7
3
1
1
0
-2
1
1
2
5
1
1
-2
1
1
0
2
-4
1
0
0
0
1
0
0
5
1
-1
0
2
1
0.5
-0.125
0.007
2
3
2
100000000
2
1000000000000
3
```

Screen:

```text
== Quadratic solver ==
ax^2 + bx + c = 0, solved exactly - and the same algebra in floats losing its small root.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> -3
c> 2
x^2 - 3x + 2 = 0
Discriminant b^2 - 4ac = 1.
It is the square of a fraction: two rational roots, x = 1 and x = 2.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 2
b> -7
c> 3
2x^2 - 7x + 3 = 0
Discriminant b^2 - 4ac = 25.
It is the square of a fraction: two rational roots, x = 1/2 and x = 3.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> 0
c> -2
x^2 - 2 = 0
Discriminant b^2 - 4ac = 8.
It is positive and not the square of a fraction: two irrational roots, x = (-b - sqrt(D)) / 2a and (-b + sqrt(D)) / 2a:
  x = -1.414213562373
  x = 1.414213562373
to 12 places, each digit worked out exactly.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> 2
c> 5
x^2 + 2x + 5 = 0
Discriminant b^2 - 4ac = -16.
It is negative: two complex roots, x = (-b / 2a) - i (sqrt(-D) / 2a) and its conjugate:
  x = -1 - 2.000000000000i
  x = -1 + 2.000000000000i
the imaginary part to 12 places.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 1
b> -2
c> 1
x^2 - 2x + 1 = 0
Discriminant b^2 - 4ac = 0.
It is zero: one double root, x = 1.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 0
b> 2
c> -4
2x - 4 = 0
With a = 0 it is not quadratic: one root, x = 2.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 0
b> 0
c> 0
0 = 0
Every x solves it: all three coefficients are zero.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 0
b> 0
c> 5
5 = 0
No x solves it: it says a number other than zero equals zero.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> -1
b> 0
c> 2
-x^2 + 2 = 0
Discriminant b^2 - 4ac = 8.
It is positive and not the square of a fraction: two irrational roots, x = (-b - sqrt(D)) / 2a and (-b + sqrt(D)) / 2a:
  x = -1.414213562373
  x = 1.414213562373
to 12 places, each digit worked out exactly.

1) solve  2) floats and cancellation  3) quit
choice> 1
a> 0.5
b> -0.125
c> 0.007
0.5x^2 - 0.125x + 0.007 = 0
Discriminant b^2 - 4ac = 0.001625.
It is positive and not the square of a fraction: two irrational roots, x = (-b - sqrt(D)) / 2a and (-b + sqrt(D)) / 2a:
  x = 0.084688711259
  x = 0.165311288741
to 12 places, each digit worked out exactly.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> 3
x^2 + 3x + 1 = 0, its root nearest zero:
  exact, 20 digits:  -3.8196601125010515179e-1
  textbook formula:  -0.3819660112501051   (right to about 15 significant digits)
  stable formula:    -0.38196601125010515   (right to 17 significant digits or more)
The textbook formula takes -b + sqrt(b^2 - 4) when the two are nearly equal, and the digits they share cancel. The stable one works out the big root, where nothing cancels, and divides: the product of the roots is c / a = 1.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> 100000000
x^2 + 100000000x + 1 = 0, its root nearest zero:
  exact, 20 digits:  -1.0000000000000001000e-8
  textbook formula:  -1.4901161193847656e-08   (off by 49%)
  stable formula:    -1.0000000000000002e-08   (right to about 16 significant digits)
The textbook formula takes -b + sqrt(b^2 - 4) when the two are nearly equal, and the digits they share cancel. The stable one works out the big root, where nothing cancels, and divides: the product of the roots is c / a = 1.

1) solve  2) floats and cancellation  3) quit
choice> 2
b, a whole number from 3 to 1000000000000> 1000000000000
x^2 + 1000000000000x + 1 = 0, its root nearest zero:
  exact, 20 digits:  -1.0000000000000000000e-12
  textbook formula:  0.0   (off by 100%)
  stable formula:    -1e-12   (right to 17 significant digits or more)
The textbook formula takes -b + sqrt(b^2 - 4) when the two are nearly equal, and the digits they share cancel. The stable one works out the big root, where nothing cancels, and divides: the product of the roots is c / a = 1.

1) solve  2) floats and cancellation  3) quit
choice> 3
Bye.
```

## Modules

### main.eml

```eml
# P114 quadratic solver: exact roots of ax^2 + bx + c = 0, real or complex,
# and the small root of x^2 + bx + 1 in floats, by the textbook formula and
# the stable one, against its exact digits.
import quad

def plain(s):
    0 => i
    len(s) => j
    while i < j and s[i] == " ":
        i + 1 => i
    while j > i and s[j - 1] == " ":
        j - 1 => j
    return s[i:j]

def number(s, most):
    if len(s) == 0 or len(s) > most:
        return -1
    0 => n
    for c in s:
        if not (c >= "0" and c <= "9"):
            return -1
        n * 10 + int(c) => n
    return n

def ask_coef(name):
    # [coefficient x 1000, text], [] to cancel, ["!"] after a complaint
    plain(input(name + "> ")) => t
    if t == "":
        return []
    quad.parse_coef(t) => c
    if len(c) == 0:
        "Type a number with at most 9 digits before the point and 3 after, such as -2 or 0.125." ^0
        return ["!"]
    return [c[0], quad.trimmed(quad.scaled_text(c[0], 3))]

def equation(a, b, c):
    "" => out
    for part in [[a, "x^2"], [b, "x"], [c, ""]]:
        part[0] => v
        if v[0] != 0:
            v[1] => t
            "+" => sign
            if t[0] == "-":
                "-" => sign
                t[1:len(t)] => t
            if (t == "1") and part[1] != "":
                "" => t
            if out == "":
                if sign == "-":
                    "-" => out
                out + t + part[1] => out
            else:
                out + " " + sign + " " + t + part[1] => out
    if out == "":
        "0" => out
    return out + " = 0"

def solve_menu():
    ask_coef("a") => a
    if len(a) == 0:
        "Cancelled." ^0
        return
    if a[0] == "!":
        return
    ask_coef("b") => b
    if len(b) == 0:
        "Cancelled." ^0
        return
    if b[0] == "!":
        return
    ask_coef("c") => c
    if len(c) == 0:
        "Cancelled." ^0
        return
    if c[0] == "!":
        return
    equation(a, b, c) ^0
    quad.solve(a[0], b[0], c[0]) => r
    if r[0] == "all":
        "Every x solves it: all three coefficients are zero." ^0
        return
    if r[0] == "none":
        "No x solves it: it says a number other than zero equals zero." ^0
        return
    if r[0] == "linear":
        ("With a = 0 it is not quadratic: one root, x = " + r[1] + ".") ^0
        return
    b[0] * b[0] - 4 * a[0] * c[0] => D
    ("Discriminant b^2 - 4ac = " + quad.trimmed(quad.scaled_text(D, 6)) + ".") ^0
    if r[0] == "double":
        ("It is zero: one double root, x = " + r[1] + ".") ^0
    elif r[0] == "two rational":
        ("It is the square of a fraction: two rational roots, x = " + r[1] + " and x = " + r[2] + ".") ^0
    elif r[0] == "two real":
        "It is positive and not the square of a fraction: two irrational roots, x = (-b - sqrt(D)) / 2a and (-b + sqrt(D)) / 2a:" ^0
        ("  x = " + r[2]) ^0
        ("  x = " + r[3]) ^0
        "to 12 places, each digit worked out exactly." ^0
    else:
        "It is negative: two complex roots, x = (-b / 2a) - i (sqrt(-D) / 2a) and its conjugate:" ^0
        ("  x = " + r[1] + " - " + r[2] + "i") ^0
        ("  x = " + r[1] + " + " + r[2] + "i") ^0
        "the imaginary part to 12 places." ^0

def printed_value(x):
    # [M, e] with the printed float equal to M x 10^e exactly, M >= 0
    str(x) => t
    if t[0] == "-":
        t[1:len(t)] => t
    0 => e
    -1 => at
    for i in [0:len(t) - 1]:
        if at == -1 and t[i] == "e":
            i => at
    if at >= 0:
        t[at + 1:len(t)] => es
        1 => sign
        if es[0] == "-":
            0 - 1 => sign
            es[1:len(es)] => es
        elif es[0] == "+":
            es[1:len(es)] => es
        sign * int(es) => e
        t[0:at] => t
    "" => ds
    0 => after
    False => dot
    for c in t:
        if c == ".":
            True => dot
        else:
            ds + c => ds
            if dot:
                after + 1 => after
    return [int(ds), e - after]

def ten_to(k):
    "1" => t
    for i in [1:k]:
        t + "0" => t
    return int(t)

def quality(x, ex):
    # how close the printed float x is to the exact root, -digits x
    # 10^(lead - 19): [k, -1] with k the largest up to 17 whose relative
    # error is at most 10^-k, or [0, whole percent] when it is over 10%
    printed_value(x) => p
    int(ex[0]) => m
    ex[1] - 19 => e
    p[1] => lo
    if e < lo:
        e => lo
    p[0] * ten_to(p[1] - lo) => xv
    if x < 0.0:
        0 - xv => xv
    0 - m * ten_to(e - lo) => ev
    xv - ev => diff
    if diff < 0:
        0 - diff => diff
    0 - ev => size
    0 => k
    while k < 17 and diff * ten_to(k + 1) <= size:
        k + 1 => k
    if k >= 1:
        return [k, 0 - 1]
    return [0, quad.quotient_digits(diff * 100, size, 0)]

def verdict(q):
    if q[1] >= 0:
        return "off by " + str(q[1]) + "%"
    if q[0] >= 17:
        return "right to 17 significant digits or more"
    return "right to about " + str(q[0]) + " significant digits"

def cancellation():
    plain(input("b, a whole number from 3 to 1000000000000> ")) => t
    if t == "":
        "Cancelled." ^0
        return
    number(t, 13) => b
    if b < 3 or b > 1000000000000:
        "Type a whole number from 3 to 1000000000000." ^0
        return
    quad.small_root(b) => r
    quad.exact_small(b, 20) => ex
    ("x^2 + " + str(b) + "x + 1 = 0, its root nearest zero:") ^0
    ("  exact, 20 digits:  -" + ex[0][0:1] + "." + ex[0][1:len(ex[0])] + "e" + str(ex[1])) ^0
    ("  textbook formula:  " + str(r[0]) + "   (" + verdict(quality(r[0], ex)) + ")") ^0
    ("  stable formula:    " + str(r[1]) + "   (" + verdict(quality(r[1], ex)) + ")") ^0
    "The textbook formula takes -b + sqrt(b^2 - 4) when the two are nearly equal, and the digits they share cancel. The stable one works out the big root, where nothing cancels, and divides: the product of the roots is c / a = 1." ^0

"== Quadratic solver ==" ^0
"ax^2 + bx + c = 0, solved exactly - and the same algebra in floats losing its small root." ^0
True => running
while running:
    "" ^0
    "1) solve  2) floats and cancellation  3) quit" ^0
    plain(input("choice> ")) => choice
    if choice == "1":
        solve_menu()
    elif choice == "2":
        cancellation()
    elif choice == "3":
        False => running
    else:
        "Pick a number from 1 to 3." ^0
"Bye." ^0
```

Python projection of main.eml:

```python
import quad

def plain(s):
    i = 0
    j = len(s)
    while i < j and s[i] == " ":
        i = i + 1
    while j > i and s[j - 1] == " ":
        j = j - 1
    return s[i:j]

def number(s, most):
    if len(s) == 0 or len(s) > most:
        return -1
    n = 0
    for c in s:
        if not (c >= "0" and c <= "9"):
            return -1
        n = n * 10 + int(c)
    return n

def ask_coef(name):
    t = plain(input(name + "> "))
    if t == "":
        return []
    c = quad.parse_coef(t)
    if len(c) == 0:
        print("Type a number with at most 9 digits before the point and 3 after, such as -2 or 0.125.")
        return ["!"]
    return [c[0], quad.trimmed(quad.scaled_text(c[0], 3))]

def equation(a, b, c):
    out = ""
    for part in [[a, "x^2"], [b, "x"], [c, ""]]:
        v = part[0]
        if v[0] != 0:
            t = v[1]
            sign = "+"
            if t[0] == "-":
                sign = "-"
                t = t[1:len(t)]
            if t == "1" and part[1] != "":
                t = ""
            if out == "":
                if sign == "-":
                    out = "-"
                out = out + t + part[1]
            else:
                out = out + " " + sign + " " + t + part[1]
    if out == "":
        out = "0"
    return out + " = 0"

def solve_menu():
    a = ask_coef("a")
    if len(a) == 0:
        print("Cancelled.")
        return
    if a[0] == "!":
        return
    b = ask_coef("b")
    if len(b) == 0:
        print("Cancelled.")
        return
    if b[0] == "!":
        return
    c = ask_coef("c")
    if len(c) == 0:
        print("Cancelled.")
        return
    if c[0] == "!":
        return
    print(equation(a, b, c))
    r = quad.solve(a[0], b[0], c[0])
    if r[0] == "all":
        print("Every x solves it: all three coefficients are zero.")
        return
    if r[0] == "none":
        print("No x solves it: it says a number other than zero equals zero.")
        return
    if r[0] == "linear":
        print("With a = 0 it is not quadratic: one root, x = " + r[1] + ".")
        return
    D = b[0] * b[0] - 4 * a[0] * c[0]
    print("Discriminant b^2 - 4ac = " + quad.trimmed(quad.scaled_text(D, 6)) + ".")
    if r[0] == "double":
        print("It is zero: one double root, x = " + r[1] + ".")
    elif r[0] == "two rational":
        print("It is the square of a fraction: two rational roots, x = " + r[1] + " and x = " + r[2] + ".")
    elif r[0] == "two real":
        print("It is positive and not the square of a fraction: two irrational roots, x = (-b - sqrt(D)) / 2a and (-b + sqrt(D)) / 2a:")
        print("  x = " + r[2])
        print("  x = " + r[3])
        print("to 12 places, each digit worked out exactly.")
    else:
        print("It is negative: two complex roots, x = (-b / 2a) - i (sqrt(-D) / 2a) and its conjugate:")
        print("  x = " + r[1] + " - " + r[2] + "i")
        print("  x = " + r[1] + " + " + r[2] + "i")
        print("the imaginary part to 12 places.")

def printed_value(x):
    t = str(x)
    if t[0] == "-":
        t = t[1:len(t)]
    e = 0
    at = -1
    for i in range(0, len(t)):
        if at == -1 and t[i] == "e":
            at = i
    if at >= 0:
        es = t[at + 1:len(t)]
        sign = 1
        if es[0] == "-":
            sign = 0 - 1
            es = es[1:len(es)]
        elif es[0] == "+":
            es = es[1:len(es)]
        e = sign * int(es)
        t = t[0:at]
    ds = ""
    after = 0
    dot = False
    for c in t:
        if c == ".":
            dot = True
        else:
            ds = ds + c
            if dot:
                after = after + 1
    return [int(ds), e - after]

def ten_to(k):
    t = "1"
    for i in range(1, k+1):
        t = t + "0"
    return int(t)

def quality(x, ex):
    p = printed_value(x)
    m = int(ex[0])
    e = ex[1] - 19
    lo = p[1]
    if e < lo:
        lo = e
    xv = p[0] * ten_to(p[1] - lo)
    if x < 0.0:
        xv = 0 - xv
    ev = 0 - m * ten_to(e - lo)
    diff = xv - ev
    if diff < 0:
        diff = 0 - diff
    size = 0 - ev
    k = 0
    while k < 17 and diff * ten_to(k + 1) <= size:
        k = k + 1
    if k >= 1:
        return [k, 0 - 1]
    return [0, quad.quotient_digits(diff * 100, size, 0)]

def verdict(q):
    if q[1] >= 0:
        return "off by " + str(q[1]) + "%"
    if q[0] >= 17:
        return "right to 17 significant digits or more"
    return "right to about " + str(q[0]) + " significant digits"

def cancellation():
    t = plain(input("b, a whole number from 3 to 1000000000000> "))
    if t == "":
        print("Cancelled.")
        return
    b = number(t, 13)
    if b < 3 or b > 1000000000000:
        print("Type a whole number from 3 to 1000000000000.")
        return
    r = quad.small_root(b)
    ex = quad.exact_small(b, 20)
    print("x^2 + " + str(b) + "x + 1 = 0, its root nearest zero:")
    print("  exact, 20 digits:  -" + ex[0][0:1] + "." + ex[0][1:len(ex[0])] + "e" + str(ex[1]))
    print("  textbook formula:  " + str(r[0]) + "   (" + verdict(quality(r[0], ex)) + ")")
    print("  stable formula:    " + str(r[1]) + "   (" + verdict(quality(r[1], ex)) + ")")
    print("The textbook formula takes -b + sqrt(b^2 - 4) when the two are nearly equal, and the digits they share cancel. The stable one works out the big root, where nothing cancels, and divides: the product of the roots is c / a = 1.")

print("== Quadratic solver ==")
print("ax^2 + bx + c = 0, solved exactly - and the same algebra in floats losing its small root.")
running = True
while running:
    print("")
    print("1) solve  2) floats and cancellation  3) quit")
    choice = plain(input("choice> "))
    if choice == "1":
        solve_menu()
    elif choice == "2":
        cancellation()
    elif choice == "3":
        running = False
    else:
        print("Pick a number from 1 to 3.")
print("Bye.")
```

### quad.eml

```eml
# P114 quadratic solver - ax^2 + bx + c = 0 with coefficients of at most
# three decimals, scaled by 1000 to whole numbers so the discriminant
# b^2 - 4ac is exact (the corpus case quadratic-solver, here with complex
# roots too). Square roots and divisions are done digit by digit on whole
# numbers - adding, subtracting, multiplying and comparing only - so every
# printed digit is exact. And as in the corpus case
# catastrophic-cancellation, the textbook formula in floats loses the small
# root of x^2 + bx + 1 when b is large, while the stable one keeps it.

def gcd(a, b):
    if a < 0:
        0 - a => a
    if b < 0:
        0 - b => b
    while b != 0:
        a % b => r
        b => a
        r => b
    return a

def parse_coef(t):
    # the coefficient times 1000, or [] (at most 9 digits before the point
    # and 3 after)
    if len(t) == 0:
        return []
    1 => sign
    t => s
    if s[0] == "-":
        0 - 1 => sign
        s[1:len(s)] => s
    "" => whole
    "" => frac
    False => dot
    for c in s:
        if c == ".":
            if dot:
                return []
            True => dot
        elif c >= "0" and c <= "9":
            if dot:
                frac + c => frac
            else:
                whole + c => whole
        else:
            return []
    if whole == "" or len(frac) > 3 or (dot and frac == "") or len(whole) > 9:
        return []
    while len(frac) < 3:
        frac + "0" => frac
    return [sign * (int(whole) * 1000 + int(frac))]

def root_digits(n, k):
    # floor(sqrt(n) x 10^k) for n >= 0, by the digit-pair method
    str(n) => s
    for i in [1:k]:
        s + "00" => s
    if len(s) % 2 == 1:
        "0" + s => s
    0 => p
    0 => r
    0 => i
    while i < len(s):
        r * 100 + int(s[i:i + 2]) => r
        0 => d
        while (20 * p + d + 1) * (d + 1) <= r:
            d + 1 => d
        r - (20 * p + d) * d => r
        p * 10 + d => p
        i + 2 => i
    return p

def is_square(n):
    root_digits(n, 0) => s
    return s * s == n

def quotient_digits(num, den, k):
    # floor(num / den x 10^k) for num >= 0, den > 0, by long division
    str(num) => s
    for i in [1:k]:
        s + "0" => s
    0 => q
    0 => r
    for c in s:
        r * 10 + int(c) => r
        0 => d
        while r >= den:
            r - den => r
            d + 1 => d
        q * 10 + d => q
    return q

def scaled_text(n, k):
    # n / 10^k as a decimal with k places
    "" => sign
    if n < 0:
        "-" => sign
        0 - n => n
    str(n) => s
    while len(s) <= k:
        "0" + s => s
    if k == 0:
        return sign + s
    return sign + s[0:len(s) - k] + "." + s[len(s) - k:len(s)]

def trimmed(t):
    # drop trailing zeros after the point, and the point if nothing follows
    t => s
    False => dot
    for c in s:
        if c == ".":
            True => dot
    if not dot:
        return s
    while s[len(s) - 1] == "0":
        s[0:len(s) - 1] => s
    if s[len(s) - 1] == ".":
        s[0:len(s) - 1] => s
    return s

def abs_(n):
    if n < 0:
        return 0 - n
    return n

def fraction_text(num, den):
    # num / den reduced, the sign in front; den is not 0
    if den < 0:
        0 - num => num
        0 - den => den
    gcd(num, den) => g
    if g == 0:
        1 => g
    quotient_digits(abs_(num), g, 0) => n
    quotient_digits(den, g, 0) => d
    "" => sign
    if num < 0:
        "-" => sign
    if d == 1:
        return sign + str(n)
    return sign + str(n) + "/" + str(d)

def decimal12(top, root_sign, d, den):
    # (top + root_sign x sqrt(d)) / den to 12 places, rounded half away from
    # zero; den > 0 and d >= 0. Worked to 18 places first: the floor of the
    # root is off by less than one unit there, far below the 12th place.
    root_digits(d, 18) => s
    "1" => t
    for i in [1:18]:
        t + "0" => t
    top * int(t) + root_sign * s => num
    "" => sign
    if num < 0:
        "-" => sign
        0 - num => num
    str(quotient_digits(num, den, 0)) => q
    while len(q) < 7:
        "0" + q => q
    int(q[0:len(q) - 6]) => whole
    if int(q[len(q) - 6:len(q)]) >= 500000:
        whole + 1 => whole
    if whole == 0:
        "" => sign
    return sign + scaled_text(whole, 12)

def solve(A, B, C):
    # [kind, ...] for whole coefficients (already scaled):
    #   ["all"] or ["none"] when A = B = 0; ["linear", root]
    #   ["double", root]; ["two rational", root, root]
    #   ["two real", D, root, root] with D not a square, roots to 12 places
    #   ["complex", real part, imaginary part]
    if A == 0:
        if B == 0:
            if C == 0:
                return ["all"]
            return ["none"]
        return ["linear", fraction_text(0 - C, B)]
    if A < 0:
        # the same roots, with a positive leading coefficient
        0 - A => A
        0 - B => B
        0 - C => C
    B * B - 4 * A * C => D
    2 * A => den
    if D == 0:
        return ["double", fraction_text(0 - B, den)]
    if D > 0 and is_square(D):
        root_digits(D, 0) => s
        return ["two rational", fraction_text(0 - B - s, den), fraction_text(0 - B + s, den)]
    if D > 0:
        return ["two real", D, decimal12(0 - B, 0 - 1, D, den), decimal12(0 - B, 1, D, den)]
    return ["complex", fraction_text(0 - B, den), decimal12(0, 1, 0 - D, den)]

def fsqrt(x):
    # a float square root by Newton's method from above, using only
    # adding, dividing and comparing - the same in every Python
    if x <= 0.0:
        return 0.0
    x => r
    if r < 1.0:
        1.0 => r
    for i in [1:200]:
        (r + x / r) / 2.0 => y
        if y >= r:
            return r
        y => r
    return r

def small_root(b):
    # x^2 + bx + 1 = 0 for whole b >= 3: [textbook, stable] in floats,
    # the root nearer zero
    b * 1.0 => fb
    fsqrt(fb * fb - 4.0) => s
    (0.0 - fb + s) / 2.0 => naive
    (0.0 - fb - s) / 2.0 => big
    1.0 / big => stable
    return [naive, stable]

def exact_small(b, n):
    # the root of x^2 + bx + 1 nearest zero is -2 / (b + sqrt(b^2 - 4));
    # [its first n significant digits, the power of ten of the first one]
    n + len(str(b)) + 5 => k
    "1" => t
    for i in [1:k]:
        t + "0" => t
    int(t) => ten
    root_digits(b * b - 4, k) => s
    str(quotient_digits(2 * ten * ten, b * ten + s, 0)) => q
    return [q[0:n], len(q) - k - 1]
```

Python projection of quad.eml:

```python
def gcd(a, b):
    if a < 0:
        a = 0 - a
    if b < 0:
        b = 0 - b
    while b != 0:
        r = a % b
        a = b
        b = r
    return a

def parse_coef(t):
    if len(t) == 0:
        return []
    sign = 1
    s = t
    if s[0] == "-":
        sign = 0 - 1
        s = s[1:len(s)]
    whole = ""
    frac = ""
    dot = False
    for c in s:
        if c == ".":
            if dot:
                return []
            dot = True
        elif c >= "0" and c <= "9":
            if dot:
                frac = frac + c
            else:
                whole = whole + c
        else:
            return []
    if whole == "" or len(frac) > 3 or dot and frac == "" or len(whole) > 9:
        return []
    while len(frac) < 3:
        frac = frac + "0"
    return [sign * (int(whole) * 1000 + int(frac))]

def root_digits(n, k):
    s = str(n)
    for i in range(1, k+1):
        s = s + "00"
    if len(s) % 2 == 1:
        s = "0" + s
    p = 0
    r = 0
    i = 0
    while i < len(s):
        r = r * 100 + int(s[i:i + 2])
        d = 0
        while (20 * p + d + 1) * (d + 1) <= r:
            d = d + 1
        r = r - (20 * p + d) * d
        p = p * 10 + d
        i = i + 2
    return p

def is_square(n):
    s = root_digits(n, 0)
    return s * s == n

def quotient_digits(num, den, k):
    s = str(num)
    for i in range(1, k+1):
        s = s + "0"
    q = 0
    r = 0
    for c in s:
        r = r * 10 + int(c)
        d = 0
        while r >= den:
            r = r - den
            d = d + 1
        q = q * 10 + d
    return q

def scaled_text(n, k):
    sign = ""
    if n < 0:
        sign = "-"
        n = 0 - n
    s = str(n)
    while len(s) <= k:
        s = "0" + s
    if k == 0:
        return sign + s
    return sign + s[0:len(s) - k] + "." + s[len(s) - k:len(s)]

def trimmed(t):
    s = t
    dot = False
    for c in s:
        if c == ".":
            dot = True
    if not dot:
        return s
    while s[len(s) - 1] == "0":
        s = s[0:len(s) - 1]
    if s[len(s) - 1] == ".":
        s = s[0:len(s) - 1]
    return s

def abs_(n):
    if n < 0:
        return 0 - n
    return n

def fraction_text(num, den):
    if den < 0:
        num = 0 - num
        den = 0 - den
    g = gcd(num, den)
    if g == 0:
        g = 1
    n = quotient_digits(abs_(num), g, 0)
    d = quotient_digits(den, g, 0)
    sign = ""
    if num < 0:
        sign = "-"
    if d == 1:
        return sign + str(n)
    return sign + str(n) + "/" + str(d)

def decimal12(top, root_sign, d, den):
    s = root_digits(d, 18)
    t = "1"
    for i in range(1, 19):
        t = t + "0"
    num = top * int(t) + root_sign * s
    sign = ""
    if num < 0:
        sign = "-"
        num = 0 - num
    q = str(quotient_digits(num, den, 0))
    while len(q) < 7:
        q = "0" + q
    whole = int(q[0:len(q) - 6])
    if int(q[len(q) - 6:len(q)]) >= 500000:
        whole = whole + 1
    if whole == 0:
        sign = ""
    return sign + scaled_text(whole, 12)

def solve(A, B, C):
    if A == 0:
        if B == 0:
            if C == 0:
                return ["all"]
            return ["none"]
        return ["linear", fraction_text(0 - C, B)]
    if A < 0:
        A = 0 - A
        B = 0 - B
        C = 0 - C
    D = B * B - 4 * A * C
    den = 2 * A
    if D == 0:
        return ["double", fraction_text(0 - B, den)]
    if D > 0 and is_square(D):
        s = root_digits(D, 0)
        return ["two rational", fraction_text(0 - B - s, den), fraction_text(0 - B + s, den)]
    if D > 0:
        return ["two real", D, decimal12(0 - B, 0 - 1, D, den), decimal12(0 - B, 1, D, den)]
    return ["complex", fraction_text(0 - B, den), decimal12(0, 1, 0 - D, den)]

def fsqrt(x):
    if x <= 0.0:
        return 0.0
    r = x
    if r < 1.0:
        r = 1.0
    for i in range(1, 201):
        y = (r + x / r) / 2.0
        if y >= r:
            return r
        r = y
    return r

def small_root(b):
    fb = b * 1.0
    s = fsqrt(fb * fb - 4.0)
    naive = (0.0 - fb + s) / 2.0
    big = (0.0 - fb - s) / 2.0
    stable = 1.0 / big
    return [naive, stable]

def exact_small(b, n):
    k = n + len(str(b)) + 5
    t = "1"
    for i in range(1, k+1):
        t = t + "0"
    ten = int(t)
    s = root_digits(b * b - 4, k)
    q = str(quotient_digits(2 * ten * ten, b * ten + s, 0))
    return [q[0:n], len(q) - k - 1]
```

## README

# P114 - Quadratic solver

ax² + bx + c = 0, solved exactly: rational, irrational and complex roots.
The program also solves x² + bx + 1 in floats, by the textbook formula and
by a stable one, against the root's exact digits.

- `main.eml` - the menu, the equation and its roots, and the float
  comparison
- `quad.eml` - reading coefficients, square roots and divisions digit by
  digit, the solver, a float square root, and the small root two ways and
  exactly

How each part works:

- Coefficients have at most three decimals. They are scaled by 1000 to
  whole numbers, which leaves the roots as they are, so the discriminant
  b² - 4ac is exact.
- The kind of root, as in the corpus case `quadratic-solver`, here with
  complex roots too:
  - a discriminant of 0 gives one double root;
  - the square of a fraction gives two rational roots;
  - any other positive one gives two irrational roots;
  - a negative one gives complex roots.

  A leading coefficient below zero is turned positive by changing every
  sign, which leaves the roots as they are.
- Square roots use the digit-pair method and divisions use long division,
  both on whole numbers. Only adding, subtracting, multiplying and
  comparing are used, so the 12 decimal places are exact: 18 are worked
  out, and then rounded half away from zero.
- Floats, as in the corpus case `catastrophic-cancellation`: for x² + bx +
  1 with b large, the textbook formula takes -b + √(b² - 4) when the two
  are nearly equal, and the digits they share cancel. The stable formula
  works out the big root, where nothing cancels, and divides: the two
  roots multiply to 1.
  The float square root is Newton's method from above, so it is the same
  in every Python. The program measures each printed float's relative
  error against 20 exact digits.

What is checked: menu choices 1 to 3; coefficients with at most 9 digits
before the point and 3 after; b from 3 to 10¹² for the float comparison.
An empty answer cancels.

Sessions: `sessions/basic.in` gives:
- x² - 3x + 2 and 2x² - 7x + 3, with rational roots;
- x² - 2, with irrational roots, and x² + 2x + 5, with complex roots;
- x² - 2x + 1, a double root;
- 2x - 4, 0 = 0 and 5 = 0, where a, or everything, is 0;
- -x² + 2, and 0.5x² - 0.125x + 0.007;
- the float comparison for b = 3, 10⁸ (the textbook root off by 49%) and
  10¹² (off by 100%: it comes out as 0).

`sessions/bad-input.in` gives:
- menu choices 0 and x;
- coefficients: a cancel, abc, four decimals, two later cancels, a comma
  and ten digits;
- the float comparison with a cancel, 2, 10¹² + 1 and x.

Built on the verified corpus cases `quadratic-solver` (the quadratic
formula and its discriminant) and `catastrophic-cancellation` (the same
algebra, and one arrangement keeps its digits).