Theory Sqrt_Babylonian_Auxiliary

(*  Title:       Computing Square Roots using the Babylonian Method
Author:      René Thiemann       <rene.thiemann@uibk.ac.at>
Maintainer:  René Thiemann
*)

(*

This file is part of IsaFoR/CeTA.

IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the
Foundation, either version 3 of the License, or (at your option) any later
version.

IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along
with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>.
*)
section ‹Auxiliary lemmas which might be moved into the Isabelle distribution.›

theory Sqrt_Babylonian_Auxiliary
imports
Complex_Main
begin

lemma mod_div_equality_int: "(n :: int) div x * x = n - n mod x"
using div_mult_mod_eq[of n x] by arith

lemma div_is_floor_divide_rat: "n div y = ⌊rat_of_int n / rat_of_int y⌋"
unfolding Fract_of_int_quotient[symmetric] floor_Fract by simp

lemma div_is_floor_divide_real: "n div y = ⌊real_of_int n / of_int y⌋"
unfolding div_is_floor_divide_rat[of n y]
by (metis Ratreal_def of_rat_divide of_rat_of_int_eq real_floor_code)

lemma floor_div_pos_int:
fixes r :: "'a :: floor_ceiling"
assumes n: "n > 0"
shows "⌊r / of_int n⌋ = ⌊r⌋ div n" (is "?l = ?r")
proof -
let ?of_int = "of_int :: int ⇒ 'a"
define rhs where "rhs = ⌊r⌋ div n"
let ?n = "?of_int n"
define m where "m = ⌊r⌋ mod n"
let ?m = "?of_int m"
from div_mult_mod_eq[of "floor r" n] have dm: "rhs * n + m = ⌊r⌋" unfolding rhs_def m_def by simp
have mn: "m < n" and m0: "m ≥ 0" using n m_def by auto
define e where "e = r - ?of_int ⌊r⌋"
have e0: "e ≥ 0" unfolding e_def
by (metis diff_self eq_iff floor_diff_of_int zero_le_floor)
have e1: "e < 1" unfolding e_def
by (metis diff_self dual_order.refl floor_diff_of_int floor_le_zero)
have "r = ?of_int ⌊r⌋ + e" unfolding e_def by simp
also have "⌊r⌋ = rhs * n + m" using dm by simp
finally have "r = ?of_int (rhs * n + m) + e" .
hence "r / ?n = ?of_int (rhs * n) / ?n + (e + ?m) / ?n" using n by (simp add: field_simps)
also have "?of_int (rhs * n) / ?n = ?of_int rhs" using n by auto
finally have *: "r / ?of_int n = (e + ?of_int m) / ?of_int n + ?of_int rhs" by simp
have "?l = rhs + floor ((e + ?m) / ?n)" unfolding * by simp
also have "floor ((e + ?m) / ?n) = 0"
proof (rule floor_unique)
show "?of_int 0 ≤ (e + ?m) / ?n" using e0 m0 n
by (metis add_increasing2 divide_nonneg_pos of_int_0 of_int_0_le_iff of_int_0_less_iff)
show "(e + ?m) / ?n < ?of_int 0 + 1"
proof (rule ccontr)
from n have n': "?n > 0" "?n ≥ 0" by simp_all
assume "¬ ?thesis"
hence "(e + ?m) / ?n ≥ 1" by auto
from mult_right_mono[OF this n'(2)]
have "?n ≤ e + ?m" using n'(1) by simp
also have "?m ≤ ?n - 1" using mn
by (metis of_int_1 of_int_diff of_int_le_iff zle_diff1_eq)
finally have "?n ≤ e + ?n - 1" by auto
with e1 show False by arith
qed
qed
finally show ?thesis unfolding rhs_def by simp
qed

lemma floor_div_neg_int:
fixes r :: "'a :: floor_ceiling"
assumes n: "n < 0"
shows "⌊r / of_int n⌋ = ⌈r⌉ div n"
proof -
from n have n': "- n > 0" by auto
have "⌊r / of_int n⌋ = ⌊ - r / of_int (- n)⌋" using n
by (metis floor_of_int floor_zero less_int_code(1) minus_divide_left minus_minus nonzero_minus_divide_right of_int_minus)
also have "… = ⌊ - r ⌋ div (- n)" by (rule floor_div_pos_int[OF n'])
also have "… = ⌈ r ⌉ div n" using n
by (metis ceiling_def div_minus_right)
finally show ?thesis .
qed

lemma divide_less_floor1: "n / y < of_int (floor (n / y)) + 1"

context linordered_idom
begin

lemma sgn_int_pow_if [simp]:
"sgn x ^ p = (if even p then 1 else sgn x)" if "x ≠ 0"
using that by (induct p) simp_all

lemma compare_pow_le_iff: "p > 0 ⟹ (x :: 'a) ≥ 0 ⟹ y ≥ 0 ⟹ (x ^ p ≤ y ^ p) = (x ≤ y)"
by (rule power_mono_iff)

lemma compare_pow_less_iff: "p > 0 ⟹ (x :: 'a) ≥ 0 ⟹ y ≥ 0 ⟹ (x ^ p < y ^ p) = (x < y)"
using compare_pow_le_iff [of p x y]
using local.dual_order.order_iff_strict local.power_strict_mono by blast

end

lemma quotient_of_int[simp]: "quotient_of (of_int i) = (i,1)"
by (metis Rat.of_int_def quotient_of_int)

lemma quotient_of_nat[simp]: "quotient_of (of_nat i) = (int i,1)"
by (metis Rat.of_int_def Rat.quotient_of_int of_int_of_nat_eq)

lemma square_lesseq_square: "⋀ x y. 0 ≤ (x :: 'a :: linordered_field) ⟹ 0 ≤ y ⟹ (x * x ≤ y * y) = (x ≤ y)"
by (metis mult_mono mult_strict_mono' not_less)

lemma square_less_square: "⋀ x y. 0 ≤ (x :: 'a :: linordered_field) ⟹ 0 ≤ y ⟹ (x * x < y * y) = (x < y)"
by (metis mult_mono mult_strict_mono' not_less)

lemma sqrt_sqrt[simp]: "x ≥ 0 ⟹ sqrt x * sqrt x = x"
by (metis real_sqrt_pow2 power2_eq_square)

lemma abs_lesseq_square: "abs (x :: real) ≤ abs y ⟷ x * x ≤ y * y"
using square_lesseq_square[of "abs x" "abs y"] by auto

end


Theory Log_Impl

(*  Title:       Computing Square Roots using the Babylonian Method
Author:      René Thiemann       <rene.thiemann@uibk.ac.at>
Maintainer:  René Thiemann
*)
section ‹A Fast Logarithm Algorithm›

theory Log_Impl
imports
Sqrt_Babylonian_Auxiliary
begin

text ‹We implement the discrete logarithm function in a manner similar to
a repeated squaring exponentiation algorithm.›

text ‹In order to prove termination of the algorithm without intermediate checks
we need to ensure that we only use proper bases,
i.e., values of at least 2. This will be encoded into a separate type.›

typedef proper_base = "{x :: int. x ≥ 2}" by auto

setup_lifting type_definition_proper_base

lift_definition get_base :: "proper_base ⇒ int" is "λ x. x" .

lift_definition square_base :: "proper_base ⇒ proper_base" is "λ x. x * x"
proof -
fix i :: int
assume i: "2 ≤ i"
have "2 * 2 ≤ i * i"
by (rule mult_mono[OF i i], insert i, auto)
thus "2 ≤ i * i" by auto
qed

lift_definition into_base :: "int ⇒ proper_base" is "λ x. if x ≥ 2 then x else 2" by auto

lemma square_base: "get_base (square_base b) = get_base b * get_base b"
by (transfer, auto)

lemma get_base_2: "get_base b ≥ 2"
by (transfer, auto)

lemma b_less_square_base_b: "get_base b < get_base (square_base b)"
unfolding square_base using get_base_2[of b] by simp

lemma b_less_div_base_b: assumes xb: "¬ x < get_base b"
shows "x div get_base b < x"
proof -
from get_base_2[of b] have b: "get_base b ≥ 2" .
with xb have x2: "x ≥ 2" by auto
with b int_div_less_self[of x "(get_base b)"]
show ?thesis by auto
qed

text ‹We now state the main algorithm.›

function log_main :: "proper_base ⇒ int ⇒ nat × int" where
"log_main b x = (if x < get_base b then (0,1) else
case log_main (square_base b) x of
(z, bz) ⇒
let l = 2 * z; bz1 = bz * get_base b
in if x < bz1 then (l,bz) else (Suc l,bz1))"
by pat_completeness auto

termination by (relation "measure (λ (b,x). nat (1 + x - get_base b))",
insert b_less_square_base_b, auto)

lemma log_main: "x > 0 ⟹ log_main b x = (y,by) ⟹ by = (get_base b)^y ∧ (get_base b)^y ≤ x ∧ x < (get_base b)^(Suc y)"
proof (induct b x arbitrary: y "by" rule: log_main.induct)
case (1 b x y "by")
note x = 1(2)
note y = 1(3)
note IH = 1(1)
let ?b = "get_base b"
show ?case
proof (cases "x < ?b")
case True
with x y show ?thesis by auto
next
case False
obtain z bz where zz: "log_main (square_base b) x = (z,bz)"
by (cases "log_main (square_base b) x", auto)
have id: "get_base (square_base b) ^ k = ?b ^ (2 * k)" for k unfolding square_base
from IH[OF False x zz, unfolded id]
have z: "?b ^ (2 * z) ≤ x" "x < ?b ^ (2 * Suc z)" and bz: "bz = get_base b ^ (2 * z)" by auto
from y[unfolded log_main.simps[of b x] Let_def zz split] bz False
have yy: "(if x < bz * ?b then (2 * z, bz) else (Suc (2 * z), bz * ?b)) =
(y, by)" by auto
show ?thesis
proof (cases "x < bz * ?b")
case True
with yy have yz: "y = 2 * z" "by = bz" by auto
from True z(1) bz show ?thesis unfolding yz by (auto simp: ac_simps)
next
case False
with yy have yz: "y = Suc (2 * z)" "by = ?b * bz" by auto
from False have "?b ^ Suc (2 * z) ≤ x" by (auto simp: bz ac_simps)
with z(2) bz show ?thesis unfolding yz by auto
qed
qed
qed

text ‹We then derive the floor- and ceiling-log functions.›

definition log_floor :: "int ⇒ int ⇒ nat" where
"log_floor b x = fst (log_main (into_base b) x)"

definition log_ceiling :: "int ⇒ int ⇒ nat" where
"log_ceiling b x = (case log_main (into_base b) x of
(y,by) ⇒ if x = by then y else Suc y)"

lemma log_floor_sound: assumes "b > 1" "x > 0" "log_floor b x = y"
shows "b^y ≤ x" "x < b^(Suc y)"
proof -
from assms(1,3) have id: "get_base (into_base b) = b" by transfer auto
obtain yy bb where log: "log_main (into_base b) x = (yy,bb)"
by (cases "log_main (into_base b) x", auto)
from log_main[OF assms(2) log] assms(3)[unfolded log_floor_def log] id
show "b^y ≤ x" "x < b^(Suc y)" by auto
qed

lemma log_ceiling_sound: assumes "b > 1" "x > 0" "log_ceiling b x = y"
shows "x ≤ b^y" "y ≠ 0 ⟹ b^(y - 1) < x"
proof -
from assms(1,3) have id: "get_base (into_base b) = b" by transfer auto
obtain yy bb where log: "log_main (into_base b) x = (yy,bb)"
by (cases "log_main (into_base b) x", auto)
from log_main[OF assms(2) log, unfolded id] assms(3)[unfolded log_ceiling_def log split]
have bnd: "b ^ yy ≤ x" "x < b ^ Suc yy" and
y: "y = (if x = b ^ yy then yy else Suc yy)" by auto
have "x ≤ b^y ∧ (y ≠ 0 ⟶ b^(y - 1) < x)"
proof (cases "x = b ^ yy")
case True
with y bnd assms(1) show ?thesis by (cases yy, auto)
next
case False
with y bnd show ?thesis by auto
qed
thus "x ≤ b^y" "y ≠ 0 ⟹ b^(y - 1) < x" by auto
qed

text ‹Finally, we connect it to the @{const log} function working on real numbers.›

lemma log_floor[simp]: assumes b: "b > 1" and x: "x > 0"
shows "log_floor b x = ⌊log b x⌋"
proof -
obtain y where y: "log_floor b x = y" by auto
note main = log_floor_sound[OF assms y]
from b x have *: "1 < real_of_int b" "0 < real_of_int (b ^ y)" "0 < real_of_int x"
and **: "1 < real_of_int b" "0 < real_of_int x" "0 < real_of_int (b ^ Suc y)"
by auto
show ?thesis unfolding y
proof (rule sym, rule floor_unique)
show "real_of_int (int y) ≤ log (real_of_int b) (real_of_int x)"
using main(1)[folded log_le_cancel_iff[OF *, unfolded of_int_le_iff]]
using log_pow_cancel[of b y] b by auto
show "log (real_of_int b) (real_of_int x) < real_of_int (int y) + 1"
using main(2)[folded log_less_cancel_iff[OF **, unfolded of_int_less_iff]]
using log_pow_cancel[of b "Suc y"] b by auto
qed
qed

lemma log_ceiling[simp]: assumes b: "b > 1" and x: "x > 0"
shows "log_ceiling b x = ⌈log b x⌉"
proof -
obtain y where y: "log_ceiling b x = y" by auto
note main = log_ceiling_sound[OF assms y]
from b x have *: "1 < real_of_int b" "0 < real_of_int (b ^ (y - 1))" "0 < real_of_int x"
and **: "1 < real_of_int b" "0 < real_of_int x" "0 < real_of_int (b ^ y)"
by auto
show ?thesis unfolding y
proof (rule sym, rule ceiling_unique)
show "log (real_of_int b) (real_of_int x) ≤ real_of_int (int y)"
using main(1)[folded log_le_cancel_iff[OF **, unfolded of_int_le_iff]]
using log_pow_cancel[of b y] b by auto
from x have x: "x ≥ 1" by auto
show "real_of_int (int y) - 1 < log (real_of_int b) (real_of_int x)"
proof (cases "y = 0")
case False
thus ?thesis
using main(2)[folded log_less_cancel_iff[OF *, unfolded of_int_less_iff]]
using log_pow_cancel[of b "y - 1"] b x by auto
next
case True
have "real_of_int (int y) - 1 = log b (1/b)" using True b
by (subst log_divide, auto)
also have "… < log b 1"
by (subst log_less_cancel_iff, insert b, auto)
also have "… ≤ log b x"
by (subst log_le_cancel_iff, insert b x, auto)
finally show "real_of_int (int y) - 1 < log (real_of_int b) (real_of_int x)" .
qed
qed
qed

end


Theory NthRoot_Impl

(*  Title:       Computing Square Roots using the Babylonian Method
Author:      René Thiemann       <rene.thiemann@uibk.ac.at>
Maintainer:  René Thiemann
*)

(*

This file is part of IsaFoR/CeTA.

IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the
Foundation, either version 3 of the License, or (at your option) any later
version.

IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along
with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>.
*)
section ‹Executable algorithms for $p$-th roots›

theory NthRoot_Impl
imports
Log_Impl
Cauchy.CauchysMeanTheorem
begin

text ‹
We implemented algorithms to decide $\sqrt[p]{n} \in \rats$ and to compute $\lfloor \sqrt[p]{n} \rfloor$.
To this end, we use a variant of Newton iteration which works with integer division instead of floating
point or rational division. To get suitable starting values for the Newton iteration, we also implemented
a function to approximate logarithms.
›

subsection ‹Logarithm›

text ‹For computing the $p$-th root of a number $n$, we must choose a starting value
in the iteration. Here, we use @{term "2 ^ (nat ⌈of_int ⌈log 2 n⌉ / p⌉)"}.
›

text ‹We use a partial efficient algorithm, which does not terminate on
corner-cases, like $b = 0$ or $p = 1$, and invoke it properly afterwards.
Then there is a second algorithm which terminates on these corner-cases by additional
guards and on which we can perform induction.
›

subsection ‹Computing the $p$-th root of an integer number›

text ‹Using the logarithm, we can define an executable version of the
intended  starting value. Its main property is the inequality
@{term "(start_value x p) ^ p ≥ x"}, i.e., the start value is larger
than the p-th root. This property is essential, since our algorithm will abort
as soon as we fall below the p-th root.›

definition start_value :: "int ⇒ nat ⇒ int" where
"start_value n p = 2 ^ (nat ⌈of_nat (log_ceiling 2 n) / rat_of_nat p⌉)"

lemma start_value_main: assumes x: "x ≥ 0" and p: "p > 0"
shows "x ≤ (start_value x p)^p ∧ start_value x p ≥ 0"
proof (cases "x = 0")
case True
with p show ?thesis unfolding start_value_def True by simp
next
case False
with x have x: "x > 0" by auto
define l2x where "l2x = ⌈log 2 x⌉"
define pow where "pow = nat ⌈rat_of_int l2x / of_nat p⌉"
have "root p x = x powr (1 / p)" by (rule root_powr_inverse, insert x p, auto)
also have "… = (2 powr (log 2 x)) powr (1 / p)" using powr_log_cancel[of 2 x] x by auto
also have "… = 2 powr (log 2 x * (1 / p))" by (rule powr_powr)
also have "log 2 x * (1 / p) = log 2 x / p" using p by auto
finally have r: "root p x = 2 powr (log 2 x / p)" .
have lp: "log 2 x ≥ 0" using x by auto
hence l2pos: "l2x ≥ 0" by (auto simp: l2x_def)
have "log 2 x / p ≤ l2x / p" using x p unfolding l2x_def
by (metis divide_right_mono le_of_int_ceiling of_nat_0_le_iff)
also have "… ≤ ⌈l2x / (p :: real)⌉" by (simp add: ceiling_correct)
also have "l2x / real p = l2x / real_of_rat (of_nat p)"
by (metis of_rat_of_nat_eq)
also have "of_int l2x = real_of_rat (of_int l2x)"
by (metis of_rat_of_int_eq)
also have "real_of_rat (of_int l2x) / real_of_rat (of_nat p) = real_of_rat (rat_of_int l2x / of_nat p)"
by (metis of_rat_divide)
also have "⌈real_of_rat (rat_of_int l2x / rat_of_nat p)⌉ = ⌈rat_of_int l2x / of_nat p⌉"
by simp
also have "⌈rat_of_int l2x / of_nat p⌉ ≤ real pow" unfolding pow_def by auto
finally have le: "log 2 x / p ≤ pow" .
from powr_mono[OF le, of 2, folded r]
have "root p x ≤ 2 powr pow" by auto
also have "… = 2 ^ pow" by (rule powr_realpow, auto)
also have "… = of_int ((2 :: int) ^ pow)" by simp
also have "pow = (nat ⌈of_int (log_ceiling 2 x) / rat_of_nat p⌉)"
unfolding pow_def l2x_def using x by simp
also have "real_of_int ((2 :: int) ^ … ) = start_value x p" unfolding start_value_def  by simp
finally have less: "root p x ≤ start_value x p" .
have "0 ≤ root p x" using p x by auto
also have "… ≤ start_value x p" by (rule less)
finally have start: "0 ≤ start_value x p" by simp
from power_mono[OF less, of p] have "root p (of_int x) ^ p ≤ of_int (start_value x p) ^ p" using p x by auto
also have "… = start_value x p ^ p" by simp
also have "root p (of_int x) ^ p = x" using p x by force
finally have "x ≤ (start_value x p) ^ p" by presburger
with start show ?thesis by auto
qed

lemma start_value: assumes x: "x ≥ 0" and p: "p > 0" shows "x ≤ (start_value x p) ^ p" "start_value x p ≥ 0"
using start_value_main[OF x p] by auto

text ‹We now define the Newton iteration to compute the $p$-th root. We are working on the integers,
where every @{term "(/)"} is replaced by @{term "(div)"}. We are proving several things within
a locale which ensures that $p > 0$, and where $pm = p - 1$.
›

locale fixed_root =
fixes p pm :: nat
assumes p: "p = Suc pm"
begin

function root_newton_int_main :: "int ⇒ int ⇒ int × bool" where
"root_newton_int_main x n = (if (x < 0 ∨ n < 0) then (0,False) else (if x ^ p ≤ n then (x, x ^ p = n)
else root_newton_int_main ((n div (x ^ pm) + x * int pm) div (int p)) n))"
by pat_completeness auto
end

text ‹For the executable algorithm we omit the guard and use a let-construction›

partial_function (tailrec) root_int_main' :: "nat ⇒ int ⇒ int ⇒ int ⇒ int ⇒ int × bool" where
[code]: "root_int_main' pm ipm ip x n = (let xpm = x^pm; xp = xpm * x in if xp ≤ n then (x, xp = n)
else root_int_main' pm ipm ip ((n div xpm + x * ipm) div ip) n)"

text ‹In the following algorithm, we
start the iteration.
It will compute @{term "⌊root p n⌋"} and a boolean to indicate whether the root is exact.›

definition root_int_main :: "nat ⇒ int ⇒ int × bool" where
"root_int_main p n ≡ if p = 0 then (1,n = 1) else
let pm = p - 1
in root_int_main' pm (int pm) (int p) (start_value n p) n"

text ‹Once we have proven soundness of @{const fixed_root.root_newton_int_main} and equivalence
to @{const root_int_main}, it
is easy to assemble the following algorithm which computes all roots for arbitrary integers.›

definition root_int :: "nat ⇒ int ⇒ int list" where
"root_int p x ≡ if p = 0 then [] else
if x = 0 then [0] else
let e = even p; s = sgn x; x' = abs x
in if x < 0 ∧ e then [] else case root_int_main p x' of (y,True) ⇒ if e then [y,-y] else [s * y] | _ ⇒ []"

context fixed_root
begin
lemma iteration_mono_eq: assumes xn: "x ^ p = (n :: int)"
shows "(n div x ^ pm + x * int pm) div int p = x"
proof -
have [simp]: "⋀ n. (x + x * n) = x * (1 + n)" by (auto simp: field_simps)
show ?thesis unfolding xn[symmetric] p by simp
qed

lemma p0: "p ≠ 0" unfolding p by auto

text ‹The following property is the essential property for
proving termination of @{const "root_newton_int_main"}.
›
lemma iteration_mono_less: assumes x: "x ≥ 0"
and n: "n ≥ 0"
and xn: "x ^ p > (n :: int)"
shows "(n div x ^ pm + x * int pm) div int p < x"
proof -
let ?sx = "(n div x ^ pm + x * int pm) div int p"
from xn have xn_le: "x ^ p ≥ n" by auto
from xn x n have x0: "x > 0"
using not_le p by fastforce
from p have xp: "x ^ p = x * x ^ pm" by auto
from x n have "n div x ^ pm * x ^ pm ≤ n"
by (auto simp add: minus_mod_eq_div_mult [symmetric] mod_int_pos_iff not_less power_le_zero_eq)
also have "… ≤ x ^ p" using xn by auto
finally have le: "n div x ^ pm ≤ x" using x x0 unfolding xp by simp
have "?sx ≤ (x^p div x ^ pm + x * int pm) div int p"
by (rule zdiv_mono1, insert le p0, unfold xp, auto)
also have "x^p div x ^ pm = x" unfolding xp by auto
also have "x + x * int pm = x * int p" unfolding p by (auto simp: field_simps)
also have "x * int p div int p = x" using p by force
finally have le: "?sx ≤ x" .
{
assume "?sx = x"
from arg_cong[OF this, of "λ x. x * int p"]
have "x * int p ≤ (n div x ^ pm + x * int pm) div (int p) * int p" using p0 by simp
also have "… ≤ n div x ^ pm + x * int pm"
unfolding mod_div_equality_int using p by auto
finally have "n div x^pm ≥ x" by (auto simp: p field_simps)
from mult_right_mono[OF this, of "x ^ pm"]
have ge: "n div x^pm * x^pm ≥ x^p" unfolding xp using x by auto
from div_mult_mod_eq[of n "x^pm"] have "n div x^pm * x^pm = n - n mod x^pm" by arith
from ge[unfolded this]
have le: "x^p ≤ n - n mod x^pm" .
from x n have ge: "n mod x ^ pm ≥ 0"
by (auto simp add: mod_int_pos_iff not_less power_le_zero_eq)
from le ge
have "n ≥ x^p" by auto
with xn have False by auto
}
with le show ?thesis unfolding p by fastforce
qed

lemma iteration_mono_lesseq: assumes x: "x ≥ 0" and n: "n ≥ 0" and xn: "x ^ p ≥ (n :: int)"
shows "(n div x ^ pm + x * int pm) div int p ≤ x"
proof (cases "x ^ p = n")
case True
from iteration_mono_eq[OF this] show ?thesis by simp
next
case False
with assms have "x ^ p > n" by auto
from iteration_mono_less[OF x n this]
show ?thesis by simp
qed

termination (* of root_newton_int_main *)
proof -
let ?mm = "λ x  n :: int. nat x"
let ?m1 = "λ (x,n). ?mm x n"
let ?m = "measures [?m1]"
show ?thesis
proof (relation ?m)
fix x n :: int
assume "¬ x ^ p ≤ n"
hence x: "x ^ p > n" by auto
assume "¬ (x < 0 ∨ n < 0)"
hence x_n: "x ≥ 0" "n ≥ 0" by auto
from x x_n have x0: "x > 0" using p by (cases "x = 0", auto)
from iteration_mono_less[OF x_n x] x0
show "(((n div x ^ pm + x * int pm) div int p, n), x, n) ∈ ?m" by auto
qed auto
qed

text ‹We next prove that @{const root_int_main'} is a correct implementation of @{const root_newton_int_main}.
We additionally prove that the result is always positive, a lower bound, and that the returned boolean indicates
whether the result has a root or not. We prove all these results in one go, so that we can share the
inductive proof.
›

abbreviation root_main' where "root_main' ≡ root_int_main' pm (int pm) (int p)"

lemmas root_main'_simps = root_int_main'.simps[of pm "int pm" "int p"]

lemma root_main'_newton_pos: "x ≥ 0 ⟹ n ≥ 0 ⟹
root_main' x n = root_newton_int_main x n ∧ (root_main' x n = (y,b) ⟶ y ≥ 0 ∧ y^p ≤ n ∧ b = (y^p = n))"
proof (induct x n rule: root_newton_int_main.induct)
case (1 x n)
have pm_x[simp]: "x ^ pm * x = x ^ p" unfolding p by simp
from 1 have id: "(x < 0 ∨ n < 0) = False" by auto
note d = root_main'_simps[of x n] root_newton_int_main.simps[of x n] id if_False Let_def
show ?case
proof (cases "x ^ p ≤ n")
case True
thus ?thesis unfolding d using 1(2) by auto
next
case False
hence id: "(x ^ p ≤ n) = False" by simp
from 1(3) 1(2) have not: "¬ (x < 0 ∨ n < 0)" by auto
then have x: "x > 0 ∨ x = 0"
by auto
with ‹0 ≤ n› have "0 ≤ (n div x ^ pm + x * int pm) div int p"
by (auto simp add: p algebra_simps pos_imp_zdiv_nonneg_iff power_0_left)
then show ?thesis unfolding d id pm_x
by (rule 1(1)[OF not False _ 1(3)])
qed
qed

lemma root_main': "x ≥ 0 ⟹ n ≥ 0 ⟹ root_main' x n = root_newton_int_main x n"
using root_main'_newton_pos by blast

lemma root_main'_pos: "x ≥ 0 ⟹ n ≥ 0 ⟹ root_main' x n = (y,b) ⟹ y ≥ 0"
using root_main'_newton_pos by blast

lemma root_main'_sound: "x ≥ 0 ⟹ n ≥ 0 ⟹ root_main' x n = (y,b) ⟹ b = (y ^ p = n)"
using root_main'_newton_pos by blast

text ‹In order to prove completeness of the algorithms, we provide sharp upper and lower bounds
for @{const root_main'}. For the upper bounds, we use Cauchy's mean theorem where we added
the non-strict variant to Porter's formalization of this theorem.›

lemma root_main'_lower: "x ≥ 0 ⟹ n ≥ 0 ⟹ root_main' x n = (y,b) ⟹ y ^ p ≤ n"
using root_main'_newton_pos by blast

lemma root_newton_int_main_upper:
shows "y ^ p ≥ n ⟹ y ≥ 0 ⟹ n ≥ 0 ⟹ root_newton_int_main y n = (x,b) ⟹ n < (x + 1) ^ p"
proof (induct y n rule: root_newton_int_main.induct)
case (1 y n)
from 1(3) have y0: "y ≥ 0" .
then have "y > 0 ∨ y = 0"
by auto
from 1(4) have n0: "n ≥ 0" .
define y' where "y' = (n div (y ^ pm) + y * int pm) div (int p)"
from ‹y > 0 ∨ y = 0› ‹n ≥ 0› have y'0: "y' ≥ 0"
by (auto simp add: y'_def p algebra_simps pos_imp_zdiv_nonneg_iff power_0_left)
let ?rt = "root_newton_int_main"
from 1(5) have rt: "?rt y n = (x,b)" by auto
from y0 n0 have not: "¬ (y < 0 ∨ n < 0)" "(y < 0 ∨ n < 0) = False" by auto
note rt = rt[unfolded root_newton_int_main.simps[of y n] not(2) if_False, folded y'_def]
note IH = 1(1)[folded y'_def, OF not(1) _ _ y'0 n0]
show ?case
proof (cases "y ^ p ≤ n")
case False note yyn = this
with rt have rt: "?rt y' n = (x,b)" by simp
show ?thesis
proof (cases "n ≤ y' ^ p")
case True
show ?thesis
by (rule IH[OF False True rt])
next
case False
with rt have x: "x = y'" unfolding root_newton_int_main.simps[of y' n]
using n0 y'0 by simp
from yyn have yyn: "y^p > n" by simp
from False have yyn': "n > y' ^ p" by auto
{
assume pm: "pm = 0"
have y': "y' = n" unfolding y'_def p pm by simp
with yyn' have False unfolding p pm by auto
}
hence pm0: "pm > 0" by auto
show ?thesis
proof (cases "n = 0")
case True
thus ?thesis unfolding p
by (metis False y'0 zero_le_power)
next
case False note n00 = this
let ?y = "of_int y :: real"
let ?n = "of_int n :: real"
from yyn n0 have y00: "y ≠ 0" unfolding p by auto
from y00 y0 have y0: "?y > 0" by auto
from n0 False have n0: "?n > 0" by auto
define Y where "Y = ?y * of_int pm"
define NY where "NY = ?n / ?y ^ pm"
note pos_intro = divide_nonneg_pos add_nonneg_nonneg mult_nonneg_nonneg
have NY0: "NY > 0" unfolding NY_def using y0 n0
by (metis NY_def zero_less_divide_iff zero_less_power)
let ?ls = "NY # replicate pm ?y"
have prod: "∏:replicate pm ?y = ?y ^ pm "
by (induct pm, auto)
have sum: "∑:replicate pm ?y = Y" unfolding Y_def
by (induct pm, auto simp: field_simps)
have pos: "pos ?ls" unfolding pos_def using NY0 y0 by auto
have "root p ?n = gmean ?ls" unfolding gmean_def using y0
by (auto simp: p NY_def prod)
also have "… < mean ?ls"
proof (rule CauchysMeanTheorem_Less[OF pos het_gt_0I])
show "NY ∈ set ?ls" by simp
from pm0 show "?y ∈ set ?ls" by simp
have "NY < ?y"
proof -
from yyn have less: "?n < ?y ^ Suc pm" unfolding p
by (metis of_int_less_iff of_int_power)
have "NY < ?y ^ Suc pm / ?y ^ pm" unfolding NY_def
by (rule divide_strict_right_mono[OF less], insert y0, auto)
thus ?thesis using y0 by auto
qed
thus "NY ≠ ?y" by blast
qed
also have "… = (NY + Y) / real p"
by (simp add: mean_def sum p)
finally have *: "root p ?n < (NY + Y) / real p" .
have "?n = (root p ?n)^p" using n0
by (metis neq0_conv p0 real_root_pow_pos)
also have "… < ((NY + Y) / real p)^p"
by (rule power_strict_mono[OF *], insert n0 p, auto)
finally have ineq1: "?n < ((NY + Y) / real p)^p" by auto
{
define s where "s = n div y ^ pm + y * int pm"
define S where "S = NY + Y"
have Y0: "Y ≥ 0" using y0 unfolding Y_def
by (metis "1.prems"(2) mult_nonneg_nonneg of_int_0_le_iff of_nat_0_le_iff)
have S0: "S > 0" using NY0 Y0 unfolding S_def by auto
from p have p0: "p > 0" by auto
have "?n / ?y ^ pm  < of_int (floor (?n / ?y^pm)) + 1"
by (rule divide_less_floor1)
also have "floor (?n / ?y ^ pm) = n div y^pm"
unfolding div_is_floor_divide_real by (metis of_int_power)
finally have "NY < of_int (n div y ^ pm) + 1" unfolding NY_def by simp
hence less: "S < of_int s + 1" unfolding Y_def s_def S_def by simp
{ (* by sledgehammer *)
have f1: "∀x⇩0. rat_of_int ⌊rat_of_nat x⇩0⌋ = rat_of_nat x⇩0"
using of_int_of_nat_eq by simp
have f2: "∀x⇩0. real_of_int ⌊rat_of_nat x⇩0⌋ = real x⇩0"
using of_int_of_nat_eq by auto
have f3: "∀x⇩0 x⇩1. ⌊rat_of_int x⇩0 / rat_of_int x⇩1⌋ = ⌊real_of_int x⇩0 / real_of_int x⇩1⌋"
using div_is_floor_divide_rat div_is_floor_divide_real by simp
have f4: "0 < ⌊rat_of_nat p⌋"
using p by simp
have "⌊S⌋ ≤ s" using less floor_le_iff by auto
hence "⌊rat_of_int ⌊S⌋ / rat_of_nat p⌋ ≤ ⌊rat_of_int s / rat_of_nat p⌋"
using f1 f3 f4 by (metis div_is_floor_divide_real zdiv_mono1)
hence "⌊S / real p⌋ ≤ ⌊rat_of_int s / rat_of_nat p⌋"
using f1 f2 f3 f4 by (metis div_is_floor_divide_real floor_div_pos_int)
hence "S / real p ≤ real_of_int (s div int p) + 1"
using f1 f3 by (metis div_is_floor_divide_real floor_le_iff floor_of_nat less_eq_real_def)
}
hence "S / real p ≤ of_int(s div p) + 1" .
note this[unfolded S_def s_def]
}
hence ge: "of_int y' + 1 ≥ (NY + Y) / p" unfolding y'_def
by simp
have pos1: "(NY + Y) / p ≥ 0" unfolding Y_def NY_def
insert y0 n0 p0, auto)
have pos2: "of_int y' + (1 :: rat) ≥ 0" using y'0 by auto
have ineq2: "(of_int y' + 1) ^ p ≥ ((NY + Y) / p) ^ p"
by (rule power_mono[OF ge pos1])
from order.strict_trans2[OF ineq1 ineq2]
have "?n < of_int ((x + 1) ^ p)" unfolding x
thus "n < (x + 1) ^ p" using of_int_less_iff by blast
qed
qed
next
case True
with rt have x: "x = y" by simp
with 1(2) True have n: "n = y ^ p" by auto
show ?thesis unfolding n x using y0 unfolding p
qed
qed

lemma root_main'_upper:
"x ^ p ≥ n ⟹ x ≥ 0 ⟹ n ≥ 0 ⟹ root_main' x n = (y,b) ⟹ n < (y + 1) ^ p"
using root_newton_int_main_upper[of n x y b] root_main'[of x n] by auto
end

text ‹Now we can prove all the nice properties of @{const root_int_main}.›

lemma root_int_main_all: assumes n: "n ≥ 0"
and rm: "root_int_main p n = (y,b)"
shows "y ≥ 0 ∧ b = (y ^ p = n) ∧ (p > 0 ⟶ y ^ p ≤ n ∧ n < (y + 1)^p)
∧ (p > 0 ⟶ x ≥ 0 ⟶ x ^ p = n ⟶ y = x ∧ b)"
proof (cases "p = 0")
case True
with rm[unfolded root_int_main_def]
have y: "y = 1" and b: "b = (n = 1)" by auto
show ?thesis unfolding True y b using n by auto
next
case False
from False have p_0: "p > 0" by auto
from False have "(p = 0) = False" by simp
from rm[unfolded root_int_main_def this Let_def]
have rm: "root_int_main' (p - 1) (int (p - 1)) (int p) (start_value n p) n = (y,b)" by simp
from start_value[OF n p_0] have start: "n ≤ (start_value n p)^p" "0 ≤ start_value n p" by auto
interpret fixed_root p "p - 1"
by (unfold_locales, insert False, auto)
from root_main'_pos[OF start(2) n rm] have y: "y ≥ 0" .
from root_main'_sound[OF start(2) n rm] have b: "b = (y ^ p = n)" .
from root_main'_lower[OF start(2) n rm] have low: "y ^ p ≤ n" .
from root_main'_upper[OF start n rm] have up: "n < (y + 1) ^ p" .
{
assume n: "x ^ p = n" and x: "x ≥ 0"
with low up have low: "y ^ p ≤ x ^ p" and up: "x ^ p < (y+1) ^ p" by auto
from power_strict_mono[of x y, OF _ x p_0] low have x: "x ≥ y" by arith
from power_mono[of "(y + 1)" x p] y up have y: "y ≥ x" by arith
from x y have "x = y" by auto
with b n
have "y = x ∧ b" by auto
}
thus ?thesis using b low up y by auto
qed

lemma root_int_main: assumes n: "n ≥ 0"
and rm: "root_int_main p n = (y,b)"
shows "y ≥ 0" "b = (y ^ p = n)" "p > 0 ⟹ y ^ p ≤ n" "p > 0 ⟹ n < (y + 1)^p"
"p > 0 ⟹ x ≥ 0 ⟹ x ^ p = n ⟹ y = x ∧ b"
using root_int_main_all[OF n rm, of x] by blast+

lemma root_int[simp]: assumes p: "p ≠ 0 ∨ x ≠ 1"
shows "set (root_int p x) = {y . y ^ p = x}"
proof (cases "p = 0")
case True
with p have "x ≠ 1" by auto
thus ?thesis unfolding root_int_def True by auto
next
case False
hence p: "(p = 0) = False" and p0: "p > 0" by auto
note d = root_int_def p if_False Let_def
show ?thesis
proof (cases "x = 0")
case True
thus ?thesis unfolding d using p0 by auto
next
case False
hence x: "(x = 0) = False" by auto
show ?thesis
proof (cases "x < 0 ∧ even p")
case True
hence left: "set (root_int p x) = {}" unfolding d by auto
{
fix y
assume x: "y ^ p = x"
with True have "y ^ p < 0 ∧ even p" by auto
hence False by presburger
}
with left show ?thesis by auto
next
case False
with x p have cond: "(x = 0) = False" "(x < 0 ∧ even p) = False" by auto
obtain y b where rt: "root_int_main p ¦x¦ = (y,b)" by force
have "abs x ≥ 0" by auto
note rm = root_int_main[OF this rt]
have "?thesis =
(set (case root_int_main p ¦x¦ of (y, True) ⇒ if even p then [y, - y] else [sgn x * y] | (y, False) ⇒ []) =
{y. y ^ p = x})" unfolding d cond by blast
also have "(case root_int_main p ¦x¦ of (y, True) ⇒ if even p then [y, - y] else [sgn x * y] | (y, False) ⇒ [])
= (if b then if even p then [y, - y] else [sgn x * y] else [])" (is "_ = ?lhs")
unfolding rt by auto
also have "set ?lhs = {y. y ^ p = x}" (is "_ = ?rhs")
proof -
{
fix z
assume idx: "z ^ p = x"
hence eq: "(abs z) ^ p = abs x" by (metis power_abs)
from idx x p0 have z: "z ≠ 0" unfolding p by auto
have "(y, b) = (¦z¦, True)"
using rm(5)[OF p0 _ eq] by auto
hence id: "y = abs z" "b = True" by auto
have "z ∈ set ?lhs" unfolding id using z by (auto simp: idx[symmetric], cases "z < 0", auto)
}
moreover
{
fix z
assume z: "z ∈ set ?lhs"
hence b: "b = True" by (cases b, auto)
note z = z[unfolded b if_True]
from rm(2) b have yx: "y ^ p = ¦x¦" by auto
from rm(1) have y: "y ≥ 0" .
from False have "odd p ∨ even p ∧ x ≥ 0" by auto
hence "z ∈ ?rhs"
proof
assume odd: "odd p"
with z have "z = sgn x * y" by auto
hence "z ^ p = (sgn x * y) ^ p" by auto
also have "… = sgn x ^ p * y ^ p" unfolding power_mult_distrib by auto
also have "… = sgn x ^ p * abs x" unfolding yx by simp
also have "sgn x ^ p = sgn x" using x odd by auto
also have "sgn x * abs x = x" by (rule mult_sgn_abs)
finally show "z ∈ ?rhs" by auto
next
assume even: "even p ∧ x ≥ 0"
from z even have "z = y ∨ z = -y" by auto
hence id: "abs z = y" using y by auto
with yx x even have z: "z ≠ 0" using p0 by (cases "y = 0", auto)
have "z ^ p = (sgn z * abs z) ^ p" by (simp add: mult_sgn_abs)
also have "… = (sgn z * y) ^ p" using id by auto
also have "… = (sgn z)^p * y ^ p" unfolding power_mult_distrib  by simp
also have "… = sgn z ^ p * x" unfolding yx using even by auto
also have "sgn z ^ p = 1" using even z by (auto)
finally show "z ∈ ?rhs" by auto
qed
}
ultimately show ?thesis by blast
qed
finally show ?thesis by auto
qed
qed
qed

lemma root_int_pos: assumes x: "x ≥ 0" and ri: "root_int p x = y # ys"
shows "y ≥ 0"
proof -
from x have abs: "abs x = x" by auto
note ri = ri[unfolded root_int_def Let_def abs]
from ri have p: "(p = 0) = False" by (cases p, auto)
note ri = ri[unfolded p if_False]
show ?thesis
proof (cases "x = 0")
case True
with ri show ?thesis by auto
next
case False
hence "(x = 0) = False" "(x < 0 ∧ even p) = False" using x by auto
note ri = ri[unfolded this if_False]
obtain y' b' where r: "root_int_main p x = (y',b')" by force
note ri = ri[unfolded this]
hence y: "y = (if even p then y' else sgn x * y')" by (cases b', auto)
from root_int_main(1)[OF x r] have y': "0 ≤ y'" .
thus ?thesis unfolding y using x False by auto
qed
qed

subsection ‹Floor and ceiling of roots›

text ‹Using the bounds for @{const root_int_main} we can easily design
algorithms which compute @{term "floor (root p x)"} and @{term "ceiling (root p x)"}.
To this end, we first develop algorithms for non-negative @{term x}, and later on
these are used for the general case.›

definition "root_int_floor_pos p x = (if p = 0 then 0 else fst (root_int_main p x))"
definition "root_int_ceiling_pos p x = (if p = 0 then 0 else (case root_int_main p x of (y,b) ⇒ if b then y else y + 1))"

lemma root_int_floor_pos_lower: assumes p0: "p ≠ 0" and x: "x ≥ 0"
shows "root_int_floor_pos p x ^ p ≤ x"
using root_int_main(3)[OF x, of p] p0 unfolding root_int_floor_pos_def
by (cases "root_int_main p x", auto)

lemma root_int_floor_pos_pos: assumes x: "x ≥ 0"
shows "root_int_floor_pos p x ≥ 0"
using root_int_main(1)[OF x, of p]
unfolding root_int_floor_pos_def
by (cases "root_int_main p x", auto)

lemma root_int_floor_pos_upper: assumes p0: "p ≠ 0" and x: "x ≥ 0"
shows "(root_int_floor_pos p x + 1) ^ p > x"
using root_int_main(4)[OF x, of p] p0 unfolding root_int_floor_pos_def
by (cases "root_int_main p x", auto)

lemma root_int_floor_pos: assumes x: "x ≥ 0"
shows "root_int_floor_pos p x = floor (root p (of_int x))"
proof (cases "p = 0")
case True
thus ?thesis by (simp add: root_int_floor_pos_def)
next
case False
hence p: "p > 0" by auto
let ?s1 = "real_of_int (root_int_floor_pos p x)"
let ?s2 = "root p (of_int x)"
from x have s1: "?s1 ≥ 0"
by (metis of_int_0_le_iff root_int_floor_pos_pos)
from x have s2: "?s2 ≥ 0"
by (metis of_int_0_le_iff real_root_pos_pos_le)
from s1 have s11: "?s1 + 1 ≥ 0" by auto
have id: "?s2 ^ p = of_int x" using x
by (metis p of_int_0_le_iff real_root_pow_pos2)
show ?thesis
proof (rule floor_unique[symmetric])
show "?s1 ≤ ?s2"
unfolding compare_pow_le_iff[OF p s1 s2, symmetric]
unfolding id
using root_int_floor_pos_lower[OF False x]
by (metis of_int_le_iff of_int_power)
show "?s2 < ?s1 + 1"
unfolding compare_pow_less_iff[OF p s2 s11, symmetric]
unfolding id
using root_int_floor_pos_upper[OF False x]
by (metis of_int_add of_int_less_iff of_int_power of_int_1)
qed
qed

lemma root_int_ceiling_pos: assumes x: "x ≥ 0"
shows "root_int_ceiling_pos p x = ceiling (root p (of_int x))"
proof (cases "p = 0")
case True
thus ?thesis by (simp add: root_int_ceiling_pos_def)
next
case False
hence p: "p > 0" by auto
obtain y b where s: "root_int_main p x = (y,b)" by force
note rm = root_int_main[OF x s]
note rm = rm(1-2) rm(3-5)[OF p]
from rm(1) have y: "y ≥ 0" by simp
let ?s = "root_int_ceiling_pos p x"
let ?sx = "root p (of_int x)"
note d = root_int_ceiling_pos_def
show ?thesis
proof (cases b)
case True
hence id: "?s = y" unfolding s d using p by auto
from rm(2) True have xy: "x = y ^ p" by auto
show ?thesis unfolding id unfolding xy using y
next
case False
hence id: "?s = root_int_floor_pos p x + 1" unfolding d root_int_floor_pos_def
using s p by simp
from False have x0: "x ≠ 0" using rm(5)[of 0] using s unfolding root_int_main_def Let_def using p
by (cases "x = 0", auto)
show ?thesis unfolding id root_int_floor_pos[OF x]
proof (rule ceiling_unique[symmetric])
show "?sx ≤ real_of_int (⌊root p (of_int x)⌋ + 1)"
let ?l = "real_of_int (⌊root p (of_int x)⌋ + 1) - 1"
let ?m = "real_of_int ⌊root p (of_int x)⌋"
have "?l = ?m" by simp
also have "… < ?sx"
proof -
have le: "?m ≤ ?sx" by (rule of_int_floor_le)
have neq: "?m ≠ ?sx"
proof
assume "?m = ?sx"
hence "?m ^ p = ?sx ^ p" by auto
also have "… = of_int x" using x False
by (metis p real_root_ge_0_iff real_root_pow_pos2 root_int_floor_pos root_int_floor_pos_pos zero_le_floor zero_less_Suc)
finally have xs: "x = ⌊root p (of_int x)⌋ ^ p"
by (metis floor_power floor_of_int)
hence "⌊root p (of_int x)⌋ ∈ set (root_int p x)" using p by simp
hence "root_int p x ≠ []" by force
with s False ‹p ≠ 0› x x0 show False unfolding root_int_def
by (cases p, auto)
qed
from le neq show ?thesis by arith
qed
finally show "?l < ?sx" .
qed
qed
qed

definition "root_int_floor p x = (if x ≥ 0 then root_int_floor_pos p x else - root_int_ceiling_pos p (- x))"
definition "root_int_ceiling p x = (if x ≥ 0 then root_int_ceiling_pos p x else - root_int_floor_pos p (- x))"

lemma root_int_floor[simp]: "root_int_floor p x = floor (root p (of_int x))"
proof -
note d = root_int_floor_def
show ?thesis
proof (cases "x ≥ 0")
case True
with root_int_floor_pos[OF True, of p] show ?thesis unfolding d by simp
next
case False
hence "- x ≥ 0" by auto
from False root_int_ceiling_pos[OF this] show ?thesis unfolding d
qed
qed

lemma root_int_ceiling[simp]: "root_int_ceiling p x = ceiling (root p (of_int x))"
proof -
note d = root_int_ceiling_def
show ?thesis
proof (cases "x ≥ 0")
case True
with root_int_ceiling_pos[OF True] show ?thesis unfolding d by simp
next
case False
hence "- x ≥ 0" by auto
from False root_int_floor_pos[OF this, of p] show ?thesis unfolding d
qed
qed

subsection ‹Downgrading algorithms to the naturals›

definition root_nat_floor :: "nat ⇒ nat ⇒ int" where
"root_nat_floor p x = root_int_floor_pos p (int x)"

definition root_nat_ceiling :: "nat ⇒ nat ⇒ int" where
"root_nat_ceiling p x = root_int_ceiling_pos p (int x)"

definition root_nat :: "nat ⇒ nat ⇒ nat list" where
"root_nat p x = map nat (take 1 (root_int p x))"

lemma root_nat_floor [simp]: "root_nat_floor p x = floor (root p (real x))"
unfolding root_nat_floor_def using root_int_floor_pos[of "int x" p]
by auto

lemma root_nat_floor_lower: assumes p0: "p ≠ 0"
shows "root_nat_floor p x ^ p ≤ x"
using root_int_floor_pos_lower[OF p0, of x] unfolding root_nat_floor_def by auto

lemma root_nat_floor_upper: assumes p0: "p ≠ 0"
shows "(root_nat_floor p x + 1) ^ p > x"
using root_int_floor_pos_upper[OF p0, of x] unfolding root_nat_floor_def by auto

lemma root_nat_ceiling [simp]: "root_nat_ceiling p x = ceiling (root p x)"
unfolding root_nat_ceiling_def using root_int_ceiling_pos[of x p]
by auto

lemma root_nat: assumes p0: "p ≠ 0 ∨ x ≠ 1"
shows "set (root_nat p x) = { y. y ^ p = x}"
proof -
{
fix y
assume "y ∈ set (root_nat p x)"
note y = this[unfolded root_nat_def]
then obtain yi ys where ri: "root_int p x = yi # ys" by (cases "root_int p x", auto)
with y have y: "y = nat yi" by auto
from root_int_pos[OF _ ri] have yi: "0 ≤ yi" by auto
from root_int[of p "int x"] p0 ri have "yi ^ p = x" by auto
from arg_cong[OF this, of nat] yi have "nat yi ^ p = x"
by (metis nat_int nat_power_eq)
hence "y ∈ {y. y ^ p = x}" using y by auto
}
moreover
{
fix y
assume yx: "y ^ p = x"
hence y: "int y ^ p = int x"
by (metis of_nat_power)
hence "set (root_int p (int x)) ≠ {}" using root_int[of p "int x"] p0
by (metis (mono_tags) One_nat_def ‹y ^ p = x› empty_Collect_eq nat_power_eq_Suc_0_iff)
then obtain yi ys where ri: "root_int p (int x) = yi # ys"
by (cases "root_int p (int x)", auto)
from root_int_pos[OF _ this] have yip: "yi ≥ 0" by auto
from root_int[of p "int x", unfolded ri] p0 have yi: "yi ^ p = int x" by auto
with y have "int y ^ p = yi ^ p" by auto
from arg_cong[OF this, of nat] have id: "y ^ p = nat yi ^ p"
by (metis ‹y ^ p = x› nat_int nat_power_eq yi yip)
{
assume p: "p ≠ 0"
hence p0: "p > 0" by auto
obtain yy b where rm: "root_int_main p (int x) = (yy,b)" by force
from root_int_main(5)[OF _ rm p0 _ y] have "yy = int y" and "b = True" by auto
note rm = rm[unfolded this]
hence "y ∈ set (root_nat p x)"
unfolding root_nat_def p root_int_def using p0 p yx
by auto
}
moreover
{
assume p: "p = 0"
with p0 have "x ≠ 1" by auto
with y p have False by auto
}
ultimately have "y ∈ set (root_nat p x)" by auto
}
ultimately show ?thesis by blast
qed

subsection ‹Upgrading algorithms to the rationals›

text ‹The main observation to lift everything from the integers to the rationals is the fact, that one
can reformulate $\frac{a}{b}^{1/p}$ as $\frac{(ab^{p-1})^{1/p}}b$.›

definition root_rat_floor :: "nat ⇒ rat ⇒ int" where
"root_rat_floor p x ≡ case quotient_of x of (a,b) ⇒ root_int_floor p (a * b^(p - 1)) div b"

definition root_rat_ceiling :: "nat ⇒ rat ⇒ int" where
"root_rat_ceiling p x ≡ - (root_rat_floor p (-x))"

definition root_rat :: "nat ⇒ rat ⇒ rat list" where
"root_rat p x ≡ case quotient_of x of (a,b) ⇒ concat
(map (λ rb. map (λ ra. of_int ra / rat_of_int rb) (root_int p a)) (take 1 (root_int p b)))"

lemma root_rat_reform: assumes q: "quotient_of x = (a,b)"
shows "root p (real_of_rat x) = root p (of_int (a * b ^ (p - 1))) / of_int b"
proof (cases "p = 0")
case False
from quotient_of_denom_pos[OF q] have b: "0 < b" by auto
hence b: "0 < real_of_int b" by auto
from quotient_of_div[OF q] have x: "root p (real_of_rat x) = root p (a / b)"
by (metis of_rat_divide of_rat_of_int_eq)
also have "a / b = a * real_of_int b ^ (p - 1) / of_int b ^ p" using b False
by (cases p, auto simp: field_simps)
also have "root p … = root p (a * real_of_int b ^ (p - 1)) / root p (of_int b ^ p)" by (rule real_root_divide)
also have "root p (of_int b ^ p) = of_int b" using False b
by (metis neq0_conv real_root_pow_pos real_root_power)
also have "a * real_of_int b ^ (p - 1) = of_int (a * b ^ (p - 1))"
by (metis of_int_mult of_int_power)
finally show ?thesis .
qed auto

lemma root_rat_floor [simp]: "root_rat_floor p x = floor (root p (of_rat x))"
proof -
obtain a b where q: "quotient_of x = (a,b)" by force
from quotient_of_denom_pos[OF q] have b: "b > 0" .
show ?thesis
unfolding root_rat_floor_def q split root_int_floor
unfolding root_rat_reform[OF q] floor_div_pos_int[OF b] ..
qed

lemma root_rat_ceiling [simp]: "root_rat_ceiling p x = ceiling (root p (of_rat x))"
unfolding
root_rat_ceiling_def
ceiling_def
real_root_minus
root_rat_floor
of_rat_minus
..

lemma root_rat[simp]: assumes p: "p ≠ 0 ∨ x ≠ 1"
shows "set (root_rat p x) = { y. y ^ p = x}"
proof (cases "p = 0")
case False
note p = this
obtain a b where q: "quotient_of x = (a,b)" by force
note x = quotient_of_div[OF q]
have b: "b > 0" by (rule quotient_of_denom_pos[OF q])
note d = root_rat_def q split set_concat set_map
{
fix q
assume "q ∈ set (root_rat p x)"
note mem = this[unfolded d]
from mem obtain rb xs where rb: "root_int p b = Cons rb xs" by (cases "root_int p b", auto)
note mem = mem[unfolded this]
from mem obtain ra where ra: "ra ∈ set (root_int p a)" and q: "q = of_int ra / of_int rb"
by (cases "root_int p a", auto)
from rb have "rb ∈ set (root_int p b)" by auto
with ra p have rb: "b = rb ^ p" and ra: "a = ra ^ p" by auto
have "q ∈ {y. y ^ p = x}" unfolding q x ra rb
by (auto simp: power_divide)
}
moreover
{
fix q
assume "q ∈ {y. y ^ p = x}"
hence "q ^ p = of_int a / of_int b" unfolding x by auto
hence eq: "of_int b * q ^ p = of_int a" using b by auto
obtain z n where quo: "quotient_of q = (z,n)" by force
note qzn = quotient_of_div[OF quo]
have n: "n > 0" using quotient_of_denom_pos[OF quo] .
from eq[unfolded qzn] have "rat_of_int b * of_int z^p / of_int n^p = of_int a"
unfolding power_divide by simp
from arg_cong[OF this, of "λ x. x * of_int n^p"] n have "rat_of_int b * of_int z^p = of_int a * of_int n ^ p"
by auto
also have "rat_of_int b * of_int z^p = rat_of_int (b * z^p)" unfolding of_int_mult of_int_power ..
also have "of_int a * rat_of_int n ^ p = of_int (a * n ^ p)" unfolding of_int_mult of_int_power ..
finally have id: "a * n ^ p = b * z ^ p" by linarith
from quotient_of_coprime[OF quo] have cop: "coprime (z ^ p) (n ^ p)"
by simp
from coprime_crossproduct_int[OF quotient_of_coprime[OF q] this] arg_cong[OF id, of abs]
have "¦n ^ p¦ = ¦b¦"
with n b have bnp: "b = n ^ p" by auto
hence rn: "n ∈ set (root_int p b)" using p by auto
then obtain rb rs where rb: "root_int p b = Cons rb rs" by (cases "root_int p b", auto)
from id[folded bnp] b have "a = z ^ p" by auto
hence a: "z ∈ set (root_int p a)" using p by auto
from root_int_pos[OF _ rb] b have rb0: "rb ≥ 0" by auto
from root_int[OF disjI1[OF p], of b] rb have "rb ^ p = b" by auto
with bnp have id: "rb ^ p = n ^ p" by auto
have "rb = n" by (rule power_eq_imp_eq_base[OF id], insert n rb0 p, auto)
with rb have b: "n ∈ set (take 1 (root_int p b))" by auto
have "q ∈ set (root_rat p x)" unfolding d qzn using b a by auto
}
ultimately show ?thesis by blast
next
case True
with p have x: "x ≠ 1" by auto
obtain a b where q: "quotient_of x = (a,b)" by force
show ?thesis unfolding True root_rat_def q split root_int_def using x
by auto
qed

end


Theory Sqrt_Babylonian

(*  Title:       Computing Square Roots using the Babylonian Method
Author:      René Thiemann       <rene.thiemann@uibk.ac.at>
Maintainer:  René Thiemann
*)

(*

This file is part of IsaFoR/CeTA.

IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the
Foundation, either version 3 of the License, or (at your option) any later
version.

IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along
with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>.
*)

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"

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)"

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)"

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) ⌋"

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) ⌉"

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]

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
have "er y = 1 / 4 * (n / (x * x) - 2  + x * x / n)" unfolding er_def y_def using x n
also have "… = 1 / 4 * er x * er x / (1 + er x)" unfolding er_def using x n
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)
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
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
`