Theory NthRoot_Impl

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

(*
Copyright 2009-2014 René Thiemann

This file is part of IsaFoR/CeTA.

IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the
terms of the GNU Lesser General Public License as published by the Free Software
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] | _ ⇒ []"

text {* We start with proving termination of @{const fixed_root.root_newton_int_main}. *}

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
          by (intro divide_nonneg_pos add_nonneg_nonneg mult_nonneg_nonneg,
          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
          by (metis of_int_1 of_int_add of_int_power)
        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
      by (metis add_le_less_mono add_less_cancel_left lessI less_add_one add.right_neutral le_iff_add power_strict_mono)
  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
      by (simp add: p real_root_power_cancel)
  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)"
        by (metis of_int_add real_of_int_floor_add_one_ge of_int_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
      by (simp add: real_root_minus ceiling_minus)
  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
      by (simp add: real_root_minus floor_minus)
  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¦"
      by (simp add: field_simps abs_mult)
    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