# 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 
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 = "
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 = " 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 (z3) 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
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
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:]"
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