<!-- canonical: https://efficientnewlanguage.org/eml-p/projects/P082-root-finder/ | updated: 2026-10-15 -->

# P082 Root finder

Bisection and Newton's method side by side on square roots and on roots of polynomials, in exact fractions, with the correct decimals of every step counted exactly - one bit a step against digits that double.

EML-P project `projects/root-finder` in the EML language repo: 2 module(s), entry `main.eml`, terminal UI. There, `eml project run projects/root-finder` runs it and `eml project verify projects/root-finder` 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: bisection-root-finder (https://efficientnewlanguage.org/cases/110-bisection-root-finder/), newton-sqrt (https://efficientnewlanguage.org/cases/120-newton-sqrt/).

## Sessions

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

Input:

```text
0
x
1
1
1000001
abc
16
2
0 1
1
1 2 3 4 5 6 7
1 x
1 0 -4
2 1
-200 0
3 4
2 5
3
```

Screen:

```text
== Root finder ==
Bisection and Newton's method, step by step, with the right digits counted exactly.

1) square root  2) root of a polynomial  3) quit
choice> 0
Pick a number from 1 to 3.

1) square root  2) root of a polynomial  3) quit
choice> x
Pick a number from 1 to 3.

1) square root  2) root of a polynomial  3) quit
choice> 1
number (2 to 1000000)> 1
Type a whole number from 2 to 1000000.
number (2 to 1000000)> 1000001
Type a whole number from 2 to 1000000.
number (2 to 1000000)> abc
Type a whole number from 2 to 1000000.
number (2 to 1000000)> 16
16 is 4 squared; there is nothing to find.

1) square root  2) root of a polynomial  3) quit
choice> 2
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 0 1
Type 2 to 6 whole numbers from -1000 to 1000, the first not 0.
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1
Type 2 to 6 whole numbers from -1000 to 1000, the first not 0.
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1 2 3 4 5 6 7
Type 2 to 6 whole numbers from -1000 to 1000, the first not 0.
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1 x
Type 2 to 6 whole numbers from -1000 to 1000, the first not 0.
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1 0 -4
an interval where it changes sign (like 1 2)> 2 1
Type two whole numbers from -100 to 100, the smaller first.
an interval where it changes sign (like 1 2)> -200 0
Type two whole numbers from -100 to 100, the smaller first.
an interval where it changes sign (like 1 2)> 3 4
x^2 - 4 has the same sign at 3 and 4, so bisection cannot start there.
an interval where it changes sign (like 1 2)> 2 5
2 is a root of x^2 - 4 exactly.

1) square root  2) root of a polynomial  3) quit
choice> 3
Bye.
```

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

Input:

```text
1
2
1
10
2
1 0 -1 -2
1 2
2
1 -6 11 -6
0 5
2
1 0 -2
-2 -1
3
```

Screen:

```text
== Root finder ==
Bisection and Newton's method, step by step, with the right digits counted exactly.

1) square root  2) root of a polynomial  3) quit
choice> 1
number (2 to 1000000)> 2
The root of x^2 - 2 is 1.414213562373095048801688724210...
  step  bisection          right  Newton             right
     1  1.500000000000000      1  1.500000000000000      1
     2  1.250000000000000      0  1.416666666666666      2
     3  1.375000000000000      1  1.414215686274509      5
     4  1.437500000000000      1  1.414213562374689     11
     5  1.406250000000000      2  1.414213562373095     24
     6  1.421875000000000      2  1.414213562373095     48
     7  1.414062500000000      3  1.414213562373095     60
     8  1.417968750000000      2  1.414213562373095     60
Bisection has 15 decimals right after 48 steps: it halves the interval, about a third of a digit a step.
Newton has them after 5 steps: near the root the number of right digits about doubles each step.
(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)
(Bisection starts from 1 and 2, Newton from 2.)

1) square root  2) root of a polynomial  3) quit
choice> 1
number (2 to 1000000)> 10
The root of x^2 - 10 is 3.162277660168379331998893544433...
  step  bisection          right  Newton             right
     1  3.500000000000000      0  5.500000000000000      0
     2  3.250000000000000      1  3.659090909090909      0
     3  3.125000000000000      1  3.196005081874647      1
     4  3.187500000000000      1  3.162455622803890      3
     5  3.156250000000000      2  3.162277665175674      8
     6  3.171875000000000      2  3.162277660168379     17
     7  3.164062500000000      2  3.162277660168379     35
     8  3.160156250000000      2  3.162277660168379     60
Bisection has 15 decimals right after 49 steps: it halves the interval, about a third of a digit a step.
Newton has them after 6 steps: near the root the number of right digits about doubles each step.
(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)
(Bisection starts from 3 and 4, Newton from 10.)

1) square root  2) root of a polynomial  3) quit
choice> 2
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1 0 -1 -2
an interval where it changes sign (like 1 2)> 1 2
The root of x^3 - x - 2 is 1.521379706804567569604080832254...
  step  bisection          right  Newton             right
     1  1.500000000000000      1  1.521739130434782      3
     2  1.750000000000000      0  1.521379805964786      7
     3  1.625000000000000      0  1.521379706804575     14
     4  1.562500000000000      1  1.521379706804567     28
     5  1.531250000000000      2  1.521379706804567     56
     6  1.515625000000000      2  1.521379706804567     60
     7  1.523437500000000      2  1.521379706804567     60
     8  1.519531250000000      2  1.521379706804567     60
Bisection has 15 decimals right after 49 steps: it halves the interval, about a third of a digit a step.
Newton has them after 4 steps: near the root the number of right digits about doubles each step.
(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)
(Both start from the interval; Newton from its middle.)

1) square root  2) root of a polynomial  3) quit
choice> 2
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1 -6 11 -6
an interval where it changes sign (like 1 2)> 0 5
The root of x^3 - 6x^2 + 11x - 6 is 3.000000000000000000000000000000...
  step  bisection          right  Newton             right
     1  2.500000000000000      0  1.000000000000000      0
     2  3.750000000000000      0  1.000000000000000      0
     3  3.125000000000000      0  1.000000000000000      0
     4  2.812500000000000      0  1.000000000000000      0
     5  2.968750000000000      1  1.000000000000000      0
     6  3.046875000000000      1  1.000000000000000      0
     7  3.007812500000000      2  1.000000000000000      0
     8  2.988281250000000      1  1.000000000000000      0
Bisection has 15 decimals right after 51 steps: it halves the interval, about a third of a digit a step.
Newton settled on a different root, 1.000000000000000: from its start the slope sent it there, not to the root bisection keeps.
(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)
(Both start from the interval; Newton from its middle.)

1) square root  2) root of a polynomial  3) quit
choice> 2
coefficients, highest power first (like 1 0 -2 for x^2 - 2)> 1 0 -2
an interval where it changes sign (like 1 2)> -2 -1
The root of x^2 - 2 is -1.414213562373095048801688724210...
  step  bisection           right  Newton              right
     1  -1.500000000000000      1  -1.416666666666666      2
     2  -1.250000000000000      0  -1.414215686274509      5
     3  -1.375000000000000      1  -1.414213562374689     11
     4  -1.437500000000000      1  -1.414213562373095     24
     5  -1.406250000000000      2  -1.414213562373095     48
     6  -1.421875000000000      2  -1.414213562373095     60
     7  -1.414062500000000      3  -1.414213562373095     60
     8  -1.417968750000000      2  -1.414213562373095     60
Bisection has 15 decimals right after 48 steps: it halves the interval, about a third of a digit a step.
Newton has them after 4 steps: near the root the number of right digits about doubles each step.
(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)
(Both start from the interval; Newton from its middle.)

1) square root  2) root of a polynomial  3) quit
choice> 3
Bye.
```

## Modules

### main.eml

```eml
# P082 root finder: bisection and Newton's method side by side, on square
# roots and on roots of polynomials, with every step's correct digits
# counted exactly.
import roots

def trim(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 integer(s):
    # A whole number, maybe negative, of up to 7 digits, as [ok, value].
    s => t
    1 => sg
    if len(t) > 1 and t[0] == "-":
        t[1:len(t)] => t
        -1 => sg
    if len(t) == 0 or len(t) > 7:
        return [False, 0]
    0 => n
    for c in t:
        if not (c >= "0" and c <= "9"):
            return [False, 0]
        n * 10 + int(c) => n
    return [True, sg * n]

def numbers(s):
    # "1 0 -2" -> [1, 0, -2], or [] if any is not a whole number.
    [] => out
    "" => w
    for c in s + " ":
        if c == " ":
            if w != "":
                integer(w) => v
                if not v[0]:
                    return []
                out + [v[1]] => out
            "" => w
        else:
            w + c => w
    return out

def lpad(s, n):
    while len(s) < n:
        " " + s => s
    return s

def pad(s, n):
    while len(s) < n:
        s + " " => s
    return s

def poly_text(cs):
    len(cs) - 1 => n
    "" => s
    for i in [0:n]:
        cs[i] => c
        n - i => p
        if c != 0:
            c => a
            if a < 0:
                0 - a => a
            if s == "":
                if c < 0:
                    "-" => s
            elif c < 0:
                s + " - " => s
            else:
                s + " + " => s
            if a != 1 or p == 0:
                s + str(a) => s
            if p == 1:
                s + "x" => s
            elif p > 1:
                s + "x^" + str(p) => s
    return s

def race(cs, a, b, x0, label):
    # The two methods on the root of cs in [a, b], step by step.
    roots.reference(cs, a, b) => ref
    roots.bisect_steps(cs, a, b, 70) => bi
    roots.newton_steps(cs, x0, 30, 60) => nw
    ("The root of " + label + " is " + roots.decimal(roots.rounded(ref, 30), 30) + "...") ^0
    [] => rows
    10 => wb
    10 => wn
    for k in [0:7]:
        roots.decimal(bi[0][k], 15) => bs
        str(roots.correct(bi[0][k], ref, 60)) => bd
        "-" => ns
        "" => nd
        if k < len(nw[0]):
            roots.decimal(nw[0][k], 15) => ns
            str(roots.correct(nw[0][k], ref, 60)) => nd
        if len(bs) > wb:
            len(bs) => wb
        if len(ns) > wn:
            len(ns) => wn
        rows + [[bs, bd, ns, nd]] => rows
    ("  step  " + pad("bisection", wb) + "  right  " + pad("Newton", wn) + "  right") ^0
    for k in [0:7]:
        rows[k] => r
        ("  " + lpad(str(k + 1), 4) + "  " + pad(r[0], wb) + lpad(r[1], 7) + "  " + pad(r[2], wn) + lpad(r[3], 7)) ^0
    -1 => need_b
    for k in [0:69]:
        if need_b == -1 and roots.correct(bi[0][k], ref, 60) >= 15:
            k + 1 => need_b
    -1 => need_n
    for k in [0:len(nw[0]) - 1]:
        if need_n == -1 and roots.correct(nw[0][k], ref, 60) >= 15:
            k + 1 => need_n
    if need_b > 0:
        ("Bisection has 15 decimals right after " + str(need_b) + " steps: it halves the interval, about a third of a digit a step.") ^0
    if need_n > 0:
        ("Newton has them after " + str(need_n) + " steps: near the root the number of right digits about doubles each step.") ^0
    elif nw[1] != "":
        ("Newton stopped after " + str(len(nw[0])) + " steps: " + nw[1] + ".") ^0
    elif len(nw[0]) >= 2 and nw[0][len(nw[0]) - 1] == nw[0][len(nw[0]) - 2]:
        ("Newton settled on a different root, " + roots.decimal(nw[0][len(nw[0]) - 1], 15) + ": from its start the slope sent it there, not to the root bisection keeps.") ^0
    else:
        "Newton did not get 15 decimals of this root in 30 steps: it wandered." ^0
    "(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)" ^0

def square_root():
    while True:
        trim(input("number (2 to 1000000)> ")) => s
        if s == "":
            "Cancelled." ^0
            return
        integer(s) => v
        if v[0] and v[1] >= 2 and v[1] <= 1000000:
            v[1] => n
            1 => k
            while (k + 1) * (k + 1) <= n:
                k + 1 => k
            if k * k == n:
                (str(n) + " is " + str(k) + " squared; there is nothing to find.") ^0
                return
            race([1, 0, 0 - n], [k, 1], [k + 1, 1], [n, 1], "x^2 - " + str(n))
            ("(Bisection starts from " + str(k) + " and " + str(k + 1) + ", Newton from " + str(n) + ".)") ^0
            return
        "Type a whole number from 2 to 1000000." ^0

def poly_root():
    while True:
        trim(input("coefficients, highest power first (like 1 0 -2 for x^2 - 2)> ")) => s
        if s == "":
            "Cancelled." ^0
            return
        numbers(s) => cs
        True => ok
        for c in cs:
            if c < -1000 or c > 1000:
                False => ok
        if ok and len(cs) >= 2 and len(cs) <= 6 and cs[0] != 0:
            poly_text(cs) => label
            while True:
                trim(input("an interval where it changes sign (like 1 2)> ")) => t
                if t == "":
                    "Cancelled." ^0
                    return
                numbers(t) => ab
                if len(ab) == 2 and ab[0] < ab[1] and ab[0] >= -100 and ab[1] <= 100:
                    [ab[0], 1] => a
                    [ab[1], 1] => b
                    roots.sign(roots.value(cs, a)) => sa
                    roots.sign(roots.value(cs, b)) => sb
                    if sa == 0:
                        (str(ab[0]) + " is a root of " + label + " exactly.") ^0
                        return
                    if sb == 0:
                        (str(ab[1]) + " is a root of " + label + " exactly.") ^0
                        return
                    if sa == sb:
                        (label + " has the same sign at " + str(ab[0]) + " and " + str(ab[1]) + ", so bisection cannot start there.") ^0
                    else:
                        race(cs, a, b, roots.mul(roots.add(a, b), [1, 2]), label)
                        "(Both start from the interval; Newton from its middle.)" ^0
                        return
                else:
                    "Type two whole numbers from -100 to 100, the smaller first." ^0
        else:
            "Type 2 to 6 whole numbers from -1000 to 1000, the first not 0." ^0

"== Root finder ==" ^0
"Bisection and Newton's method, step by step, with the right digits counted exactly." ^0
True => running
while running:
    "" ^0
    "1) square root  2) root of a polynomial  3) quit" ^0
    trim(input("choice> ")) => choice
    if choice == "1":
        square_root()
    elif choice == "2":
        poly_root()
    elif choice == "3":
        False => running
    else:
        "Pick a number from 1 to 3." ^0
"Bye." ^0
```

Python projection of main.eml:

```python
import roots

def trim(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 integer(s):
    t = s
    sg = 1
    if len(t) > 1 and t[0] == "-":
        t = t[1:len(t)]
        sg = -1
    if len(t) == 0 or len(t) > 7:
        return [False, 0]
    n = 0
    for c in t:
        if not (c >= "0" and c <= "9"):
            return [False, 0]
        n = n * 10 + int(c)
    return [True, sg * n]

def numbers(s):
    out = []
    w = ""
    for c in s + " ":
        if c == " ":
            if w != "":
                v = integer(w)
                if not v[0]:
                    return []
                out = out + [v[1]]
            w = ""
        else:
            w = w + c
    return out

def lpad(s, n):
    while len(s) < n:
        s = " " + s
    return s

def pad(s, n):
    while len(s) < n:
        s = s + " "
    return s

def poly_text(cs):
    n = len(cs) - 1
    s = ""
    for i in range(0, n+1):
        c = cs[i]
        p = n - i
        if c != 0:
            a = c
            if a < 0:
                a = 0 - a
            if s == "":
                if c < 0:
                    s = "-"
            elif c < 0:
                s = s + " - "
            else:
                s = s + " + "
            if a != 1 or p == 0:
                s = s + str(a)
            if p == 1:
                s = s + "x"
            elif p > 1:
                s = s + "x^" + str(p)
    return s

def race(cs, a, b, x0, label):
    ref = roots.reference(cs, a, b)
    bi = roots.bisect_steps(cs, a, b, 70)
    nw = roots.newton_steps(cs, x0, 30, 60)
    print("The root of " + label + " is " + roots.decimal(roots.rounded(ref, 30), 30) + "...")
    rows = []
    wb = 10
    wn = 10
    for k in range(0, 8):
        bs = roots.decimal(bi[0][k], 15)
        bd = str(roots.correct(bi[0][k], ref, 60))
        ns = "-"
        nd = ""
        if k < len(nw[0]):
            ns = roots.decimal(nw[0][k], 15)
            nd = str(roots.correct(nw[0][k], ref, 60))
        if len(bs) > wb:
            wb = len(bs)
        if len(ns) > wn:
            wn = len(ns)
        rows = rows + [[bs, bd, ns, nd]]
    print("  step  " + pad("bisection", wb) + "  right  " + pad("Newton", wn) + "  right")
    for k in range(0, 8):
        r = rows[k]
        print("  " + lpad(str(k + 1), 4) + "  " + pad(r[0], wb) + lpad(r[1], 7) + "  " + pad(r[2], wn) + lpad(r[3], 7))
    need_b = -1
    for k in range(0, 70):
        if need_b == -1 and roots.correct(bi[0][k], ref, 60) >= 15:
            need_b = k + 1
    need_n = -1
    for k in range(0, len(nw[0])):
        if need_n == -1 and roots.correct(nw[0][k], ref, 60) >= 15:
            need_n = k + 1
    if need_b > 0:
        print("Bisection has 15 decimals right after " + str(need_b) + " steps: it halves the interval, about a third of a digit a step.")
    if need_n > 0:
        print("Newton has them after " + str(need_n) + " steps: near the root the number of right digits about doubles each step.")
    elif nw[1] != "":
        print("Newton stopped after " + str(len(nw[0])) + " steps: " + nw[1] + ".")
    elif len(nw[0]) >= 2 and nw[0][len(nw[0]) - 1] == nw[0][len(nw[0]) - 2]:
        print("Newton settled on a different root, " + roots.decimal(nw[0][len(nw[0]) - 1], 15) + ": from its start the slope sent it there, not to the root bisection keeps.")
    else:
        print("Newton did not get 15 decimals of this root in 30 steps: it wandered.")
    print("(Bisection's midpoints are exact fractions; Newton's values are kept to 60 decimals, as a float keeps a fixed number of bits.)")

def square_root():
    while True:
        s = trim(input("number (2 to 1000000)> "))
        if s == "":
            print("Cancelled.")
            return
        v = integer(s)
        if v[0] and v[1] >= 2 and v[1] <= 1000000:
            n = v[1]
            k = 1
            while (k + 1) * (k + 1) <= n:
                k = k + 1
            if k * k == n:
                print(str(n) + " is " + str(k) + " squared; there is nothing to find.")
                return
            race([1, 0, 0 - n], [k, 1], [k + 1, 1], [n, 1], "x^2 - " + str(n))
            print("(Bisection starts from " + str(k) + " and " + str(k + 1) + ", Newton from " + str(n) + ".)")
            return
        print("Type a whole number from 2 to 1000000.")

def poly_root():
    while True:
        s = trim(input("coefficients, highest power first (like 1 0 -2 for x^2 - 2)> "))
        if s == "":
            print("Cancelled.")
            return
        cs = numbers(s)
        ok = True
        for c in cs:
            if c < -1000 or c > 1000:
                ok = False
        if ok and len(cs) >= 2 and len(cs) <= 6 and cs[0] != 0:
            label = poly_text(cs)
            while True:
                t = trim(input("an interval where it changes sign (like 1 2)> "))
                if t == "":
                    print("Cancelled.")
                    return
                ab = numbers(t)
                if len(ab) == 2 and ab[0] < ab[1] and ab[0] >= -100 and ab[1] <= 100:
                    a = [ab[0], 1]
                    b = [ab[1], 1]
                    sa = roots.sign(roots.value(cs, a))
                    sb = roots.sign(roots.value(cs, b))
                    if sa == 0:
                        print(str(ab[0]) + " is a root of " + label + " exactly.")
                        return
                    if sb == 0:
                        print(str(ab[1]) + " is a root of " + label + " exactly.")
                        return
                    if sa == sb:
                        print(label + " has the same sign at " + str(ab[0]) + " and " + str(ab[1]) + ", so bisection cannot start there.")
                    else:
                        race(cs, a, b, roots.mul(roots.add(a, b), [1, 2]), label)
                        print("(Both start from the interval; Newton from its middle.)")
                        return
                else:
                    print("Type two whole numbers from -100 to 100, the smaller first.")
        else:
            print("Type 2 to 6 whole numbers from -1000 to 1000, the first not 0.")

print("== Root finder ==")
print("Bisection and Newton's method, step by step, with the right digits counted exactly.")
running = True
while running:
    print("")
    print("1) square root  2) root of a polynomial  3) quit")
    choice = trim(input("choice> "))
    if choice == "1":
        square_root()
    elif choice == "2":
        poly_root()
    elif choice == "3":
        running = False
    else:
        print("Pick a number from 1 to 3.")
print("Bye.")
```

### roots.eml

```eml
# P082 root finder - bisection and Newton's method on square roots and on
# polynomials with whole-number coefficients, in exact fractions [n, d]
# (d > 0), so how many digits each step has right can be counted exactly.
# The corpus cases bisection-root-finder and newton-sqrt do the same in
# floats, where both stop at about 16 digits.

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

def frac(n, d):
    if d < 0:
        0 - n => n
        0 - d => d
    gcd(n, d) => g
    if g > 1:
        divide(n, g) => n
        divide(d, g) => d
    return [n, d]

def divide(a, b):
    # a / b for b > 0 that divides a exactly (a may be negative), by long
    # division digit by digit - int(a / b) would go through a float.
    if a < 0:
        return 0 - divide(0 - a, b)
    0 => q
    0 => r
    for ch in str(a):
        r * 10 + int(ch) => r
        0 => k
        while r >= b:
            r - b => r
            k + 1 => k
        q * 10 + k => q
    return q

def add(x, y):
    return frac(x[0] * y[1] + y[0] * x[1], x[1] * y[1])

def sub(x, y):
    return frac(x[0] * y[1] - y[0] * x[1], x[1] * y[1])

def mul(x, y):
    return frac(x[0] * y[0], x[1] * y[1])

def div(x, y):
    return frac(x[0] * y[1], x[1] * y[0])

def value(coeffs, x):
    # Horner's rule; coeffs are whole numbers, highest power first.
    [0, 1] => v
    for c in coeffs:
        add(mul(v, x), [c, 1]) => v
    return v

def slope(coeffs):
    # The derivative's coefficients.
    len(coeffs) - 1 => n
    [] => out
    for i in [0:n - 1]:
        out + [coeffs[i] * (n - i)] => out
    return out

def sign(x):
    if x[0] > 0:
        return 1
    if x[0] < 0:
        return -1
    return 0

def rounded(x, places):
    # x to the nearest multiple of 10^-places (halves up): Newton's values
    # are kept to this many places, as a float keeps a fixed number of bits.
    1 => s
    for i in [1:places]:
        s * 10 => s
    x[0] * s * 2 + x[1] => n
    2 * x[1] => d
    # floor(n / d), with n possibly negative
    if n >= 0:
        divide(n - n % d, d) => q
    else:
        0 - divide((0 - n) + (d - (0 - n) % d) % d, d) => q
    return frac(q, s)

def bisect_steps(coeffs, a, b, steps):
    # The midpoints of `steps` halvings of [a, b], where coeffs changes
    # sign: each time keep the half whose ends still differ in sign.
    # Returns [midpoints, a, b].
    sign(value(coeffs, a)) => sa
    [] => mids
    for k in [1:steps]:
        mul(add(a, b), [1, 2]) => m
        mids + [m] => mids
        sign(value(coeffs, m)) => sm
        if sm == 0:
            m => a
            m => b
        elif sm == sa:
            m => a
        else:
            m => b
    return [mids, a, b]

def newton_steps(coeffs, x, steps, places):
    # x <- x - f(x)/f'(x), each value kept to `places` decimals. Returns
    # [values, "" or why it stopped].
    slope(coeffs) => d
    [] => xs
    for k in [1:steps]:
        value(d, x) => fd
        if fd[0] == 0:
            return [xs, "the slope is 0 there"]
        rounded(sub(x, div(value(coeffs, x), fd)), places) => x
        xs + [x] => xs
    return [xs, ""]

def reference(coeffs, a, b):
    # The root in [a, b] by 220 halvings: within (b - a) / 2^221 of it, which
    # for an interval up to 200 wide is below 10^-63 - past the 60 decimals
    # the program ever counts.
    bisect_steps(coeffs, a, b, 220) => r
    return mul(add(r[1], r[2]), [1, 2])

def correct(x, ref, most):
    # How many decimals x has right: the largest d <= most with
    # |x - ref| < 10^-d.
    sub(x, ref) => e
    e[0] => n
    if n < 0:
        0 - n => n
    0 => d
    1 => p
    while d < most and n * p * 10 < e[1]:
        d + 1 => d
        p * 10 => p
    return d

def decimal(x, places):
    # x to `places` decimals, cut off (not rounded), with its sign.
    x[0] => n
    "" => sg
    if n < 0:
        "-" => sg
        0 - n => n
    divide(n - n % x[1], x[1]) => w
    n % x[1] => r
    "" => digits
    for i in [1:places]:
        r * 10 => r
        0 => k
        while r >= x[1]:
            r - x[1] => r
            k + 1 => k
        digits + str(k) => digits
    return sg + str(w) + "." + digits
```

Python projection of roots.eml:

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

def frac(n, d):
    if d < 0:
        n = 0 - n
        d = 0 - d
    g = gcd(n, d)
    if g > 1:
        n = divide(n, g)
        d = divide(d, g)
    return [n, d]

def divide(a, b):
    if a < 0:
        return 0 - divide(0 - a, b)
    q = 0
    r = 0
    for ch in str(a):
        r = r * 10 + int(ch)
        k = 0
        while r >= b:
            r = r - b
            k = k + 1
        q = q * 10 + k
    return q

def add(x, y):
    return frac(x[0] * y[1] + y[0] * x[1], x[1] * y[1])

def sub(x, y):
    return frac(x[0] * y[1] - y[0] * x[1], x[1] * y[1])

def mul(x, y):
    return frac(x[0] * y[0], x[1] * y[1])

def div(x, y):
    return frac(x[0] * y[1], x[1] * y[0])

def value(coeffs, x):
    v = [0, 1]
    for c in coeffs:
        v = add(mul(v, x), [c, 1])
    return v

def slope(coeffs):
    n = len(coeffs) - 1
    out = []
    for i in range(0, n):
        out = out + [coeffs[i] * (n - i)]
    return out

def sign(x):
    if x[0] > 0:
        return 1
    if x[0] < 0:
        return -1
    return 0

def rounded(x, places):
    s = 1
    for i in range(1, places+1):
        s = s * 10
    n = x[0] * s * 2 + x[1]
    d = 2 * x[1]
    if n >= 0:
        q = divide(n - n % d, d)
    else:
        q = 0 - divide(0 - n + (d - (0 - n) % d) % d, d)
    return frac(q, s)

def bisect_steps(coeffs, a, b, steps):
    sa = sign(value(coeffs, a))
    mids = []
    for k in range(1, steps+1):
        m = mul(add(a, b), [1, 2])
        mids = mids + [m]
        sm = sign(value(coeffs, m))
        if sm == 0:
            a = m
            b = m
        elif sm == sa:
            a = m
        else:
            b = m
    return [mids, a, b]

def newton_steps(coeffs, x, steps, places):
    d = slope(coeffs)
    xs = []
    for k in range(1, steps+1):
        fd = value(d, x)
        if fd[0] == 0:
            return [xs, "the slope is 0 there"]
        x = rounded(sub(x, div(value(coeffs, x), fd)), places)
        xs = xs + [x]
    return [xs, ""]

def reference(coeffs, a, b):
    r = bisect_steps(coeffs, a, b, 220)
    return mul(add(r[1], r[2]), [1, 2])

def correct(x, ref, most):
    e = sub(x, ref)
    n = e[0]
    if n < 0:
        n = 0 - n
    d = 0
    p = 1
    while d < most and n * p * 10 < e[1]:
        d = d + 1
        p = p * 10
    return d

def decimal(x, places):
    n = x[0]
    sg = ""
    if n < 0:
        sg = "-"
        n = 0 - n
    w = divide(n - n % x[1], x[1])
    r = n % x[1]
    digits = ""
    for i in range(1, places+1):
        r = r * 10
        k = 0
        while r >= x[1]:
            r = r - x[1]
            k = k + 1
        digits = digits + str(k)
    return sg + str(w) + "." + digits
```

## README

# P082 - Root finder

Bisection and Newton's method side by side, on square roots and on the
roots of polynomials. Every number is an exact fraction, so the correct
decimals of every step can be counted exactly - one bit a step against
digits that double.

- `main.eml` - the menu, reading numbers and polynomials, and the race
  table
- `roots.eml` - fractions, Horner's rule and the derivative, bisection,
  Newton's method, the reference root, and counting correct decimals

How each part works:

- Bisection, as in the corpus case `bisection-root-finder`: an interval
  where the function changes sign holds a root. Halve it and keep the half
  whose ends still differ in sign. Each step gains one bit, about a third
  of a decimal digit; it needs only the sign, and it cannot fail.
- Newton's method, as in the corpus case `newton-sqrt`: replace x by
  x - f(x)/f'(x), the point where the tangent crosses zero. Near a simple
  root the number of correct digits about doubles each step. It needs the
  derivative, and from a bad start it can settle on a different root.
- The corpus cases work in floats, where both stop at about 16 digits.
  Here bisection's midpoints are exact fractions. Newton's values are
  rounded to 60 decimals each step, as a float keeps a fixed number of
  bits; without that, its fractions would double in length every step.
- Correct decimals: the root is first found by 220 halvings, which is
  within (b - a) / 2^221 of it - below 10^-63 for any interval allowed.
  A value then has d correct decimals when it is within 10^-d of that
  root.

What is checked: menu choices 1 to 3; a whole number from 2 to 1,000,000
for a square root (a perfect square is reported at once); 2 to 6 whole
coefficients from -1000 to 1000, the first not 0; an interval of two whole
numbers from -100 to 100, the smaller first, where the polynomial changes
sign. An empty answer cancels.

Sessions: `sessions/basic.in` gives:
- the square roots of 2 and 10: bisection needs 48 and 49 steps for 15
  decimals, Newton 5 and 6;
- the root of x^3 - x - 2 between 1 and 2: Newton in 4 steps;
- x^3 - 6x^2 + 11x - 6 between 0 and 5: bisection finds 3, while Newton,
  started at 2.5, lands on the root 1 in one step;
- the negative root of x^2 - 2.

`sessions/bad-input.in` gives:
- menu choices 0 and x;
- 1, 1000001 and abc, then 16, a perfect square;
- coefficients "0 1", "1", seven numbers and "1 x";
- for x^2 - 4: intervals "2 1" and "-200 0", "3 4" with no sign change,
  and "2 5", where 2 is a root exactly.

Built on the verified corpus cases `bisection-root-finder` (halving an
interval with a sign change) and `newton-sqrt` (the tangent step, and its
digits that double).