Session Complex_Geometry

Theory More_Transcendental

(* ---------------------------------------------------------------------------- *)
section ‹Introduction›
(* ---------------------------------------------------------------------------- *)

text ‹The complex plane or some of its parts (e.g., the unit disc or the upper half plane) are often
taken as the domain in which models of various geometries (both Euclidean and non-Euclidean ones)
are formalized. The complex plane gives simpler and more compact formulas than the Cartesian plane.
Within complex plane is easier to describe geometric objects and perform the calculations (usually
shedding some new light on the subject). We give a formalization of the extended complex
plane (given both as a complex projective space and as the Riemann sphere), its objects (points,
circles and lines), and its transformations (Möbius transformations).›

(* ---------------------------------------------------------------------------- *)
section ‹Related work›
(* ---------------------------------------------------------------------------- *)

text‹During the last decade, there have been many results in formalizing
geometry in proof-assistants. Parts of Hilbert’s seminal book
,,Foundations of Geometry'' \cite{hilbert} have been formalized both
in Coq and Isabelle/Isar.  Formalization of first two groups of axioms
in Coq, in an intuitionistic setting was done by Dehlinger et
al. \cite{hilbert-coq}. First formalization in Isabelle/HOL was done
by Fleuriot and Meikele \cite{hilbert-isabelle}, and some further
developments were made in master thesis of Scott \cite{hilbert-scott}.
Large fragments of Tarski's geometry \cite{tarski} have been
formalized in Coq by Narboux et al. \cite{narboux-tarski}. Within Coq,
there are also formalizations of von Plato’s constructive geometry by
Kahn \cite{vonPlato,von-plato-formalization}, French high school
geometry by Guilhot \cite{guilhot} and ruler and compass geometry by
Duprat \cite{duprat2008}, etc.

In our previous work \cite{petrovic2012formalizing}, we have already
formally investigated a Cartesian model of Euclidean geometry. 
›

(* ---------------------------------------------------------------------------- *)
section ‹Background theories› 
(* ---------------------------------------------------------------------------- *)

text ‹In this section we introduce some basic mathematical notions and prove some lemmas needed in the rest of our
formalization. We describe:

     trigonometric functions,

     complex numbers, 

     systems of two and three linear equations with two unknowns (over arbitrary fields), 

     quadratic equations (over real and complex numbers), systems of quadratic and real
      equations, and systems of two quadratic equations,

     two-dimensional vectors and matrices over complex numbers.
›

(* -------------------------------------------------------------------------- *)
subsection ‹Library Additions for Trigonometric Functions›
(* -------------------------------------------------------------------------- *)

theory More_Transcendental
  imports Complex_Main "HOL-Library.Periodic_Fun"
begin

text ‹Additional properties of @{term sin} and @{term cos} functions that are later used in proving
conjectures for argument of complex number.›

text ‹Sign of trigonometric functions on some characteristic intervals.›

lemma cos_lt_zero_on_pi2_pi [simp]:
  assumes "x > pi/2" and "x  pi"
  shows "cos x < 0"
  using cos_gt_zero_pi[of "pi - x"] assms
  by simp

text ‹Value of trigonometric functions in points $k\pi$ and $\frac{\pi}{2} + k\pi$.›

lemma sin_kpi [simp]:
  fixes k::int
  shows "sin (k * pi) = 0"
  by (simp add: sin_zero_iff_int2)

lemma cos_odd_kpi [simp]:
  fixes k::int
  assumes "odd k"
  shows "cos (k * pi) = -1"
  by (simp add: assms mult.commute)

lemma cos_even_kpi [simp]:
  fixes k::int
  assumes "even k"
  shows "cos (k * pi) = 1"
  by (simp add: assms mult.commute)

lemma sin_pi2_plus_odd_kpi [simp]:
  fixes k::int
  assumes "odd k"
  shows "sin (pi / 2 + k * pi) = -1"
  using assms
  by (simp add: sin_add)

lemma sin_pi2_plus_even_kpi [simp]:
  fixes k::int
  assumes "even k"
  shows "sin (pi / 2 + k * pi) = 1"
  using assms
  by (simp add: sin_add)

text ‹Solving trigonometric equations and systems with special values (0, 1, or -1) of sine and cosine functions›

lemma cos_0_iff_canon:
  assumes "cos φ = 0" and "-pi < φ" and "φ  pi"
  shows "φ = pi/2  φ = -pi/2"
  by (smt (verit, best) arccos_0 arccos_cos assms cos_minus divide_minus_left)

lemma sin_0_iff_canon:
  assumes "sin φ = 0" and "-pi < φ" and "φ  pi"
  shows "φ = 0  φ = pi"
  using assms sin_eq_0_pi by force

lemma cos0_sin1:
  assumes "sin φ = 1"
  shows " k::int. φ = pi/2 + 2*k*pi"
  by (smt (verit, ccfv_threshold) assms cos_diff cos_one_2pi_int cos_pi_half mult_cancel_right1 sin_pi_half sin_plus_pi)

(* TODO: add lemmas for cos = -1, sin = 0 and cos = 0, sin = -1 *)


text ‹Sine is injective on $[-\frac{\pi}{2}, \frac{\pi}{2}]$›

lemma sin_inj:
  assumes "-pi/2  α  α  pi/2" and "-pi/2  α'  α'  pi/2"
  assumes "α  α'"
  shows "sin α  sin α'"
  by (metis assms divide_minus_left sin_inj_pi)

text ‹Periodicity of trigonometric functions›

text ‹The following are available in HOL-Decision\_Procs.Approximation\_Bounds, but we want to avoid
that dependency›

lemma sin_periodic_nat [simp]: 
  fixes n :: nat
  shows "sin (x + n * (2 * pi)) = sin x"
  by (metis (no_types, hide_lams) add.commute add.left_neutral cos_2npi cos_one_2pi_int mult.assoc mult.commute mult.left_neutral mult_zero_left sin_add sin_int_2pin)

lemma sin_periodic_int [simp]:
  fixes i :: int
  shows "sin (x + i * (2 * pi)) = sin x"
  by (metis add.right_neutral cos_int_2pin mult.commute mult.right_neutral mult_zero_right sin_add sin_int_2pin)

lemma cos_periodic_nat [simp]: 
  fixes n :: nat
  shows "cos (x + n * (2 * pi)) = cos x"
  by (metis add.left_neutral cos_2npi cos_add cos_periodic mult.assoc mult_2 mult_2_right of_nat_numeral sin_periodic sin_periodic_nat)

lemma cos_periodic_int [simp]:
  fixes i :: int
  shows "cos (x + i * (2 * pi)) = cos x"
  by (metis cos_add cos_int_2pin diff_zero mult.commute mult.right_neutral mult_zero_right sin_int_2pin)

text ‹Values of both sine and cosine are repeated only after multiples of $2\cdot \pi$›

lemma sin_cos_eq:
  fixes a b :: real
  assumes "cos a = cos b" and "sin a = sin b"
  shows " k::int. a - b = 2*k*pi"
  by (metis assms cos_diff cos_one_2pi_int mult.commute sin_cos_squared_add3)

text ‹The following two lemmas are consequences of surjectivity of cosine for the range $[-1, 1]$.›

lemma ex_cos_eq:
  assumes "-pi/2  α  α  pi/2"
  assumes "a  0" and "a < 1"
  shows " α'. -pi/2  α'  α'  pi/2  α'  α  cos (α - α') = a"
proof-
  have "arccos a > 0" "arccos a  pi/2"
    using a  0 a < 1
    using arccos_lt_bounded arccos_le_pi2
    by auto

  show ?thesis
  proof (cases "α - arccos a  - pi/2")
    case True
    thus ?thesis
      using assms ‹arccos a > 0 ‹arccos a  pi/2
      by (rule_tac x = "α - arccos a" in exI) auto
  next
    case False
    thus ?thesis
      using assms ‹arccos a > 0 ‹arccos a  pi/2
      by (rule_tac x = "α + arccos a" in exI) auto
  qed
qed

lemma ex_cos_gt:
  assumes "-pi/2  α  α  pi/2"
  assumes "a < 1"
  shows " α'. -pi/2  α'  α'  pi/2  α'  α  cos (α - α') > a"
proof-
  obtain a' where "a'  0" "a' > a" "a' < 1"
    by (metis assms(2) dense_le_bounded linear not_one_le_zero)
  thus ?thesis
    using ex_cos_eq[of α a'] assms
    by auto
qed

text ‹The function @{term atan2} is a generalization of @{term arctan} that takes a pair of coordinates
of non-zero points returns its angle in the range $[-\pi, \pi)$.›

definition atan2 where
  "atan2 y x =
    (if x > 0 then arctan (y/x)
     else if x < 0 then
          if y > 0 then arctan (y/x) + pi else arctan (y/x) - pi
     else
          if y > 0 then pi/2 else if y < 0 then -pi/2 else 0)"

lemma atan2_bounded: 
  shows "-pi  atan2 y x  atan2 y x < pi"
  using arctan_bounded[of "y/x"] zero_le_arctan_iff[of "y/x"] arctan_le_zero_iff[of "y/x"] zero_less_arctan_iff[of "y/x"] arctan_less_zero_iff[of "y/x"]
  using divide_neg_neg[of y x] divide_neg_pos[of y x] divide_pos_pos[of y x]  divide_pos_neg[of y x]
  unfolding atan2_def
  by (simp (no_asm_simp)) auto

end

Theory Canonical_Angle

(* -------------------------------------------------------------------------- *)
subsection ‹Canonical angle› 
(* -------------------------------------------------------------------------- *)

text ‹Canonize any angle to $(-\pi, \pi]$ (taking account of $2\pi$ periodicity of @{term sin} and
@{term cos}). With this function, for example, multiplicative properties of @{term arg} for complex
numbers can easily be expressed and proved.›

theory Canonical_Angle
imports More_Transcendental
begin


abbreviation canon_ang_P where
 "canon_ang_P α α'  (-pi < α'  α'  pi)  ( k::int. α - α' = 2*k*pi)"

definition canon_ang :: "real  real" ("_") where
  "α = (THE α'. canon_ang_P α α')"

text ‹There is a canonical angle for every angle.›
lemma canon_ang_ex:
  shows " α'. canon_ang_P α α'"
proof-
  have ***: " α::real.  α'. 0 < α'  α'  1  ( k::int. α' = α - k)"
  proof
    fix α::real
    show "α'>0. α'  1  (k::int. α' = α - k)"
    proof (cases "α = floor α")
      case True
      thus ?thesis
        by (rule_tac x="α - floor α + 1" in exI, auto) (rule_tac x="floor α - 1" in exI, auto)
    next
      case False
      thus ?thesis
        using real_of_int_floor_ge_diff_one[of "α"]
        using of_int_floor_le[of "α"]
        by (rule_tac x="α - floor α" in exI) smt
    qed
  qed

  have **: " α::real.  α'. 0 < α'  α'  2  ( k::int. α - α' = 2*k - 1)"
  proof
    fix α::real
    from ***[rule_format, of "(α + 1) /2"]
    obtain α' and k::int where "0 < α'" "α'  1" "α' = (α + 1)/2 - k"
      by force
    hence "0 < α'" "α'  1" "α' = α/2 - k + 1/2"
      by auto
    thus "α'>0. α'  2  (k::int. α - α' = real_of_int (2 * k - 1))"
      by (rule_tac x="2*α'" in exI) auto
  qed
  have *: " α::real.  α'. -1 < α'  α'  1  ( k::int. α - α' = 2*k)"
  proof
    fix α::real
    from ** obtain α' and k :: int where
      "0 < α'  α'  2  α - α' = 2*k - 1"
      by force
    thus "α'>-1. α'  1  (k. α - α' = real_of_int (2 * (k::int)))"
      by (rule_tac x="α' - 1" in exI) (auto simp add: field_simps)
  qed
  obtain α' k where 1: "α' >- 1  α'  1" and 2: "α / pi - α' = real_of_int (2 * k)"
    using *[rule_format, of "α / pi"]
    by auto
  have "α'*pi > -pi  α'*pi  pi" 
    using 1
    by (smt mult.commute mult_le_cancel_left1 mult_minus_right pi_gt_zero)
  moreover
  have "α - α'*pi = 2 * real_of_int k * pi"
    using 2
    by (auto simp add: field_simps)
  ultimately
  show ?thesis
    by auto
qed

text ‹Canonical angle of any angle is unique.›
lemma canon_ang_unique:
  assumes "canon_ang_P α α1" and "canon_ang_P α α2"
  shows "α1 = α2"
proof-
  obtain k1::int where "α - α1 = 2*k1*pi"
    using assms(1)
    by auto
  obtain k2::int where "α - α2 = 2*k2*pi"
    using assms(2)
    by auto
  hence *: "-α1 + α2 = 2*(k1 - k2)*pi"
    using α - α1 = 2*k1*pi›
    by (simp add:field_simps)
  moreover
  have "-α1 + α2 < 2 * pi" "-α1 + α2 > -2*pi"
    using assms
    by auto
  ultimately
  have "-α1 + α2 = 0"
    using mult_less_cancel_right[of "-2" pi "real_of_int(2 * (k1 - k2))"]
    by auto
  thus ?thesis
    by auto
qed

text ‹Canonical angle is always in $(-\pi, \pi]$ and differs from the starting angle by $2k\pi$.›
lemma canon_ang:
  shows "-pi < α" and "α  pi" and " k::int. α - α = 2*k*pi"
proof-
  obtain α' where "canon_ang_P α α'"
    using canon_ang_ex[of α]
    by auto
  have "canon_ang_P α α"
    unfolding canon_ang_def
  proof (rule theI[where a="α'"])
    show "canon_ang_P α α'"
      by fact
  next
    fix α''
    assume "canon_ang_P α α''"
    thus "α'' = α'"
      using ‹canon_ang_P α α'
      using canon_ang_unique[of α' α α'']
      by simp
  qed
  thus "-pi < α" "α  pi" " k::int. α - α = 2*k*pi"
    by auto
qed

text ‹Angles in $(-\pi, \pi]$ are already canonical.›
lemma canon_ang_id:
  assumes  "-pi < α  α  pi"
  shows "α = α"
  using assms
  using canon_ang_unique[of "canon_ang α" α α] canon_ang[of α]
  by auto

text ‹Angles that differ by $2k\pi$ have equal canonical angles.›
lemma canon_ang_eq:
  assumes " k::int. α1 - α2 = 2*k*pi"
  shows "α1 = α2"
proof-
  obtain k'::int where *: "- pi < α1" "α1  pi" "α1 - α1 = 2 * k' * pi"
    using canon_ang[of α1]
    by auto

  obtain k''::int where **: "- pi < α2" "α2  pi" "α2 - α2 = 2 * k'' * pi"
    using canon_ang[of α2]
    by auto

  obtain k::int where ***: "α1 - α2 = 2*k*pi"
    using assms
    by auto

  have "m::int. α1 - α2 = 2 * m * pi"
    using **(3) ***
    by (rule_tac x="k+k''" in exI) (auto simp add: field_simps)

  thus ?thesis
    using canon_ang_unique[of "α1" α1 "α2"] * **
    by auto
qed

text ‹Introduction and elimination rules›
lemma canon_ang_eqI:
  assumes "k::int. α' - α = 2 * k * pi" and "- pi < α'  α'  pi"
  shows "α = α'"
  using assms
  using canon_ang_eq[of α' α]
  using canon_ang_id[of α']
  by auto

lemma canon_ang_eqE:
  assumes "α1 = α2"
  shows " (k::int). α1 - α2 = 2 *k * pi"
proof-
  obtain k1 k2 :: int where
    "α1 - α1 = 2 * k1 * pi"
    "α2 - α2 = 2 * k2 * pi"
    using canon_ang[of α1] canon_ang[of α2]
    by auto
  thus ?thesis
    using assms
    by (rule_tac x="k1 - k2" in exI) (auto simp add: field_simps)
qed

text ‹Canonical angle of opposite angle›

lemma canon_ang_uminus:
  assumes "α  pi"
  shows "-α = -α"
proof (rule canon_ang_eqI)
  show "x::int. - α - - α = 2 * x * pi"
    using canon_ang(3)[of α]
    by (metis minus_diff_eq minus_diff_minus)
next
  show "- pi < - α  - α  pi"
    using canon_ang(1)[of α] canon_ang(2)[of α] assms
    by auto
qed

lemma canon_ang_uminus_pi:
  assumes "α = pi"
  shows "-α = α"
proof (rule canon_ang_eqI)
  obtain k::int where "α - α = 2 * k * pi"
    using canon_ang(3)[of α]
    by auto
  thus "x::int. α - - α = 2 * x * pi"
    using assms
    by (rule_tac x="k+(1::int)" in exI) (auto simp add: field_simps)
next
  show "- pi < α  α  pi"
    using assms
    by auto
qed

text ‹Canonical angle of difference of two angles›
lemma canon_ang_diff:
  shows "α - β = α - β"
proof (rule canon_ang_eq)
  show "x::int. α - β - (α - β) = 2 * x * pi"
  proof-
    obtain k1::int where "α - α = 2*k1*pi"
      using canon_ang(3)
      by auto
    moreover
    obtain k2::int where "β - β = 2*k2*pi"
      using canon_ang(3)
      by auto
    ultimately
    show ?thesis
      by (rule_tac x="k1 - k2" in exI) (auto simp add: field_simps)
  qed
qed

text ‹Canonical angle of sum of two angles›
lemma canon_ang_sum:
  shows "α + β = α + β"
proof (rule canon_ang_eq)
  show "x::int. α + β - (α + β) = 2 * x * pi"
  proof-
    obtain k1::int where "α - α = 2*k1*pi"
      using canon_ang(3)
      by auto
    moreover
    obtain k2::int where "β - β = 2*k2*pi"
      using canon_ang(3)
      by auto
    ultimately
    show ?thesis
      by (rule_tac x="k1 + k2" in exI) (auto simp add: field_simps)
  qed
qed

text ‹Canonical angle of angle from $(0, 2\pi]$ shifted by $\pi$›

lemma canon_ang_plus_pi1:
  assumes "0 < α" and "α  2*pi"
  shows "α + pi = α - pi"
proof (rule canon_ang_eqI)
  show " x::int. α - pi - (α + pi) = 2 * x * pi"
    by (rule_tac x="-1" in exI) auto
next
  show "- pi < α - pi  α - pi  pi"
    using assms
    by auto
qed

lemma canon_ang_minus_pi1:
  assumes "0 < α" and "α  2*pi"
  shows "α - pi = α - pi"
proof (rule canon_ang_id)
  show "- pi < α - pi  α - pi  pi"
    using assms
    by auto
qed

text ‹Canonical angle of angles from $(-2\pi, 0]$ shifted by $\pi$›

lemma canon_ang_plus_pi2:
  assumes "-2*pi < α" and "α  0"
  shows "α + pi = α + pi"
proof (rule canon_ang_id)
  show "- pi < α + pi  α + pi  pi"
    using assms
    by auto
qed

lemma canon_ang_minus_pi2:
  assumes "-2*pi < α" and "α  0"
  shows "α - pi = α + pi"
proof (rule canon_ang_eqI)
  show " x::int. α + pi - (α - pi) = 2 * x * pi"
    by (rule_tac x="1" in exI) auto
next
  show "- pi < α + pi  α + pi  pi"
    using assms
    by auto
qed

text ‹Canonical angle of angle in $(\pi, 3\pi]$.›
lemma canon_ang_pi_3pi: 
  assumes "pi < α" and "α  3 * pi"
  shows "α = α - 2*pi"
proof-
  have "x. - pi = pi * real_of_int x"
    by (rule_tac x="-1" in exI, simp)
  thus ?thesis
    using assms canon_ang_eqI[of "α - 2*pi" "α"]
    by auto
qed

text ‹Canonical angle of angle in $(-3\pi, -\pi]$.›
lemma canon_ang_minus_3pi_minus_pi: 
  assumes "-3*pi < α" and "α  -pi"
  shows "α = α + 2*pi"
proof-
  have "x. pi = pi * real_of_int x"
    by (rule_tac x="1" in exI, simp)
  thus ?thesis
    using assms canon_ang_eqI[of "α + 2*pi" "α"]
    by auto
qed

text ‹Canonical angles for some special angles›

lemma zero_canonical [simp]:
  shows "0 = 0"
  using canon_ang_eqI[of 0 0]
  by simp

lemma pi_canonical [simp]:
  shows "pi = pi"
  by (simp add: canon_ang_id)

lemma two_pi_canonical [simp]:
  shows "2 * pi = 0"
  using canon_ang_plus_pi1[of "pi"]
  by simp

text ‹Canonization preserves sine and cosine›
lemma canon_ang_sin [simp]:
  shows "sin α = sin α"
proof-
  obtain x::int where "α = α + pi * (x * 2)"
    using canon_ang(3)[of α]
    by (auto simp add: field_simps)
  thus ?thesis
    using sin_periodic_int[of "α" x]
    by (simp add: field_simps)
qed

lemma canon_ang_cos [simp]:
  shows "cos α = cos α"
proof-
  obtain x::int where "α = α + pi * (x * 2)"
    using canon_ang(3)[of α]
    by (auto simp add: field_simps)
  thus ?thesis
    using cos_periodic_int[of "α" x]
    by (simp add: field_simps)
qed

end

Theory More_Complex

(* -------------------------------------------------------------------------- *)
subsection ‹Library Additions for Complex Numbers›
(* -------------------------------------------------------------------------- *)

text ‹Some additional lemmas about complex numbers.›

theory More_Complex
  imports Complex_Main More_Transcendental Canonical_Angle
begin

text ‹Conjugation and @{term cis}
          
declare cis_cnj[simp] 

lemma rcis_cnj: 
  shows "cnj a = rcis (cmod a) (- arg a)"
  by (subst rcis_cmod_arg[of a, symmetric]) (simp add: rcis_def)

lemmas complex_cnj = complex_cnj_diff complex_cnj_mult complex_cnj_add complex_cnj_divide complex_cnj_minus

text ‹Some properties for @{term complex_of_real}. Also, since it is often used in our
formalization we abbreviate it to @{term cor}.›

abbreviation cor :: "real  complex" where
  "cor  complex_of_real"

lemma cmod_cis [simp]:
  assumes "a  0"
  shows "cor (cmod a) * cis (arg a) = a"
  using assms
  by (metis rcis_cmod_arg rcis_def)

lemma cis_cmod [simp]:
  assumes "a  0"
  shows "cis (arg a) * cor (cmod a) = a"
  using assms cmod_cis[of a]
  by (simp add: field_simps)

lemma cor_squared:
  shows "(cor x)2 = cor (x2)"
  by (simp add: power2_eq_square)

lemma cor_sqrt_mult_cor_sqrt [simp]:
  shows "cor (sqrt A) * cor (sqrt A) = cor ¦A¦"
  by (metis of_real_mult real_sqrt_mult_self)

lemma cor_eq_0: "cor x + 𝗂 * cor y = 0  x = 0  y = 0"
 by (metis Complex_eq Im_complex_of_real Im_i_times Re_complex_of_real add_cancel_left_left of_real_eq_0_iff plus_complex.sel(2) zero_complex.code)

lemma one_plus_square_neq_zero [simp]:
  shows "1 + (cor x)2  0"
  by (metis (hide_lams, no_types) of_real_1 of_real_add of_real_eq_0_iff of_real_power power_one sum_power2_eq_zero_iff zero_neq_one)

text ‹Additional lemmas about @{term Complex} constructor. Following newer versions of Isabelle,
these should be deprecated.›

lemma complex_real_two [simp]:
  shows "Complex 2 0 = 2"
  by (simp add: Complex_eq)

lemma complex_double [simp]:
  shows "(Complex a b) * 2 = Complex (2*a) (2*b)"
  by (simp add: Complex_eq)

lemma complex_half [simp]: 
  shows "(Complex a b) / 2 = Complex (a/2) (b/2)"
  by (subst complex_eq_iff) auto

lemma Complex_scale1:
  shows "Complex (a * b) (a * c) = cor a * Complex b c"
  unfolding complex_of_real_def
  unfolding Complex_eq
  by (auto simp add: field_simps)

lemma Complex_scale2: 
  shows "Complex (a * c) (b * c) = Complex a b * cor c"
  unfolding complex_of_real_def
  unfolding Complex_eq
  by (auto simp add: field_simps)

lemma Complex_scale3: 
  shows "Complex (a / b) (a / c) = cor a * Complex (1 / b) (1 / c)"
  unfolding complex_of_real_def
  unfolding Complex_eq
  by (auto simp add: field_simps)

lemma Complex_scale4:
  shows "c  0  Complex (a / c) (b / c) = Complex a b / cor c"
  unfolding complex_of_real_def
  unfolding Complex_eq
  by (auto simp add: field_simps power2_eq_square)

lemma Complex_Re_express_cnj:
  shows "Complex (Re z) 0 = (z + cnj z) / 2"
  by (cases z) (simp add: Complex_eq)

lemma Complex_Im_express_cnj:
  shows "Complex 0 (Im z) = (z - cnj z)/2"
  by (cases z) (simp add: Complex_eq)

text ‹Additional properties of @{term cmod}.›

lemma complex_mult_cnj_cmod:
  shows "z * cnj z = cor ((cmod z)2)"
  using complex_norm_square
  by auto

lemma cmod_square: 
  shows "(cmod z)2 = Re (z * cnj z)"
  using complex_mult_cnj_cmod[of z]
  by (simp add: power2_eq_square)

lemma cor_cmod_power_4 [simp]:
  shows "cor (cmod z) ^ 4 = (z * cnj z)2"
  by (simp add: complex_mult_cnj_cmod)

lemma cnjE:
  assumes "x  0"
  shows "cnj x = cor ((cmod x)2) / x"
  using complex_mult_cnj_cmod[of x] assms
  by (auto simp add: field_simps)

lemma cmod_cor_divide [simp]:
  shows "cmod (z / cor k) = cmod z / ¦k¦"
  by (simp add: norm_divide)

lemma cmod_mult_minus_left_distrib [simp]:
  shows "cmod (z*z1 - z*z2) = cmod z * cmod(z1 - z2)"
  by (metis norm_mult right_diff_distrib)

lemma cmod_eqI:
  assumes "z1 * cnj z1 = z2 * cnj z2"
  shows "cmod z1 = cmod z2"
  using assms
  by (subst complex_mod_sqrt_Re_mult_cnj)+ auto

lemma cmod_eqE:
  assumes "cmod z1 = cmod z2"
  shows "z1 * cnj z1 = z2 * cnj z2"
  by (simp add: assms complex_mult_cnj_cmod)

lemma cmod_eq_one [simp]:
  shows "cmod a = 1  a*cnj a = 1"
  by (metis cmod_eqE cmod_eqI complex_cnj_one monoid_mult_class.mult.left_neutral norm_one)

text ‹We introduce @{term is_real} (the imaginary part of complex number is zero) and @{term is_imag}
(real part of complex number is zero) operators and prove some of their properties.›

abbreviation is_real where
  "is_real z  Im z = 0"

abbreviation is_imag where
  "is_imag z  Re z = 0"

lemma real_imag_0:
  assumes "is_real a" "is_imag a" 
  shows "a = 0"
  using assms
  by (simp add: complex.expand)

lemma complex_eq_if_Re_eq:
  assumes "is_real z1" and "is_real z2"
  shows "z1 = z2  Re z1 = Re z2"
  using assms
  by (cases z1, cases z2) auto

lemma mult_reals [simp]:
  assumes "is_real a" and "is_real b"
  shows "is_real (a * b)"
  using assms
  by auto

lemma div_reals [simp]:
  assumes "is_real a" and "is_real b"
  shows "is_real (a / b)"
  using assms
  by (simp add: complex_is_Real_iff)

lemma complex_of_real_Re [simp]:
  assumes "is_real k"
  shows "cor (Re k) = k"
  using assms
  by (cases k) (auto simp add: complex_of_real_def)

lemma cor_cmod_real:
  assumes "is_real a"
  shows "cor (cmod a) = a  cor (cmod a) = -a"
  using assms
  unfolding cmod_def
  by (cases "Re a > 0") auto

lemma eq_cnj_iff_real:
  shows "cnj z = z  is_real z"
  by (cases z) (simp add: Complex_eq)

lemma eq_minus_cnj_iff_imag:
  shows "cnj z = -z  is_imag z"
  by (cases z) (simp add: Complex_eq)

lemma Re_divide_real:
  assumes "is_real b" and "b  0"
  shows "Re (a / b) = (Re a) / (Re b)"
  using assms
  by (simp add: complex_is_Real_iff)

lemma Re_mult_real:
  assumes "is_real a"
  shows "Re (a * b) = (Re a) * (Re b)"
  using assms
  by simp

lemma Im_mult_real:
  assumes "is_real a"
  shows "Im (a * b) = (Re a) * (Im b)"
  using assms
  by simp

lemma Im_divide_real:
  assumes "is_real b" and "b  0"
  shows "Im (a / b) = (Im a) / (Re b)"
  using assms
  by (simp add: complex_is_Real_iff)

lemma Re_sgn:
  assumes "is_real R"
  shows "Re (sgn R) = sgn (Re R)"
  using assms
  by (metis Re_sgn complex_of_real_Re norm_of_real real_sgn_eq)

lemma is_real_div:
  assumes "b  0"
  shows "is_real (a / b)  a*cnj b = b*cnj a"
  using assms
  by (metis complex_cnj_divide complex_cnj_zero_iff eq_cnj_iff_real frac_eq_eq mult.commute)

lemma is_real_mult_real:
  assumes "is_real a" and "a  0"
  shows "is_real b  is_real (a * b)"
  using assms
  by (cases a, auto simp add: Complex_eq)

lemma Im_express_cnj:
  shows "Im z = (z - cnj z) / (2 * 𝗂)"
  by (simp add: complex_diff_cnj field_simps)

lemma Re_express_cnj: 
  shows "Re z = (z + cnj z) / 2"
  by (simp add: complex_add_cnj)

text ‹Rotation of complex number for 90 degrees in the positive direction.›

abbreviation rot90 where
  "rot90 z  Complex (-Im z) (Re z)"

lemma rot90_ii: 
  shows "rot90 z = z * 𝗂"
  by (metis Complex_mult_i complex_surj)

text ‹With @{term cnj_mix} we introduce scalar product between complex vectors. This operation shows
to be useful to succinctly express some conditions.›

abbreviation cnj_mix where
  "cnj_mix z1 z2  cnj z1 * z2 + z1 * cnj z2"

abbreviation scalprod where
  "scalprod z1 z2  cnj_mix z1 z2 / 2"

lemma cnj_mix_minus:
  shows "cnj z1*z2 - z1*cnj z2 = 𝗂 * cnj_mix (rot90 z1) z2"
  by (cases z1, cases z2) (simp add: Complex_eq field_simps)

lemma cnj_mix_minus':
  shows "cnj z1*z2 - z1*cnj z2 = rot90 (cnj_mix (rot90 z1) z2)"
  by (cases z1, cases z2) (simp add: Complex_eq field_simps)

lemma cnj_mix_real [simp]:
  shows "is_real (cnj_mix z1 z2)"
  by (cases z1, cases z2) simp

lemma scalprod_real [simp]:
  shows "is_real (scalprod z1 z2)"
  using cnj_mix_real
  by simp

text ‹Additional properties of @{term cis} function.›

lemma cis_minus_pi2 [simp]:
  shows "cis (-pi/2) = -𝗂"
  by (simp add: cis_inverse[symmetric])

lemma cis_pi2_minus_x [simp]:
  shows "cis (pi/2 - x) = 𝗂 * cis(-x)"
  using cis_divide[of "pi/2" x, symmetric]
  using cis_divide[of 0 x, symmetric]
  by simp

lemma cis_pm_pi [simp]: 
  shows "cis (x - pi) = - cis x" and  "cis (x + pi) = - cis x"
  by (simp add: cis.ctr complex_minus)+


lemma cis_times_cis_opposite [simp]: 
  shows "cis φ * cis (- φ) = 1"
  by (simp add: cis_mult)

text @{term cis} repeats only after $2k\pi$›
lemma cis_eq:
  assumes "cis a = cis b"
  shows " k::int. a - b = 2 * k * pi"
  using assms sin_cos_eq[of a b]
  using cis.sel[of a] cis.sel[of b]
  by (cases "cis a", cases "cis b") auto

text @{term cis} is injective on $(-\pi, \pi]$.›
lemma cis_inj:
  assumes "-pi < α" and "α  pi" and "-pi < α'" and "α'  pi"
  assumes "cis α = cis α'"
  shows "α = α'"
  using assms
  by (metis arg_unique sgn_cis)

text @{term cis} of an angle combined with @{term cis} of the opposite angle›

lemma cis_diff_cis_opposite [simp]: 
  shows "cis φ - cis (- φ) = 2 * 𝗂 * sin φ"
  using Im_express_cnj[of "cis φ"]
  by simp

lemma cis_opposite_diff_cis [simp]:
  shows "cis (-φ) - cis (φ) = - 2 * 𝗂 * sin φ"
  using cis_diff_cis_opposite[of "-φ"]
  by simp

lemma cis_add_cis_opposite [simp]: 
  shows "cis φ + cis (-φ) = 2 * cos φ"
  by (metis cis.sel(1) cis_cnj complex_add_cnj)

text @{term cis} equal to 1 or -1›
lemma cis_one [simp]:
  assumes "sin φ = 0" and "cos φ = 1"
  shows "cis φ = 1"
  using assms
  by (auto simp add: cis.ctr one_complex.code)

lemma cis_minus_one [simp]:
  assumes "sin φ = 0" and "cos φ = -1"
  shows "cis φ = -1"
  using assms
  by (auto simp add: cis.ctr Complex_eq_neg_1)

(* -------------------------------------------------------------------------- *)
subsubsection ‹Additional properties of complex number argument›
(* -------------------------------------------------------------------------- *)

text @{term arg} of real numbers›

lemma is_real_arg1:
  assumes "arg z = 0  arg z = pi"
  shows "is_real z"
  using assms
  using rcis_cmod_arg[of z] Im_rcis[of "cmod z" "arg z"]
  by auto

lemma is_real_arg2:
  assumes "is_real z"
  shows "arg z = 0  arg z = pi"
proof (cases "z = 0")
  case False
  thus ?thesis
    using arg_bounded[of z]
    by (smt (verit, best) Im_sgn assms cis.simps(2) cis_arg div_0 sin_zero_pi_iff)
qed (auto simp add: arg_zero)

lemma arg_complex_of_real_positive [simp]:
  assumes "k > 0"
  shows "arg (cor k) = 0"
proof-
  have "cos (arg (Complex k 0)) > 0"
    using assms
    using rcis_cmod_arg[of "Complex k 0"] Re_rcis[of "cmod (Complex k 0)" "arg (Complex k 0)"]
    using cmod_eq_Re by force
  thus ?thesis
    using assms is_real_arg2[of "cor k"]
    unfolding complex_of_real_def
    by auto
qed

lemma arg_complex_of_real_negative [simp]:
  assumes "k < 0"
  shows "arg (cor k) = pi"
proof-
  have "cos (arg (Complex k 0)) < 0"
    using k < 0 rcis_cmod_arg[of "Complex k 0"] Re_rcis[of "cmod (Complex k 0)" "arg (Complex k 0)"]
    by (metis complex.sel(1) mult_less_0_iff norm_not_less_zero)
  thus ?thesis
    using assms is_real_arg2[of "cor k"]
    unfolding complex_of_real_def
    by auto
qed

lemma arg_0_iff:
  shows "z  0  arg z = 0  is_real z  Re z > 0"
  by (smt arg_complex_of_real_negative arg_complex_of_real_positive arg_zero complex_of_real_Re is_real_arg1 pi_gt_zero zero_complex.simps)

lemma arg_pi_iff:
  shows "arg z = pi  is_real z  Re z < 0"
  by (smt arg_complex_of_real_negative arg_complex_of_real_positive arg_zero complex_of_real_Re is_real_arg1 pi_gt_zero zero_complex.simps)


text @{term arg} of imaginary numbers›

lemma is_imag_arg1:
  assumes "arg z = pi/2  arg z = -pi/2"
  shows "is_imag z"
  using assms
  using rcis_cmod_arg[of z] Re_rcis[of "cmod z" "arg z"]
  by (metis cos_minus cos_pi_half minus_divide_left mult_eq_0_iff)

lemma is_imag_arg2:
  assumes "is_imag z" and "z  0"
  shows "arg z = pi/2  arg z = -pi/2"
  using arg_bounded assms cos_0_iff_canon cos_arg_i_mult_zero by presburger

lemma arg_complex_of_real_times_i_positive [simp]:
  assumes "k > 0"
  shows "arg (cor k * 𝗂) = pi / 2"
proof-
  have "sin (arg (Complex 0 k)) > 0"
    using k > 0 rcis_cmod_arg[of "Complex 0 k"] Im_rcis[of "cmod (Complex 0 k)" "arg (Complex 0 k)"]
    by (smt complex.sel(2) mult_nonneg_nonpos norm_ge_zero)
  thus ?thesis
    using assms is_imag_arg2[of "cor k * 𝗂"]
    using arg_zero complex_of_real_i
    by force
qed

lemma arg_complex_of_real_times_i_negative [simp]:
  assumes "k < 0"
  shows "arg (cor k * 𝗂) = - pi / 2"
proof-
  have "sin (arg (Complex 0 k)) < 0"
    using k < 0 rcis_cmod_arg[of "Complex 0 k"] Im_rcis[of "cmod (Complex 0 k)" "arg (Complex 0 k)"]
    by (metis complex.sel(2) mult_less_0_iff norm_not_less_zero)
  thus ?thesis
    using assms is_imag_arg2[of "cor k * 𝗂"]
    using arg_zero complex_of_real_i[of k]
    by (smt complex.sel(1) sin_pi_half sin_zero)
qed

lemma arg_pi2_iff:
  shows "z  0  arg z = pi / 2  is_imag z  Im z > 0"
  by (smt Im_rcis Re_i_times Re_rcis arcsin_minus_1 cos_pi_half divide_minus_left mult.commute mult_cancel_right1 rcis_cmod_arg is_imag_arg2 sin_arcsin sin_pi_half zero_less_mult_pos zero_less_norm_iff)

lemma arg_minus_pi2_iff:
  shows "z  0  arg z = - pi / 2  is_imag z  Im z < 0"
  by (smt arg_pi2_iff complex.expand divide_cancel_right pi_neq_zero is_imag_arg1 is_imag_arg2 zero_complex.simps(1) zero_complex.simps(2))

lemma arg_ii [simp]:
  shows "arg 𝗂 = pi/2"
  by (metis arg_pi2_iff imaginary_unit.sel zero_less_one)

lemma arg_minus_ii [simp]: 
  shows "arg (-𝗂) = -pi/2"
proof-
  have "-𝗂 = cis (arg (- 𝗂))"
    using rcis_cmod_arg[of "-𝗂"]
    by (simp add: rcis_def)
  hence "cos (arg (-𝗂)) = 0" "sin (arg (-𝗂)) = -1"
    using cis.simps[of "arg (-𝗂)"]
    by auto
  thus ?thesis
    using cos_0_iff_canon[of "arg (-𝗂)"] arg_bounded[of "-𝗂"]
    by fastforce
qed

text ‹Argument is a canonical angle›

lemma canon_ang_arg:
  shows "arg z = arg z"
  using canon_ang_id[of "arg z"] arg_bounded
  by simp

lemma arg_cis:
  shows "arg (cis φ) = φ"
  using arg_unique canon_ang canon_ang_cos canon_ang_sin cis.ctr sgn_cis by presburger

text ‹Cosine and sine of @{term arg}

lemma cos_arg:
  assumes "z  0"
  shows "cos (arg z) = Re z / cmod z"
  by (metis Complex.Re_sgn cis.simps(1) assms cis_arg)

lemma sin_arg:
  assumes "z  0"
  shows "sin (arg z) = Im z / cmod z"
  by (metis Complex.Im_sgn cis.simps(2) assms cis_arg)

text ‹Argument of product›

lemma cis_arg_mult:
  assumes "z1 * z2  0"
  shows "cis (arg (z1 * z2)) = cis (arg z1 + arg z2)"
  by (metis assms cis_arg cis_mult mult_eq_0_iff sgn_mult)

lemma arg_mult_2kpi:
  assumes "z1 * z2  0"
  shows " k::int. arg (z1 * z2) = arg z1 + arg z2 + 2*k*pi"
proof-
  have "cis (arg (z1*z2)) = cis (arg z1 + arg z2)"
    by (rule cis_arg_mult[OF assms])
  thus ?thesis
    using cis_eq[of "arg (z1*z2)" "arg z1 + arg z2"]
    by (auto simp add: field_simps)
qed

lemma arg_mult:
  assumes "z1 * z2  0"
  shows "arg(z1 * z2) = arg z1 + arg z2"
proof-
  obtain k::int where "arg(z1 * z2) = arg z1 + arg z2 + 2*k*pi"
    using arg_mult_2kpi[of z1 z2]
    using assms
    by auto
  hence "arg(z1 * z2) = arg z1 + arg z2"
    using canon_ang_eq
    by(simp add:field_simps)
  thus ?thesis
    using canon_ang_arg[of "z1*z2"]
    by auto
qed

lemma arg_mult_real_positive [simp]:
  assumes "k > 0"
  shows "arg (cor k * z) = arg z"
proof (cases "z = 0")
  case False
  thus ?thesis
    using arg_mult assms canon_ang_arg by force
qed (auto simp: arg_zero)

lemma arg_mult_real_negative [simp]:
  assumes "k < 0"
  shows "arg (cor k * z) = arg (-z)"
proof (cases "z = 0")
  case False
  thus ?thesis
    using assms
    by (metis arg_mult_real_positive minus_mult_commute neg_0_less_iff_less of_real_minus minus_minus)
qed (auto simp: arg_zero)

lemma arg_div_real_positive [simp]:
  assumes "k > 0"
  shows "arg (z / cor k) = arg z"
proof(cases "z = 0")
  case True
  thus ?thesis
    by auto
next
  case False
  thus ?thesis
    using assms
    using arg_mult_real_positive[of "1/k" z]
    by auto
qed

lemma arg_div_real_negative [simp]:
  assumes "k < 0"
  shows "arg (z / cor k) = arg (-z)"
proof(cases "z = 0")
  case True
  thus ?thesis
    by auto
next
  case False
  thus ?thesis
    using assms
    using arg_mult_real_negative[of "1/k" z]
    by auto
qed

lemma arg_mult_eq:
  assumes "z * z1  0" and "z * z2  0"
  assumes "arg (z * z1) = arg (z * z2)"
  shows "arg z1 = arg z2"
  by (metis (no_types, lifting) arg_cis assms canon_ang_arg cis_arg mult_eq_0_iff nonzero_mult_div_cancel_left sgn_divide)

text ‹Argument of conjugate›

lemma arg_cnj_pi:
  assumes "arg z = pi"
  shows "arg (cnj z) = pi"
  using arg_pi_iff assms by auto

lemma arg_cnj_not_pi:
  assumes "arg z  pi"
  shows "arg (cnj z) = -arg z"
proof(cases "arg z = 0")
  case True
  thus ?thesis
    using eq_cnj_iff_real[of z] is_real_arg1[of z] by force
next
  case False
  have "arg (cnj z) = arg z  arg(cnj z) = -arg z"
    using arg_bounded[of z] arg_bounded[of "cnj z"]
    by (smt (verit, best) arccos_cos arccos_cos2 cnj.sel(1) complex_cnj_zero_iff complex_mod_cnj cos_arg)
  moreover
  have "arg (cnj z)  arg z"
    using sin_0_iff_canon[of "arg (cnj z)"] arg_bounded False assms
    by (metis complex_mod_cnj eq_cnj_iff_real is_real_arg2 rcis_cmod_arg)
  ultimately
  show ?thesis
    by auto
qed

text ‹Argument of reciprocal›

lemma arg_inv_not_pi:
  assumes "z  0" and "arg z  pi"
  shows "arg (1 / z) = - arg z"
proof-
  have "1/z = cnj z / cor ((cmod z)2 )"
    using z  0 complex_mult_cnj_cmod[of z]
    by (auto simp add:field_simps)
  thus ?thesis
    using arg_div_real_positive[of "(cmod z)2" "cnj z"] z  0
    using arg_cnj_not_pi[of z] ‹arg z  pi›
    by auto
qed

lemma arg_inv_pi:
  assumes "z  0" and "arg z = pi"
  shows "arg (1 / z) = pi"
proof-
  have "1/z = cnj z / cor ((cmod z)2 )"
    using z  0 complex_mult_cnj_cmod[of z]
    by (auto simp add:field_simps)
  thus ?thesis
    using arg_div_real_positive[of "(cmod z)2" "cnj z"] z  0
    using arg_cnj_pi[of z] ‹arg z = pi›
    by auto
qed

lemma arg_inv_2kpi:
  assumes "z  0"
  shows " k::int. arg (1 / z) = - arg z + 2*k*pi"
  using arg_inv_pi[OF assms]
  using arg_inv_not_pi[OF assms]
  by (cases "arg z = pi") (rule_tac x="1" in exI, simp, rule_tac x="0" in exI, simp)

lemma arg_inv:
  assumes "z  0"
  shows "arg (1 / z) = - arg z"
  by (metis arg_inv_not_pi arg_inv_pi assms canon_ang_arg canon_ang_uminus_pi)

text ‹Argument of quotient›

lemma arg_div_2kpi:
  assumes "z1  0" and "z2  0"
  shows " k::int. arg (z1 / z2) = arg z1 - arg z2 + 2*k*pi"
proof-
  obtain x1 where "arg (z1 * (1 / z2)) = arg z1 + arg (1 / z2) + 2 * real_of_int x1 * pi"
    using assms arg_mult_2kpi[of z1 "1/z2"]
    by auto
  moreover
  obtain x2 where "arg (1 / z2) = - arg z2 + 2 * real_of_int x2 * pi"
    using assms arg_inv_2kpi[of z2]
    by auto
  ultimately
  show ?thesis
    by (rule_tac x="x1 + x2" in exI, simp add: field_simps)
qed

lemma arg_div:
  assumes "z1  0" and "z2  0"
  shows "arg(z1 / z2) = arg z1 - arg z2"
proof-
  obtain k::int where "arg(z1 / z2) = arg z1 - arg z2 + 2*k*pi"
    using arg_div_2kpi[of z1 z2]
    using assms
    by auto
  hence "canon_ang(arg(z1 / z2)) = canon_ang(arg z1 - arg z2)"
    using canon_ang_eq
    by(simp add:field_simps)
  thus ?thesis
    using canon_ang_arg[of "z1/z2"]
    by auto
qed

text ‹Argument of opposite›

lemma arg_uminus:
  assumes "z  0"
  shows "arg (-z) = arg z + pi"
  using assms
  using arg_mult[of "-1" z]
  using arg_complex_of_real_negative[of "-1"]
  by (auto simp add: field_simps)

lemma arg_uminus_opposite_sign:
  assumes "z  0"
  shows "arg z > 0  ¬ arg (-z) > 0"
proof (cases "arg z = 0")
  case True
  thus ?thesis
    using assms
    by (simp add: arg_uminus)
next
  case False
  show ?thesis
  proof (cases "arg z > 0")
    case True
    thus ?thesis
      using assms
      using arg_bounded[of z]
      using canon_ang_plus_pi1[of "arg z"]
      by (simp add: arg_uminus)
  next
    case False
    thus ?thesis
      using ‹arg z  0
      using assms
      using arg_bounded[of z]
      using canon_ang_plus_pi2[of "arg z"]
      by (simp add: arg_uminus)
  qed
qed

text ‹Sign of argument is the same as the sign of the Imaginary part›

lemma arg_Im_sgn:
  assumes "¬ is_real z"
  shows "sgn (arg z) = sgn (Im z)"
proof-
  have "z  0"
    using assms
    by auto
  then obtain r φ where polar: "z = cor r * cis φ" "φ = arg z" "r > 0"
    by (smt cmod_cis mult_eq_0_iff norm_ge_zero of_real_0)
  hence "Im z = r * sin φ"
    by (metis Im_mult_real Re_complex_of_real cis.simps(2) Im_complex_of_real)
  hence  "Im z > 0  sin φ > 0" "Im z < 0  sin φ < 0"
    using r > 0
    using mult_pos_pos mult_nonneg_nonneg zero_less_mult_pos mult_less_cancel_left
    by smt+
  moreover
  have "φ  pi" "φ  0"
    using ¬ is_real z polar cis_pi
    by force+
  hence "sin φ > 0  φ > 0" "φ < 0  sin φ < 0"
    using φ = arg z φ  0 φ  pi›
    using arg_bounded[of z]
    by (smt sin_gt_zero sin_le_zero sin_pi_minus sin_0_iff_canon sin_ge_zero)+
  ultimately
  show ?thesis
    using φ = arg z
    by auto
qed


subsubsection ‹Complex square root›

definition
  "ccsqrt z = rcis (sqrt (cmod z)) (arg z / 2)"

lemma square_ccsqrt [simp]:
  shows "(ccsqrt x)2 = x"
  unfolding ccsqrt_def
  by (subst DeMoivre2) (simp add: rcis_cmod_arg)

lemma ex_complex_sqrt:
  shows " s::complex. s*s = z"
  unfolding power2_eq_square[symmetric]
  by (rule_tac x="csqrt z" in exI) simp

lemma ccsqrt:
  assumes "s * s = z"
  shows "s = ccsqrt z  s = -ccsqrt z"
proof (cases "s = 0")
  case True
  thus ?thesis
    using assms
    unfolding ccsqrt_def
    by simp
next
  case False
  then obtain k::int where "cmod s * cmod s = cmod z" "2 * arg s - arg z = 2*k*pi"
    using assms
    using rcis_cmod_arg[of z] rcis_cmod_arg[of s]
    using arg_mult[of s s]
    using canon_ang(3)[of "2*arg s"]
    by (auto simp add: norm_mult arg_mult)
  have *: "sqrt (cmod z) = cmod s"
    using ‹cmod s * cmod s = cmod z
    by (smt norm_not_less_zero real_sqrt_abs2)

  have **: "arg z / 2 = arg s - k*pi"
    using 2 * arg s - arg z = 2*k*pi›
    by simp

  have "cis (arg s - k*pi) = cis (arg s)  cis (arg s - k*pi) = -cis (arg s)"
  proof (cases "even k")
    case True
    hence "cis (arg s - k*pi) = cis (arg s)"
      by (simp add: cis_def complex.corec cos_diff sin_diff)
    thus ?thesis
      by simp
  next
    case False
    hence "cis (arg s - k*pi) = -cis (arg s)"
      by (simp add: cis_def complex.corec Complex_eq cos_diff sin_diff)
    thus ?thesis
      by simp
  qed
  thus ?thesis
  proof
    assume ***: "cis (arg s - k * pi) = cis (arg s)"
    hence "s = ccsqrt z"
      using rcis_cmod_arg[of s]
      unfolding ccsqrt_def rcis_def
      by (subst *, subst **, subst ***, simp)
    thus ?thesis
      by simp
  next
    assume ***: "cis (arg s - k * pi) = -cis (arg s)"
    hence "s = - ccsqrt z"
      using rcis_cmod_arg[of s]
      unfolding ccsqrt_def rcis_def
      by (subst *, subst **, subst ***, simp)
    thus ?thesis
      by simp
  qed
qed

lemma null_ccsqrt [simp]:
  shows "ccsqrt x = 0  x = 0"
  unfolding ccsqrt_def
  by auto

lemma ccsqrt_mult:
  shows "ccsqrt (a * b) = ccsqrt a * ccsqrt b 
         ccsqrt (a * b) = - ccsqrt a * ccsqrt b"
proof (cases "a = 0  b = 0")
  case True
  thus ?thesis
    by auto
next
  case False
  obtain k::int where "arg a + arg b - arg a + arg b = 2 * real_of_int k * pi"
    using canon_ang(3)[of "arg a + arg b"]
    by auto
  hence *: "arg a + arg b = arg a + arg b - 2 * (real_of_int k) * pi"
    by (auto simp add: field_simps)

  have "cis (arg a + arg b / 2) = cis (arg a / 2 + arg b / 2)  cis (arg a + arg b / 2) = - cis (arg a / 2 + arg b / 2)"
    using cos_even_kpi[of k] cos_odd_kpi[of k]
    by ((subst *)+, (subst diff_divide_distrib)+, (subst add_divide_distrib)+)
       (cases "even k", auto simp add: cis_def complex.corec Complex_eq cos_diff sin_diff)
  thus ?thesis
    using False
    unfolding ccsqrt_def
    by (smt (verit, best) arg_mult mult_minus_left mult_minus_right no_zero_divisors norm_mult rcis_def rcis_mult real_sqrt_mult)
qed

lemma csqrt_real:
  assumes "is_real x"
  shows "(Re x  0  ccsqrt x = cor (sqrt (Re x))) 
         (Re x < 0  ccsqrt x = 𝗂 * cor (sqrt (- (Re x))))"
proof (cases "x = 0")
  case True
  thus ?thesis
    by auto
next
  case False
  show ?thesis
  proof (cases "Re x > 0")
    case True
    hence "arg x = 0"
      using ‹is_real x
      by (metis arg_complex_of_real_positive complex_of_real_Re)
    thus ?thesis
      using ‹Re x > 0 ‹is_real x
      unfolding ccsqrt_def
      by (simp add: cmod_eq_Re)
  next
    case False
    hence "Re x < 0"
      using x  0 ‹is_real x
      using complex_eq_if_Re_eq by auto
    hence "arg x = pi"
      using ‹is_real x
      by (metis arg_complex_of_real_negative complex_of_real_Re)
    thus ?thesis
      using ‹Re x < 0 ‹is_real x
      unfolding ccsqrt_def rcis_def
      by (simp add: cis_def complex.corec Complex_eq cmod_eq_Re)
  qed
qed


text ‹Rotation of complex vector to x-axis.›

lemma is_real_rot_to_x_axis:
  assumes "z  0"
  shows "is_real (cis (-arg z) * z)"
proof (cases "arg z = pi")
  case True
  thus ?thesis
    using is_real_arg1[of z]
    by auto
next
  case False
  hence "- arg z = - arg z"
    using canon_ang_eqI[of "- arg z" "-arg z"]
    using arg_bounded[of z]
    by (auto simp add: field_simps)
  hence "arg (cis (- (arg z)) * z) = 0"
    using arg_mult[of "cis (- (arg z))" z] z  0
    using arg_cis[of "- arg z"]
    by simp
  thus ?thesis
    using is_real_arg1[of "cis (- arg z) * z"]
    by auto
qed

lemma positive_rot_to_x_axis:
  assumes "z  0"
  shows "Re (cis (-arg z) * z) > 0"
  using assms
  by (smt Re_complex_of_real cis_rcis_eq mult_cancel_right1 rcis_cmod_arg rcis_mult rcis_zero_arg zero_less_norm_iff)

text ‹Inequalities involving @{term cmod}.›

lemma cmod_1_plus_mult_le:
  shows "cmod (1 + z*w)  sqrt((1 + (cmod z)2) * (1 + (cmod w)2))"
proof-
  have "Re ((1+z*w)*(1+cnj z*cnj w))  Re (1+z*cnj z)* Re (1+w*cnj w)"
  proof-
    have "Re ((w - cnj z)*cnj(w - cnj z))  0"
      by (subst complex_mult_cnj_cmod) (simp add: power2_eq_square)
    hence "Re (z*w + cnj z * cnj w)  Re (w*cnj w) + Re(z*cnj z)"
      by (simp add: field_simps)
    thus ?thesis
      by (simp add: field_simps)
  qed
  hence "(cmod (1 + z * w))2  (1 + (cmod z)2) * (1 + (cmod w)2)"
    by (subst cmod_square)+ simp
  thus ?thesis
    by (metis abs_norm_cancel real_sqrt_abs real_sqrt_le_iff)
qed

lemma cmod_diff_ge: 
  shows "cmod (b - c)  sqrt (1 + (cmod b)2) - sqrt (1 + (cmod c)2)"
proof-
  have "(cmod (b - c))2 + (1/2*Im(b*cnj c - c*cnj b))2  0"
    by simp
  hence "(cmod (b - c))2  - (1/2*Im(b*cnj c - c*cnj b))2"
    by simp
  hence "(cmod (b - c))2  (1/2*Re(b*cnj c + c*cnj b))2 - Re(b*cnj b*c*cnj c) "
    by (auto simp add: power2_eq_square field_simps)
  hence "Re ((b - c)*(cnj b - cnj c))  (1/2*Re(b*cnj c + c*cnj b))2 - Re(b*cnj b*c*cnj c)"
    by (subst (asm) cmod_square) simp
  moreover
  have "(1 + (cmod b)2) * (1 + (cmod c)2) = 1 + Re(b*cnj b) + Re(c*cnj c) + Re(b*cnj b*c*cnj c)"
    by (subst cmod_square)+ (simp add: field_simps power2_eq_square)
  moreover
  have "(1 + Re (scalprod b c))2 = 1 + 2*Re(scalprod b c) + ((Re (scalprod b c))2)"
    by (subst power2_sum) simp
  hence "(1 + Re (scalprod b c))2 = 1 + Re(b*cnj c + c*cnj b) + (1/2 * Re (b*cnj c + c*cnj b))2"
    by simp
  ultimately
  have "(1 + (cmod b)2) * (1 + (cmod c)2)  (1 + Re (scalprod b c))2"
    by (simp add: field_simps)
  moreover
  have "sqrt((1 + (cmod b)2) * (1 + (cmod c)2))  0"
    by (metis one_power2 real_sqrt_sum_squares_mult_ge_zero)
  ultimately
  have "sqrt((1 + (cmod b)2) * (1 + (cmod c)2))  1 + Re (scalprod b c)"
    by (metis power2_le_imp_le real_sqrt_ge_0_iff real_sqrt_pow2_iff)
  hence "Re ((b - c) * (cnj b - cnj c))  1 + Re (c*cnj c) + 1 + Re (b*cnj b) - 2*sqrt((1 + (cmod b)2) * (1 + (cmod c)2))"
    by (simp add: field_simps)
  hence *: "(cmod (b - c))2  (sqrt (1 + (cmod b)2) - sqrt (1 + (cmod c)2))2"
    apply (subst cmod_square)+
    apply (subst (asm) cmod_square)+
    apply (subst power2_diff)
    apply (subst real_sqrt_pow2, simp)
    apply (subst real_sqrt_pow2, simp)
    apply (simp add: real_sqrt_mult)
    done
  thus ?thesis
  proof (cases "sqrt (1 + (cmod b)2) - sqrt (1 + (cmod c)2) > 0")
    case True
    thus ?thesis
      using power2_le_imp_le[OF *]
      by simp
  next
    case False
    hence "0  sqrt (1 + (cmod b)2) - sqrt (1 + (cmod c)2)"
      by (metis less_eq_real_def linorder_neqE_linordered_idom)
    moreover
    have "cmod (b - c)  0"
      by simp
    ultimately
    show ?thesis
      by (metis add_increasing monoid_add_class.add.right_neutral)
  qed
qed

lemma cmod_diff_le:
  shows "cmod (b - c)  sqrt (1 + (cmod b)2) + sqrt (1 + (cmod c)2)"
proof-
  have "(cmod (b + c))2 + (1/2*Im(b*cnj c - c*cnj b))2  0"
    by simp
  hence "(cmod (b + c))2  - (1/2*Im(b*cnj c - c*cnj b))2"
    by simp
  hence "(cmod (b + c))2  (1/2*Re(b*cnj c + c*cnj b))2 - Re(b*cnj b*c*cnj c) "
    by (auto simp add: power2_eq_square field_simps)
  hence "Re ((b + c)*(cnj b + cnj c))  (1/2*Re(b*cnj c + c*cnj b))2 - Re(b*cnj b*c*cnj c)"
    by (subst (asm) cmod_square) simp
  moreover
  have "(1 + (cmod b)2) * (1 + (cmod c)2) = 1 + Re(b*cnj b) + Re(c*cnj c) + Re(b*cnj b*c*cnj c)"
    by (subst cmod_square)+ (simp add: field_simps power2_eq_square)
  moreover
  have ++: "2*Re(scalprod b c) = Re(b*cnj c + c*cnj b)"
    by simp
  have "(1 - Re (scalprod b c))2 = 1 - 2*Re(scalprod b c) + ((Re (scalprod b c))2)"
    by (subst power2_diff) simp
  hence "(1 - Re (scalprod b c))2 = 1 - Re(b*cnj c + c*cnj b) + (1/2 * Re (b*cnj c + c*cnj b))2"
    by (subst ++[symmetric]) simp
  ultimately
  have "(1 + (cmod b)2) * (1 + (cmod c)2)  (1 - Re (scalprod b c))2"
    by (simp add: field_simps)
  moreover
  have "sqrt((1 + (cmod b)2) * (1 + (cmod c)2))  0"
    by (metis one_power2 real_sqrt_sum_squares_mult_ge_zero)
  ultimately
  have "sqrt((1 + (cmod b)2) * (1 + (cmod c)2))  1 - Re (scalprod b c)"
    by (metis power2_le_imp_le real_sqrt_ge_0_iff real_sqrt_pow2_iff)
  hence "Re ((b - c) * (cnj b - cnj c))  1 + Re (c*cnj c) + 1 + Re (b*cnj b) + 2*sqrt((1 + (cmod b)2) * (1 + (cmod c)2))"
    by (simp add: field_simps)
  hence *: "(cmod (b - c))2  (sqrt (1 + (cmod b)2) + sqrt (1 + (cmod c)2))2"
    apply (subst cmod_square)+
    apply (subst (asm) cmod_square)+
    apply (subst power2_sum)
    apply (subst real_sqrt_pow2, simp)
    apply (subst real_sqrt_pow2, simp)
    apply (simp add: real_sqrt_mult)
    done
  thus ?thesis
    using power2_le_imp_le[OF *]
    by simp
qed


text ‹Definition of Euclidean distance between two complex numbers.›

definition cdist where
  [simp]: "cdist z1 z2  cmod (z2 - z1)"

text ‹Misc. properties of complex numbers.›

lemma ex_complex_to_complex [simp]:
  fixes z1 z2 :: complex
  assumes "z1  0" and "z2  0"
  shows "k. k  0  z2 = k * z1"
  using assms
  by (rule_tac x="z2/z1" in exI) simp

lemma ex_complex_to_one [simp]:
  fixes z::complex
  assumes "z  0"
  shows "k. k  0  k * z = 1"
  using assms
  by (rule_tac x="1/z" in exI) simp

lemma ex_complex_to_complex2 [simp]:
  fixes z::complex
  shows "k. k  0  k * z = z"
  by (rule_tac x="1" in exI) simp

lemma complex_sqrt_1:
  fixes z::complex
  assumes "z  0"
  shows "z = 1 / z  z = 1  z = -1"
  using assms
  using nonzero_eq_divide_eq square_eq_iff
  by fastforce

end

Theory Angles

(* ---------------------------------------------------------------------------- *)
subsection ‹Angle between two vectors›
(* ---------------------------------------------------------------------------- *)

text ‹In this section we introduce different measures of angle between two vectors (represented by complex numbers).›

theory Angles
imports More_Transcendental Canonical_Angle More_Complex
begin

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Oriented angle›
(* ---------------------------------------------------------------------------- *)

text ‹Oriented angle between two vectors (it is always in the interval $(-\pi, \pi]$).›
definition ang_vec ("") where
  [simp]: " z1 z2  arg z2 - arg z1"

lemma ang_vec_bounded:
  shows "-pi <  z1 z2   z1 z2  pi"
  by (simp add: canon_ang(1) canon_ang(2))

lemma ang_vec_sym:
  assumes " z1 z2  pi"
  shows " z1 z2 = -  z2 z1"
  using assms
  unfolding ang_vec_def
  using canon_ang_uminus[of "arg z2 - arg z1"]
  by simp

lemma ang_vec_sym_pi:
  assumes " z1 z2 = pi"
  shows " z1 z2 =  z2 z1"
  using assms
  unfolding ang_vec_def
  using canon_ang_uminus_pi[of "arg z2 - arg z1"]
  by simp

lemma ang_vec_plus_pi1:
  assumes " z1 z2 > 0"
  shows " z1 z2 + pi =  z1 z2 - pi"
proof (rule canon_ang_eqI)
  show " x::int.  z1 z2 - pi - ( z1 z2 + pi) = 2 * real_of_int x * pi"
    by (rule_tac x="-1" in exI) auto
next
  show "- pi <  z1 z2 - pi   z1 z2 - pi  pi"
    using assms
    unfolding ang_vec_def
    using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
    by auto
qed

lemma ang_vec_plus_pi2:
  assumes " z1 z2  0"
  shows " z1 z2 + pi =  z1 z2 + pi"
proof (rule canon_ang_id)
  show "- pi <  z1 z2 + pi   z1 z2 + pi  pi"
    using assms
    unfolding ang_vec_def
    using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
    by auto
qed

lemma ang_vec_opposite1:
  assumes "z1  0"
  shows " (-z1) z2 =  z1 z2 - pi"
proof-
  have " (-z1) z2 = arg z2 - (arg z1 + pi)"
    unfolding ang_vec_def
    using arg_uminus[OF assms] 
    using canon_ang_arg[of z2, symmetric]
    using canon_ang_diff[of "arg z2" "arg z1 + pi", symmetric]
    by simp
  moreover
  have " z1 z2 - pi = arg z2 - arg z1 - pi"
    using canon_ang_id[of pi, symmetric]
    using canon_ang_diff[of "arg z2 - arg z1" "pi", symmetric]
    by simp_all
  ultimately
  show ?thesis
    by (simp add: field_simps)
qed

lemma ang_vec_opposite2:
  assumes "z2  0"
  shows " z1 (-z2) =  z1 z2 + pi"
  unfolding ang_vec_def
  using arg_mult[of "-1" "z2"] assms
  using arg_complex_of_real_negative[of "-1"]
  using canon_ang_diff[of "arg (-1) + arg z2" "arg z1", symmetric]
  using canon_ang_sum[of "arg z2 - arg z1" "pi", symmetric]
  using canon_ang_id[of pi] canon_ang_arg[of z1]
  by (auto simp: algebra_simps)
  

lemma ang_vec_opposite_opposite:
  assumes "z1  0" and "z2  0"
  shows " (-z1) (-z2) =  z1 z2"
proof-
  have " (-z1) (-z2) =  z1 z2 + pi - pi"
    using ang_vec_opposite1[OF assms(1)]
    using ang_vec_opposite2[OF assms(2)]
    using canon_ang_id[of pi, symmetric]
    by simp_all
  also have "... =  z1 z2"
    by (subst canon_ang_diff[symmetric], simp)
  finally
  show ?thesis
    by (metis ang_vec_def canon_ang(1) canon_ang(2) canon_ang_id)
qed

lemma ang_vec_opposite_opposite':
  assumes "z1  z" and "z2  z"
  shows " (z - z1) (z - z2) =  (z1 - z) (z2 - z)"
using ang_vec_opposite_opposite[of "z - z1" "z - z2"] assms
by (simp add: field_simps del: ang_vec_def)

text ‹Cosine, scalar product and the law of cosines›

lemma cos_cmod_scalprod:
  shows "cmod z1 * cmod z2 * (cos ( z1 z2)) = Re (scalprod z1 z2)"
proof (cases "z1 = 0  z2 = 0")
  case True
  thus ?thesis
    by auto
next
  case False
  thus ?thesis
    by (simp add: cos_diff cos_arg sin_arg field_simps)
qed

lemma cos0_scalprod0:
  assumes "z1  0" and "z2  0"
  shows "cos ( z1 z2) = 0  scalprod z1 z2 = 0"
  using assms
  using cnj_mix_real[of z1 z2]
  using cos_cmod_scalprod[of z1 z2]
  by (auto simp add: complex_eq_if_Re_eq)

lemma ortho_scalprod0:
  assumes "z1  0" and "z2  0"
  shows " z1 z2 = pi/2   z1 z2 = -pi/2  scalprod z1 z2 = 0"
  using cos0_scalprod0[OF assms]
  using ang_vec_bounded[of z1 z2]
  using cos_0_iff_canon[of " z1 z2"]
  by (metis cos_minus cos_pi_half divide_minus_left)

lemma law_of_cosines:
  shows "(cdist B C)2 = (cdist A C)2 + (cdist A B)2 - 2*(cdist A C)*(cdist A B)*(cos ( (C-A) (B-A)))"
proof-
  let ?a = "C-B" and ?b = "C-A" and ?c = "B-A"
  have "?a = ?b - ?c"
    by simp
  hence "(cmod ?a)2 = (cmod (?b - ?c))2"
    by metis
  also have "... = Re (scalprod (?b-?c) (?b-?c))"
    by (simp add: cmod_square)
  also have "... = (cmod ?b)2 + (cmod ?c)2 - 2*Re (scalprod ?b ?c)"
    by (simp add: cmod_square field_simps)
  finally
  show ?thesis
    using cos_cmod_scalprod[of ?b ?c]
    by simp
qed

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Unoriented angle›
(* ---------------------------------------------------------------------------- *)

text ‹Convex unoriented angle between two vectors (it is always in the interval $[0, pi]$).›
definition ang_vec_c ("∠c") where
  [simp]:"∠c z1 z2  abs ( z1 z2)"

lemma ang_vec_c_sym:
  shows "∠c z1 z2 = ∠c z2 z1"
  unfolding ang_vec_c_def
  using ang_vec_sym_pi[of z1 z2] ang_vec_sym[of z1 z2]
  by (cases " z1 z2 = pi") auto

lemma ang_vec_c_bounded: "0  ∠c z1 z2  ∠c z1 z2  pi"
  using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
  by auto

text ‹Cosine and scalar product›

lemma cos_c_: "cos (∠c z1 z2) = cos ( z1 z2)"
  unfolding ang_vec_c_def
  by (smt cos_minus)

lemma ortho_c_scalprod0:
  assumes "z1  0" and "z2  0"
  shows "∠c z1 z2 = pi/2  scalprod z1 z2 = 0"
proof-
  have " z1 z2 = pi / 2   z1 z2 = - pi / 2  ∠c z1 z2 = pi/2"
    unfolding ang_vec_c_def
    using arctan 
    by force
  thus ?thesis
    using ortho_scalprod0[OF assms]
    by simp
qed

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Acute angle›
(* ---------------------------------------------------------------------------- *)

text ‹Acute or right angle (non-obtuse) between two vectors (it is always in the interval $[0, \frac{\pi}{2}$]).
We will use this to measure angle between two circles, since it can always be acute (or right).›

definition acute_ang where
  [simp]: "acute_ang α = (if α > pi / 2 then pi - α else α)"

definition ang_vec_a ("∠a") where
  [simp]: "∠a z1 z2  acute_ang (∠c z1 z2)"

lemma ang_vec_a_sym:
  "∠a z1 z2 = ∠a z2 z1"
  unfolding ang_vec_a_def
  using ang_vec_c_sym
  by auto

lemma ang_vec_a_opposite2:
  "∠a z1 z2 = ∠a z1 (-z2)"
proof(cases "z2  = 0")
  case True
  thus ?thesis
    by (metis minus_zero)
next
  case False
  thus ?thesis
  proof(cases " z1 z2 < -pi / 2")
    case True
    hence " z1 z2 < 0"
      using pi_not_less_zero
      by linarith
    have "∠a z1 z2 = pi +  z1 z2"
      using True  z1 z2 < 0
      unfolding ang_vec_a_def ang_vec_c_def ang_vec_a_def abs_real_def
      by auto
    moreover
    have "∠a z1 (-z2) = pi +  z1 z2"
      unfolding ang_vec_a_def ang_vec_c_def abs_real_def
      using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
      using ang_vec_plus_pi2[of z1 z2] True  z1 z2 < 0 z2  0
      using ang_vec_opposite2[of z2 z1]
      by auto
    ultimately
    show ?thesis
      by auto
  next
    case False
    show ?thesis
    proof (cases " z1 z2  0")
      case True
      have "∠a z1 z2 = -  z1 z2"
        using ¬  z1 z2 < - pi / 2 True
        unfolding ang_vec_a_def ang_vec_c_def ang_vec_a_def abs_real_def
        by auto
      moreover
      have "∠a z1 (-z2) = -  z1 z2"
        using ¬  z1 z2 < - pi / 2 True
        unfolding ang_vec_a_def ang_vec_c_def abs_real_def
        using ang_vec_plus_pi2[of z1 z2]
        using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
        using z2  0 ang_vec_opposite2[of z2 z1]
        by auto
      ultimately
      show ?thesis
        by simp
    next
      case False
      show ?thesis
      proof (cases " z1 z2 < pi / 2")
        case True
        have "∠a z1 z2 =  z1 z2"
          using ¬  z1 z2  0 True
          unfolding ang_vec_a_def ang_vec_c_def ang_vec_a_def abs_real_def
          by auto
        moreover
        have "∠a z1 (-z2) =  z1 z2"
          using ¬  z1 z2  0 True
          unfolding ang_vec_a_def ang_vec_c_def abs_real_def
          using ang_vec_plus_pi1[of z1 z2]
          using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
          using z2  0 ang_vec_opposite2[of z2 z1]
          by auto
        ultimately
        show ?thesis
          by simp
      next
        case False
        have " z1 z2 > 0"
          using False
          by (metis less_linear less_trans pi_half_gt_zero)
        have "∠a z1 z2 = pi -  z1 z2"
          using False  z1 z2 > 0
          unfolding ang_vec_a_def ang_vec_c_def ang_vec_a_def abs_real_def
          by auto
        moreover
        have "∠a z1 (-z2) = pi -  z1 z2"
          unfolding ang_vec_a_def ang_vec_c_def abs_real_def
          using False  z1 z2 > 0
          using ang_vec_plus_pi1[of z1 z2]
          using canon_ang(1)[of "arg z2 - arg z1"] canon_ang(2)[of "arg z2 - arg z1"]
          using z2  0 ang_vec_opposite2[of z2 z1]
          by auto
        ultimately
        show ?thesis
          by auto
      qed
    qed
  qed
qed

lemma ang_vec_a_opposite1:
  shows "∠a z1 z2 = ∠a (-z1) z2"
  using ang_vec_a_sym[of "-z1" z2] ang_vec_a_opposite2[of z2 z1] ang_vec_a_sym[of z2 z1]
  by auto

lemma ang_vec_a_scale1:
  assumes "k  0"
  shows "∠a (cor k * z1) z2 = ∠a z1 z2"
proof (cases "k > 0")
  case True
  thus ?thesis
    unfolding ang_vec_a_def ang_vec_c_def ang_vec_def
    using arg_mult_real_positive[of k z1]
    by auto
next
  case False
  hence "k < 0"
    using assms
    by auto
  thus ?thesis
    using arg_mult_real_negative[of k z1]
    using ang_vec_a_opposite1[of z1 z2]
    unfolding ang_vec_a_def ang_vec_c_def ang_vec_def
    by simp
qed

lemma ang_vec_a_scale2:
  assumes "k  0"
  shows "∠a z1 (cor k * z2) = ∠a z1 z2"
  using ang_vec_a_sym[of z1 "complex_of_real k * z2"]
  using ang_vec_a_scale1[OF assms, of z2 z1]
  using ang_vec_a_sym[of z1 z2]
  by auto

lemma ang_vec_a_scale:
  assumes "k1  0" and "k2  0"
  shows "∠a (cor k1 * z1) (cor k2 * z2) = ∠a z1 z2"
  using ang_vec_a_scale1[OF assms(1)] ang_vec_a_scale2[OF assms(2)]
  by auto

lemma ang_a_cnj_cnj:
  shows "∠a z1 z2 = ∠a (cnj z1) (cnj z2)"
unfolding ang_vec_a_def ang_vec_c_def ang_vec_def
proof(cases "arg z1  pi  arg z2  pi")
  case True
  thus "acute_ang ¦arg z2 - arg z1¦ = acute_ang ¦arg (cnj z2) - arg (cnj z1)¦"
    using arg_cnj_not_pi[of z1] arg_cnj_not_pi[of z2]
    apply (auto simp del:acute_ang_def)
    proof(cases "arg z2 - arg z1 = pi")
      case True
      thus "acute_ang ¦arg z2 - arg z1¦ = acute_ang ¦arg z1 - arg z2¦"
        using  canon_ang_uminus_pi[of "arg z2 - arg z1"]
        by (auto simp add:field_simps)
    next
      case False
      thus "acute_ang ¦arg z2 - arg z1¦ = acute_ang ¦arg z1 - arg z2¦"
        using  canon_ang_uminus[of "arg z2 - arg z1"]
        by (auto simp add:field_simps)
    qed
  next
    case False
    thus "acute_ang ¦arg z2 - arg z1¦ = acute_ang ¦arg (cnj z2) - arg (cnj z1)¦"
    proof(cases "arg z1 = pi")
      case False
      hence "arg z2 = pi"
        using ¬ (arg z1  pi  arg z2  pi)
        by auto
      thus ?thesis
        using False
        using arg_cnj_not_pi[of z1] arg_cnj_pi[of z2]
        apply (auto simp del:acute_ang_def)
      proof(cases "arg z1 > 0")
          case True
          hence "-arg z1  0"
            by auto
          thus "acute_ang ¦pi - arg z1¦ = acute_ang ¦pi + arg z1¦"
            using True canon_ang_plus_pi1[of "arg z1"]
            using arg_bounded[of z1] canon_ang_plus_pi2[of "-arg z1"]
            by (auto simp add:field_simps)
        next
          case False
          hence "-arg z1  0"
             by simp
          thus "acute_ang ¦pi - arg z1¦ = acute_ang ¦pi + arg z1¦"
          proof(cases "arg z1 = 0")
            case True
            thus ?thesis
              by (auto simp del:acute_ang_def)
          next
            case False
            hence "-arg z1 > 0"
              using -arg z1  0
              by auto
            thus ?thesis
            using False canon_ang_plus_pi1[of "-arg z1"]
            using arg_bounded[of z1] canon_ang_plus_pi2[of "arg z1"]
            by (auto simp add:field_simps)
        qed
      qed
    next
      case True
      thus ?thesis
        using arg_cnj_pi[of z1]
        apply (auto simp del:acute_ang_def)
      proof(cases "arg z2 = pi")
        case True
        thus "acute_ang ¦arg z2 - pi¦ = acute_ang ¦arg (cnj z2) - pi¦"
          using arg_cnj_pi[of z2]
          by auto
      next
        case False
        thus "acute_ang ¦arg z2 - pi¦ = acute_ang ¦arg (cnj z2) - pi¦"
          using arg_cnj_not_pi[of z2]
          apply (auto simp del:acute_ang_def)
        proof(cases "arg z2 > 0")
          case True
          hence "-arg z2  0"
            by auto
          thus "acute_ang ¦arg z2 - pi¦ = acute_ang ¦- arg z2 - pi¦"
            using True canon_ang_minus_pi1[of "arg z2"]
            using arg_bounded[of z2] canon_ang_minus_pi2[of "-arg z2"]
            by (auto simp add: field_simps)
        next
          case False
          hence "-arg z2  0"
             by simp
          thus "acute_ang ¦arg z2 - pi¦ = acute_ang ¦- arg z2 - pi¦"
          proof(cases "arg z2 = 0")
            case True
            thus ?thesis
              by (auto simp del:acute_ang_def)
          next
            case False
            hence "-arg z2 > 0"
              using -arg z2  0
              by auto
            thus ?thesis
            using False canon_ang_minus_pi1[of "-arg z2"]
            using arg_bounded[of z2] canon_ang_minus_pi2[of "arg z2"]
            by (auto simp add:field_simps)
        qed
      qed
    qed
  qed
qed

text ‹Cosine and scalar product›

lemma ortho_a_scalprod0:
  assumes "z1  0" and "z2  0"
  shows "∠a z1 z2 = pi/2  scalprod z1 z2 = 0"
  unfolding ang_vec_a_def
  using assms ortho_c_scalprod0[of z1 z2]
  by auto

declare ang_vec_c_def[simp del]

lemma cos_a_c: "cos (∠a z1 z2) = abs (cos (∠c z1 z2))"
proof-
  have "0  ∠c z1 z2" "∠c z1 z2  pi"
    using ang_vec_c_bounded[of z1 z2]
    by auto
  show ?thesis
  proof (cases "∠c z1 z2 = pi/2")
    case True
    thus ?thesis
      unfolding ang_vec_a_def acute_ang_def
      by (smt cos_pi_half pi_def pi_half)
  next
    case False
    show ?thesis
    proof (cases "∠c z1 z2 < pi / 2")
      case True
      thus ?thesis
        using 0  ∠c z1 z2
        using cos_gt_zero_pi[of "∠c z1 z2"]
        unfolding ang_vec_a_def
        by simp
    next
      case False
      hence "∠c z1 z2 > pi/2"
        using ∠c z1 z2  pi/2
        by simp
      hence "cos (∠c z1 z2) < 0"
        using ∠c z1 z2  pi›
        using cos_lt_zero_on_pi2_pi[of "∠c z1 z2"] 
        by simp
      thus ?thesis
        using ∠c z1 z2 > pi/2
        unfolding ang_vec_a_def
        by simp
    qed
  qed
qed

end

Theory More_Set

(* ---------------------------------------------------------------------------- *)
subsection ‹Library Aditions for Set Cardinality›
(* ---------------------------------------------------------------------------- *)

text ‹In this section some additional simple lemmas about set cardinality are proved.›

theory More_Set
imports Main
begin

text ‹Every infinite set has at least two different elements›
lemma infinite_contains_2_elems:
  assumes "infinite A"
  shows " x y. x  y  x  A  y  A"
  by (metis assms finite.simps is_singletonI' is_singleton_def)

text ‹Every infinite set has at least three different elements›
lemma infinite_contains_3_elems:
  assumes "infinite A"
  shows " x y z. x  y  x  z  y  z  x  A  y  A  z  A"
  by (metis Diff_iff assms infinite_contains_2_elems infinite_remove insertI1)

text ‹Every set with cardinality greater than 1 has at least two different elements›
lemma card_geq_2_iff_contains_2_elems:
  shows "card A  2  finite A  ( x y. x  y  x  A  y  A)"
proof (intro iffI conjI)
  assume *: "finite A  ( x y. x  y  x  A  y  A)"
  thus "card A  2"
    by (metis card_0_eq card_Suc_eq empty_iff leI less_2_cases singletonD)
next
  assume *: "2  card A"
  then show "finite A"
    using card.infinite by force
  show " x y. x  y  x  A  y  A"
    by (meson "*" card_2_iff' in_mono obtain_subset_with_card_n)
qed

text ‹Set cardinality is at least 3 if and only if it contains three different elements›
lemma card_geq_3_iff_contains_3_elems:
  shows "card A  3  finite A  ( x y z. x  y  x  z  y  z  x  A  y  A  z  A)"
proof (intro iffI conjI)
  assume *: "card A  3"
  then show "finite A"
    using card.infinite by force
  show " x y z. x  y  x  z  y  z  x  A  y  A  z  A"
    by (smt (verit, best) "*" card_2_iff' card_geq_2_iff_contains_2_elems le_cases3 not_less_eq_eq numeral_2_eq_2 numeral_3_eq_3)
next
  assume *: "finite A  ( x y z. x  y  x  z  y  z  x  A  y  A  z  A)"
  thus "card A  3"
    by (metis One_nat_def Suc_le_eq card_2_iff' card_le_Suc0_iff_eq leI numeral_3_eq_3 one_add_one order_class.order.eq_iff plus_1_eq_Suc)
qed

text ‹Set cardinality of A is equal to 2 if and only if A={x, y} for two different elements x and y›
lemma card_eq_2_iff_doubleton: "card A = 2  ( x y. x  y  A = {x, y})"
  using card_geq_2_iff_contains_2_elems[of A]
  using card_geq_3_iff_contains_3_elems[of A]
  by auto (rule_tac x=x in exI, rule_tac x=y in exI, auto)

lemma card_eq_2_doubleton:
  assumes "card A = 2" and "x  y" and "x  A" and "y  A"
  shows "A = {x, y}"
  using assms card_eq_2_iff_doubleton[of A]
  by auto

text ‹Bijections map singleton to singleton sets›

lemma bij_image_singleton:
  shows "f ` A = {b}; f a = b; bij f  A = {a}"
  by (metis bij_betw_def image_empty image_insert inj_image_eq_iff)

end

Theory Linear_Systems

(* ---------------------------------------------------------------------------- *)
subsection ‹Systems of linear equations›
(* ---------------------------------------------------------------------------- *)
(* TODO: merge with matrices *)

text ‹In this section some simple properties of systems of linear equations with two or three unknowns are derived.
Existence and uniqueness of solutions of regular and singular homogenous and non-homogenous systems is characterized.›

theory Linear_Systems
imports Main
begin

text ‹Determinant of 2x2 matrix›
definition det2 :: "('a::field)  'a  'a  'a  'a" where
  [simp]: "det2 a11 a12 a21 a22  a11*a22 - a12*a21"

text ‹Regular homogenous system has only trivial solution›
lemma regular_homogenous_system:
  fixes a11 a12 a21 a22 x1 x2 :: "'a::field"
  assumes "det2 a11 a12 a21 a22  0"
  assumes "a11*x1 + a12*x2 = 0" and
          "a21*x1 + a22*x2 = 0"
  shows "x1 = 0  x2 = 0"
proof (cases "a11 = 0")
  case True
  with assms(1) have "a12  0" "a21  0"
    by auto
  thus ?thesis
    using a11 = 0 assms(2) assms(3)
    by auto
next
  case False
  hence "x1 = - a12*x2 / a11"
    using assms(2)
    by (metis eq_neg_iff_add_eq_0 mult_minus_left nonzero_mult_div_cancel_left)
  hence "a21 * (- a12 * x2 / a11) + a22 * x2 = 0"
    using assms(3)
    by simp
  hence "a21 * (- a12 * x2) + a22 * x2 * a11 = 0"
    using  a11  0
    by auto
  hence "(a11*a22 - a12*a21)*x2 = 0"
    by (simp add: field_simps)
  thus ?thesis
    using assms(1) assms(2) a11  0
    by auto
qed

text ‹Regular system has a unique solution›
lemma regular_system:
  fixes a11 a12 a21 a22 b1 b2 :: "'a::field"
  assumes "det2 a11 a12 a21 a22  0"
  shows "∃! x. a11*(fst x) + a12*(snd x) = b1 
               a21*(fst x) + a22*(snd x) = b2"
proof
  let ?d = "a11*a22 - a12*a21" and ?d1 = "b1*a22 - b2*a12" and ?d2 = "b2*a11 - b1*a21"
  let ?x = "(?d1 / ?d, ?d2 / ?d)"
  have "a11 * ?d1 + a12 * ?d2 = b1*?d" "a21 * ?d1 + a22 * ?d2 = b2*?d"
    by (auto simp add: field_simps)
  thus "a11 * fst ?x + a12 * snd ?x = b1  a21 * fst ?x + a22 * snd ?x = b2"
    using assms
    by (metis (hide_lams, no_types) det2_def add_divide_distrib eq_divide_imp fst_eqD snd_eqD times_divide_eq_right)

  fix x'
  assume "a11 * fst x' + a12 * snd x' = b1  a21 * fst x' + a22 * snd x' = b2"
  with a11 * fst ?x + a12 * snd ?x = b1  a21 * fst ?x + a22 * snd ?x = b2
  have "a11 * (fst x' - fst ?x) + a12 * (snd x' - snd ?x) = 0  a21 * (fst x' - fst ?x) + a22 * (snd x' - snd ?x) = 0"
    by (auto simp add: field_simps)
  thus "x' = ?x"
    using regular_homogenous_system[OF assms, of "fst x' - fst ?x" "snd x' - snd ?x"]
    by (cases x') auto
qed

text ‹Singular system does not have a unique solution›
lemma singular_system:
  fixes a11 a12 a21 a22 ::"'a::field"
  assumes "det2 a11 a12 a21 a22 = 0" and "a11  0  a12  0"
  assumes x0: "a11*fst x0 + a12*snd x0 = b1"
              "a21*fst x0 + a22*snd x0 = b2"
  assumes x: "a11*fst x + a12*snd x = b1"
  shows "a21*fst x + a22*snd x = b2"
proof (cases "a11 = 0")
  case True
  with assms have "a21 = 0" "a12  0"
    by auto
  let ?k = "a22 / a12"
  have "b2 = ?k * b1"
    using x0 a11 = 0 a21 = 0 a12  0
    by auto
  thus ?thesis
    using a11 = 0 a21 = 0 a12  0 x
    by auto
next
  case False
  let ?k = "a21 / a11"
  from x
  have "?k * a11 * fst x + ?k * a12 * snd x = ?k * b1"
    using a11  0
    by (auto simp add: field_simps)
  moreover
  have "a21 = ?k * a11" "a22 = ?k * a12" "b2 = ?k * b1"
    using assms(1) x0 a11  0
    by (auto simp add: field_simps)
  ultimately
  show ?thesis
    by auto
qed

text ‹All solutions of a homogenous system of 2 equations with 3 unknows are proportional›
lemma linear_system_homogenous_3_2:
  fixes a11 a12 a13 a21 a22 a23 x1 y1 z1 x2 y2 z2 :: "'a::field"
  assumes "f1 = (λ x y z. a11 * x + a12 * y + a13 * z)"
  assumes "f2 = (λ x y z. a21 * x + a22 * y + a23 * z)"
  assumes "f1 x1 y1 z1 = 0" and "f2 x1 y1 z1 = 0"
  assumes "f1 x2 y2 z2 = 0" and "f2 x2 y2 z2 = 0"
  assumes "x2  0  y2  0  z2  0"
  assumes "det2 a11 a12 a21 a22  0  det2 a11 a13 a21 a23  0  det2 a12 a13 a22 a23  0"
  shows " k. x1 = k * x2  y1 = k * y2  z1 = k * z2"
proof-
  let ?Dz = "det2 a11 a12 a21 a22"
  let ?Dy = "det2 a11 a13 a21 a23"
  let ?Dx = "det2 a12 a13 a22 a23"

  have "a21 * (f1 x1 y1 z1) - a11 * (f2 x1 y1 z1) = 0"
    using assms
    by simp
  hence yz1: "?Dz*y1 + ?Dy*z1 = 0"
    using assms
    by (simp add: field_simps)

  have "a21 * (f1 x2 y2 z2) - a11 * (f2 x2 y2 z2) = 0"
    using assms
    by simp
  hence yz2: "?Dz*y2 + ?Dy*z2 = 0"
    using assms
    by (simp add: field_simps)
                                     
  have "a22 * (f1 x1 y1 z1) - a12 * (f2 x1 y1 z1) = 0"
    using assms
    by simp                          
  hence xz1: "-?Dz*x1 + ?Dx*z1 = 0"
    using assms
    by (simp add: field_simps)

  have "a22 * (f1 x2 y2 z2) - a12 * (f2 x2 y2 z2) = 0"
    using assms
    by simp                          
  hence xz2: "-?Dz*x2 + ?Dx*z2 = 0"
    using assms
    by (simp add: field_simps)

  have "a23 * (f1 x1 y1 z1) - a13 * (f2 x1 y1 z1) = 0"
    using assms
    by simp                          
  hence xy1: "?Dy*x1 + ?Dx*y1 = 0"
    using assms
    by (simp add: field_simps)

  have "a23 * (f1 x2 y2 z2) - a13 * (f2 x2 y2 z2) = 0"
    using assms
    by simp                          
  hence xy2: "?Dy*x2 + ?Dx*y2 = 0"
    using assms
    by (simp add: field_simps)

  show ?thesis
    using ?Dz  0  ?Dy  0  ?Dx  0
  proof safe
    assume "?Dz  0"
    
    hence *:
      "x2 = (?Dx / ?Dz) * z2" "y2 = - (?Dy / ?Dz) * z2"
      "x1 = (?Dx / ?Dz) * z1" "y1 = - (?Dy / ?Dz) * z1"
      using xy2 xz2 xy1 xz1 yz1 yz2
      by (simp_all add: field_simps)     

    hence "z2  0"
      using x2  0  y2  0  z2  0
      by auto

    thus ?thesis
      using * ?Dz  0
      by (rule_tac x="z1/z2" in exI) auto
  next
    assume "?Dy  0"
    hence *:
      "x2 = - (?Dx / ?Dy) * y2" "z2 = - (?Dz / ?Dy) * y2"
      "x1 = - (?Dx / ?Dy) * y1" "z1 = - (?Dz / ?Dy) * y1"
      using xy2 xz2 xy1 xz1 yz1 yz2
      by (simp_all add: field_simps)     

    hence "y2  0"
      using x2  0  y2  0  z2  0
      by auto

    thus ?thesis
      using * ?Dy  0
      by (rule_tac x="y1/y2" in exI) auto
  next
    assume "?Dx  0"
    hence *:
      "y2 = - (?Dy / ?Dx) * x2" "z2 = (?Dz / ?Dx) * x2"
      "y1 = - (?Dy / ?Dx) * x1" "z1 = (?Dz / ?Dx) * x1"
      using xy2 xz2 xy1 xz1 yz1 yz2
      by (simp_all add: field_simps)     

    hence "x2  0"
      using x2  0  y2  0  z2  0
      by auto

    thus ?thesis
      using * ?Dx  0
      by (rule_tac x="x1/x2" in exI) auto
  qed
qed

end

Theory Quadratic

(* ----------------------------------------------------------------- *)
subsection ‹Quadratic equations›
(* ----------------------------------------------------------------- *)

text ‹In this section some simple properties of quadratic equations and their roots are derived.
Quadratic equations over reals and over complex numbers, but also systems of quadratic equations and
systems of quadratic and linear equations are analysed.›

theory Quadratic
  imports More_Complex "HOL-Library.Quadratic_Discriminant"
begin

(* ----------------------------------------------------------------- *)
subsubsection ‹Real quadratic equations, Viette rules›
(* ----------------------------------------------------------------- *)

lemma viette2_monic:
  fixes b c ξ1 ξ2 :: real
  assumes "b2 - 4*c  0" and "ξ12 + b*ξ1 + c = 0" and "ξ22 + b*ξ2 + c = 0" and "ξ1  ξ2"
  shows "ξ1*ξ2 = c"
  using assms
  by algebra

lemma viette2:
  fixes a b c ξ1 ξ2 :: real
  assumes "a  0" and "b2 - 4*a*c  0" and "a*ξ12 + b*ξ1 + c = 0" and "a*ξ22 + b*ξ2 + c = 0" and "ξ1  ξ2"
  shows "ξ1*ξ2 = c/a"
proof (rule viette2_monic[of "b/a" "c/a" ξ1 ξ2])
  have "(b / a)2 - 4 * (c / a) = (b2 - 4*a*c) / a2"
    using a  0
    by (auto simp add: power2_eq_square field_simps)
  thus "0  (b / a)2 - 4 * (c / a)"
    using b2 - 4*a*c  0
    by simp
next
  show "ξ12 + b / a * ξ1 + c / a = 0" "ξ22 + b / a * ξ2 + c / a = 0"
    using assms
    by (auto simp add: power2_eq_square field_simps)
next
  show "ξ1  ξ2"
    by fact
qed

lemma viette2'_monic:
  fixes b c ξ :: real
  assumes "b2 - 4*c = 0" and "ξ2 + b*ξ + c = 0"
  shows "ξ*ξ = c"
  using assms
  by algebra

lemma viette2':
  fixes a b c ξ :: real
  assumes "a  0" and "b2 - 4*a*c = 0" and "a*ξ2 + b*ξ + c = 0"
  shows "ξ*ξ = c/a"
proof (rule viette2'_monic)
  have "(b / a)2 - 4 * (c / a) = (b2 - 4*a*c) / a2"
    using a  0
    by (auto simp add: power2_eq_square field_simps)
  thus "(b / a)2 - 4 * (c / a) = 0"
    using b2 - 4*a*c = 0
    by simp
next
  show "ξ2 + b / a * ξ + c / a = 0"
    using assms
    by (auto simp add: power2_eq_square field_simps)
qed

(* ----------------------------------------------------------------- *)
subsubsection ‹Complex quadratic equations›
(* ----------------------------------------------------------------- *)

lemma complex_quadratic_equation_monic_only_two_roots:
  fixes ξ :: complex
  assumes "ξ2 + b * ξ + c = 0"
  shows "ξ = (-b + ccsqrt(b2 - 4*c)) / 2  ξ = (-b - ccsqrt(b2 - 4*c)) / 2"
using assms
proof-
  from assms have "(2 * (ξ + b/2))2 = b2 - 4*c"
    by (simp add: power2_eq_square field_simps)
       (metis (no_types, lifting) distrib_right_numeral mult.assoc mult_zero_left)
  hence "2 * (ξ + b/2) = ccsqrt (b2 - 4*c)  2 * (ξ + b/2) = - ccsqrt (b2 - 4*c)"
    using ccsqrt[of "(2 * (ξ + b / 2))" "b2 - 4 * c"]
    by (simp add: power2_eq_square)
  thus ?thesis
    using mult_cancel_right[of "b + ξ * 2" 2 "ccsqrt (b2 - 4*c)"]
    using mult_cancel_right[of "b + ξ * 2" 2 "-ccsqrt (b2 - 4*c)"]
    by (auto simp add: field_simps) (metis add_diff_cancel diff_minus_eq_add minus_diff_eq)
qed

lemma complex_quadratic_equation_monic_roots:
  fixes ξ :: complex
  assumes "ξ = (-b + ccsqrt(b2 - 4*c)) / 2 
           ξ = (-b - ccsqrt(b2 - 4*c)) / 2"
  shows  "ξ2 + b * ξ + c = 0"
using assms
proof
  assume *: "ξ = (- b + ccsqrt (b2 - 4 * c)) / 2"
  show ?thesis
    by ((subst *)+) (subst power_divide, subst power2_sum, simp add: field_simps, simp add: power2_eq_square)
next
  assume *: "ξ = (- b - ccsqrt (b2 - 4 * c)) / 2"
  show ?thesis
    by ((subst *)+, subst power_divide, subst power2_diff, simp add: field_simps, simp add: power2_eq_square)
qed

lemma complex_quadratic_equation_monic_distinct_roots:
  fixes b c :: complex
  assumes "b2 - 4*c  0"
  shows " k1 k2. k1  k2  k12 + b*k1 + c = 0  k22 + b*k2 + c = 0"
proof-
  let ?ξ1 = "(-b + ccsqrt(b2 - 4*c)) / 2"
  let ?ξ2 = "(-b - ccsqrt(b2 - 4*c)) / 2"
  show ?thesis
    apply (rule_tac x="?ξ1" in exI)
    apply (rule_tac x="?ξ2" in exI)
    using assms 
    using complex_quadratic_equation_monic_roots[of ?ξ1 b c]
    using complex_quadratic_equation_monic_roots[of ?ξ2 b c]
    by simp
qed

lemma complex_quadratic_equation_two_roots:
  fixes ξ :: complex
  assumes "a  0" and "a*ξ2 + b * ξ + c = 0"
  shows "ξ = (-b + ccsqrt(b2 - 4*a*c)) / (2*a) 
         ξ = (-b - ccsqrt(b2 - 4*a*c)) / (2*a)"
proof-
  from assms have "ξ2 + (b/a) * ξ + (c/a) = 0"
    by (simp add: field_simps)
  hence "ξ = (-(b/a) + ccsqrt((b/a)2 - 4*(c/a))) / 2  ξ = (-(b/a) - ccsqrt((b/a)2 - 4*(c/a))) / 2"
    using complex_quadratic_equation_monic_only_two_roots[of ξ "b/a" "c/a"]
    by simp
  hence " k. ξ = (-(b/a) + (-1)^k * ccsqrt((b/a)2 - 4*(c/a))) / 2"
    by safe (rule_tac x="2" in exI, simp, rule_tac x="1" in exI, simp)
  then obtain k1 where "ξ = (-(b/a) + (-1)^k1 * ccsqrt((b/a)2 - 4*(c/a))) / 2"
    by auto
  moreover
  have "(b / a)2 - 4 * (c / a) = (b2 - 4 * a * c) * (1 / a2)"
    using a  0
    by (simp add: field_simps power2_eq_square)
  hence "ccsqrt ((b / a)2 - 4 * (c / a)) = ccsqrt (b2 - 4 * a * c) * ccsqrt (1/a2) 
         ccsqrt ((b / a)2 - 4 * (c / a)) = - ccsqrt (b2 - 4 * a * c) * ccsqrt (1/a2)"
    using ccsqrt_mult[of "b2 - 4 * a * c" "1/a2"]
    by auto
  hence " k.  ccsqrt ((b / a)2 - 4 * (c / a)) = (-1)^k * ccsqrt (b2 - 4 * a * c) * ccsqrt (1 / a2)"
    by safe (rule_tac x="2" in exI, simp, rule_tac x="1" in exI, simp)
  then obtain k2 where "ccsqrt ((b / a)2 - 4 * (c / a)) = (-1)^k2 * ccsqrt (b2 - 4 * a * c) * ccsqrt (1 / a2)"
    by auto
  moreover
  have "ccsqrt (1 / a2) = 1/a  ccsqrt (1 / a2) = -1/a"
    using ccsqrt[of "1/a" "1 / a2"]
    by (auto simp add: power2_eq_square)
  hence " k. ccsqrt (1 / a2) = (-1)^k * 1/a"
    by safe (rule_tac x="2" in exI, simp, rule_tac x="1" in exI, simp)
  then obtain k3 where "ccsqrt (1 / a2) = (-1)^k3 * 1/a"
    by auto
  ultimately
  have "ξ = (- (b / a) + ((-1) ^ k1 * (-1) ^ k2 * (-1) ^ k3) * ccsqrt (b2 - 4 * a * c) * 1/a) / 2"
    by simp
  moreover
  have "(-(1::complex)) ^ k1 * (-1) ^ k2 * (-1) ^ k3 = 1  (-(1::complex)) ^ k1 * (-1) ^ k2 * (-1) ^ k3 = -1"
    using neg_one_even_power[of "k1 + k2 + k3"]
    using neg_one_odd_power[of "k1 + k2 + k3"]
    by (smt power_add)
  ultimately
  have "ξ = (- (b / a) + ccsqrt (b2 - 4 * a * c) * 1 / a) / 2  ξ = (- (b / a) - ccsqrt (b2 - 4 * a * c) * 1 / a) / 2"
    by auto
  thus ?thesis
    using a  0
    by (simp add: field_simps)
qed

lemma complex_quadratic_equation_only_two_roots:
  fixes x :: complex
  assumes "a  0"
  assumes "qf = (λ x. a*x2 + b*x + c)"
          "qf x1 = 0" and "qf x2 = 0" and "x1  x2"
          "qf x = 0"
  shows "x = x1  x = x2"
  using assms
  using complex_quadratic_equation_two_roots
  by blast


(* ----------------------------------------------------------------- *)
subsubsection ‹Intersections of linear and quadratic forms›
(* ----------------------------------------------------------------- *)
(* These lemmas are not used *)

lemma quadratic_linear_at_most_2_intersections_help:
  fixes x y :: complex
  assumes "(a11, a12, a22)  (0, 0, 0)" and "k2  0"
          "qf = (λ x y. a11*x2 + 2*a12*x*y + a22*y2 + b1*x + b2*y + c)" and "lf = (λ x y. k1*x + k2*y + n)"
          "qf x y = 0" and "lf x y = 0"
          "pf = (λ x. (a11 - 2*a12*k1/k2 + a22*k12/k22)*x2 + (-2*a12*n/k2  + b1 + a22*2*n*k1/k22 - b2*k1/k2)*x + a22*n2/k22 - b2*n/k2 + c)"
          "yf = (λ x. (-n - k1*x) / k2)"
  shows "pf x = 0" and "y = yf x"
proof -
  show "y = yf x"
    using assms
    by (simp add:field_simps eq_neg_iff_add_eq_0)
next
  have "2*a12*x*(-n - k1*x)/k2 = (-2*a12*n/k2)*x - (2*a12*k1/k2)*x2"
    by algebra
  have "a22*((-n - k1*x)/k2)2 = a22*n2/k22 + (a22*2*n*k1/k22)*x + (a22*k12/k22)*x2"
    by (simp add: power_divide) algebra
  have "2*a12*x*(-n - k1*x)/k2 = (-2*a12*n/k2)*x - (2*a12*k1/k2)*x2"
    by algebra
  have "b2*(-n - k1*x)/k2 = -b2*n/k2 - (b2*k1/k2)*x"
    by algebra

  have *: "y = (-n - k1*x)/k2"
    using assms(2, 4, 6)
    by (simp add:field_simps eq_neg_iff_add_eq_0)

  have "0 = a11*x2 + 2*a12*x*y + a22*y2 + b1*x + b2*y + c"
    using assms
    by simp
  hence "0 = a11*x2 + 2*a12*x*(-n - k1*x)/k2 + a22*((-n - k1*x)/k2)2 + b1*x + b2*(-n - k1*x)/k2 + c"
    by (subst (asm) *, subst (asm) *, subst (asm) *) auto
  also have "... = (a11 - 2*a12*k1/k2 + a22*k12/k22)*x2 + (-2*a12*n/k2  + b1 + a22*2*n*k1/k22 - b2*k1/k2)*x + a22*n2/k22 -b2*n/k2 + c"
    using 2*a12*x*(-n - k1*x)/k2 = (-2*a12*n/k2)*x - (2*a12*k1/k2)*x2
    using a22*((-n - k1*x)/k2)2 = a22*n2/k22 + (a22*2*n*k1/k22)*x + (a22*k12/k22)*x2
    using b2*(-n - k1*x)/k2 = -b2*n/k2 - (b2*k1/k2)*x
    by (simp add:field_simps)
  finally show "pf x = 0"
    using assms(7)
    by auto
qed

lemma quadratic_linear_at_most_2_intersections_help':
  fixes x y :: complex
  assumes "qf = (λ x y. a11*x2 + 2*a12*x*y + a22*y2 + b1*x + b2*y + c)"
          "x = -n/k1" and "k1  0" and "qf x y = 0"
          "yf = (λ y. k12*a22*y2 + (-2*a12*n*k1 + b2*k12)*y + a11*n2 - b1*n*k1 + c*k12)"
  shows "yf y = 0"
proof-
  have "0 = a11*n2/k12 - 2*a12*n*y/k1 + a22*y2 - b1*n/k1 + b2*y + c"
    using assms(1, 2, 4)
    by (simp add: power_divide)
  hence "0 = a11*n2 - 2*a12*n*k1*y + a22*y2*k12 - b1*n*k1 + b2*y*k12 + c*k12"
    using assms(3)
    apply (simp add:field_simps power2_eq_square)
    by algebra
  thus ?thesis
    using assms(1, 4, 5)
    by (simp add:field_simps)
qed

lemma quadratic_linear_at_most_2_intersections:
  fixes x y x1 y1 x2 y2 :: complex
  assumes "(a11, a12, a22)  (0, 0, 0)" and "(k1, k2)  (0, 0)"
  assumes "a11*k22 - 2*a12*k1*k2 + a22*k12  0"
  assumes "qf = (λ x y. a11*x2 + 2*a12*x*y + a22*y2 + b1*x + b2*y + c)" and "lf = (λ x y. k1*x + k2*y + n)"
          "qf x1 y1 = 0" and "lf x1 y1 = 0"
          "qf x2 y2 = 0" and "lf x2 y2 = 0"
          "(x1, y1)  (x2, y2)"
          "qf x y = 0" and "lf x y = 0"
  shows "(x, y) = (x1, y1)  (x, y) = (x2, y2)"
proof(cases "k2 = 0")
  case True
  hence "k1  0"
    using assms(2)
    by simp

  have "a22*k12  0"
    using assms(3) True
    by auto

  have "x1 = -n/k1"
    using k1  0 assms(5, 7) True
    by (metis add.right_neutral add_eq_0_iff2 mult_zero_left nonzero_mult_div_cancel_left)
  have "x2 = -n/k1"
    using k1  0 assms(5, 9) True
    by (metis add.right_neutral add_eq_0_iff2 mult_zero_left nonzero_mult_div_cancel_left)
  have "x = -n/k1"
    using k1  0 assms(5, 12) True
    by (metis add.right_neutral add_eq_0_iff2 mult_zero_left nonzero_mult_div_cancel_left)

  let ?yf =  "(λ y. k12*a22*y2 + (-2*a12*n*k1 + b2*k12)*y + a11*n2 - b1*n*k1 + c*k12)"

  have "?yf y = 0"
    using quadratic_linear_at_most_2_intersections_help'[of qf a11 a12 a22 b1 b2 c x n k1 y ?yf]
    using assms(4, 11) k1  0 x = -n/k1
    by auto
  have "?yf y1 = 0"
    using quadratic_linear_at_most_2_intersections_help'[of qf a11 a12 a22 b1 b2 c x1 n k1 y1 ?yf]
    using assms(4, 6) k1  0 x1 = -n/k1
    by auto
  have "?yf y2 = 0"
    using quadratic_linear_at_most_2_intersections_help'[of qf a11 a12 a22 b1 b2 c x2 n k1 y2 ?yf]
    using assms(4, 8) k1  0 x2 = -n/k1
    by auto

  have "y1  y2"
    using assms(10) x1 = -n/k1 x2 = -n/k1
    by blast

  have "y = y1  y = y2"
    using complex_quadratic_equation_only_two_roots[of "a22*k12" ?yf "-2*a12*n*k1 + b2*k12" "a11*n2 - b1*n*k1 + c*k12"
                                                        y1 y2 y]
    using a22*k12  0 ?yf y1 = 0 y1  y2 ?yf y2 = 0 ?yf y = 0
    by fastforce

  thus ?thesis
    using x1 = -n/k1 x2 = -n/k1  x = -n/k1
    by auto
next
  case False

  let ?py = "(λ x. (-n - k1*x)/k2)"
  let ?pf = "(λ x. (a11 - 2*a12*k1/k2 + a22*k12/k22)*x2 + (-2*a12*n/k2  + b1 + a22*2*n*k1/k22 - b2*k1/k2)*x + a22*n2/k22 -b2*n/k2 + c)"
  have "?pf x1 = 0" "y1 = ?py x1"
    using quadratic_linear_at_most_2_intersections_help[of a11 a12 a22 k2 qf b1 b2 c lf k1 n x1 y1]
    using assms(1, 4, 5, 6, 7) False
    by auto
  have "?pf x2 = 0" "y2 = ?py x2"
    using quadratic_linear_at_most_2_intersections_help[of a11 a12 a22 k2 qf b1 b2 c lf k1 n x2 y2]
    using assms(1, 4, 5, 8, 9) False
    by auto
  have "?pf x = 0" "y = ?py x"
    using quadratic_linear_at_most_2_intersections_help[of a11 a12 a22 k2 qf b1 b2 c lf k1 n x y]
    using assms(1, 4, 5, 11, 12) False
    by auto

  have "x1  x2"
    using assms(10) y1 = ?py x1 y2 = ?py x2
    by auto

  have "a11 - 2*a12*k1/k2 + a22*k12/k22 = (a11 * k22 - 2 * a12 * k1 * k2 + a22 * k12)/k22"
    by (simp add: False power2_eq_square add_divide_distrib diff_divide_distrib)
  also have "...   0"
    using False assms(3)
    by simp
  finally have "a11 - 2*a12*k1/k2 + a22*k12/k22  0"
    .

  have "x = x1  x = x2"
    using complex_quadratic_equation_only_two_roots[of "a11 - 2*a12*k1/k2 + a22*k12/k22" ?pf
                                                       "(- 2 * a12 * n / k2 + b1 + a22 * 2 * n * k1 / k22 - b2 * k1 / k2)"
                                                       "a22 * n2 / k22 - b2 * n / k2 + c" x1 x2 x]
    using ?pf x2 = 0 ?pf x1 = 0 ?pf x = 0
    using a11 - 2 * a12 * k1 / k2 + a22 * k12 / k22  0
    using x1  x2
    by fastforce

  thus ?thesis
    using y = ?py x y1 = ?py x1 y2 = ?py x2
    by (cases "x = x1", auto)
qed

lemma quadratic_quadratic_at_most_2_intersections':
  fixes x y x1 y1 x2 y2 :: complex
  assumes "b2  B2  b1  B1"
          "(b2 - B2)2 + (b1 - B1)2  0"
  assumes "qf1 = (λ x y. x2 + y2 + b1*x + b2*y + c)"
          "qf2 = (λ x y. x2 + y2 + B1*x + B2*y + C)"
          "qf1 x1 y1 = 0" "qf2 x1 y1 = 0"
          "qf1 x2 y2 = 0" "qf2 x2 y2 = 0"
          "(x1, y1)  (x2, y2)"
          "qf1 x y = 0" "qf2 x y = 0"
  shows "(x, y) = (x1, y1)  (x, y) = (x2, y2)"
proof-
  have "x2 + y2 + b1*x + b2*y + c = 0"
    using assms by auto
  have "x2 + y2 + B1*x + B2*y + C = 0"
    using assms by auto
  hence "0 = x2 + y2 + b1*x + b2*y + c - (x2 + y2 + B1*x + B2*y + C)"
    using x2 + y2 + b1*x + b2*y + c = 0
    by auto
  hence "0 = (b1 - B1)*x + (b2 - B2)*y + c - C"
    by (simp  add:field_simps) 
  
  have "x12 + y12 + b1*x1 + b2*y1 + c = 0"
    using assms  by auto
  have "x12 + y12 + B1*x1 + B2*y1 + C = 0"
    using assms by auto
  hence "0 = x12 + y12 + b1*x1 + b2*y1 + c - (x12 + y12 + B1*x1 + B2*y1 + C)"
    using x12 + y12 + b1*x1 + b2*y1 + c = 0
    by auto
  hence "0 = (b1 - B1)*x1 + (b2 - B2)*y1 + c - C"
    by (simp  add:field_simps) 

  have "x22 + y22 + b1*x2 + b2*y2 + c = 0"
    using assms  by auto
  have "x22 + y22 + B1*x2 + B2*y2 + C = 0"
    using assms  by auto
  hence "0 = x22 + y22 + b1*x2 + b2*y2 + c - (x22 + y22 + B1*x2 + B2*y2 + C)"
    using x22 + y22 + b1*x2 + b2*y2 + c = 0
    by auto
  hence "0 = (b1 - B1)*x2 + (b2 - B2)*y2 + c - C"
    by (simp  add:field_simps)  

  have "(b1 - B1, b2 - B2)  (0, 0)"
    using assms(1) by auto

  let ?lf = "(λ x y. (b1 - B1)*x + (b2 - B2)*y + c - C)"

  have "?lf x y = 0" "?lf x1 y1 = 0" "?lf x2 y2 = 0"
    using 0 = (b1 - B1)*x2 + (b2 - B2)*y2 + c - C
          0 = (b1 - B1)*x1 + (b2 - B2)*y1 + c - C
          0 = (b1 - B1)*x + (b2 - B2)*y + c - C
    by auto

  thus ?thesis
    using quadratic_linear_at_most_2_intersections[of 1 0 1 "b1 - B1" "b2 - B2" qf1 b1 b2 c ?lf "c - C" x1 y1 x2 y2 x y]
    using (b1 - B1, b2 - B2)  (0, 0)
    using assms (b1 - B1, b2 - B2)  (0, 0)
    using (b1 - B1) * x + (b2 - B2) * y + c - C = 0 (b1 - B1) * x1 + (b2 - B2) * y1 + c - C = 0
    by (simp add: add_diff_eq)
qed

lemma quadratic_change_coefficients:
  fixes x y :: complex
  assumes "A1  0" 
  assumes "qf = (λ x y. A1*x2 + A1*y2 + b1*x + b2*y + c)"
          "qf x y = 0"
          "qf_1 = (λ x y. x2 + y2 + (b1/A1)*x + (b2/A1)*y + c/A1)"
  shows "qf_1 x y = 0"
proof-
  have "0 = A1*x2 + A1*y2 + b1*x + b2*y + c"
    using assms by auto
  hence "0/A1 = (A1*x2 + A1*y2 + b1*x + b2*y + c)/A1"
    using assms(1) by auto
  also have "... = A1*x2/A1 + A1*y2/A1 + b1*x/A1 + b2*y/A1 + c/A1"
    by (simp add: add_divide_distrib)
  also have "... = x2 + y2 + (b1/A1)*x + (b2/A1)*y + c/A1"
    using assms(1)
    by (simp add:field_simps)
  finally have "0 = x2 + y2 + (b1/A1)*x + (b2/A1)*y + c/A1"
    by simp
  thus ?thesis
    using assms
    by simp
qed

lemma quadratic_quadratic_at_most_2_intersections:
  fixes x y x1 y1 x2 y2 :: complex
  assumes "A1  0" and "A2  0"
  assumes "qf1 = (λ x y. A1*x2 + A1*y2 + b1*x + b2*y + c)" and
          "qf2 = (λ x y. A2*x2 + A2*y2 + B1*x + B2*y + C)" and
          "qf1 x1 y1 = 0" and "qf2 x1 y1 = 0" and
          "qf1 x2 y2 = 0" and "qf2 x2 y2 = 0" and
          "(x1, y1)  (x2, y2)" and
          "qf1 x y = 0" and "qf2 x y = 0"
  assumes "(b2*A2 - B2*A1)2 + (b1*A2 - B1*A1)2  0" and
          "b2*A2  B2*A1  b1*A2  B1*A1"
  shows "(x, y) = (x1, y1)  (x, y) = (x2, y2)"
proof-
  have *: "b2 / A1  B2 / A2  b1 / A1  B1 / A2"
    using assms(1, 2) assms(13)
    by (simp add:field_simps)
  have **: "(b2 / A1 - B2 / A2)2 + (b1 / A1 - B1 / A2)2  0"
    using assms(1, 2) assms(12)
    by (simp add:field_simps)

  let ?qf_1 = "(λ x y. x2 + y2 + (b1/A1)*x + (b2/A1)*y + c/A1)"
  let ?qf_2 = "(λ x y. x2 + y2 + (B1/A2)*x + (B2/A2)*y + C/A2)"

  have "?qf_1 x1 y1 = 0" "?qf_1 x2 y2 = 0" "?qf_1 x y = 0"
       "?qf_2 x1 y1 = 0" "?qf_2 x2 y2 = 0" "?qf_2 x y = 0"
    using assms quadratic_change_coefficients[of A1 qf1 b1 b2 c x2 y2 ?qf_1]
          quadratic_change_coefficients[of A1 qf1 b1 b2 c x1 y1 ?qf_1]
          quadratic_change_coefficients[of A2 qf2 B1 B2 C x1 y1 ?qf_2]
          quadratic_change_coefficients[of A2 qf2 B1 B2 C x2 y2 ?qf_2]
          quadratic_change_coefficients[of A1 qf1 b1 b2 c x y ?qf_1]
          quadratic_change_coefficients[of A2 qf2 B1 B2 C x y ?qf_2]
    by auto

  thus ?thesis
    using quadratic_quadratic_at_most_2_intersections'
              [of "b2 / A1" "B2 / A2" "b1 / A1" "B1 / A2" ?qf_1 "c / A1" ?qf_2 "C / A2" x1 y1 x2 y2 x y]
    using * ** (x1, y1)  (x2, y2)
    by fastforce
qed

end

Theory Matrices

(* ---------------------------------------------------------------------------- *)
subsection ‹Vectors and Matrices in $\mathbb{C}^2$›
(* ---------------------------------------------------------------------------- *)

text ‹Representing vectors and matrices of arbitrary dimensions pose a challenge in formal theorem
proving \cite{harrison05}, but we only need to consider finite dimension spaces $\mathbb{C}^2$ and
$\mathbb{R}^3$.›

theory Matrices
imports More_Complex Linear_Systems Quadratic
begin

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Vectors in $\mathbb{C}^2$›
(* ---------------------------------------------------------------------------- *)

text ‹Type of complex vector›

type_synonym complex_vec = "complex × complex"

definition vec_zero :: "complex_vec" where
  [simp]: "vec_zero = (0, 0)"

text ‹Vector scalar multiplication›

fun mult_sv :: "complex  complex_vec  complex_vec" (infixl "*sv" 100) where
  "k *sv (x, y) = (k*x, k*y)"

lemma fst_mult_sv [simp]: 
  shows "fst (k *sv v) = k * fst v"
  by (cases v) simp

lemma snd_mult_sv [simp]:
  shows "snd (k *sv v) = k * snd v"
  by (cases v) simp

lemma mult_sv_mult_sv [simp]: 
  shows "k1 *sv (k2 *sv v) = (k1*k2) *sv v"
  by (cases v) simp

lemma one_mult_sv [simp]:
  shows "1 *sv v =  v"
  by (cases v) simp

lemma mult_sv_ex_id1 [simp]:
  shows " k::complex. k  0  k *sv v = v"
  by (rule_tac x=1 in exI, simp)

lemma mult_sv_ex_id2 [simp]:
  shows " k::complex. k  0  v = k *sv v"
  by (rule_tac x=1 in exI, simp)

text ‹Scalar product of two vectors›

fun mult_vv :: "complex × complex  complex × complex  complex" (infixl "*vv" 100) where
 "(x, y) *vv (a, b) = x*a + y*b"

lemma mult_vv_commute:
  shows "v1 *vv v2 = v2 *vv v1"
  by (cases v1, cases v2) auto

lemma mult_vv_scale_sv1:
  shows "(k *sv v1) *vv v2 = k * (v1 *vv v2)"
  by (cases v1, cases v2) (auto simp add: field_simps)

lemma mult_vv_scale_sv2:
  shows "v1 *vv (k *sv v2) = k * (v1 *vv v2)"
  by (cases v1, cases v2) (auto simp add: field_simps)

text ‹Conjugate vector›

fun vec_map where
 "vec_map f (x, y) = (f x, f y)"

definition vec_cnj where
  "vec_cnj = vec_map cnj"

lemma vec_cnj_vec_cnj [simp]:
  shows "vec_cnj (vec_cnj v) = v"
  by (cases v) (simp add: vec_cnj_def)

lemma cnj_mult_vv:
  shows "cnj (v1 *vv v2) = (vec_cnj v1) *vv (vec_cnj v2)"
  by (cases v1, cases v2) (simp add: vec_cnj_def)

lemma vec_cnj_sv [simp]:
  shows "vec_cnj (k *sv A) = cnj k *sv vec_cnj A"
  by (cases A) (auto simp add: vec_cnj_def)

lemma scalsquare_vv_zero:
  shows "(vec_cnj v) *vv v = 0  v = vec_zero"
  apply (cases v)
  apply (auto simp add: vec_cnj_def field_simps complex_mult_cnj_cmod power2_eq_square)
  apply (metis (no_types) norm_eq_zero of_real_0 of_real_add of_real_eq_iff of_real_mult sum_squares_eq_zero_iff)+
  done

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Matrices in $\mathbb{C}^2$›
(* ---------------------------------------------------------------------------- *)

text ‹Type of complex matrices›

type_synonym complex_mat = "complex × complex × complex × complex"

text ‹Matrix scalar multiplication›

fun mult_sm :: "complex  complex_mat  complex_mat" (infixl "*sm" 100) where
  "k *sm (a, b, c, d) = (k*a, k*b, k*c, k*d)"

lemma mult_sm_distribution [simp]:
  shows "k1 *sm (k2 *sm A) = (k1*k2) *sm A"
  by (cases A) auto

lemma mult_sm_neutral [simp]:
  shows "1 *sm A = A"
  by (cases A) auto

lemma mult_sm_inv_l:
  assumes "k  0" and "k *sm A = B"
  shows "A = (1/k) *sm B"
  using assms
  by auto

lemma mult_sm_ex_id1 [simp]:
  shows " k::complex. k  0  k *sm M = M"
  by (rule_tac x=1 in exI, simp)

lemma mult_sm_ex_id2 [simp]:
  shows " k::complex. k  0  M = k *sm M"
  by (rule_tac x=1 in exI, simp)

text ‹Matrix addition and subtraction›

definition mat_zero :: "complex_mat" where [simp]: "mat_zero = (0, 0, 0, 0)"

fun mat_plus :: "complex_mat  complex_mat  complex_mat" (infixl "+mm" 100) where
  "mat_plus (a1, b1, c1, d1) (a2, b2, c2, d2) = (a1+a2, b1+b2, c1+c2, d1+d2)"

fun mat_minus :: "complex_mat  complex_mat  complex_mat" (infixl "-mm" 100) where
  "mat_minus (a1, b1, c1, d1) (a2, b2, c2, d2) = (a1-a2, b1-b2, c1-c2, d1-d2)"

fun mat_uminus :: "complex_mat  complex_mat" where
  "mat_uminus (a, b, c, d) = (-a, -b, -c, -d)"

lemma nonzero_mult_real:
  assumes "A  mat_zero" and "k  0"
  shows "k *sm A  mat_zero"
  using assms
  by (cases A) simp

text ‹Matrix multiplication.›

fun mult_mm :: "complex_mat  complex_mat  complex_mat" (infixl "*mm" 100) where
  "(a1, b1, c1, d1) *mm (a2, b2, c2, d2) =
   (a1*a2 + b1*c2, a1*b2 + b1*d2, c1*a2+d1*c2, c1*b2+d1*d2)"

lemma mult_mm_assoc:
  shows "A *mm (B *mm C) = (A *mm B) *mm C"
  by (cases A, cases B, cases C) (auto simp add: field_simps)

lemma mult_assoc_5:
  shows "A *mm (B *mm C *mm D) *mm E = (A *mm B) *mm C *mm (D *mm E)"
  by (simp only: mult_mm_assoc)

lemma mat_zero_r [simp]:
  shows "A *mm mat_zero = mat_zero"
  by (cases A) simp

lemma mat_zero_l [simp]:
  shows "mat_zero *mm A = mat_zero"
  by (cases A) simp

definition eye :: "complex_mat" where
  [simp]: "eye = (1, 0, 0, 1)"

lemma mat_eye_l:
  shows "eye *mm A = A"
  by (cases A) auto

lemma mat_eye_r:
  shows "A *mm eye = A"
  by (cases A) auto

lemma mult_mm_sm [simp]:
  shows "A *mm (k *sm B) = k *sm (A *mm B)"
  by (cases A, cases B) (simp add: field_simps)

lemma mult_sm_mm [simp]:
  shows "(k *sm A) *mm B = k *sm (A *mm B)"
  by (cases A, cases B) (simp add: field_simps)

lemma mult_sm_eye_mm [simp]:
  shows "k *sm eye *mm A = k *sm A"
  by (cases A) simp

text ‹Matrix determinant›

fun mat_det where "mat_det (a, b, c, d) = a*d - b*c"

lemma mat_det_mult [simp]:
  shows "mat_det (A *mm B) = mat_det A * mat_det B"
  by (cases A, cases B) (auto simp add: field_simps)

lemma mat_det_mult_sm [simp]:
  shows "mat_det (k *sm A) = (k*k) * mat_det A"
  by (cases A) (auto simp add: field_simps)

text ‹Matrix inverse›

fun mat_inv :: "complex_mat  complex_mat" where
  "mat_inv (a, b, c, d) = (1/(a*d - b*c)) *sm (d, -b, -c, a)"

lemma mat_inv_r:
  assumes "mat_det A  0"
  shows "A *mm (mat_inv A) = eye"
  using assms
proof (cases A, auto simp add: field_simps)
  fix a b c d :: complex
  assume "a * (a * (d * d)) + b * (b * (c * c)) = a * (b * (c * (d * 2)))"
  hence "(a*d - b*c)*(a*d - b*c) = 0"
    by (auto simp add: field_simps)
  hence *: "a*d - b*c = 0"
    by auto
  assume "a*d  b*c"
  with * show False
    by auto
qed

lemma mat_inv_l:
  assumes "mat_det A  0"
  shows "(mat_inv A) *mm A  = eye"
  using assms
proof (cases A, auto simp add: field_simps)
  fix a b c d :: complex
  assume "a * (a * (d * d)) + b * (b * (c * c)) = a * (b * (c * (d * 2)))"
  hence "(a*d - b*c)*(a*d - b*c) = 0"
    by (auto simp add: field_simps)
  hence *: "a*d - b*c = 0"
    by auto
  assume "a*d  b*c"
  with * show False
    by auto
qed

lemma mat_det_inv:
  assumes "mat_det A  0"
  shows "mat_det (mat_inv A) = 1 / mat_det A"
proof-
  have "mat_det eye = mat_det A * mat_det (mat_inv A)"
    using mat_inv_l[OF assms, symmetric]
    by simp
  thus ?thesis
    using assms
    by (simp add: field_simps)
qed

lemma mult_mm_inv_l:
  assumes "mat_det A  0" and "A *mm B = C"
  shows "B = mat_inv A *mm C"
  using assms mat_eye_l[of B]
  by (auto simp add: mult_mm_assoc mat_inv_l)

lemma mult_mm_inv_r:
  assumes "mat_det B  0" and "A *mm B = C"
  shows "A = C *mm mat_inv B"
  using assms mat_eye_r[of A]
  by (auto simp add: mult_mm_assoc[symmetric] mat_inv_r)

lemma mult_mm_non_zero_l:
  assumes "mat_det A  0" and "B  mat_zero"
  shows "A *mm B  mat_zero"
  using assms mat_zero_r
  using mult_mm_inv_l[OF assms(1), of B mat_zero]
  by auto

lemma mat_inv_mult_mm:
  assumes "mat_det A  0" and "mat_det B  0"
  shows "mat_inv (A *mm B) = mat_inv B *mm mat_inv A"
  using assms
proof-
  have "(A *mm B) *mm (mat_inv B *mm mat_inv A) = eye"
    using assms
    by (metis mat_inv_r mult_mm_assoc mult_mm_inv_r)
  thus ?thesis
    using mult_mm_inv_l[of "A *mm B" "mat_inv B *mm mat_inv A" eye] assms mat_eye_r
    by simp
qed

lemma mult_mm_cancel_l:
  assumes "mat_det M  0"  "M *mm A = M *mm B"
  shows "A = B"
  using assms
  by (metis mult_mm_inv_l)

lemma mult_mm_cancel_r:
  assumes "mat_det M  0"  "A *mm M = B *mm M"
  shows "A = B"
  using assms
  by (metis mult_mm_inv_r)

lemma mult_mm_non_zero_r:
  assumes "A  mat_zero" and "mat_det B  0"
  shows "A *mm B  mat_zero"
  using assms mat_zero_l
  using mult_mm_inv_r[OF assms(2), of A mat_zero]
  by auto

lemma mat_inv_mult_sm:
  assumes "k  0"
  shows "mat_inv (k *sm A) = (1 / k) *sm mat_inv A"
proof-
  obtain a b c d where "A = (a, b, c, d)"
    by (cases A) auto
  thus ?thesis
    using assms
    by auto (subst mult.assoc[of k a "k*d"], subst mult.assoc[of k b "k*c"], subst right_diff_distrib[of k "a*(k*d)" "b*(k*c)", symmetric], simp, simp add: field_simps)+
qed

lemma mat_inv_inv [simp]:
  assumes "mat_det M  0"
  shows "mat_inv (mat_inv M) = M"
proof-
  have "mat_inv M *mm M = eye"
    using mat_inv_l[OF assms]
    by simp
  thus ?thesis
    using assms mat_det_inv[of M]
    using mult_mm_inv_l[of "mat_inv M" M eye] mat_eye_r
    by (auto simp del: eye_def)
qed

text ‹Matrix transpose›

fun mat_transpose where
  "mat_transpose (a, b, c, d) = (a, c, b, d)"

lemma mat_t_mat_t [simp]:
  shows "mat_transpose (mat_transpose A) = A"
  by (cases A) auto

lemma mat_t_mult_sm [simp]:
  shows "mat_transpose (k *sm A) = k *sm (mat_transpose A)"
  by (cases A) simp

lemma mat_t_mult_mm [simp]:
  shows "mat_transpose (A *mm B) = mat_transpose B *mm mat_transpose A"
  by (cases A, cases B) auto

lemma mat_inv_transpose:
  shows "mat_transpose (mat_inv M) = mat_inv (mat_transpose M)"
  by (cases M) auto

lemma mat_det_transpose [simp]:
  fixes M :: "complex_mat"
  shows "mat_det (mat_transpose M) = mat_det M"
  by (cases M) auto

text ‹Diagonal matrices definition›

fun mat_diagonal where
 "mat_diagonal (A, B, C, D) = (B = 0  C = 0)"

text ‹Matrix conjugate›

fun mat_map where
 "mat_map f (a, b, c, d) = (f a, f b, f c, f d)"

definition mat_cnj where
  "mat_cnj = mat_map cnj"

lemma mat_cnj_cnj [simp]:
  shows "mat_cnj (mat_cnj A) = A"
  unfolding mat_cnj_def
  by (cases A) auto

lemma mat_cnj_sm [simp]:
  shows "mat_cnj (k *sm A) = cnj k *sm (mat_cnj A)"
  by (cases A) (simp add: mat_cnj_def)

lemma mat_det_cnj [simp]: 
  shows "mat_det (mat_cnj A) = cnj (mat_det A)"
  by (cases A) (simp add: mat_cnj_def)

lemma nonzero_mat_cnj:
  shows "mat_cnj A = mat_zero  A = mat_zero"
  by (cases A) (auto simp add: mat_cnj_def)

lemma mat_inv_cnj:
  shows "mat_cnj (mat_inv M) = mat_inv (mat_cnj M)"
  unfolding mat_cnj_def
  by (cases M) auto

text ‹Matrix adjoint - the conjugate traspose matrix ($A^* = \overline{A^t}$)›

definition mat_adj where
  "mat_adj A = mat_cnj (mat_transpose A)"

lemma mat_adj_mult_mm [simp]:
  shows "mat_adj (A *mm B) = mat_adj B *mm mat_adj A"
  by (cases A, cases B) (auto simp add: mat_adj_def mat_cnj_def)

lemma mat_adj_mult_sm [simp]:
  shows "mat_adj (k *sm A) = cnj k *sm mat_adj A"
  by (cases A) (auto simp add: mat_adj_def mat_cnj_def)

lemma mat_det_adj: 
  shows "mat_det (mat_adj A) = cnj (mat_det A)"
  by (cases A) (auto simp add: mat_adj_def mat_cnj_def)

lemma mat_adj_inv:
  assumes "mat_det M  0"
  shows "mat_adj (mat_inv M) = mat_inv (mat_adj M)"
  by (cases M) (auto simp add: mat_adj_def mat_cnj_def)

lemma mat_transpose_mat_cnj:
  shows "mat_transpose (mat_cnj A) = mat_adj A"
  by (cases A)  (auto simp add: mat_adj_def mat_cnj_def)

lemma mat_adj_adj [simp]:
  shows "mat_adj (mat_adj A) = A"
  unfolding mat_adj_def
  by (subst mat_transpose_mat_cnj) (simp add: mat_adj_def)

lemma mat_adj_eye [simp]:
  shows "mat_adj eye = eye"
  by (auto simp add: mat_adj_def mat_cnj_def)

text ‹Matrix trace›

fun mat_trace where
  "mat_trace (a, b, c, d) = a + d"

text ‹Multiplication of matrix and a vector›

fun mult_mv :: "complex_mat  complex_vec  complex_vec" (infixl "*mv" 100)  where
  "(a, b, c, d) *mv (x, y) = (x*a + y*b, x*c + y*d)"

fun mult_vm :: "complex_vec  complex_mat  complex_vec" (infixl "*vm" 100) where
  "(x, y) *vm (a, b, c, d)  = (x*a + y*c, x*b + y*d)"

lemma eye_mv_l [simp]:
  shows "eye *mv v = v"
  by (cases v) simp

lemma mult_mv_mv [simp]: 
  shows "B *mv (A *mv v) = (B *mm A) *mv v"
  by (cases v, cases A, cases B) (auto simp add: field_simps)

lemma mult_vm_vm [simp]:
  shows "(v *vm A) *vm B = v *vm (A *mm B)"
  by (cases v, cases A, cases B) (auto simp add: field_simps)

lemma mult_mv_inv:
  assumes "x =  A *mv y" and "mat_det A  0"
  shows "y = (mat_inv A) *mv x"
  using assms
  by (cases y) (simp add: mat_inv_l)

lemma mult_vm_inv:
  assumes "x =  y *vm A" and "mat_det A  0"
  shows "y = x *vm (mat_inv A) "
  using assms
  by (cases y) (simp add: mat_inv_r)

lemma mult_mv_cancel_l:
  assumes "mat_det A  0" and "A *mv v = A *mv v'"
  shows "v = v'"
  using assms
  using mult_mv_inv
  by blast

lemma mult_vm_cancel_r:
  assumes "mat_det A  0" and "v *vm A = v' *vm A"
  shows "v = v'"
  using assms
  using mult_vm_inv
  by blast

lemma vec_zero_l [simp]:
  shows "A *mv vec_zero = vec_zero"
  by (cases A) simp

lemma vec_zero_r [simp]:
  shows "vec_zero *vm A = vec_zero"
  by (cases A) simp

lemma mult_mv_nonzero:
  assumes "v  vec_zero" and "mat_det A  0"
  shows "A *mv v  vec_zero"
  apply (rule ccontr)
  using assms mult_mv_inv[of vec_zero A v] mat_inv_l vec_zero_l
  by auto

lemma mult_vm_nonzero:
  assumes "v  vec_zero" and "mat_det A  0"
  shows "v *vm A  vec_zero"
  apply (rule ccontr)
  using assms mult_vm_inv[of vec_zero v A] mat_inv_r vec_zero_r
  by auto

lemma mult_sv_mv:
  shows "k *sv (A *mv v) = (A *mv (k *sv v))"
  by (cases A, cases v) (simp add: field_simps)

lemma mult_mv_mult_vm: 
  shows "A *mv x = x *vm (mat_transpose A)"
  by (cases A, cases x) auto

lemma mult_mv_vv:
  shows "A *mv v1 *vv v2 = v1 *vv (mat_transpose A *mv v2)"
  by (cases v1, cases v2, cases A) (auto simp add: field_simps)

lemma mult_vv_mv:
  shows "x *vv (A *mv y)  = (x *vm A) *vv y"
  by (cases x, cases y, cases A) (auto simp add: field_simps)

lemma vec_cnj_mult_mv:
  shows "vec_cnj (A *mv x) =  (mat_cnj A) *mv (vec_cnj x)"
  by (cases A, cases x) (auto simp add: vec_cnj_def mat_cnj_def)

lemma vec_cnj_mult_vm:
  shows "vec_cnj (v *vm A) = vec_cnj v *vm mat_cnj A"
  unfolding vec_cnj_def mat_cnj_def
  by (cases A, cases v, auto)

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Eigenvalues and eigenvectors›
(* ---------------------------------------------------------------------------- *)

definition eigenpair where
  [simp]: "eigenpair k v H  v  vec_zero  H *mv v = k *sv v"

definition eigenval where
  [simp]: "eigenval k H  ( v. v  vec_zero  H *mv v = k *sv v)"

lemma eigen_equation:
  shows "eigenval k H  k2 - mat_trace H * k + mat_det H = 0" (is "?lhs  ?rhs")
proof-
  obtain A B C D where HH: "H = (A, B, C, D)"
    by (cases H) auto
  show ?thesis
  proof
    assume ?lhs
    then obtain v where "v  vec_zero" "H *mv v = k *sv v"
      unfolding eigenval_def
      by blast
    obtain v1 v2 where vv: "v = (v1, v2)"
      by (cases v) auto
    from H *mv v = k *sv v have "(H -mm (k *sm eye)) *mv v = vec_zero"
      using HH vv
      by (auto simp add: field_simps)
    hence "mat_det (H -mm (k *sm eye)) = 0"
      using v  vec_zero› vv HH
      using regular_homogenous_system[of "A - k" B C "D - k" v1 v2]
      unfolding det2_def
      by (auto simp add: field_simps)
    thus ?rhs
      using HH
      by (auto simp add: power2_eq_square field_simps)
  next
    assume ?rhs
    hence *: "mat_det (H -mm (k *sm eye)) = 0"
      using HH
      by (auto simp add: field_simps power2_eq_square)
    show ?lhs
    proof (cases "H -mm (k *sm eye) = mat_zero")
      case True
      thus ?thesis
        using HH
        by (auto) (rule_tac x=1 in exI, simp)
    next
      case False
      hence "(A - k  0  B  0)  (D - k  0  C  0)"
        using HH
        by auto
      thus ?thesis
      proof
        assume "A - k  0  B  0"
        hence "C * B + (D - k) * (k - A) = 0"
          using * singular_system[of "A-k" "D-k" B C "(0, 0)" 0 0  "(B, k-A)"] HH
          by (auto simp add: field_simps)
        hence  "(B, k-A)  vec_zero" "(H -mm (k *sm eye)) *mv (B, k-A) = vec_zero"
          using HH A - k  0  B  0
          by (auto simp add: field_simps)
        then obtain v where "v  vec_zero  (H -mm (k *sm eye)) *mv v = vec_zero"
          by blast
        thus ?thesis
          using HH
          unfolding eigenval_def
          by (rule_tac x="v" in exI) (case_tac v, simp add: field_simps)
      next
        assume "D - k  0  C  0"
        hence "C * B + (D - k) * (k - A) = 0"
          using * singular_system[of "D-k" "A-k" C B "(0, 0)" 0 0  "(C, k-D)"] HH
          by (auto simp add: field_simps)
        hence  "(k-D, C)  vec_zero" "(H -mm (k *sm eye)) *mv (k-D, C) = vec_zero"
          using HH D - k  0  C  0
          by (auto simp add: field_simps)
        then obtain v where "v  vec_zero  (H -mm (k *sm eye)) *mv v = vec_zero"
          by blast
        thus ?thesis
          using HH
          unfolding eigenval_def
          by (rule_tac x="v" in exI) (case_tac v, simp add: field_simps)
      qed
    qed
  qed
qed

(* ---------------------------------------------------------------------------- *)
subsubsection ‹Bilinear and Quadratic forms, Congruence, and Similarity›
(* ---------------------------------------------------------------------------- *)

text ‹Bilinear forms›

definition bilinear_form where
  [simp]: "bilinear_form v1 v2 H = (vec_cnj v1) *vm H *vv v2"

lemma bilinear_form_scale_m:
  shows "bilinear_form v1 v2 (k *sm H) = k * bilinear_form v1 v2 H"
  by (cases v1, cases v2, cases H) (simp add: vec_cnj_def field_simps)

lemma bilinear_form_scale_v1:
  shows "bilinear_form (k *sv v1) v2 H = cnj k * bilinear_form v1 v2 H"
  by (cases v1, cases v2, cases H) (simp add: vec_cnj_def field_simps)

lemma bilinear_form_scale_v2:
  shows "bilinear_form  v1 (k *sv v2) H = k * bilinear_form v1 v2 H"
  by (cases v1, cases v2, cases H) (simp add: vec_cnj_def field_simps)

text ‹Quadratic forms›

definition quad_form where
  [simp]: "quad_form v H = (vec_cnj v) *vm H *vv v"

lemma quad_form_bilinear_form: 
  shows "quad_form v H = bilinear_form v v H"
  by simp

lemma quad_form_scale_v:
  shows "quad_form (k *sv v) H = cor ((cmod k)2) * quad_form v H"
  using bilinear_form_scale_v1 bilinear_form_scale_v2
  by (simp add: complex_mult_cnj_cmod field_simps)

lemma quad_form_scale_m:
  shows "quad_form v (k *sm H) = k * quad_form v H"
  using bilinear_form_scale_m
  by simp

lemma cnj_quad_form [simp]:
  shows "cnj (quad_form z H) = quad_form z (mat_adj H)"
  by (cases H, cases z) (auto simp add: mat_adj_def mat_cnj_def vec_cnj_def field_simps)

text ‹Matrix congruence›

text ‹Two matrices are congruent iff they represent the same quadratic form with respect to different
bases (for example if one circline can be transformed to another by a Möbius trasformation).›

definition congruence where
  [simp]: "congruence M H  mat_adj M *mm H *mm M"

lemma congruence_nonzero:
  assumes "H  mat_zero" and "mat_det M  0"
  shows "congruence M H  mat_zero"
  using assms
  unfolding congruence_def
  by (subst mult_mm_non_zero_r, subst mult_mm_non_zero_l) (auto simp add: mat_det_adj)

lemma congruence_congruence:
  shows "congruence M1 (congruence M2 H) = congruence (M2 *mm M1) H"
  unfolding congruence_def
  apply (subst mult_mm_assoc)
  apply (subst mult_mm_assoc)
  apply (subst mat_adj_mult_mm)
  apply (subst mult_mm_assoc)
  by simp

lemma congruence_eye [simp]: 
  shows "congruence eye H = H"
  by (cases H) (simp add: mat_adj_def mat_cnj_def)

lemma congruence_congruence_inv [simp]:
  assumes "mat_det M  0"
  shows "congruence M (congruence (mat_inv M) H) = H"
  using assms congruence_congruence[of M "mat_inv M" H]
  using mat_inv_l[of M] mat_eye_l mat_eye_r
  unfolding congruence_def
  by (simp del: eye_def)

lemma congruence_inv:
  assumes "mat_det M  0" and "congruence M H = H'"
  shows "congruence (mat_inv M) H' = H"
  using assms
  using ‹mat_det M  0 mult_mm_inv_l[of "mat_adj M" "H *mm M" "H'"]
  using mult_mm_inv_r[of M "H" "mat_inv (mat_adj M) *mm H'"]
  by (simp add: mat_det_adj mult_mm_assoc mat_adj_inv)

lemma congruence_scale_m [simp]:
  shows "congruence M (k *sm H) = k *sm (congruence M H)"
  by (cases M, cases H) (auto simp add: mat_adj_def mat_cnj_def field_simps)

lemma inj_congruence:
  assumes "mat_det M  0" and "congruence M H = congruence M H'"
  shows "H = H'"
proof-
  have "H *mm M = H' *mm M "
    using assms
    using mult_mm_cancel_l[of "mat_adj M" "H *mm M" "H' *mm M"]
    by (simp add: mat_det_adj mult_mm_assoc)
  thus ?thesis
    using assms
    using mult_mm_cancel_r[of "M" "H" "H'"]
    by simp
qed

lemma mat_det_congruence [simp]:
  "mat_det (congruence M H) = (cor ((cmod (mat_det M))2)) * mat_det H"
  using complex_mult_cnj_cmod[of "mat_det M"]
  by (auto simp add: mat_det_adj field_simps)

lemma det_sgn_congruence [simp]:
  assumes "mat_det M  0"
  shows "sgn (mat_det (congruence M H)) = sgn (mat_det H)"
  using assms
  by (subst mat_det_congruence, auto simp add: sgn_mult power2_eq_square) (simp add: sgn_of_real)

lemma Re_det_sgn_congruence [simp]:
  assumes "mat_det M  0"
  shows "sgn (Re (mat_det (congruence M H))) = sgn (Re (mat_det H))"
proof-
  have *: "Re (mat_det (congruence M H)) = (cmod (mat_det M))2 * Re (mat_det H)"
    by (subst mat_det_congruence, subst Re_mult_real, rule Im_complex_of_real) (subst Re_complex_of_real, simp)
  show ?thesis
    using assms
    by (subst *) (auto simp add: sgn_mult)
qed

text ‹Transforming a matrix $H$ by a regular matrix $M$ preserves its bilinear and quadratic forms.›

lemma bilinear_form_congruence [simp]:
  assumes "mat_det M  0"
  shows "bilinear_form (M *mv v1) (M *mv v2) (congruence (mat_inv M) H) =
         bilinear_form v1 v2 H"
proof-
  have "mat_det (mat_adj M)  0"
    using assms
    by (simp add: mat_det_adj)
  show ?thesis
    unfolding bilinear_form_def congruence_def
    apply (subst mult_mv_mult_vm)
    apply (subst vec_cnj_mult_vm)
    apply (subst mat_adj_def[symmetric])
    apply (subst mult_vm_vm)
    apply (subst mult_vv_mv)
    apply (subst mult_vm_vm)
    apply (subst mat_adj_inv[OF ‹mat_det M  0])
    apply (subst mult_assoc_5)
    apply (subst mat_inv_r[OF ‹mat_det (mat_adj M)  0])
    apply (subst mat_inv_l[OF ‹mat_det M  0])
    apply (subst mat_eye_l, subst mat_eye_r)
    by simp
qed

lemma quad_form_congruence [simp]:
  assumes "mat_det M  0"
  shows "quad_form (M *mv z) (congruence (mat_inv M) H) = quad_form z H"
  using bilinear_form_congruence[OF assms]
  by simp


text ‹Similar matrices›

text ‹Two matrices are similar iff they represent the same linear operator with respect to (possibly)
different bases (e.g., if they represent the same Möbius transformation after changing the
coordinate system)›

definition similarity where
  "similarity A M = mat_inv A *mm M *mm A"

lemma mat_det_similarity [simp]:
  assumes "mat_det A  0"
  shows "mat_det (similarity A M) = mat_det M"
  using assms
  unfolding similarity_def
  by (simp add: mat_det_inv)

lemma mat_trace_similarity [simp]:
  assumes "mat_det A  0"
  shows "mat_trace (similarity A M) = mat_trace M"
proof-
  obtain a b c d where AA: "A = (a, b, c, d)"
    by (cases A) auto
  obtain mA mB mC mD where MM: "M = (mA, mB, mC, mD)"
    by (cases M) auto
  have "mA * (a * d) / (a * d - b * c) + mD * (a * d) / (a * d - b * c) =
        mA + mD + mA * (b * c) / (a * d - b * c) + mD * (b * c) / (a * d - b * c)"
    using assms AA
    by (simp add: field_simps)
  thus ?thesis
    using AA MM
    by (simp add: field_simps similarity_def)
qed

lemma similarity_eye [simp]:
  shows "similarity eye M = M"
  unfolding similarity_def
  using mat_eye_l mat_eye_r
  by auto


lemma similarity_eye' [simp]:
  shows "similarity (1, 0, 0, 1) M = M"
  unfolding eye_def[symmetric]
  by (simp del: eye_def)

lemma similarity_comp [simp]:
  assumes "mat_det A1  0" and "mat_det A2  0"
  shows "similarity A1 (similarity A2 M) = similarity (A2*mmA1) M"
  using assms
  unfolding similarity_def
  by (simp add: mult_mm_assoc mat_inv_mult_mm)

lemma similarity_inv:
  assumes "similarity A M1 = M2" and "mat_det A  0"
  shows "similarity (mat_inv A) M2 = M1"
  using assms
  unfolding similarity_def
  by (metis mat_det_mult mult_mm_assoc mult_mm_inv_l mult_mm_inv_r mult_zero_left)

end

Theory Unitary_Matrices

(* -------------------------------------------------------------------------- *)
subsection ‹Generalized Unitary Matrices›
(* -------------------------------------------------------------------------- *)

theory Unitary_Matrices
imports Matrices More_Complex
begin

text ‹In this section (generalized) $2\times 2$ unitary matrices are introduced.›

text ‹Unitary matrices›
definition unitary where
  "unitary M  mat_adj M *mm M = eye"

text ‹Generalized unitary matrices›
definition unitary_gen where
  "unitary_gen M 
   ( k::complex. k  0  mat_adj M *mm M