Theory Complex_Roots

section ‹$n$-th roots of complex numbers›

theory Complex_Roots
  imports 
    Complex_Geometry.More_Complex
    Algebraic_Numbers.Complex_Algebraic_Numbers
    Factor_Algebraic_Polynomial.Roots_via_IA
    "HOL-Library.Product_Lexorder"
begin

subsection ‹An algorithm to compute all complex roots of (algebraic) complex numbers›

definition all_croots :: "nat ⇒ complex ⇒ complex list" where
  "all_croots n x = (if n = 0 then [] else
     if algebraic x then
       (let p = min_int_poly x;
         q = poly_nth_root n p;
         xs = complex_roots_of_int_poly q
         in filter (λ y. y^n = x) xs)
     else (SOME ys. set ys = {y. y^n = x}))" 

lemma all_croots: assumes n0: "n ≠ 0" shows "set (all_croots n x) = {y. y^n = x}" 
proof (cases "algebraic x")
  case True
  hence id: "(if n = 0 then y else if algebraic x then z else u) = z" 
    for y z u :: "complex list" using n0 by auto
  define p where "p = poly_nth_root n (min_int_poly x)" 
  show ?thesis unfolding Let_def p_def[symmetric] all_croots_def id
  proof (standard, force, standard, simp)
    fix y
    assume y: "y ^n  = x"     
    have "min_int_poly x represents x" using True by auto
    from represents_nth_root[OF n0 y this]
    have "p represents y" unfolding p_def by auto
    thus "y ∈ set (complex_roots_of_int_poly p)" 
      by (subst complex_roots_of_int_poly, auto)
  qed
next
  case False
  hence id: "(if n = 0 then y else if algebraic x then z else u) = u" 
    for y z u :: "complex list" using n0 by auto
  show ?thesis unfolding Let_def all_croots_def id
    by (rule someI_ex, rule finite_list, insert n0, blast)
qed

text ‹TODO: One might change @{const complex_roots_of_int_poly} to @{const complex_roots_of_int_poly3}
  in order to avoid an unnecessary factorization of an integer polynomial. However, then
  this change already needs to be performed within the definition of @{const all_croots}.› 
lift_definition all_croots_part1 :: "nat ⇒ complex ⇒ complex genuine_roots_aux" is
  "λ n x. if n = 0 ∨ x = 0 ∨ ¬ algebraic x then (1,[],0, filter_fun_complex 1) 
         else let p = min_int_poly x;
              q = poly_nth_root n p;
              zeros = complex_roots_of_int_poly q;
              r = Polynomial.monom 1 n - [:x:]
         in (r,zeros, n, filter_fun_complex r)"
  subgoal for n x
  proof (cases "n = 0 ∨ x = 0 ∨ ¬ algebraic x")
    case True
    thus ?thesis by (simp add: filter_fun_complex)
  next
    case False
    hence *: "algebraic x" "n ≠ 0" "x ≠ 0" by auto
    {
      fix z
      assume zn: "z^n = x" 
      from *(1) have repr: "min_int_poly x represents x" by auto
      from represents_nth_root[OF *(2) zn repr]
      have "poly_nth_root n (min_int_poly x) represents z" .
    }
    moreover have "card {z. z ^ n = x} = n"
      by (rule card_nth_roots) (use * in auto)
    ultimately show ?thesis using * 
      by (auto simp: Let_def complex_roots_of_int_poly filter_fun_complex poly_monom)
  qed
  done

lemma all_croots_code[code]: 
  "all_croots n x = (if n = 0 then [] else if x = 0 then [0]
     else if algebraic x then genuine_roots_impl (all_croots_part1 n x)
     else Code.abort (STR ''all_croots invoked on non-algebraic number'') (λ _. all_croots n x))"
proof (cases "n = 0")
  case True
  thus ?thesis unfolding all_croots_def by simp
next
  case n: False
  show ?thesis
  proof (cases "x = 0")
    case x: False
    show ?thesis
    proof (cases "algebraic x")
      case False
      with n x show ?thesis by simp
    next
      case True
      define t where "t = ?thesis" 
      have "t ⟷ filter (λy. y ^ n = x)
                (complex_roots_of_int_poly (poly_nth_root n (min_int_poly x)))
            = genuine_roots_impl (all_croots_part1 n x)" 
        unfolding t_def
        by (subst all_croots_def[of n x], unfold Let_def, insert n x True, auto)
      also have "…" using n x True unfolding genuine_roots_impl_def
        by (transfer, simp add: Let_def genuine_roots_def poly_monom)
      finally show ?thesis unfolding t_def by simp
    qed
  next
    case x: True
    have "set (all_croots n 0) = {0}" unfolding all_croots[OF n] using n by simp
    moreover have "distinct (all_croots n 0)" unfolding all_croots_def using n
      by (auto intro!: distinct_filter complex_roots_of_int_poly)
    ultimately have "all_croots n 0 = [0]"
      by (smt (verit, del_insts) distinct.simps(2) distinct_singleton insert_ident list.set_cases list.set_intros(1) list.simps(15) mem_Collect_eq set_empty singleton_conv)
    moreover have "?thesis ⟷ all_croots n 0 = [0]" using n x by simp
    ultimately show ?thesis by auto
  qed
qed


subsection ‹A definition of \emph{the} complex root of a complex number›

text ‹While the definition of the complex root is quite natural and easy,
  the main task is a criterion to determine which of all possible roots of a
  complex number is the chosen one.›

definition croot :: "nat ⇒ complex ⇒ complex" where
  "croot n x = (rcis (root n (cmod x)) (Arg x / of_nat n))" 

lemma croot_0[simp]: "croot n 0 = 0" "croot 0 x = 0"
  unfolding croot_def by auto

lemma croot_power: assumes n: "n ≠ 0" 
  shows "(croot n x) ^ n = x" 
  unfolding croot_def DeMoivre2
  by (subst real_root_pow_pos2, insert n, auto simp: rcis_cmod_Arg)

lemma Arg_of_real: "Arg (of_real x) = 
  (if x < 0 then pi else 0)"
proof (cases "x = 0")
  case False
  hence "x < 0 ∨ x > 0" by auto
  thus ?thesis by (intro cis_Arg_unique, auto
      simp: complex_sgn_def scaleR_complex.ctr complex_eq_iff)
qed (auto simp: Arg_def)


lemma Arg_rcis_cis[simp]: assumes "x > 0" 
  shows "Arg (rcis x y) = Arg (cis y)" 
  using assms unfolding rcis_def by simp

lemma cis_Arg_1[simp]: "cis (Arg 1) = 1" 
  using Arg_of_real[of 1] by simp

lemma cis_Arg_power[simp]: assumes "x ≠ 0" 
  shows "cis (Arg (x ^ n)) = cis (Arg x * real n)" 
proof (induct n)
  case (Suc n)
  show ?case unfolding power.simps
  proof (subst cis_arg_mult)    
    show "cis (Arg x + Arg (x ^ n)) = cis (Arg x * real (Suc n))" 
      unfolding mult.commute[of "Arg x"] DeMoivre[symmetric]
      unfolding power.simps using Suc
      by (metis DeMoivre cis_mult mult.commute)
    show "x * x ^ n ≠ 0" using assms by auto
  qed
qed simp

lemma Arg_croot[simp]: "Arg (croot n x) = Arg x / real n" 
proof (cases "n = 0 ∨ x = 0")
  case True
  thus ?thesis by (auto simp: Arg_def)
next
  case False
  hence n: "n ≠ 0" and x: "x ≠ 0" by auto
  let ?root = "croot n x" 
  from n have n1: "real n ≥ 1" "real n > 0" "real n ≠ 0" by auto
  have bounded: "- pi < Arg x / real n ∧ Arg x / real n ≤ pi" 
  proof (cases "Arg x < 0")
    case True
    from Arg_bounded[of x] have "- pi < Arg x" by auto
    also have "… ≤ Arg x / real n" using n1 True
      by (smt (verit) div_by_1 divide_minus_left frac_le)
    finally have one: "- pi < Arg x / real n" .
    have "Arg x / real n ≤ 0" using True n1 
      by (smt (verit) divide_less_0_iff)
    also have "… ≤ pi" by simp
    finally show ?thesis using one by auto
  next
    case False
    hence ax: "Arg x ≥ 0" by auto
    have "Arg x / real n ≤ Arg x" using n1 ax
      by (smt (verit) div_by_1 frac_le)
    also have "… ≤ pi" using Arg_bounded[of x] by simp
    finally have one: "Arg x / real n ≤ pi" .
    have "-pi < 0" by simp
    also have "… ≤ Arg x / real n" using ax n1 by simp
    finally show ?thesis using one by auto
  qed
  have "Arg ?root = Arg (cis (Arg x / real n))" 
    unfolding croot_def using x n by simp 
  also have "… = Arg x / real n" 
    by (rule cis_Arg_unique, force, insert bounded, auto)
  finally show ?thesis .
qed

lemma cos_abs[simp]: "cos (abs x :: real) = cos x" 
proof (cases "x < 0")
  case True
  hence abs: "abs x = - x" by simp
  show ?thesis unfolding abs by simp
qed simp

lemma cos_mono_le: assumes "abs x ≤ pi" 
  and "abs y ≤ pi" 
shows "cos x ≤ cos y ⟷ abs y ≤ abs x" 
proof -
  have "cos x ≤ cos y ⟷ cos (abs x) ≤ cos (abs y)" by simp
  also have "… ⟷ abs y ≤ abs x" 
    by (subst cos_mono_le_eq, insert assms, auto)
  finally show ?thesis .
qed

lemma abs_add_2_mult_bound: fixes x :: "'a :: linordered_idom" 
  assumes xy: "¦x¦ ≤ y" 
  shows "¦x¦ ≤ ¦x + 2 * of_int i * y¦" 
proof (cases "i = 0")
  case i: False
  let ?oi = "of_int :: int ⇒ 'a" 
  from xy have y: "y ≥ 0" by auto
  consider (pp) "x ≥ 0" "i ≥ 0" 
    | (nn) "x ≤ 0" "i ≤ 0" 
    | (pn) "x ≥ 0" "i ≤ 0" 
    | (np) "x ≤ 0" "i ≥ 0" 
    by linarith
  thus ?thesis
  proof cases
    case pp
    thus ?thesis using y by simp
  next
    case nn 
    have "x ≥ x + 2 * ?oi i * y" 
      using nn y by (simp add: mult_nonneg_nonpos2)
    with nn show ?thesis by linarith
  next
    case pn
    with i have "0 ≤ x" "i < 0" by auto
    define j where "j = nat (-i) - 1" 
    define z where "z = x - 2 * y" 
    define u where "u = 2 * ?oi (nat j) * y" 
    have u: "u ≥ 0" unfolding u_def using y by auto 
    have i: "i = - int (Suc j)" 
      using ‹i < 0› unfolding j_def by simp
    have id: "x + 2 * ?oi i * y = z - u" 
      unfolding i z_def u_def by (simp add: field_simps)
    have z: "z ≤ 0" "abs z ≥ x" using xy y pn(1)
      unfolding z_def by auto
    show ?thesis unfolding id using pn(1) z u by simp
  next
    case np
    with i have "0 ≥ x" "i > 0" by auto
    define j where "j = nat i - 1" 
    have i: "i = int (Suc j)" 
      using ‹i > 0› unfolding j_def by simp
    define u where "u = 2 * ?oi (nat j) * y" 
    have u: "u ≥ 0" unfolding u_def using y by auto 
    define z where "z = - x - 2 * y" 
    have id: "x + 2 * ?oi i * y = - z + u" 
      unfolding i z_def u_def by (simp add: field_simps)
    have z: "z ≤ 0" "abs z ≥ - x" using xy y np(1)
      unfolding z_def by auto
    show ?thesis unfolding id using np(1) z u by simp
  qed
qed simp

lemma abs_eq_add_2_mult: fixes y :: "'a :: linordered_idom" 
  assumes abs_id: "¦x¦ = ¦x + 2 * of_int i * y¦" 
  and xy: "- y < x" "x ≤ y" 
  and i: "i ≠ 0" 
shows "x = y ∧ i = -1" 
proof -  
  let ?oi = "of_int :: int ⇒ 'a" 
  from xy have y: "y > 0" by auto
  consider (pp) "x ≥ 0" "i ≥ 0" 
    | (nn) "x < 0" "i ≤ 0" 
    | (pn) "x ≥ 0" "i ≤ 0" 
    | (np) "x < 0" "i ≥ 0" 
    by linarith
  hence "?thesis ∨ x = ?oi (- i) * y" 
  proof cases
    case pp
    thus ?thesis using y abs_id xy i by simp
  next
    case nn
    hence "¦x + 2 * ?oi i * y¦ = 
      - (x + 2 * ?oi i * y)"
      using y nn
      by (intro abs_of_nonpos add_nonpos_nonpos, 
          force, simp, intro mult_nonneg_nonpos, auto)
    thus ?thesis using y abs_id xy i nn
      by auto
  next
    case pn
    with i have "0 ≤ x" "i < 0" by auto
    define j where "j = nat (-i) - 1" 
    define z where "z = x - 2 * y" 
    define u where "u = 2 * ?oi (nat j) * y" 
    have u: "u ≥ 0" unfolding u_def using y by auto 
    have i: "i = - int (Suc j)" 
      using ‹i < 0› unfolding j_def by simp
    have id: "x + 2 * ?oi i * y = z - u" 
      unfolding i z_def u_def by (simp add: field_simps)
    have z: "z ≤ 0" "abs z ≥ x" using xy y pn(1)
      unfolding z_def by auto
    from abs_id[unfolded id] have "z - u = -x " 
      using z u pn by auto
    from this[folded id] have "x = of_int (-i) * y" 
      by auto
    thus ?thesis by auto
  next
    case np
    with i have "0 ≥ x" "i > 0" by auto
    define j where "j = nat i - 1" 
    have i: "i = int (Suc j)" 
      using ‹i > 0› unfolding j_def by simp
    define u where "u = 2 * ?oi (nat j) * y" 
    have u: "u ≥ 0" unfolding u_def using y by auto 
    define z where "z = - x - 2 * y" 
    have id: "x + 2 * ?oi i * y = - z + u" 
      unfolding i z_def u_def by (simp add: field_simps)
    have z: "z ≤ 0" 
      using xy y np(1) unfolding z_def by auto
    from abs_id[unfolded id] have "- z + u = - x" 
      using u z np by auto
    from this[folded id] have "x = of_int (- i) * y" 
      by auto
    thus ?thesis by auto
  qed
  thus ?thesis
  proof
    assume "x = ?oi (- i) * y" 
    with xy i y
    show ?thesis
      by (smt (verit, ccfv_SIG) less_le minus_less_iff mult_le_cancel_right2 mult_minus1_right mult_minus_left mult_of_int_commute of_int_hom.hom_one of_int_le_1_iff of_int_minus)
  qed
qed

text ‹This is the core lemma. It tells us that @{const croot} will choose the 
  principal root, i.e. the root with
  largest real part and if there are two roots with identical real part, then the largest imaginary part.
  This criterion will be crucial for implementing @{const croot}.›

lemma croot_principal: assumes n: "n ≠ 0"
  and y: "y ^ n = x"
  and neq: "y ≠ croot n x"
shows "Re y < Re (croot n x) ∨ Re y = Re (croot n x) ∧ Im y < Im (croot n x)" 
proof (cases "x = 0")
  case True
  with neq y have False by auto
  thus ?thesis ..
next
  case x: False
  let ?root = "croot n x" 
  from n have n1: "real n ≥ 1" "real n > 0" "real n ≠ 0" by auto
  from x y n have y0: "y ≠ 0" by auto
  from croot_power[OF n, of x] y
  have id: "?root ^ n = y ^ n" by simp
  hence "cmod (?root ^ n) = cmod (y ^ n)" by simp
  hence norm_eq: "cmod ?root = cmod y" using n unfolding norm_power
    by (meson gr_zeroI norm_ge_zero power_eq_imp_eq_base)
  have "cis (Arg y * real n) = cis (Arg (y^n))" by (subst cis_Arg_power[OF y0], simp) 
  also have "… = cis (Arg x)" using y by simp
  finally have ciseq: "cis (Arg y * real n) = cis (Arg x)" by simp
  from cis_eq[OF ciseq] obtain i where
    "Arg y * real n - Arg x = 2 * real_of_int i * pi" 
    by auto
  hence "Arg y * real n = Arg x + 2 * real_of_int i * pi" by auto
  from arg_cong[OF this, of "λ x. x / real n"] n1
  have Argy: "Arg y = Arg ?root + 2 * real_of_int i * pi / real n" 
    by (auto simp: field_simps)
  have i0: "i ≠ 0" 
  proof
    assume "i = 0" 
    hence "Arg y = Arg ?root" unfolding Argy by simp
    with norm_eq have "?root = y" by (metis rcis_cmod_Arg)
    with neq show False by simp
  qed
  from y0 have cy0: "cmod y > 0" by auto
  from Arg_bounded[of x] have abs_pi: "abs (Arg x) ≤ pi" by auto
  have "Re y ≤ Re ?root ⟷ Re y / cmod y ≤ Re ?root / cmod y"
    using cy0 unfolding divide_le_cancel by simp
  also have cosy: "Re y / cmod y = cos (Arg y)" unfolding cos_arg[OF y0] ..
  also have cosrt: "Re ?root / cmod y = cos (Arg ?root)" 
    unfolding norm_eq[symmetric] by (subst cos_arg, insert norm_eq cy0, auto)
  also have "cos (Arg y) ≤ cos (Arg ?root) ⟷ abs (Arg ?root) ≤ abs (Arg y)" 
    by (rule cos_mono_le, insert Arg_bounded[of y] Arg_bounded[of ?root], auto)
  also have "… ⟷ abs (Arg ?root) * real n ≤ abs (Arg y) * real n" 
    unfolding mult_le_cancel_right using n1 by simp
  also have "… ⟷ abs (Arg x) ≤ ¦Arg x + 2 * real_of_int i * pi¦" 
    unfolding Argy using n1 by (simp add: field_simps)
  also have "…" using abs_pi
    by (rule abs_add_2_mult_bound)
  finally have le: "Re y ≤ Re (croot n x)" .
  show ?thesis
  proof (cases "Re y = Re (croot n x)")
    case False
    with le show ?thesis by auto
  next
    case True
    hence "Re y / cmod y = Re ?root / cmod y" by simp
    hence "cos (Arg y) = cos (Arg ?root)" unfolding cosy cosrt .
    hence "cos (abs (Arg y)) = cos (abs (Arg ?root))" unfolding cos_abs .
    from cos_inj_pi[OF _ _ _ _ this]
    have "abs (Arg y) = abs (Arg ?root)" 
      using Arg_bounded[of y] Arg_bounded[of ?root] by auto
    hence "abs (Arg y) * real n = abs (Arg ?root) * real n" by simp
    hence "abs (Arg x) = ¦Arg x + 2 * real_of_int i * pi¦" unfolding Argy 
      using n1 by (simp add: field_simps)
    from abs_eq_add_2_mult[OF this _ _ ‹i ≠ 0›] Arg_bounded[of x]  
    have Argx: "Arg x = pi" and i: "i = -1" by auto
    have Argy: "Arg y = -pi / real n" 
      unfolding Argy Arg_croot i Argx by simp
    have "Im ?root > Im y ⟷ Im ?root / cmod ?root > Im y / cmod y"
      unfolding norm_eq using cy0
      by (meson divide_less_cancel divide_strict_right_mono) 
    also have "… ⟷ sin (Arg ?root) > sin (Arg y)" 
      by (subst (1 2) sin_arg, insert y0 norm_eq, auto)
    also have "… ⟷ sin (- pi / real n) < sin (pi / real n)" 
      unfolding Argy Arg_croot Argx by simp
    also have …
    proof -
      have "sin (- pi / real n) < 0" 
        using n1 by (smt (verit) Arg_bounded Argy divide_neg_pos sin_gt_zero sin_minus)
      also have "… < sin (pi / real n)" 
        using n1 calculation by fastforce
      finally show ?thesis .
    qed
    finally show ?thesis using le by auto
  qed
qed

lemma croot_unique: assumes n: "n ≠ 0"
  and y: "y ^ n = x"
  and y_max_Re_Im: "⋀ z. z ^ n = x ⟹
      Re z < Re y ∨ Re z = Re y ∧ Im z ≤ Im y" 
shows "croot n x = y" 
proof (rule ccontr)
  assume "croot n x ≠ y" 
  from croot_principal[OF n y this[symmetric]]
  have "Re y < Re (croot n x) ∨
    Re y = Re (croot n x) ∧ Im y < Im (croot n x)" .
  with y_max_Re_Im[OF croot_power[OF n]]
  show False by auto
qed

lemma csqrt_is_croot_2: "csqrt = croot 2" 
proof
  fix x
  show "csqrt x = croot 2 x" 
  proof (rule sym, rule croot_unique, force, force)
    let ?p = "[:-x,0,1:]" 
    let ?cx = "csqrt x" 
    have p: "?p = [:?cx,1:] * [:-?cx,1:]" 
      by (simp add: power2_eq_square[symmetric])
    fix y
    assume "y^2 = x" 
    hence "True ⟷ poly ?p y = 0" 
      by (auto simp: power2_eq_square)
    also have "… ⟷ y = - ?cx ∨ y = ?cx" 
      unfolding p poly_mult mult_eq_0_iff poly_root_factor by auto
    finally have "y = - ?cx ∨ y = ?cx" by simp
    thus "Re y < Re ?cx ∨ Re y = Re ?cx ∧ Im y ≤ Im ?cx"
    proof
      assume y: "y = - ?cx"
      show ?thesis
      proof (cases "Re ?cx = 0")
        case False
        with csqrt_principal[of x] have "Re ?cx > 0" by simp
        thus ?thesis unfolding y by simp
      next
        case True
        with csqrt_principal[of x] have "Im ?cx ≥ 0" by simp
        thus ?thesis unfolding y using True by auto
      qed
    qed auto
  qed
qed

lemma croot_via_root_selection: assumes roots: "set ys = { y. y^n = x}" 
  and n: "n ≠ 0" 
shows "croot n x = arg_min_list (λ y. (- Re y, - Im y)) ys" 
  (is "_ = arg_min_list ?f ys")
proof (rule croot_unique[OF n])
  let ?y = "arg_min_list ?f ys" 
  have rt: "croot n x ^ n = x" using n by (rule croot_power)
  hence "croot n x ∈ set ys" unfolding roots by auto
  hence ys: "ys ≠ []" by auto
  from arg_min_list_in[OF this] have "?y ∈ set ys" by auto
  from this[unfolded roots] 
  show "?y^n = x" by auto
  fix z
  assume "z^n = x"
  hence z: "z ∈ set ys" unfolding roots by auto
  from f_arg_min_list_f[OF ys, of ?f] z 
  have "?f ?y ≤ ?f z" by simp
  thus "Re z < Re ?y ∨ Re z = Re ?y ∧ Im z ≤ Im ?y" by auto
qed

lemma croot_impl[code]: "croot n x = (if n = 0 then 0 else 
  arg_min_list (λ y. (- Re y, - Im y)) (all_croots n x))" 
proof (cases "n = 0")
  case n0: False
  hence id: "(if n = 0 then y else z) = z" 
    for y z u :: complex by auto
  show ?thesis unfolding id Let_def 
    by (rule croot_via_root_selection[OF _ n0], rule all_croots[OF n0])
qed auto

end