Project 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.
2 modules · 2 recorded sessionstext-menu UI in the terminalupdated 2026-10-15
Every screen below was recorded under CPython. When this page was built, the EML interpreter replayed each session from the same input and printed the same bytes.
About
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 tableroots.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).
Recorded sessions
What the screen shows while someone uses the program. Each typed line appears after its prompt, the way a terminal shows it.
bad-input
interpreter: byte-equal== 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.
What was typed (18 lines)
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
basic
interpreter: byte-equal== 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.
What was typed (14 lines)
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
Modules
The program as written, entry module first. Each module transpiles to its own Python file, which is what eml project run executes.
main.eml(entry)
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 (main.py)
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 (roots.py)
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