theory Sqrt_Babylonian
imports
Sqrt_Babylonian_Auxiliary
NthRoot_Impl
begin
section {*Executable algorithms for square roots*}
text {*
This theory provides executable algorithms for computing square-roots of numbers which
are all based on the Babylonian method (which is also known as Heron's method or Newton's method).
For integers / naturals / rationals precise algorithms are given, i.e., here $sqrt\ x$ delivers
a list of all integers / naturals / rationals $y$ where $y^2 = x$.
To this end, the Babylonian method has been adapted by using integer-divisions.
In addition to the precise algorithms, we also provide approximation algorithms. One works for
arbitrary linear ordered fields, where some number $y$ is computed such that
@{term "abs(y^2 - x) < ε"}. Moreover, for the naturals, integers, and rationals we provide algorithms to compute
@{term "floor (sqrt x)"} and @{term "ceiling (sqrt x)"} which are all based
on the underlying algorithm that is used to compute the precise square-roots on integers, if these
exist.
The major motivation for developing the precise algorithms was given by \ceta{} \cite{CeTA},
a tool for certifiying termination proofs. Here, non-linear equations of the form
$(a_1x_1 + \dots a_nx_n)^2 = p$ had to be solved over the integers, where $p$ is a concrete polynomial.
For example, for the equation $(ax + by)^2 = 4x^2 - 12xy + 9y^2$ one easily figures out that
$a^2 = 4, b^2 = 9$, and $ab = -6$, which results in a possible solution $a = \sqrt 4 = 2, b = - \sqrt 9 = -3$.
*}
subsection {* The Babylonian method *}
text {*
The Babylonian method for computing $\sqrt n$ iteratively computes
\[
x_{i+1} = \frac{\frac n{x_i} + x_i}2
\]
until $x_i^2 \approx n$. Note that if $x_0^2 \geq n$, then for all $i$ we have both
$x_i^2 \geq n$ and $x_i \geq x_{i+1}$.
*}
subsection {* The Babylonian method using integer division *}
text {*
First, the algorithm is developed for the non-negative integers.
Here, the division operation $\frac xy$ is replaced by @{term "x div y = ⌊of_int x / of_int y⌋"}.
Note that replacing @{term "⌊of_int x / of_int y⌋"} by @{term "⌈of_int x / of_int y⌉"} would lead to non-termination
in the following algorithm.
We explicititly develop the algorithm on the integers and not on the naturals, as the calculations
on the integers have been much easier. For example, $y - x + x = y$ on the integers, which would require
the side-condition $y \geq x$ for the naturals. These conditions will make the reasoning much more tedious---as
we have experienced in an earlier state of this development where everything was based on naturals.
Since the elements
$x_0, x_1, x_2,\dots$ are monotone decreasing, in the main algorithm we abort as soon as $x_i^2 \leq n$. *}
text {* \textbf{Since in the meantime, all of these algorithms have been generalized to arbitrary
$p$-th roots in @{theory Sqrt_Babylonian.NthRoot_Impl}, we just instantiate the general algorithms by $p = 2$ and then provide
specialized code equations which are more efficient than the general purpose algorithms.} *}
definition sqrt_int_main' :: "int ⇒ int ⇒ int × bool" where
[simp]: "sqrt_int_main' x n = root_int_main' 1 1 2 x n"
lemma sqrt_int_main'_code[code]: "sqrt_int_main' x n = (let x2 = x * x in if x2 ≤ n then (x, x2 = n)
else sqrt_int_main' ((n div x + x) div 2) n)"
using root_int_main'.simps[of 1 1 2 x n]
unfolding Let_def by auto
definition sqrt_int_main :: "int ⇒ int × bool" where
[simp]: "sqrt_int_main x = root_int_main 2 x"
lemma sqrt_int_main_code[code]: "sqrt_int_main x = sqrt_int_main' (start_value x 2) x"
by (simp add: root_int_main_def Let_def)
definition sqrt_int :: "int ⇒ int list" where
"sqrt_int x = root_int 2 x"
lemma sqrt_int_code[code]: "sqrt_int x = (if x < 0 then [] else case sqrt_int_main x of (y,True) ⇒ if y = 0 then [0] else [y,-y] | _ ⇒ [])"
proof -
interpret fixed_root 2 1 by (unfold_locales, auto)
obtain b y where res: "root_int_main 2 x = (b,y)" by force
show ?thesis
unfolding sqrt_int_def root_int_def Let_def
using root_int_main[OF _ res]
using res
by simp
qed
lemma sqrt_int[simp]: "set (sqrt_int x) = {y. y * y = x}"
unfolding sqrt_int_def by (simp add: power2_eq_square)
lemma sqrt_int_pos: assumes res: "sqrt_int x = Cons s ms"
shows "s ≥ 0"
proof -
note res = res[unfolded sqrt_int_code Let_def, simplified]
from res have x0: "x ≥ 0" by (cases ?thesis, auto)
obtain ss b where call: "sqrt_int_main x = (ss,b)" by force
from res[unfolded call] x0 have "ss = s"
by (cases b, cases "ss = 0", auto)
from root_int_main(1)[OF x0 call[unfolded this sqrt_int_main_def]]
show ?thesis .
qed
definition [simp]: "sqrt_int_floor_pos x = root_int_floor_pos 2 x"
lemma sqrt_int_floor_pos_code[code]: "sqrt_int_floor_pos x = fst (sqrt_int_main x)"
by (simp add: root_int_floor_pos_def)
lemma sqrt_int_floor_pos: assumes x: "x ≥ 0"
shows "sqrt_int_floor_pos x = ⌊ sqrt (of_int x) ⌋"
using root_int_floor_pos[OF x, of 2] by (simp add: sqrt_def)
definition [simp]: "sqrt_int_ceiling_pos x = root_int_ceiling_pos 2 x"
lemma sqrt_int_ceiling_pos_code[code]: "sqrt_int_ceiling_pos x = (case sqrt_int_main x of (y,b) ⇒ if b then y else y + 1)"
by (simp add: root_int_ceiling_pos_def)
lemma sqrt_int_ceiling_pos: assumes x: "x ≥ 0"
shows "sqrt_int_ceiling_pos x = ⌈ sqrt (of_int x) ⌉"
using root_int_ceiling_pos[OF x, of 2] by (simp add: sqrt_def)
definition "sqrt_int_floor x = root_int_floor 2 x"
lemma sqrt_int_floor_code[code]: "sqrt_int_floor x = (if x ≥ 0 then sqrt_int_floor_pos x else - sqrt_int_ceiling_pos (- x))"
unfolding sqrt_int_floor_def root_int_floor_def by simp
lemma sqrt_int_floor[simp]: "sqrt_int_floor x = ⌊ sqrt (of_int x) ⌋"
by (simp add: sqrt_int_floor_def sqrt_def)
definition "sqrt_int_ceiling x = root_int_ceiling 2 x"
lemma sqrt_int_ceiling_code[code]: "sqrt_int_ceiling x = (if x ≥ 0 then sqrt_int_ceiling_pos x else - sqrt_int_floor_pos (- x))"
unfolding sqrt_int_ceiling_def root_int_ceiling_def by simp
lemma sqrt_int_ceiling[simp]: "sqrt_int_ceiling x = ⌈ sqrt (of_int x) ⌉"
by (simp add: sqrt_int_ceiling_def sqrt_def)
lemma sqrt_int_ceiling_bound: "0 ≤ x ⟹ x ≤ (sqrt_int_ceiling x)^2"
unfolding sqrt_int_ceiling using le_of_int_ceiling sqrt_le_D
by (metis of_int_power_le_of_int_cancel_iff)
subsection {* Square roots for the naturals *}
definition sqrt_nat :: "nat ⇒ nat list"
where "sqrt_nat x = root_nat 2 x"
lemma sqrt_nat_code[code]: "sqrt_nat x ≡ map nat (take 1 (sqrt_int (int x)))"
unfolding sqrt_nat_def root_nat_def sqrt_int_def by simp
lemma sqrt_nat[simp]: "set (sqrt_nat x) = { y. y * y = x}"
unfolding sqrt_nat_def using root_nat[of 2 x] by (simp add: power2_eq_square)
definition sqrt_nat_floor :: "nat ⇒ int" where
"sqrt_nat_floor x = root_nat_floor 2 x"
lemma sqrt_nat_floor_code[code]: "sqrt_nat_floor x = sqrt_int_floor_pos (int x)"
unfolding sqrt_nat_floor_def root_nat_floor_def by simp
lemma sqrt_nat_floor[simp]: "sqrt_nat_floor x = ⌊ sqrt (real x) ⌋"
unfolding sqrt_nat_floor_def by (simp add: sqrt_def)
definition sqrt_nat_ceiling :: "nat ⇒ int" where
"sqrt_nat_ceiling x = root_nat_ceiling 2 x"
lemma sqrt_nat_ceiling_code[code]: "sqrt_nat_ceiling x = sqrt_int_ceiling_pos (int x)"
unfolding sqrt_nat_ceiling_def root_nat_ceiling_def by simp
lemma sqrt_nat_ceiling[simp]: "sqrt_nat_ceiling x = ⌈ sqrt (real x) ⌉"
unfolding sqrt_nat_ceiling_def by (simp add: sqrt_def)
subsection {* Square roots for the rationals *}
definition sqrt_rat :: "rat ⇒ rat list" where
"sqrt_rat x = root_rat 2 x"
lemma sqrt_rat_code[code]: "sqrt_rat x = (case quotient_of x of (z,n) ⇒ (case sqrt_int n of
[] ⇒ []
| sn # xs ⇒ map (λ sz. of_int sz / of_int sn) (sqrt_int z)))"
proof -
obtain z n where q: "quotient_of x = (z,n)" by force
show ?thesis
unfolding sqrt_rat_def root_rat_def q split sqrt_int_def
by (cases "root_int 2 n", auto)
qed
lemma sqrt_rat[simp]: "set (sqrt_rat x) = { y. y * y = x}"
unfolding sqrt_rat_def using root_rat[of 2 x]
by (simp add: power2_eq_square)
lemma sqrt_rat_pos: assumes sqrt: "sqrt_rat x = Cons s ms"
shows "s ≥ 0"
proof -
obtain z n where q: "quotient_of x = (z,n)" by force
note sqrt = sqrt[unfolded sqrt_rat_code q, simplified]
let ?sz = "sqrt_int z"
let ?sn = "sqrt_int n"
from q have n: "n > 0" by (rule quotient_of_denom_pos)
from sqrt obtain sz mz where sz: "?sz = sz # mz" by (cases ?sn, auto)
from sqrt obtain sn mn where sn: "?sn = sn # mn" by (cases ?sn, auto)
from sqrt_int_pos[OF sz] sqrt_int_pos[OF sn] have pos: "0 ≤ sz" "0 ≤ sn" by auto
from sqrt sz sn have s: "s = of_int sz / of_int sn" by auto
show ?thesis unfolding s using pos
by (metis of_int_0_le_iff zero_le_divide_iff)
qed
definition sqrt_rat_floor :: "rat ⇒ int" where
"sqrt_rat_floor x = root_rat_floor 2 x"
lemma sqrt_rat_floor_code[code]: "sqrt_rat_floor x = (case quotient_of x of (a,b) ⇒ sqrt_int_floor (a * b) div b)"
unfolding sqrt_rat_floor_def root_rat_floor_def by (simp add: sqrt_def)
lemma sqrt_rat_floor[simp]: "sqrt_rat_floor x = ⌊ sqrt (of_rat x) ⌋"
unfolding sqrt_rat_floor_def by (simp add: sqrt_def)
definition sqrt_rat_ceiling :: "rat ⇒ int" where
"sqrt_rat_ceiling x = root_rat_ceiling 2 x"
lemma sqrt_rat_ceiling_code[code]: "sqrt_rat_ceiling x = - (sqrt_rat_floor (-x))"
unfolding sqrt_rat_ceiling_def sqrt_rat_floor_def root_rat_ceiling_def by simp
lemma sqrt_rat_ceiling: "sqrt_rat_ceiling x = ⌈ sqrt (of_rat x) ⌉"
unfolding sqrt_rat_ceiling_def by (simp add: sqrt_def)
lemma sqr_rat_of_int: assumes x: "x * x = rat_of_int i"
shows "∃ j :: int. j * j = i"
proof -
from x have mem: "x ∈ set (sqrt_rat (rat_of_int i))" by simp
from x have "rat_of_int i ≥ 0" by (metis zero_le_square)
hence *: "quotient_of (rat_of_int i) = (i,1)" by (metis quotient_of_int)
have 1: "sqrt_int 1 = [1,-1]" by code_simp
from mem sqrt_rat_code * split 1
have x: "x ∈ rat_of_int ` {y. y * y = i}" by auto
thus ?thesis by auto
qed
subsection {* Approximating square roots *}
text {*
The difference to the previous algorithms is that now we abort, once the distance is below
$\epsilon$.
Moreover, here we use standard division and not integer division.
This part is not yet generalized by @{theory Sqrt_Babylonian.NthRoot_Impl}.
We first provide the executable version without guard @{term "x > 0"} as partial function,
and afterwards prove termination and soundness for a similar algorithm that is defined within the upcoming
locale.
*}
partial_function (tailrec) sqrt_approx_main_impl :: "'a :: linordered_field ⇒ 'a ⇒ 'a ⇒ 'a" where
[code]: "sqrt_approx_main_impl ε n x = (if x * x - n < ε then x else sqrt_approx_main_impl ε n
((n / x + x) / 2))"
text {* We setup a locale where we ensure that we have standard assumptions: positive $\epsilon$ and
positive $n$. We require sort @{term floor_ceiling}, since @{term "⌊ x ⌋"} is used for the termination
argument. *}
locale sqrt_approximation =
fixes ε :: "'a :: {linordered_field,floor_ceiling}"
and n :: 'a
assumes ε : "ε > 0"
and n: "n > 0"
begin
function sqrt_approx_main :: "'a ⇒ 'a" where
"sqrt_approx_main x = (if x > 0 then (if x * x - n < ε then x else sqrt_approx_main
((n / x + x) / 2)) else 0)"
by pat_completeness auto
text {* Termination essentially is a proof of convergence. Here, one complication is the fact
that the limit is not always defined. E.g., if @{typ "'a"} is @{typ rat} then there is no
square root of 2. Therefore, the error-rate $\frac x{\sqrt n} - 1$ is not expressible.
Instead we use the expression $\frac{x^2}n - 1$ as error-rate which
does not require any square-root operation. *}
termination
proof -
define er where "er x = (x * x / n - 1)" for x
define c where "c = 2 * n / ε"
define m where "m x = nat ⌊ c * er x ⌋" for x
have c: "c > 0" unfolding c_def using n ε by auto
show ?thesis
proof
show "wf (measures [m])" by simp
next
fix x
assume x: "0 < x" and xe: "¬ x * x - n < ε"
define y where "y = (n / x + x) / 2"
show "((n / x + x) / 2,x) ∈ measures [m]" unfolding y_def[symmetric]
proof (rule measures_less)
from n have inv_n: "1 / n > 0" by auto
from xe have "x * x - n ≥ ε" by simp
from this[unfolded mult_le_cancel_left_pos[OF inv_n, of ε, symmetric]]
have erxen: "er x ≥ ε / n" unfolding er_def using n by (simp add: field_simps)
have en: "ε / n > 0" and ne: "n / ε > 0" using ε n by auto
from en erxen have erx: "er x > 0" by linarith
have pos: "er x * 4 + er x * (er x * 4) > 0" using erx
by (auto intro: add_pos_nonneg)
have "er y = 1 / 4 * (n / (x * x) - 2 + x * x / n)" unfolding er_def y_def using x n
by (simp add: field_simps)
also have "… = 1 / 4 * er x * er x / (1 + er x)" unfolding er_def using x n
by (simp add: field_simps)
finally have "er y = 1 / 4 * er x * er x / (1 + er x)" .
also have "… < 1 / 4 * (1 + er x) * er x / (1 + er x)" using erx erx pos
by (auto simp: field_simps)
also have "… = er x / 4" using erx by (simp add: field_simps)
finally have er_y_x: "er y ≤ er x / 4" by linarith
from erxen have "c * er x ≥ 2" unfolding c_def mult_le_cancel_left_pos[OF ne, of _ "er x", symmetric]
using n ε by (auto simp: field_simps)
hence pos: "⌊c * er x⌋ > 0" "⌊c * er x⌋ ≥ 2" by auto
show "m y < m x" unfolding m_def nat_mono_iff[OF pos(1)]
proof -
have "⌊c * er y⌋ ≤ ⌊c * (er x / 4)⌋"
by (rule floor_mono, unfold mult_le_cancel_left_pos[OF c], rule er_y_x)
also have "… < ⌊c * er x / 4 + 1⌋" by auto
also have "… ≤ ⌊c * er x⌋"
by (rule floor_mono, insert pos(2), simp add: field_simps)
finally show "⌊c * er y⌋ < ⌊c * er x⌋" .
qed
qed
qed
qed
text {* Once termination is proven, it is easy to show equivalence of
@{const sqrt_approx_main_impl} and @{const sqrt_approx_main}. *}
lemma sqrt_approx_main_impl: "x > 0 ⟹ sqrt_approx_main_impl ε n x = sqrt_approx_main x"
proof (induct x rule: sqrt_approx_main.induct)
case (1 x)
hence x: "x > 0" by auto
hence nx: "0 < (n / x + x) / 2" using n by (auto intro: pos_add_strict)
note simps = sqrt_approx_main_impl.simps[of _ _ x] sqrt_approx_main.simps[of x]
show ?case
proof (cases "x * x - n < ε")
case True
thus ?thesis unfolding simps using x by auto
next
case False
show ?thesis using 1(1)[OF x False nx] unfolding simps using x False by auto
qed
qed
text {* Also soundness is not complicated. *}
lemma sqrt_approx_main_sound: assumes x: "x > 0" and xx: "x * x > n"
shows "sqrt_approx_main x * sqrt_approx_main x > n ∧ sqrt_approx_main x * sqrt_approx_main x - n < ε"
using assms
proof (induct x rule: sqrt_approx_main.induct)
case (1 x)
from 1 have x: "x > 0" "(x > 0) = True" by auto
note simp = sqrt_approx_main.simps[of x, unfolded x if_True]
show ?case
proof (cases "x * x - n < ε")
case True
with 1 show ?thesis unfolding simp by simp
next
case False
let ?y = "(n / x + x) / 2"
from False simp have simp: "sqrt_approx_main x = sqrt_approx_main ?y" by simp
from n x have y: "?y > 0" by (auto intro: pos_add_strict)
note IH = 1(1)[OF x(1) False y]
from x have x4: "4 * x * x > 0" by (auto intro: mult_sign_intros)
show ?thesis unfolding simp
proof (rule IH)
show "n < ?y * ?y"
unfolding mult_less_cancel_left_pos[OF x4, of n, symmetric]
proof -
have id: "4 * x * x * (?y * ?y) = 4 * x * x * n + (n - x * x) * (n - x * x)" using x(1)
by (simp add: field_simps)
from 1(3) have "x * x - n > 0" by auto
from mult_pos_pos[OF this this]
show "4 * x * x * n < 4 * x * x * (?y * ?y)" unfolding id
by (simp add: field_simps)
qed
qed
qed
qed
end
text {* It remains to assemble everything into one algorithm. *}
definition sqrt_approx :: "'a :: {linordered_field,floor_ceiling} ⇒ 'a ⇒ 'a" where
"sqrt_approx ε x ≡ if ε > 0 then (if x = 0 then 0 else let xpos = abs x in sqrt_approx_main_impl ε xpos (xpos + 1)) else 0"
lemma sqrt_approx: assumes ε: "ε > 0"
shows "¦sqrt_approx ε x * sqrt_approx ε x - ¦x¦¦ < ε"
proof (cases "x = 0")
case True
with ε show ?thesis unfolding sqrt_approx_def by auto
next
case False
let ?x = "¦x¦"
let ?sqrti = "sqrt_approx_main_impl ε ?x (?x + 1)"
let ?sqrt = "sqrt_approximation.sqrt_approx_main ε ?x (?x + 1)"
define sqrt where "sqrt = ?sqrt"
from False have x: "?x > 0" "?x + 1 > 0" by auto
interpret sqrt_approximation ε ?x
by (unfold_locales, insert x ε, auto)
from False ε have "sqrt_approx ε x = ?sqrti" unfolding sqrt_approx_def by (simp add: Let_def)
also have "?sqrti = ?sqrt"
by (rule sqrt_approx_main_impl, auto)
finally have id: "sqrt_approx ε x = sqrt" unfolding sqrt_def .
have sqrt: "sqrt * sqrt > ?x ∧ sqrt * sqrt - ?x < ε" unfolding sqrt_def
by (rule sqrt_approx_main_sound[OF x(2)], insert x mult_pos_pos[OF x(1) x(1)], auto simp: field_simps)
show ?thesis unfolding id using sqrt by auto
qed
subsection {* Some tests *}
text {* Testing executabity and show that sqrt 2 is irrational *}
lemma "¬ (∃ i :: rat. i * i = 2)"
proof -
have "set (sqrt_rat 2) = {}" by eval
thus ?thesis by simp
qed
text {* Testing speed *}
lemma "¬ (∃ i :: int. i * i = 1234567890123456789012345678901234567890)"
proof -
have "set (sqrt_int 1234567890123456789012345678901234567890) = {}" by eval
thus ?thesis by simp
qed
text {* The following test *}
value "let ε = 1 / 100000000 :: rat; s = sqrt_approx ε 2 in (s, s * s - 2, ¦s * s - 2¦ < ε)"
text {* results in (1.4142135623731116, 4.738200762148612e-14, True). *}
end