Theory Chordal_Metric

(* -------------------------------------------------------------------------- *)
subsection ‹Chordal Metric›
(* -------------------------------------------------------------------------- *)

text ‹Riemann sphere can be made a metric space. We are going to introduce distance on Riemann sphere
and to prove that it is a metric space. The distance between two points on the sphere is defined as
the length of the chord that connects them. This metric can be used in formalization of elliptic
geometry.›

theory Chordal_Metric
  imports Homogeneous_Coordinates Riemann_Sphere Oriented_Circlines "HOL-Analysis.Inner_Product" "HOL-Analysis.Euclidean_Space"
begin

(* -------------------------------------------------------------------------- *)
subsubsection ‹Inner product and norm›
(* -------------------------------------------------------------------------- *)

definition inprod_cvec :: "complex_vec ⇒ complex_vec ⇒ complex" where
 [simp]: "inprod_cvec z w =
             (let (z1, z2) = z;
                  (w1, w2) = w
               in vec_cnj (z1, z2) *⇩v⇩v (w1, w2))"
syntax
  "_inprod_cvec" :: "complex_vec ⇒ complex_vec ⇒ complex"  ("⟨_,_⟩")
translations
  "⟨z,w⟩" == "CONST inprod_cvec z w"

lemma real_inprod_cvec [simp]:
  shows "is_real ⟨z,z⟩"
  by (cases z, simp add: vec_cnj_def)

lemma inprod_cvec_ge_zero [simp]:
  shows "Re ⟨z,z⟩ ≥ 0"
  by (cases z, simp add: vec_cnj_def)

lemma inprod_cvec_bilinear1 [simp]:
  assumes "z' = k *⇩s⇩v  z"
  shows "⟨z',w⟩ = cnj k * ⟨z,w⟩"
  using assms
  by (cases z, cases z', cases w) (simp add: vec_cnj_def field_simps)

lemma inprod_cvec_bilinear2 [simp]:
  assumes "z' = k *⇩s⇩v z"
  shows "⟨w, z'⟩ = k * ⟨w, z⟩"
  using assms
  by (cases z, cases z', cases w) (simp add: vec_cnj_def field_simps)

lemma inprod_cvec_g_zero [simp]:
  assumes "z ≠ vec_zero"
  shows "Re ⟨z, z⟩ > 0"
proof-
  have "∀ a b. a ≠ 0 ∨ b ≠ 0 ⟶ 0 < (Re a * Re a + Im a * Im a) + (Re b * Re b + Im b * Im b)"
    by (smt complex_eq_0 not_sum_squares_lt_zero power2_eq_square)
  thus ?thesis
    using assms
    by (cases z, simp add: vec_cnj_def)
qed

definition norm_cvec :: "complex_vec ⇒ real" where
  [simp]: "norm_cvec z = sqrt (Re ⟨z,z⟩)"
syntax
  "_norm_cvec" :: "complex_vec ⇒ complex"  ("⟨_⟩")
translations
  "⟨z⟩" == "CONST norm_cvec z"

lemma norm_cvec_square:
  shows "⟨z⟩⇧2 = Re (⟨z,z⟩)"
  by (simp del: inprod_cvec_def)

lemma norm_cvec_gt_0:
  assumes "z ≠ vec_zero"
  shows "⟨z⟩ > 0"
  using assms
  by (simp del: inprod_cvec_def)

lemma norm_cvec_scale:
  assumes "z' = k *⇩s⇩v z"
  shows "⟨z'⟩⇧2 = Re (cnj k * k) * ⟨z⟩⇧2"
  unfolding norm_cvec_square
  using inprod_cvec_bilinear1[OF assms, of z']
  using inprod_cvec_bilinear2[OF assms, of z]
  by (simp del: inprod_cvec_def add: field_simps)

lift_definition inprod_hcoords :: "complex_homo_coords ⇒ complex_homo_coords ⇒ complex" is inprod_cvec
  done

lift_definition norm_hcoords :: "complex_homo_coords ⇒ real" is norm_cvec
  done

(* -------------------------------------------------------------------------- *)
subsubsection ‹Distance in $\mathbb{C}P^1$ - defined by Fubini-Study metric.›
(* -------------------------------------------------------------------------- *)

text ‹Formula for the chordal distance between the two points can be given directly based
on the homogenous coordinates of their stereographic projections in the plane. This is
called the Fubini-Study metric.›

definition dist_fs_cvec :: "complex_vec ⇒ complex_vec ⇒ real" where [simp]:
  "dist_fs_cvec z1 z2 =
     (let (z1x, z1y) = z1;
          (z2x, z2y) = z2;
          num = (z1x*z2y - z2x*z1y) * (cnj z1x*cnj z2y - cnj z2x*cnj z1y);
          den = (z1x*cnj z1x + z1y*cnj z1y) * (z2x*cnj z2x + z2y*cnj z2y)
       in 2*sqrt(Re num / Re den))"

lemma dist_fs_cvec_iff:
  assumes "z ≠ vec_zero" and "w ≠ vec_zero"
  shows "dist_fs_cvec z w = 2*sqrt(1 - (cmod ⟨z,w⟩)⇧2 / (⟨z⟩⇧2 * ⟨w⟩⇧2))"
proof-
  obtain z1 z2 w1 w2 where *: "z = (z1, z2)" "w = (w1, w2)"
    by (cases "z", cases "w") auto
  have 1: "2*sqrt(1 - (cmod ⟨z,w⟩)⇧2 / (⟨z⟩⇧2 * ⟨w⟩⇧2)) = 2*sqrt((⟨z⟩⇧2 * ⟨w⟩⇧2 - (cmod ⟨z,w⟩)⇧2) / (⟨z⟩⇧2 * ⟨w⟩⇧2))"
    using norm_cvec_gt_0[of z] norm_cvec_gt_0[of w] assms
    by (simp add: field_simps)

  have 2: "⟨z⟩⇧2 * ⟨w⟩⇧2 = Re ((z1*cnj z1 + z2*cnj z2) * (w1*cnj w1 + w2*cnj w2))"
    using assms *
    by (simp add: vec_cnj_def)

  have 3: "⟨z⟩⇧2 * ⟨w⟩⇧2 - (cmod ⟨z,w⟩)⇧2 = Re ((z1*w2 - w1*z2) * (cnj z1*cnj w2 - cnj w1*cnj z2))"
    apply (subst cmod_square, (subst norm_cvec_square)+)
    using *
    by (simp add: vec_cnj_def field_simps)

  thus ?thesis
    using 1 2 3
    using *
    unfolding dist_fs_cvec_def Let_def
    by simp
qed

lift_definition dist_fs_hcoords :: "complex_homo_coords ⇒ complex_homo_coords ⇒ real" is dist_fs_cvec
  done

lift_definition dist_fs :: "complex_homo ⇒ complex_homo ⇒ real" is dist_fs_hcoords
proof transfer
  fix z1 z2 z1' z2' :: complex_vec
  obtain z1x z1y z2x z2y z1'x z1'y z2'x z2'y where
    zz: "z1 = (z1x, z1y)" "z2 = (z2x, z2y)" "z1' = (z1'x, z1'y)" "z2' = (z2'x, z2'y)"
    by (cases "z1", cases "z2", cases "z1'", cases "z2'") blast

  assume 1: "z1 ≠ vec_zero" "z2 ≠ vec_zero" "z1' ≠ vec_zero" "z2' ≠ vec_zero" "z1 ≈⇩v z1'" "z2 ≈⇩v z2'"
  then obtain k1 k2 where
    *: "k1 ≠ 0" "z1' = k1 *⇩s⇩v z1" and
    **: "k2 ≠ 0" "z2' = k2 *⇩s⇩v z2"
    by auto
  have "(cmod ⟨z1,z2⟩)⇧2 / (⟨z1⟩⇧2 * ⟨z2⟩⇧2) = (cmod ⟨z1',z2'⟩)⇧2 / (⟨z1'⟩⇧2 * ⟨z2'⟩⇧2)"
    using ‹k1 ≠ 0› ‹k2 ≠ 0›
    using cmod_square[symmetric, of k1] cmod_square[symmetric, of k2]
    apply (subst norm_cvec_scale[OF *(2)])
    apply (subst norm_cvec_scale[OF **(2)])
    apply (subst inprod_cvec_bilinear1[OF *(2)])
    apply (subst inprod_cvec_bilinear2[OF **(2)])
    by (simp add: power2_eq_square norm_mult)
  thus "dist_fs_cvec z1 z2 = dist_fs_cvec z1' z2'"
    using 1 dist_fs_cvec_iff
    by simp
qed

lemma dist_fs_finite:
  shows "dist_fs (of_complex z1) (of_complex z2) = 2 * cmod(z1 - z2) / (sqrt (1+(cmod z1)⇧2) * sqrt (1+(cmod z2)⇧2))"
  apply transfer
  apply transfer
  apply (subst cmod_square)+
  apply (simp add: real_sqrt_divide cmod_def power2_eq_square)
  apply (subst real_sqrt_mult[symmetric])
  apply (simp add: field_simps)
  done

lemma dist_fs_infinite1:
  shows "dist_fs (of_complex z1) ∞⇩h = 2 / sqrt (1+(cmod z1)⇧2)"
  by (transfer, transfer) (subst cmod_square, simp add: real_sqrt_divide)

lemma dist_fs_infinite2:
  shows "dist_fs ∞⇩h (of_complex z1) = 2 / sqrt (1+(cmod z1)⇧2)"
  by (transfer, transfer) (subst cmod_square, simp add: real_sqrt_divide)

lemma dist_fs_cvec_zero:
  assumes "z ≠ vec_zero" and "w ≠ vec_zero"
  shows  "dist_fs_cvec z w = 0 ⟷ (cmod ⟨z,w⟩)⇧2 = (⟨z⟩⇧2 * ⟨w⟩⇧2)"
  using assms norm_cvec_gt_0[of z]  norm_cvec_gt_0[of w]
  by (subst dist_fs_cvec_iff) auto

lemma dist_fs_zero1 [simp]:
  shows "dist_fs z z = 0"
  by (transfer, transfer)
     (subst dist_fs_cvec_zero, simp, (subst norm_cvec_square)+, subst cmod_square, simp del: inprod_cvec_def)

lemma dist_fs_zero2 [simp]:
  assumes "dist_fs z1 z2 = 0"
  shows "z1 = z2"
  using assms
proof (transfer, transfer)
  fix z w :: complex_vec
  obtain z1 z2 w1 w2 where *: "z = (z1, z2)" "w = (w1, w2)"
    by (cases "z", cases "w", auto)
  let ?x = "(z1*w2 - w1*z2) * (cnj z1*cnj w2 - cnj w1*cnj z2)"
  assume "z ≠ vec_zero" "w ≠ vec_zero" "dist_fs_cvec z w = 0"
  hence "(cmod ⟨z,w⟩)⇧2 = ⟨z⟩⇧2 * ⟨w⟩⇧2"
    by (subst (asm) dist_fs_cvec_zero, simp_all)
  hence "Re ?x = 0"
    using *
    by (subst (asm) cmod_square) ((subst (asm) norm_cvec_square)+, simp add: vec_cnj_def field_simps)
  hence "?x = 0"
    using complex_mult_cnj_cmod[of "z1*w2 - w1*z2"] zero_complex.simps
    by (subst complex_eq_if_Re_eq[of ?x 0]) (simp add: power2_eq_square, simp, linarith)
  moreover
  have "z1 * w2 - w1 * z2 = 0 ⟷ cnj z1 * cnj w2 - cnj w1 * cnj z2 = 0"
    by (metis complex_cnj_diff complex_cnj_mult complex_cnj_zero_iff)
  ultimately
  show "z ≈⇩v w"
    using * ‹z ≠ vec_zero› ‹w ≠ vec_zero›
    using complex_cvec_eq_mix[of z1 z2 w1 w2]
    by auto
qed

lemma dist_fs_sym:
  shows "dist_fs z1 z2 = dist_fs z2 z1"
  by (transfer, transfer) (simp add: split_def field_simps)

(* -------------------------------------------------------------------------- *)
subsubsection ‹Triangle inequality for Fubini-Study metric›
(* -------------------------------------------------------------------------- *)

lemma dist_fs_triangle_finite:
  shows "cmod(a - b) / (sqrt (1+(cmod a)⇧2) * sqrt (1+(cmod b)⇧2)) ≤ cmod (a - c) / (sqrt (1+(cmod a)⇧2) * sqrt (1+(cmod c)⇧2)) + cmod (c - b) / (sqrt (1+(cmod b)⇧2) * sqrt (1+(cmod c)⇧2))"
proof-
  let ?cc = "1+(cmod c)⇧2" and ?bb = "1+(cmod b)⇧2" and ?aa = "1+(cmod a)⇧2"
  have "sqrt ?cc > 0" "sqrt ?aa > 0" "sqrt ?bb > 0"
    by (smt real_sqrt_gt_zero zero_compare_simps(12))+
  have "(a - b)*(1+cnj c*c) = (a-c)*(1+cnj c*b) + (c-b)*(1 + cnj c*a)"
    by (simp add: field_simps)
  moreover
  have "1 + cnj c * c = 1 + (cmod c)⇧2"
    using complex_norm_square
    by auto
  hence "cmod ((a - b)*(1+cnj c*c)) = cmod(a - b) * (1+(cmod c)⇧2)"
    by (smt norm_mult norm_of_real zero_compare_simps(12))
  ultimately
  have "cmod(a - b) * (1+(cmod c)⇧2) ≤ cmod (a-c) * cmod (1+cnj c*b) + cmod (c-b) * cmod(1 + cnj c*a)"
    using complex_mod_triangle_ineq2[of "(a-c)*(1+cnj c*b)" "(c-b)*(1 + cnj c*a)"]
    by (simp add: norm_mult)
  moreover
  have *: "⋀ a b c d b' d'. ⟦b ≤ b'; d ≤ d'; a ≥ (0::real); c ≥ 0⟧ ⟹ a*b + c*d ≤ a*b' + c*d'"
    by (simp add: add_mono_thms_linordered_semiring(1) mult_left_mono)
  have "cmod (a-c) * cmod (1+cnj c*b) + cmod (c-b) * cmod(1 + cnj c*a) ≤ cmod (a - c) * (sqrt (1+(cmod c)⇧2) * sqrt (1+(cmod b)⇧2)) + cmod (c - b) * (sqrt (1+(cmod c)⇧2) * sqrt (1+(cmod a)⇧2))"
    using *[OF cmod_1_plus_mult_le[of "cnj c" b] cmod_1_plus_mult_le[of "cnj c" a], of "cmod (a-c)" "cmod (c-b)"]
    by (simp add: field_simps real_sqrt_mult[symmetric])
  ultimately
  have "cmod(a - b) * ?cc ≤ cmod (a - c) * sqrt ?cc * sqrt ?bb + cmod (c - b) * sqrt ?cc * sqrt ?aa"
    by simp
  moreover
  hence "0 ≤ ?cc * sqrt ?aa * sqrt ?bb"
    using mult_right_mono[of 0 "sqrt ?aa"  "sqrt ?bb"]
    using mult_right_mono[of 0 "?cc" "sqrt ?aa * sqrt ?bb"]
    by simp
  moreover
  have "sqrt ?cc / ?cc = 1 / sqrt ?cc"
    using ‹sqrt ?cc > 0›
    by (simp add: field_simps)
  hence "sqrt ?cc / (?cc * sqrt ?aa) = 1 / (sqrt ?aa * sqrt ?cc)"
    using times_divide_eq_right[of "1/sqrt ?aa" "sqrt ?cc" "?cc"]
    using ‹sqrt ?aa > 0›
    by simp
  hence "cmod (a - c) * sqrt ?cc / (?cc * sqrt ?aa) = cmod (a - c) / (sqrt ?aa * sqrt ?cc)"
    using times_divide_eq_right[of "cmod (a - c)" "sqrt ?cc" "(?cc * sqrt ?aa)"]
    by simp
  moreover
  have "sqrt ?cc / ?cc = 1 / sqrt ?cc"
    using ‹sqrt ?cc > 0›
    by (simp add: field_simps)
  hence "sqrt ?cc / (?cc * sqrt ?bb) = 1 / (sqrt ?bb * sqrt ?cc)"
    using times_divide_eq_right[of "1/sqrt ?bb" "sqrt ?cc" "?cc"]
    using ‹sqrt ?bb > 0›
    by simp
  hence "cmod (c - b) * sqrt ?cc / (?cc * sqrt ?bb) = cmod (c - b) / (sqrt ?bb * sqrt ?cc)"
    using times_divide_eq_right[of "cmod (c - b)" "sqrt ?cc" "?cc * sqrt ?bb"]
    by simp
  ultimately
  show ?thesis
    using divide_right_mono[of "cmod (a - b) * ?cc" "cmod (a - c) * sqrt ?cc * sqrt ?bb + cmod (c - b) * sqrt ?cc * sqrt ?aa" "?cc * sqrt ?aa * sqrt ?bb"] ‹sqrt ?aa > 0› ‹sqrt ?bb > 0› ‹sqrt ?cc > 0›
    by (simp add: add_divide_distrib)
qed

lemma dist_fs_triangle_infinite1: 
  shows "1 / sqrt(1 + (cmod b)⇧2) ≤ 1 / sqrt(1 + (cmod c)⇧2) + cmod (b - c) / (sqrt(1 + (cmod b)⇧2) * sqrt(1 + (cmod c)⇧2))"
proof-
  let ?bb = "sqrt (1 + (cmod b)⇧2)" and ?cc = "sqrt (1 + (cmod c)⇧2)"
  have "?bb > 0" "?cc > 0"
    by (metis add_strict_increasing real_sqrt_gt_0_iff zero_le_power2 zero_less_one)+
  hence *: "?bb * ?cc ≥ 0"
    by simp
  have **: "(?cc - ?bb) / (?bb * ?cc) = 1 / ?bb - 1 / ?cc"
    using ‹sqrt (1 + (cmod b)⇧2) > 0›  ‹sqrt (1 + (cmod c)⇧2) > 0›
    by (simp add: field_simps)
  show "1 / ?bb ≤ 1 / ?cc + cmod (b - c) / (?bb * ?cc)"
    using divide_right_mono[OF cmod_diff_ge[of c b] *]
    by (subst (asm) **) (simp add: field_simps norm_minus_commute)
qed

lemma dist_fs_triangle_infinite2:
  shows "1 / sqrt(1 + (cmod a)⇧2) ≤ cmod (a - c) / (sqrt (1+(cmod a)⇧2) * sqrt (1+(cmod c)⇧2)) + 1 / sqrt(1 + (cmod c)⇧2)"
  using dist_fs_triangle_infinite1[of a c]
  by simp

lemma dist_fs_triangle_infinite3:
  shows "cmod(a - b) / (sqrt (1+(cmod a)⇧2) * sqrt (1+(cmod b)⇧2)) ≤ 1 / sqrt(1 + (cmod a)⇧2) + 1 / sqrt(1 + (cmod b)⇧2)"
proof-
  let ?aa = "sqrt (1 + (cmod a)⇧2)" and ?bb = "sqrt (1 + (cmod b)⇧2)"
  have "?aa > 0" "?bb > 0"
    by (metis add_strict_increasing real_sqrt_gt_0_iff zero_le_power2 zero_less_one)+
  hence *: "?aa * ?bb ≥ 0"
    by simp
  have **: "(?aa + ?bb) / (?aa * ?bb) = 1 / ?aa + 1 / ?bb"
    using ‹?aa > 0› ‹?bb > 0›
    by (simp add: field_simps)
  show "cmod (a - b) / (?aa * ?bb) ≤ 1 / ?aa + 1 / ?bb"
    using divide_right_mono[OF cmod_diff_le[of a b] *]
    by (subst (asm) **) (simp add: field_simps norm_minus_commute)
qed

lemma dist_fs_triangle:
  shows "dist_fs A B ≤ dist_fs A C + dist_fs C B"
proof (cases "A = ∞⇩h")
  case True
  show ?thesis
  proof (cases "B = ∞⇩h")
    case True
    show ?thesis
    proof (cases "C = ∞⇩h")
      case True
      show ?thesis
        using ‹A = ∞⇩h› ‹B = ∞⇩h› ‹C = ∞⇩h›
        by simp
    next
      case False
      then obtain c where "C = of_complex c"
        using inf_or_of_complex[of C]
        by auto
      show ?thesis
        using ‹A = ∞⇩h› ‹B = ∞⇩h› ‹C = of_complex c›
        by (simp add: dist_fs_infinite2 dist_fs_sym)
    qed
  next
    case False
    then obtain b where "B = of_complex b"
      using inf_or_of_complex[of B]
      by auto
    show ?thesis
    proof (cases "C = ∞⇩h")
      case True
      show ?thesis
        using ‹A = ∞⇩h› ‹C = ∞⇩h› ‹B = of_complex b›
        by simp
    next
      case False
      then obtain c where "C = of_complex c"
        using inf_or_of_complex[of C]
        by auto
      show ?thesis
        using ‹A = ∞⇩h› ‹B = of_complex b› ‹C = of_complex c›
        using mult_left_mono[OF dist_fs_triangle_infinite1[of b c], of 2]
        by (simp add: dist_fs_finite dist_fs_infinite1 dist_fs_infinite2 dist_fs_sym)
    qed
  qed
next
  case False
  then obtain a where "A = of_complex a"
    using inf_or_of_complex[of A]
    by auto
  show ?thesis
  proof (cases "B = ∞⇩h")
    case True
    show ?thesis
    proof (cases "C = ∞⇩h")
      case True
      show ?thesis
        using ‹B = ∞⇩h› ‹C = ∞⇩h› ‹A = of_complex a›
        by (simp add: dist_fs_infinite2)
    next
      case False
      then obtain c where "C = of_complex c"
        using inf_or_of_complex[of C]
        by auto
      show ?thesis
        using ‹B = ∞⇩h› ‹C = of_complex c› ‹A = of_complex a›
        using mult_left_mono[OF dist_fs_triangle_infinite2[of a c], of 2]
        by (simp add: dist_fs_finite dist_fs_infinite1 dist_fs_infinite2)
    qed
  next
    case False
    then obtain b where "B = of_complex b"
      using inf_or_of_complex[of B]
      by auto
    show ?thesis
    proof (cases "C = ∞⇩h")
      case True
      thus ?thesis
        using ‹C = ∞⇩h› ‹B = of_complex b› ‹A = of_complex a›
        using mult_left_mono[OF dist_fs_triangle_infinite3[of a b], of 2]
        by (simp add: dist_fs_finite dist_fs_infinite1 dist_fs_infinite2)
    next
      case False
      then obtain c where "C = of_complex c"
        using inf_or_of_complex[of C]
        by auto
      show ?thesis
        using ‹A = of_complex a› ‹B = of_complex b› ‹C = of_complex c›
        using mult_left_mono[OF dist_fs_triangle_finite[of a b c], of 2]
        by (simp add: dist_fs_finite norm_minus_commute dist_fs_sym)
    qed
  qed
qed

(* -------------------------------------------------------------------------- *)
subsubsection ‹$\mathbb{C}P^1$ with Fubini-Study metric is a metric space›
(* -------------------------------------------------------------------------- *)

text ‹Using the (already available) fact that $\mathbb{R}^3$ is a metric space (under the distance
function $\lambda\ x\ y.\ @{term norm}(x - y)$), it was not difficult to show that the type @{term
complex_homo} equipped with @{term dist_fs} is a metric space, i.e., an instantiation of the @{term
metric_space} locale.›

instantiation complex_homo :: metric_space
begin
definition "dist_complex_homo = dist_fs"
definition "(uniformity_complex_homo :: (complex_homo × complex_homo) filter) = (INF e∈{0<..}. principal {(x, y). dist_class.dist x y < e})"
definition "open_complex_homo (U :: complex_homo set) = (∀ x ∈ U. eventually (λ(x', y). x' = x ⟶ y ∈ U) uniformity)"
instance
proof
  fix x y :: complex_homo
  show "(dist_class.dist x y = 0) = (x = y)"
    unfolding dist_complex_homo_def
    using dist_fs_zero1[of x] dist_fs_zero2[of x y]
    by auto
next
  fix x y z :: complex_homo
  show "dist_class.dist x y ≤ dist_class.dist x z + dist_class.dist y z"
    unfolding dist_complex_homo_def
    using dist_fs_triangle[of x y z]
    by (simp add: dist_fs_sym)
qed (simp_all add: open_complex_homo_def uniformity_complex_homo_def)
end

(* -------------------------------------------------------------------------- *)
subsubsection ‹Chordal distance on the Riemann sphere›
(* -------------------------------------------------------------------------- *)

text ‹Distance of the two points is given by the length of the chord. We show that it corresponds to
the Fubini-Study metric in the plane.›

definition dist_riemann_sphere_r3 :: "R3 ⇒ R3 ⇒ real" where [simp]:
  "dist_riemann_sphere_r3 M1 M2 =
     (let (x1, y1, z1) = M1;
          (x2, y2, z2) = M2
       in norm (x1 - x2, y1 - y2, z1 - z2))"

lemma dist_riemann_sphere_r3_inner:
  assumes "M1 ∈ unit_sphere" and "M2 ∈ unit_sphere"
  shows  "(dist_riemann_sphere_r3 M1 M2)⇧2 = 2 - 2 * inner M1 M2"
  using assms
  apply (cases M1, cases M2)
  apply (auto simp add: norm_prod_def)
  apply (simp add: power2_eq_square field_simps)
  done


lift_definition dist_riemann_sphere' :: "riemann_sphere ⇒ riemann_sphere ⇒ real" is dist_riemann_sphere_r3
  done

lemma dist_riemann_sphere_ge_0 [simp]: 
  shows "dist_riemann_sphere' M1 M2 ≥ 0"
  apply transfer
  using norm_ge_zero
  by (simp add: split_def Let_def)

text ‹Using stereographic projection we prove the connection between chordal metric on the spehere
and Fubini-Study metric in the plane.›

lemma dist_stereographic_finite:
  assumes "stereographic M1 = of_complex m1"  and "stereographic M2 = of_complex m2"
  shows "dist_riemann_sphere' M1 M2 = 2 * cmod (m1 - m2) / (sqrt (1 + (cmod m1)⇧2) * sqrt (1 + (cmod m2)⇧2))"
  using assms
proof-
  have *: "M1 = inv_stereographic (of_complex m1)"  "M2 = inv_stereographic (of_complex m2)"
    using inv_stereographic_is_inv assms
    by (metis inv_stereographic_stereographic)+
  have "(1 + (cmod m1)⇧2) ≠ 0"  "(1 + (cmod m2)⇧2) ≠ 0"
    by (smt power2_less_0)+
  have "(1 + (cmod m1)⇧2) > 0"  "(1 + (cmod m2)⇧2) > 0"
    by (smt realpow_square_minus_le)+
  hence "(1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2) > 0"
    by (metis norm_mult_less norm_zero power2_eq_square zero_power2)
  hence ++: "sqrt ((1 + cmod m1 * cmod m1) * (1 + cmod m2 * cmod m2)) > 0"
    using real_sqrt_gt_0_iff
    by (simp add: power2_eq_square)
  hence **: "(2 * cmod (m1 - m2) / sqrt ((1 + cmod m1 * cmod m1) * (1 + cmod m2 * cmod m2))) ≥ 0 ⟷ cmod (m1 - m2) ≥ 0"
    by (metis diff_self divide_nonneg_pos mult_2 norm_ge_zero norm_triangle_ineq4 norm_zero)

  have "(dist_riemann_sphere' M1 M2)⇧2 * (1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2) = 4 * (cmod (m1 - m2))⇧2"
    using *
  proof (transfer, transfer)
    fix m1 m2 M1 M2
    assume us: "M1 ∈ unit_sphere" "M2 ∈ unit_sphere" and
        *: "M1 = inv_stereographic_cvec_r3 (of_complex_cvec m1)" "M2 = inv_stereographic_cvec_r3 (of_complex_cvec m2)"
    have "(1 + (cmod m1)⇧2) ≠ 0"  "(1 + (cmod m2)⇧2) ≠ 0"
      by (smt power2_less_0)+
    thus "(dist_riemann_sphere_r3 M1 M2)⇧2 * (1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2) =
          4 * (cmod (m1 - m2))⇧2"
      apply (subst dist_riemann_sphere_r3_inner[OF us])
      apply (subst *)+
      apply (simp add: dist_riemann_sphere_r3_inner[OF us] complex_mult_cnj_cmod)
      apply (subst left_diff_distrib[of 2])
      apply (subst left_diff_distrib[of "2*(1+(cmod m1)⇧2)"])
      apply (subst distrib_right[of _ _ "(1 + (cmod m1)⇧2)"])
      apply (subst distrib_right[of _ _ "(1 + (cmod m1)⇧2)"])
      apply simp
      apply (subst distrib_right[of _ _ "(1 + (cmod m2)⇧2)"])
      apply (subst distrib_right[of _ _ "(1 + (cmod m2)⇧2)"])
      apply (subst distrib_right[of _ _ "(1 + (cmod m2)⇧2)"])
      apply simp
      apply (subst (asm) cmod_square)+
      apply (subst cmod_square)+
      apply (simp add: field_simps)
      done
  qed
  hence "(dist_riemann_sphere' M1 M2)⇧2 = 4 * (cmod (m1 - m2))⇧2 / ((1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2))"
    using ‹(1 + (cmod m1)⇧2) ≠ 0›  ‹(1 + (cmod m2)⇧2) ≠ 0›
    using eq_divide_imp[of "(1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2)" "(dist_riemann_sphere' M1 M2)⇧2" "4 * (cmod (m1 - m2))⇧2"]
    by simp
  thus "dist_riemann_sphere' M1 M2 = 2 * cmod (m1 - m2) / (sqrt (1 + (cmod m1)⇧2) * sqrt (1 + (cmod m2)⇧2))"
    using power2_eq_iff[of "dist_riemann_sphere' M1 M2" "2 * (cmod (m1 - m2)) / sqrt ((1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2))"]
    using ‹(1 + (cmod m1)⇧2) * (1 + (cmod m2)⇧2) > 0›  ‹(1 + (cmod m1)⇧2) > 0› ‹(1 + (cmod m2)⇧2) > 0›
    apply (auto simp add: power2_eq_square real_sqrt_mult[symmetric])
    using dist_riemann_sphere_ge_0[of M1 M2] **
    using ++ divide_le_0_iff by force
qed


lemma dist_stereographic_infinite:
  assumes "stereographic M1 = ∞⇩h" and "stereographic M2 = of_complex m2"
  shows "dist_riemann_sphere' M1 M2 = 2 / sqrt (1 + (cmod m2)⇧2)"
proof-
  have *: "M1 = inv_stereographic ∞⇩h"  "M2 = inv_stereographic (of_complex m2)"
    using inv_stereographic_is_inv assms
    by (metis inv_stereographic_stereographic)+
  have "(1 + (cmod m2)⇧2) ≠ 0"
    by (smt power2_less_0)
  have "(1 + (cmod m2)⇧2) > 0"
    by (smt realpow_square_minus_le)+
  hence "sqrt (1 + cmod m2 * cmod m2) > 0"
    using real_sqrt_gt_0_iff
    by (simp add: power2_eq_square)
  hence **: "2 / sqrt (1 + cmod m2 * cmod m2) > 0"
    by simp

  have "(dist_riemann_sphere' M1 M2)⇧2 * (1 + (cmod m2)⇧2) = 4"
    using *
    apply transfer
    apply transfer
  proof-
    fix M1 M2 m2
    assume us: "M1 ∈ unit_sphere" "M2 ∈ unit_sphere" and
           *: "M1 = inv_stereographic_cvec_r3 ∞⇩v" "M2 = inv_stereographic_cvec_r3 (of_complex_cvec m2)"
    have "(1 + (cmod m2)⇧2) ≠ 0"
      by (smt power2_less_0)
    thus "(dist_riemann_sphere_r3 M1 M2)⇧2 * (1 + (cmod m2)⇧2) = 4"
      apply (subst dist_riemann_sphere_r3_inner[OF us])
      apply (subst *)+
      apply (simp add: complex_mult_cnj_cmod)
      apply (subst left_diff_distrib[of 2], simp)
      done
  qed
  hence "(dist_riemann_sphere' M1 M2)⇧2 = 4 / (1 + (cmod m2)⇧2)"
    using ‹(1 + (cmod m2)⇧2) ≠ 0›
    by (simp add: field_simps)
  thus "dist_riemann_sphere' M1 M2 = 2 / sqrt (1 + (cmod m2)⇧2)"
    using power2_eq_iff[of "dist_riemann_sphere' M1 M2" "2 / sqrt (1 + (cmod m2)⇧2)"]
    using ‹(1 + (cmod m2)⇧2) > 0›
    apply (auto simp add: power2_eq_square real_sqrt_mult[symmetric])
    using dist_riemann_sphere_ge_0[of M1 M2] **
    by simp
qed

lemma dist_rieman_sphere_zero [simp]: 
  shows "dist_riemann_sphere' M M = 0"
  by transfer auto

lemma dist_riemann_sphere_sym: 
  shows "dist_riemann_sphere' M1 M2 = dist_riemann_sphere' M2 M1"
proof transfer
  fix M1 M2 :: R3
  obtain x1 y1 z1 x2 y2 z2 where MM: "(x1, y1, z1) = M1" "(x2, y2, z2) = M2"
    by (cases "M1", cases "M2", auto)
  show "dist_riemann_sphere_r3 M1 M2 = dist_riemann_sphere_r3 M2 M1"
    using norm_minus_cancel[of "(x1 - x2, y1 - y2, z1 - z2)"] MM[symmetric]
    by simp
qed

text ‹Central theorem that connects the two metrics.›
lemma dist_stereographic:
  shows "dist_riemann_sphere' M1 M2 = dist_fs (stereographic M1) (stereographic M2)"
proof (cases "M1 = North")
  case True
  hence "stereographic M1 = ∞⇩h"
    by (simp add: stereographic_North)
  show ?thesis
  proof (cases "M2 = North")
    case True
    show ?thesis
      using ‹M1 = North› ‹M2 = North›
      by auto
  next
    case False
    hence "stereographic M2 ≠ ∞⇩h"
      using stereographic_North[of M2]
      by simp
    then obtain m2 where "stereographic M2 = of_complex m2"
      using inf_or_of_complex[of "stereographic M2"]
      by auto
    show ?thesis
      using ‹stereographic M2 = of_complex m2› ‹stereographic M1 = ∞⇩h›
      using dist_fs_infinite1 dist_stereographic_infinite
      by (simp add: dist_fs_sym)
  qed
next
  case False
  hence "stereographic M1 ≠ ∞⇩h"
    by (simp add: stereographic_North)
  then obtain m1 where "stereographic M1 = of_complex m1"
    using inf_or_of_complex[of "stereographic M1"]
    by auto
  show ?thesis
  proof (cases "M2 = North")
    case True
    hence "stereographic M2 = ∞⇩h"
      by (simp add: stereographic_North)
    show ?thesis
      using ‹stereographic M1 = of_complex m1› ‹stereographic M2 = ∞⇩h›
      using dist_fs_infinite2 dist_stereographic_infinite
      by (subst dist_riemann_sphere_sym, simp add: dist_fs_sym)
  next
    case False
    hence "stereographic M2 ≠ ∞⇩h"
      by (simp add: stereographic_North)
    then obtain m2 where "stereographic M2 = of_complex m2"
      using inf_or_of_complex[of "stereographic M2"]
      by auto
    show ?thesis
      using ‹stereographic M1 = of_complex m1› ‹stereographic M2 = of_complex m2›
      using dist_fs_finite dist_stereographic_finite
      by simp
  qed
qed

text ‹Other direction›
lemma dist_stereographic':
  shows "dist_fs A B = dist_riemann_sphere' (inv_stereographic A) (inv_stereographic B)"
  by (subst dist_stereographic) (metis stereographic_inv_stereographic)

text ‹The @{term riemann_sphere} equipped with @{term dist_riemann_sphere'} is a metric space, i.e.,
an instantiation of the @{term metric_space} locale.›

instantiation riemann_sphere :: metric_space
begin
definition "dist_riemann_sphere = dist_riemann_sphere'"
definition "(uniformity_riemann_sphere :: (riemann_sphere × riemann_sphere) filter) = (INF e∈{0<..}. principal {(x, y). dist_class.dist x y < e})"
definition "open_riemann_sphere (U :: riemann_sphere set) = (∀ x ∈ U. eventually (λ(x', y). x' = x ⟶ y ∈ U) uniformity)"
instance
proof
  fix x y :: riemann_sphere
  show "(dist_class.dist x y = 0) = (x = y)"
    unfolding dist_riemann_sphere_def
  proof transfer
    fix x y :: R3
    obtain x1 y1 z1 x2 y2 z2 where *: "(x1, y1, z1) = x" "(x2, y2, z2) = y"
      by (cases x, cases y, auto)
    assume "x ∈ unit_sphere" "y ∈ unit_sphere"
    thus "(dist_riemann_sphere_r3 x y = 0) = (x = y)"
      using norm_eq_zero[of "(x1 - y2, y1 - y2, z1 - z2)"] using *[symmetric]
      by (simp add: zero_prod_def)
  qed
next
  fix x y z :: riemann_sphere
  show "dist_class.dist x y ≤ dist_class.dist x z + dist_class.dist y z"
    unfolding dist_riemann_sphere_def
  proof transfer
    fix x y z :: R3
    obtain x1 y1 z1 x2 y2 z2 x3 y3 z3 where MM: "(x1, y1, z1) = x" "(x2, y2, z2) = y" "(x3, y3, z3) = z"
      by (cases "x", cases "y", cases "z", auto)

    assume "x ∈ unit_sphere" "y ∈ unit_sphere" "z ∈ unit_sphere"
    thus "dist_riemann_sphere_r3 x y ≤ dist_riemann_sphere_r3 x z + dist_riemann_sphere_r3 y z"
      using MM[symmetric] norm_minus_cancel[of "(x3 - x2, y3 - y2, z3 - z2)"] norm_triangle_ineq[of "(x1 - x3, y1 - y3, z1 - z3)" "(x3 - x2, y3 - y2, z3 - z2)"]
      by simp
  qed
qed (simp_all add: uniformity_riemann_sphere_def open_riemann_sphere_def)
end

text ‹The @{term riemann_sphere} metric space is perfect, i.e., it does not have isolated points.›
instantiation riemann_sphere :: perfect_space
begin
instance proof
  fix M :: riemann_sphere
  show "¬ open {M}"
    unfolding open_dist dist_riemann_sphere_def
    apply (subst dist_riemann_sphere_sym)
  proof transfer
    fix M
    assume "M ∈ unit_sphere"
    obtain x y z where MM: "M = (x, y, z)"
      by (cases "M") auto
    then obtain α β where *: "x = cos α * cos β" "y = cos α * sin β" "z = sin α" "-pi / 2 ≤ α ∧ α ≤ pi / 2"
      using ‹M ∈ unit_sphere›
      using ex_sphere_params[of x y z]
      by auto
    have "⋀ e. e > 0 ⟹ (∃y. y ∈ unit_sphere ∧ dist_riemann_sphere_r3 M y < e ∧ y ≠ M)"
    proof-
      fix e :: real
      assume "e > 0"
      then obtain α' where "1 - (e*e/2) < cos (α - α')" "α ≠ α'" "-pi/2 ≤ α'" "α' ≤ pi/2"
        using ex_cos_gt[of α "1 - (e*e/2)"] ‹- pi / 2 ≤ α ∧ α ≤ pi / 2›
        by auto
      hence "sin α ≠ sin α'"
        using ‹-pi / 2 ≤ α ∧ α ≤ pi / 2› sin_inj[of α α']
        by auto

      have "2 - 2 * cos (α - α') < e*e"
        using mult_strict_right_mono[OF ‹1 - (e*e/2) < cos (α - α')›, of 2]
        by (simp add: field_simps)
      have "2 - 2 * cos (α - α') ≥ 0"
        using cos_le_one[of "α - α'"]
        by (simp add: algebra_split_simps)
      let ?M' = "(cos α' * cos β,  cos α' * sin β, sin α')"
      have "dist_riemann_sphere_r3 M ?M' = sqrt ((cos α - cos α')⇧2 + (sin α - sin α')⇧2)"
        using MM * sphere_params_on_sphere[of _ α' β]
        using sin_cos_squared_add[of β]
        apply (simp add: dist_riemann_sphere'_def Abs_riemann_sphere_inverse norm_prod_def)
        apply (subst left_diff_distrib[symmetric])+
        apply (subst power_mult_distrib)+
        apply (subst distrib_left[symmetric])
        apply simp
        done
      also have "... = sqrt (2 - 2*cos (α - α'))"
        by (simp add: power2_eq_square field_simps cos_diff)
      finally
      have "(dist_riemann_sphere_r3 M ?M')⇧2 = 2 - 2*cos (α - α')"
        using ‹2 - 2 * cos (α - α') ≥ 0›
        by simp
      hence "(dist_riemann_sphere_r3 M ?M')⇧2 < e⇧2"
        using ‹2 - 2 * cos (α - α') < e*e›
        by (simp add: power2_eq_square)
      hence "dist_riemann_sphere_r3 M ?M' < e"
        apply (rule power2_less_imp_less)
        using ‹e > 0›
        by simp
      moreover
      have "M ≠ ?M'"
        using MM  ‹sin α ≠ sin α'› *
        by simp
      moreover
      have "?M' ∈ unit_sphere"
        using sphere_params_on_sphere by auto
      ultimately
      show "∃y. y ∈ unit_sphere ∧ dist_riemann_sphere_r3 M y < e ∧ y ≠ M"
        unfolding dist_riemann_sphere_def
        by (rule_tac x="?M'" in exI, simp)
    qed
    thus "¬ (∀x∈{M}. ∃e>0. ∀y∈{x. x ∈ unit_sphere}. dist_riemann_sphere_r3 x y < e ⟶ y ∈ {M})"
      by auto
  qed
qed
end

text ‹The @{term complex_homo} metric space is perfect, i.e., it does not have isolated points.›
instantiation complex_homo :: perfect_space
begin
instance proof
  fix x::complex_homo
  show "¬ open {x}"
    unfolding open_dist
  proof (auto)
    fix e::real
    assume "e > 0"
    thus "∃ y. dist_class.dist y x < e ∧ y ≠ x"
      using not_open_singleton[of "inv_stereographic x"]
      unfolding open_dist
      unfolding dist_complex_homo_def dist_riemann_sphere_def
      apply (subst dist_stereographic', auto)
      apply (erule_tac x=e in allE, auto)
      apply (rule_tac x="stereographic y" in exI, auto)
      done
  qed
qed

end

lemma Lim_within:
  shows "(f ⤏ l) (at a within S) ⟷
    (∀e >0. ∃d>0. ∀x ∈ S. 0 < dist_class.dist x a ∧ dist_class.dist x a  < d ⟶ dist_class.dist (f x) l < e)"
  by (auto simp: tendsto_iff eventually_at)

lemma continuous_on_iff:
  shows "continuous_on s f ⟷
    (∀x∈s. ∀e>0. ∃d>0. ∀x'∈s. dist_class.dist x' x < d ⟶ dist_class.dist (f x') (f x) < e)"
  unfolding continuous_on_def Lim_within
  by (metis dist_pos_lt dist_self)

text ‹Using the chordal metric in the extended plane, and the Euclidean metric on the sphere in
$\mathbb{R}^3$, the stereographic and inverse stereographic projections are proved to be
continuous.›

lemma "continuous_on UNIV stereographic"
unfolding continuous_on_iff
unfolding dist_complex_homo_def dist_riemann_sphere_def
by (subst dist_stereographic', auto)

lemma "continuous_on UNIV inv_stereographic"
unfolding continuous_on_iff
unfolding dist_complex_homo_def dist_riemann_sphere_def
by (subst dist_stereographic, auto)

(* -------------------------------------------------------------------------- *)
subsubsection ‹Chordal circles›
(* -------------------------------------------------------------------------- *)

text ‹Real circlines are sets of points that are equidistant from some given point in the chordal
metric. There are exactly two such points (two chordal centers). On the Riemann sphere, these two
points are obtained as intersections of the sphere and a line that goes trough center of the circle,
and its orthogonal to its plane.›

text ‹The circline for the given chordal center and radius.›
definition chordal_circle_cvec_cmat :: "complex_vec ⇒ real ⇒ complex_mat" where
 [simp]: "chordal_circle_cvec_cmat a r =
           (let (a1, a2) = a
             in ((4*a2*cnj a2 - (cor r)⇧2*(a1*cnj a1 + a2*cnj a2)), (-4*a1*cnj a2), (-4*cnj a1*a2), (4*a1*cnj a1 - (cor r)⇧2*(a1*cnj a1 + a2*cnj a2))))"

lemma chordal_circle_cmat_hermitean_nonzero [simp]:
  assumes "a ≠ vec_zero"
  shows "chordal_circle_cvec_cmat a r ∈ hermitean_nonzero"
  using assms
  by (cases a) (auto simp add: hermitean_def mat_adj_def mat_cnj_def Let_def)

lift_definition chordal_circle_hcoords_clmat :: "complex_homo_coords ⇒ real ⇒ circline_mat" is chordal_circle_cvec_cmat
  using chordal_circle_cmat_hermitean_nonzero
  by (simp del: chordal_circle_cvec_cmat_def)

lift_definition chordal_circle :: "complex_homo ⇒ real ⇒ circline" is chordal_circle_hcoords_clmat
proof transfer
  fix a b :: complex_vec and r :: real
  assume *: "a ≠ vec_zero" "b ≠ vec_zero"
  obtain a1 a2 where aa: "a = (a1, a2)"
    by (cases a, auto)
  obtain b1 b2 where bb: "b = (b1, b2)"
    by (cases b, auto)
  assume "a ≈⇩v b"
  then obtain k where "b = (k * a1, k * a2)" "k ≠ 0"
    using aa bb
    by auto
  moreover
  have "cor (Re (k * cnj k)) = k * cnj k"
    by (metis complex_In_mult_cnj_zero complex_of_real_Re)
  ultimately
  show "circline_eq_cmat (chordal_circle_cvec_cmat a r) (chordal_circle_cvec_cmat b r)"
    using * aa bb
    by simp (rule_tac x="Re (k*cnj k)" in exI, auto simp add: Let_def field_simps)
qed

lemma sqrt_1_plus_square:
  shows "sqrt (1 + a⇧2) ≠ 0"
  by (smt real_sqrt_less_mono real_sqrt_zero realpow_square_minus_le)

lemma
  assumes "dist_fs z a = r"
  shows "z ∈ circline_set (chordal_circle a r)"
proof (cases "a ≠ ∞⇩h")
  case True
  then obtain a' where  "a = of_complex a'"
    using inf_or_of_complex
    by auto
  let ?A = "4 - (cor r)⇧2 * (1 + (a'*cnj a'))" and ?B = "-4*a'" and ?C="-4*cnj a'" and ?D = "4*a'*cnj a' -  (cor r)⇧2 * (1 + (a'*cnj a'))"
  have hh: "(?A, ?B, ?C, ?D) ∈ {H. hermitean H ∧ H ≠ mat_zero}"
    by (auto simp add: hermitean_def mat_adj_def mat_cnj_def power2_eq_square)
  hence *: "chordal_circle (of_complex a') r = mk_circline ?A ?B ?C ?D"
    by (transfer, transfer, simp, rule_tac x=1 in exI, simp)

  show ?thesis
  proof (cases "z ≠ ∞⇩h")
    case True
    then obtain z' where "z = of_complex z'"
      using inf_or_of_complex[of z] inf_or_of_complex[of a]
      by auto
    have "2 * cmod (z' - a') / (sqrt (1 + (cmod z')⇧2) * sqrt (1 + (cmod a')⇧2)) = r"
      using dist_fs_finite[of z' a'] assms ‹z = of_complex z'› ‹a = of_complex a'›
      by auto
    hence "4 * (cmod (z' - a'))⇧2 / ((1 + (cmod z')⇧2) * (1 + (cmod a')⇧2)) = r⇧2 "
      by (auto simp add: field_simps)
    moreover
    have "sqrt (1 + (cmod z')⇧2) ≠ 0" "sqrt (1 + (cmod a')⇧2) ≠ 0"
      using sqrt_1_plus_square
      by simp+
    hence "(1 + (cmod z')⇧2) * (1 + (cmod a')⇧2) ≠ 0"
      by simp
    ultimately
    have "4 * (cmod (z' - a'))⇧2  = r⇧2 * ((1 + (cmod z')⇧2) * (1 + (cmod a')⇧2))"
      by (simp add: field_simps)
    hence "4 * Re ((z' - a')*cnj (z' - a')) = r⇧2 * (1 + Re (z'*cnj z')) * (1 + Re (a'*cnj a'))"
      by ((subst cmod_square[symmetric])+, simp)
    hence "4 * (Re(z'*cnj z') - Re(a'*cnj z') - Re(cnj a'*z') + Re(a'*cnj a')) = r⇧2 * (1 + Re (z'*cnj z')) * (1 + Re (a'*cnj a'))"
      by (simp add: field_simps)
    hence "Re (?A * z' * cnj z' + ?B * cnj z' + ?C * z' + ?D) = 0"
      by (simp add: power2_eq_square field_simps)
    hence "?A * z' * cnj z' + ?B * cnj z' + ?C * z' + ?D = 0"
      by (subst complex_eq_if_Re_eq) (auto simp add: power2_eq_square)
    hence "(cnj z' * ?A + ?C) * z' + (cnj z' * ?B + ?D) = 0"
      by algebra
    hence "z ∈ circline_set (mk_circline ?A ?B ?C ?D)"
      using ‹z = of_complex z'› hh
      unfolding circline_set_def
      by simp (transfer, transfer, simp add: vec_cnj_def)
    thus ?thesis
      using *
      by (subst ‹a = of_complex a'›) simp
  next
    case False
    hence "2 / sqrt (1 + (cmod a')⇧2) = r"
      using assms ‹a = of_complex a'›
      using dist_fs_infinite2[of a']
      by simp
    moreover
    have "sqrt (1 + (cmod a')⇧2) ≠ 0"
      using sqrt_1_plus_square
      by simp
    ultimately
    have "2 = r * sqrt (1 + (cmod a')⇧2)"
      by (simp add: field_simps)
    hence "4 = (r * sqrt (1 + (cmod a')⇧2))⇧2"
      by simp
    hence "4 = r⇧2 * (1 + (cmod a')⇧2)"
      by (simp add: power_mult_distrib)
    hence "Re (4 - (cor r)⇧2 * (1 + (a' * cnj a'))) = 0"
      by (subst (asm) cmod_square) (simp add: field_simps power2_eq_square)
    hence "4 - (cor r)⇧2 * (1 + (a'*cnj a')) = 0"
      by (subst complex_eq_if_Re_eq) (auto simp add: power2_eq_square)
    hence "circline_A0 (mk_circline ?A ?B ?C ?D)"
      using hh
      by (simp, transfer, transfer, simp)
    hence "z ∈ circline_set (mk_circline ?A ?B ?C ?D)"
      using inf_in_circline_set[of "mk_circline ?A ?B ?C ?D"]
      using ‹¬ z ≠ ∞⇩h›
      by simp
    thus ?thesis
      using *
      by (subst ‹a = of_complex a'›) simp
  qed
next
  case False
  let ?A = "-(cor r)⇧2" and ?B = "0" and ?C = "0" and ?D = "4 -(cor r)⇧2"
  have hh: "(?A, ?B, ?C, ?D) ∈ {H. hermitean H ∧ H ≠ mat_zero}"
    by (auto simp add: hermitean_def mat_adj_def mat_cnj_def power2_eq_square)
  hence *: "chordal_circle a r = mk_circline ?A ?B ?C ?D"
    using ‹¬ a ≠ ∞⇩h›
    by simp (transfer, transfer, simp, rule_tac x=1 in exI, simp)

  show ?thesis
  proof (cases "z = ∞⇩h")
    case True
    show ?thesis
      using assms ‹z = ∞⇩h› ‹¬ a ≠ ∞⇩h›
      using * hh
      by (simp, subst inf_in_circline_set, transfer, transfer, simp)
  next
    case False
    then obtain z' where "z = of_complex z'"
      using inf_or_of_complex[of z]
      by auto
    have "2 / sqrt (1 + (cmod z')⇧2) = r"
      using assms ‹z = of_complex z'›‹¬ a ≠ ∞⇩h›
      using dist_fs_infinite2[of z']
      by (simp add: dist_fs_sym)
    moreover
    have "sqrt (1 + (cmod z')⇧2) ≠ 0"
      using sqrt_1_plus_square
      by simp
    ultimately
    have "2 = r * sqrt (1 + (cmod z')⇧2)"
      by (simp add: field_simps)
    hence "4 = (r * sqrt (1 + (cmod z')⇧2))⇧2"
      by simp
    hence "4 = r⇧2 * (1 + (cmod z')⇧2)"
      by (simp add: power_mult_distrib)
    hence "Re (4 - (cor r)⇧2 * (1 + (z' * cnj z'))) = 0"
      by (subst (asm) cmod_square) (simp add: field_simps power2_eq_square)
    hence "- (cor r)⇧2 * z'*cnj z' + 4 - (cor r)⇧2 = 0"
      by (subst complex_eq_if_Re_eq) (auto simp add: power2_eq_square field_simps)
    hence "z ∈ circline_set (mk_circline ?A ?B ?C ?D)"
      using hh
      unfolding circline_set_def
      by (subst ‹z = of_complex z'›, simp) (transfer, transfer, auto simp add: vec_cnj_def field_simps)
    thus ?thesis
      using *
      by simp
  qed
qed

lemma
  assumes "z ∈ circline_set (chordal_circle a r)" and "r ≥ 0"
  shows "dist_fs z a = r"
proof (cases "a = ∞⇩h")
  case False
  then obtain a' where "a = of_complex a'"
    using inf_or_of_complex
    by auto

  let ?A = "4 - (cor r)⇧2 * (1 + (a'*cnj a'))" and ?B = "-4*a'" and ?C="-4*cnj a'" and ?D = "4*a'*cnj a' -  (cor r)⇧2 * (1 + (a'*cnj a'))"
  have hh: "(?A, ?B, ?C, ?D) ∈ {H. hermitean H ∧ H ≠ mat_zero}"
    by (auto simp add: hermitean_def mat_adj_def mat_cnj_def power2_eq_square)
  hence *: "chordal_circle (of_complex a') r = mk_circline ?A ?B ?C ?D"
    by (transfer, transfer, simp, rule_tac x=1 in exI, simp)

  show ?thesis
  proof (cases "z = ∞⇩h")
    case False
    then obtain z' where "z = of_complex z'"
      using inf_or_of_complex[of z] inf_or_of_complex[of a]
      by auto
    hence "z ∈ circline_set (mk_circline ?A ?B ?C ?D)"
      using assms ‹a = of_complex a'› *
      by simp
    hence "(cnj z' * ?A + ?C) * z' + (cnj z' * ?B + ?D) = 0"
      using hh
      unfolding circline_set_def
      by (subst (asm) ‹z = of_complex z'›, simp) (transfer, transfer, simp add: vec_cnj_def)
    hence "?A * z' * cnj z' + ?B * cnj z' + ?C * z' + ?D = 0"
      by algebra
    hence "Re (?A * z' * cnj z' + ?B * cnj z' +?C * z' +?D) = 0"
      by (simp add: power2_eq_square field_simps)
    hence "4 * Re ((z' - a')*cnj (z' - a')) = r⇧2 * (1 + Re (z'*cnj z')) * (1 + Re (a'*cnj a'))"
      by (simp add: field_simps power2_eq_square)
    hence "4 * (cmod (z' - a'))⇧2  = r⇧2 * ((1 + (cmod z')⇧2) * (1 + (cmod a')⇧2))"
      by (subst cmod_square)+ simp
    moreover
    have "sqrt (1 + (cmod z')⇧2) ≠ 0" "sqrt (1 + (cmod a')⇧2) ≠ 0"
      using sqrt_1_plus_square
      by simp+
    hence "(1 + (cmod z')⇧2) * (1 + (cmod a')⇧2) ≠ 0"
      by simp
    ultimately
    have "4 * (cmod (z' - a'))⇧2 / ((1 + (cmod z')⇧2) * (1 + (cmod a')⇧2)) = r⇧2 "
      by (simp add: field_simps)
    hence "2 * cmod (z' - a') / (sqrt (1 + (cmod z')⇧2) * sqrt (1 + (cmod a')⇧2)) = r"
      using ‹r ≥ 0›
      by (subst (asm) real_sqrt_eq_iff[symmetric]) (simp add: real_sqrt_mult real_sqrt_divide)
    thus ?thesis
      using ‹z = of_complex z'› ‹a = of_complex a'›
      using dist_fs_finite[of z' a']
      by simp
  next
    case True
    have "z ∈ circline_set (mk_circline ?A ?B ?C ?D)"
      using assms ‹a = of_complex a'› *
      by simp
    hence "circline_A0 (mk_circline ?A ?B ?C ?D)"
      using inf_in_circline_set[of "mk_circline ?A ?B ?C ?D"]
      using ‹z = ∞⇩h›
      by simp
    hence "4 - (cor r)⇧2 * (1 + (a'*cnj a')) = 0"
      using hh
      by (transfer, transfer, simp)
    hence "Re (4 - (cor r)⇧2 * (1 + (a' * cnj a'))) = 0"
      by simp
    hence "4 = r⇧2 * (1 + (cmod a')⇧2)"
      by (subst cmod_square) (simp add: power2_eq_square)
    hence "2 = r * sqrt (1 + (cmod a')⇧2)"
      using ‹r ≥ 0›
      by (subst (asm) real_sqrt_eq_iff[symmetric]) (simp add: real_sqrt_mult)
    moreover
    have "sqrt (1 + (cmod a')⇧2) ≠ 0"
      using sqrt_1_plus_square
      by simp
    ultimately
    have "2 / sqrt (1 + (cmod a')⇧2) = r"
      by (simp add: field_simps)
    thus ?thesis
      using ‹a = of_complex a'› ‹z = ∞⇩h›
      using dist_fs_infinite2[of a']
      by simp
  qed
next
  case True
  let ?A = "-(cor r)⇧2" and ?B = "0" and ?C = "0" and ?D = "4 -(cor r)⇧2"
  have hh: "(?A, ?B, ?C, ?D) ∈ {H. hermitean H ∧ H ≠ mat_zero}"
    by (auto simp add: hermitean_def mat_adj_def mat_cnj_def power2_eq_square)
  hence *: "chordal_circle a r = mk_circline ?A ?B ?C ?D"
    using ‹a = ∞⇩h›
    by simp (transfer, transfer, simp, rule_tac x=1 in exI, simp)

  show ?thesis
  proof (cases "z = ∞⇩h")
    case True
    thus ?thesis
      using ‹a = ∞⇩h› assms * hh
      by simp (subst (asm) inf_in_circline_set, transfer, transfer, simp)
  next
    case False
    then obtain z' where "z = of_complex z'"
      using inf_or_of_complex
      by auto
    hence "z ∈ circline_set (mk_circline ?A ?B ?C ?D)"
      using assms *
      by simp
    hence "- (cor r)⇧2 * z'*cnj z' + 4 - (cor r)⇧2 = 0"
      using hh
      unfolding circline_set_def
      apply (subst (asm) ‹z = of_complex z'›)
      by (simp, transfer, transfer, simp add: vec_cnj_def, algebra)
    hence "4 - (cor r)⇧2 * (1 + (z'*cnj z')) = 0"
      by (simp add: field_simps)
    hence "Re (4 - (cor r)⇧2 * (1 + (z' * cnj z'))) = 0"
      by simp
    hence "4 = r⇧2 * (1 + (cmod z')⇧2)"
      by (subst cmod_square) (simp add: power2_eq_square)
    hence "2 = r * sqrt (1 + (cmod z')⇧2)"
      using ‹r ≥ 0›
      by (subst (asm) real_sqrt_eq_iff[symmetric]) (simp add: real_sqrt_mult)
    moreover
    have "sqrt (1 + (cmod z')⇧2) ≠ 0"
      using sqrt_1_plus_square
      by simp
    ultimately
    have "2 / sqrt (1 + (cmod z')⇧2) = r"
      by (simp add: field_simps)
    thus ?thesis
      using ‹z = of_complex z'› ‹a = ∞⇩h›
      using dist_fs_infinite2[of z']
      by (simp add: dist_fs_sym)
  qed
qed

text ‹Two chordal centers and radii for the given circline›
definition chordal_circles_cmat :: "complex_mat ⇒ (complex × real) × (complex × real)"  where
 [simp]: "chordal_circles_cmat H =
     (let (A, B, C, D) = H;
          dsc = sqrt(Re ((D-A)⇧2 + 4 * (B*cnj B)));
          a1 = (A - D + cor dsc) / (2 * C);
          r1 = sqrt((4 - Re((-4 * a1/B) * A)) / (1 + Re (a1*cnj a1)));
          a2 = (A - D - cor dsc) / (2 * C);
          r2 = sqrt((4 - Re((-4 * a2/B) * A)) / (1 + Re (a2*cnj a2)))
       in ((a1, r1), (a2, r2)))"

lift_definition chordal_circles_clmat :: "circline_mat ⇒ (complex × real) × (complex × real)" is chordal_circles_cmat
  done

lift_definition chordal_circles :: "ocircline ⇒ (complex × real) × (complex × real)" is chordal_circles_clmat
proof transfer
  fix H1 H2 :: complex_mat
  obtain A1 B1 C1 D1 where hh1: "(A1, B1, C1, D1) = H1"
    by (cases H1) auto
  obtain A2 B2 C2 D2 where hh2: "(A2, B2, C2, D2) = H2"
    by (cases H2) auto

  assume "ocircline_eq_cmat H1 H2"
  then obtain k where *: "k > 0" "A2 = cor k * A1" "B2 = cor k * B1" "C2 = cor k * C1" "D2 = cor k * D1"
    using hh1[symmetric] hh2[symmetric]
    by auto
  let ?dsc1 = "sqrt (Re ((D1 - A1)⇧2 + 4 * (B1 * cnj B1)))" and ?dsc2 = "sqrt (Re ((D2 - A2)⇧2 + 4 * (B2 * cnj B2)))"
  let ?a11 = "(A1 - D1 + cor ?dsc1) / (2 * C1)" and ?a12 = "(A2 - D2 + cor ?dsc2) / (2 * C2)"
  let ?a21 = "(A1 - D1 - cor ?dsc1) / (2 * C1)" and ?a22 = "(A2 - D2 - cor ?dsc2) / (2 * C2)"
  let ?r11 = "sqrt((4 - Re((-4 * ?a11/B1) * A1)) / (1 + Re (?a11*cnj ?a11)))"
  let ?r12 = "sqrt((4 - Re((-4 * ?a12/B2) * A2)) / (1 + Re (?a12*cnj ?a12)))"
  let ?r21 = "sqrt((4 - Re((-4 * ?a21/B1) * A1)) / (1 + Re (?a21*cnj ?a21)))"
  let ?r22 = "sqrt((4 - Re((-4 * ?a22/B2) * A2)) / (1 + Re (?a22*cnj ?a22)))"

  have "Re ((D2 - A2)⇧2 + 4 * (B2 * cnj B2)) = k⇧2 * Re ((D1 - A1)⇧2 + 4 * (B1 * cnj B1))"
    using *
    by (simp add: power2_eq_square field_simps)
  hence "?dsc2 = k * ?dsc1"
    using ‹k > 0›
    by (simp add: real_sqrt_mult)
  hence "A2 - D2 + cor ?dsc2 = cor k * (A1 - D1 + cor ?dsc1)" "A2 - D2 - cor ?dsc2 = cor k * (A1 - D1 - cor ?dsc1)" "2*C2 = cor k * (2*C1)"
    using *
    by (auto simp add: field_simps)
  hence "?a12 = ?a11" "?a22 = ?a21"
    using ‹k > 0›
    by simp_all
  moreover
  have "Re((-4 * ?a12/B2) * A2) = Re((-4 * ?a11/B1) * A1)"
    using *
    by (subst ‹?a12 = ?a11›) (simp, simp add: field_simps)
  have "?r12 = ?r11"
    by (subst ‹Re((-4 * ?a12/B2) * A2) = Re((-4 * ?a11/B1) * A1)›, (subst ‹?a12 = ?a11›)+) simp
  moreover
  have "Re((-4 * ?a22/B2) * A2) = Re((-4 * ?a21/B1) * A1)"
    using *
    by (subst ‹?a22 = ?a21›) (simp, simp add: field_simps)
  have "?r22 = ?r21"
    by (subst ‹Re((-4 * ?a22/B2) * A2) = Re((-4 * ?a21/B1) * A1)›, (subst ‹?a22 = ?a21›)+) simp
  moreover
  have "chordal_circles_cmat H1 = ((?a11, ?r11), (?a21, ?r21))"
    using hh1[symmetric]
    unfolding chordal_circles_cmat_def Let_def
    by (simp del: times_complex.sel)
  moreover
  have "chordal_circles_cmat H2 = ((?a12, ?r12), (?a22, ?r22))"
    using hh2[symmetric]
    unfolding chordal_circles_cmat_def Let_def
    by (simp del: times_complex.sel)
  ultimately
  show "chordal_circles_cmat H1 = chordal_circles_cmat H2"
    by metis
qed

lemma chordal_circle_radius_positive:
  assumes "hermitean (A, B, C, D)" and "Re (mat_det (A, B, C, D)) ≤ 0" and "B ≠ 0" and
  "dsc = sqrt(Re ((D-A)⇧2 + 4 * (B*cnj B)))" and 
  "a1 = (A - D + cor dsc) / (2 * C)" and "a2 = (A - D - cor dsc) / (2 * C)"
  shows "Re (A*a1/B) ≥ -1 ∧ Re (A*a2/B) ≥ -1"
proof-
  from assms have "is_real A" "is_real D" "C = cnj B"
    using hermitean_elems
    by auto
  have *: "A*a1/B = ((A - D + cor dsc) / (2 * (B * cnj B))) * A"
    using ‹B ≠ 0› ‹C = cnj B› ‹a1 = (A - D + cor dsc) / (2 * C)›
    by (simp add: field_simps) algebra
  have **: "A*a2/B = ((A - D - cor dsc) / (2 * (B * cnj B))) * A"
    using ‹B ≠ 0› ‹C = cnj B› ‹a2 = (A - D - cor dsc) / (2 * C)›
    by (simp add: field_simps) algebra
  have "dsc ≥ 0"
  proof-
    have "0 ≤ Re ((D - A)⇧2) + 4 * Re ((cor (cmod B))⇧2)"
      using ‹is_real A› ‹is_real D› by simp
    thus ?thesis
      using ‹dsc = sqrt(Re ((D-A)⇧2 + 4*(B*cnj B)))›
      by (subst (asm) complex_mult_cnj_cmod) simp
  qed
  hence "Re (A - D - cor dsc) ≤ Re (A - D + cor dsc)"
    by simp
  moreover
  have "Re (2 * (B * cnj B)) > 0"
    using ‹B ≠ 0›
    by (subst complex_mult_cnj_cmod, simp add: power2_eq_square)
  ultimately
  have xxx: "Re (A - D + cor dsc) / Re (2 * (B * cnj B)) ≥ Re (A - D - cor dsc) / Re (2 * (B * cnj B))" (is "?lhs ≥ ?rhs")
    by (metis divide_right_mono less_eq_real_def)

  have "Re A * Re D ≤ Re (B*cnj B)"
    using ‹Re (mat_det (A, B, C, D)) ≤ 0› ‹C = cnj B› ‹is_real A› ‹is_real D›
    by simp


  have "(Re (2 * (B * cnj B)) / Re A) / Re (2 * B * cnj B) = 1 / Re A"
    using ‹Re (2 * (B * cnj B)) > 0›
    apply (subst divide_divide_eq_left)
    apply (subst mult.assoc)
    apply (subst nonzero_divide_mult_cancel_right)
    by simp_all

  show ?thesis
  proof (cases "Re A > 0")
    case True
    hence "Re (A*a1/B) ≥ Re (A*a2/B)"
      using * ** ‹Re (2 * (B * cnj B)) > 0› ‹B ≠ 0› ‹is_real A› ‹is_real D› xxx
      using mult_right_mono[of ?rhs ?lhs "Re A"]
      apply simp
      apply (subst Re_divide_real, simp, simp)
      apply (subst Re_divide_real, simp, simp)
      apply (subst Re_mult_real, simp)+
      apply simp
      done
    moreover
    have "Re (A*a2/B) ≥ -1"
    proof-
      from ‹Re A * Re D ≤ Re (B*cnj B)›
      have "Re (A⇧2) ≤ Re (B*cnj B) + Re ((A - D)*A)"
        using ‹Re A > 0› ‹is_real A› ‹is_real D›
        by (simp add: power2_eq_square field_simps)
      have "1 ≤ Re (B*cnj B) / Re (A⇧2) + Re (A - D) / Re A"
        using ‹Re A > 0› ‹is_real A› ‹is_real D›
        using divide_right_mono[OF ‹Re (A⇧2) ≤ Re (B*cnj B) + Re ((A - D)*A)›, of "Re (A⇧2)"]
        by (simp add: power2_eq_square add_divide_distrib)
      have "4 * Re(B*cnj B) ≤ 4 * (Re (B*cnj B))⇧2 / Re (A⇧2)  + 2*Re (A - D) / Re A * 2 * Re(B*cnj B)"
        using mult_right_mono[OF ‹1 ≤ Re (B*cnj B) / Re (A⇧2) + Re (A - D) / Re A›, of "4 * Re (B*cnj B)"]
        by (simp add: distrib_right) (simp add: power2_eq_square field_simps)
      moreover
      have "A ≠ 0"
        using ‹Re A > 0›
        by auto
      hence "4 * (Re (B*cnj B))⇧2 / Re (A⇧2) = Re (4 * (B*cnj B)⇧2 / A⇧2)"
        using Re_divide_real[of "A⇧2" "4 * (B*cnj B)⇧2"] ‹Re A > 0› ‹is_real A›
        by (auto simp add: power2_eq_square)
      moreover
      have "2*Re (A - D) / Re A * 2 * Re(B*cnj B) = Re (2 * (A - D) / A * 2 * B * cnj B)"
        using ‹is_real A› ‹is_real D› ‹A ≠ 0›
        using Re_divide_real[of "A" "(4 * A - 4 * D) * B * cnj B"]
        by (simp add: field_simps)
      ultimately
      have "Re ((D - A)⇧2 + 4 * B*cnj B) ≤ Re((A - D)⇧2 + 4 * (B*cnj B)⇧2 / A⇧2  + 2*(A - D) / A * 2 * B*cnj B)"
        by (simp add: field_simps power2_eq_square)
      hence "Re ((D - A)⇧2 + 4 * B*cnj B) ≤ Re(((A - D) +  2 * B*cnj B / A)⇧2)"
        using ‹A ≠ 0›
        by (subst power2_sum) (simp add: power2_eq_square field_simps)
      hence "dsc ≤ sqrt (Re(((A - D) +  2 * B*cnj B / A)⇧2))"
        using ‹dsc = sqrt(Re ((D-A)⇧2 + 4*(B*cnj B)))›
        by simp
      moreover
      have "Re(((A - D) +  2 * B*cnj B / A)⇧2) = (Re((A - D) +  2 * B*cnj B / A))⇧2"
        using ‹is_real A› ‹is_real D› div_reals
        by (simp add: power2_eq_square)
      ultimately
      have "dsc ≤ ¦Re (A - D + 2 * B * cnj B / A)¦"
        by simp
      moreover
      have "Re (A - D + 2 * B * cnj B / A) ≥ 0"
      proof-
        have *: "Re (A⇧2 + B*cnj B) ≥ 0"
          using ‹is_real A›
          by (simp add: power2_eq_square)
        also have "Re (A⇧2 + 2*B*cnj B - A*D) ≥ Re (A⇧2 + B*cnj B)"
          using ‹Re A * Re D ≤ Re (B*cnj B)›
          using ‹is_real A› ‹is_real D›
          by simp
        finally
        have "Re (A⇧2 + 2*B*cnj B - A*D) ≥ 0"
          by simp
        show ?thesis
          using divide_right_mono[OF ‹Re (A⇧2 + 2*B*cnj B - A*D) ≥ 0›, of "Re A"] ‹Re A > 0› ‹is_real A› ‹A ≠ 0›
          by (simp add: add_divide_distrib diff_divide_distrib) (subst Re_divide_real, auto simp add: power2_eq_square field_simps)
      qed
      ultimately
      have "dsc ≤ Re (A - D + 2 * B * cnj B / A)"
        by simp
      hence "- Re (2 * (B * cnj B) / A) ≤ Re ((A - D - cor dsc))"
        by (simp add: field_simps)
      hence *: "- (Re (2 * (B * cnj B)) / Re A) ≤ Re (A - D - cor dsc)"
        using ‹is_real A› ‹A ≠ 0›
        by (subst (asm) Re_divide_real, auto)
      from divide_right_mono[OF this, of "Re (2 * B * cnj B)"]
      have "- 1 / Re A ≤ Re (A - D - cor dsc) / Re (2 * B * cnj B)"
        using ‹Re A > 0› ‹B ≠ 0› ‹A ≠ 0› ‹0 < Re (2 * (B * cnj B))›
        using ‹(Re (2 * (B * cnj B)) / Re A) / Re (2 * B * cnj B) = 1 / Re A›
        by simp
      from mult_right_mono[OF this, of "Re A"]
      show ?thesis
        using ‹is_real A› ‹is_real D› ‹B ≠ 0› ‹Re A > 0› ‹A ≠ 0›
        apply (subst **)
        apply (subst Re_mult_real, simp)
        apply (subst Re_divide_real, simp, simp)
        apply (simp add: field_simps)
        done
    qed
    ultimately
    show ?thesis
      by simp
  next
    case False
    show ?thesis
    proof (cases "Re A < 0")
      case True
      hence "Re (A*a1/B) ≤ Re (A*a2/B)"
        using * ** ‹Re (2 * (B * cnj B)) > 0› ‹B ≠ 0› ‹is_real A› ‹is_real D› xxx
        using mult_right_mono_neg[of ?rhs ?lhs "Re A"]
        apply simp
        apply (subst Re_divide_real, simp, simp)
        apply (subst Re_divide_real, simp, simp)
        apply (subst Re_mult_real, simp)+
        apply simp
        done
      moreover
      have "Re (A*a1/B) ≥ -1"
      proof-
        from ‹Re A * Re D ≤ Re (B*cnj B)›
        have "Re (A⇧2) ≤ Re (B*cnj B) - Re ((D - A)*A)"
          using ‹Re A < 0› ‹is_real A› ‹is_real D›
          by (simp add: power2_eq_square field_simps)
        hence "1 ≤ Re (B*cnj B) / Re (A⇧2) - Re (D - A) / Re A"
          using ‹Re A < 0› ‹is_real A› ‹is_real D›
          using divide_right_mono[OF ‹Re (A⇧2) ≤ Re (B*cnj B) - Re ((D - A)*A)›, of "Re (A⇧2)"]
          by (simp add: power2_eq_square diff_divide_distrib)
        have "4 * Re(B*cnj B) ≤ 4 * (Re (B*cnj B))⇧2 / Re (A⇧2)  - 2*Re (D - A) / Re A * 2 * Re(B*cnj B)"
          using mult_right_mono[OF ‹1 ≤ Re (B*cnj B) / Re (A⇧2) - Re (D - A) / Re A›, of "4 * Re (B*cnj B)"]
          by (simp add: left_diff_distrib) (simp add: power2_eq_square field_simps)
        moreover
        have "A ≠ 0"
          using ‹Re A < 0›
          by auto
        hence "4 * (Re (B*cnj B))⇧2 / Re (A⇧2) = Re (4 * (B*cnj B)⇧2 / A⇧2)"
          using Re_divide_real[of "A⇧2" "4 * (B*cnj B)⇧2"] ‹Re A < 0› ‹is_real A›
          by (auto simp add: power2_eq_square)
        moreover
        have "2*Re (D - A) / Re A * 2 * Re(B*cnj B) = Re (2 * (D - A) / A * 2 * B * cnj B)"
          using ‹is_real A› ‹is_real D› ‹A ≠ 0›
          using Re_divide_real[of "A" "(4 * D - 4 * A) * B * cnj B"]
          by (simp add: field_simps)
        ultimately
        have "Re ((D - A)⇧2 + 4 * B*cnj B) ≤ Re((D - A)⇧2 + 4 * (B*cnj B)⇧2 / A⇧2  - 2*(D - A) / A * 2 * B*cnj B)"
          by (simp add: field_simps power2_eq_square)
        hence "Re ((D - A)⇧2 + 4 * B*cnj B) ≤ Re(((D - A) -  2 * B*cnj B / A)⇧2)"
          using ‹A ≠ 0›
          by (subst power2_diff) (simp add: power2_eq_square field_simps)
        hence "dsc ≤ sqrt (Re(((D - A) -  2 * B*cnj B / A)⇧2))"
          using ‹dsc = sqrt(Re ((D-A)⇧2 + 4*(B*cnj B)))›
          by simp
        moreover
        have "Re(((D - A) -  2 * B*cnj B / A)⇧2) = (Re((D - A) -  2 * B*cnj B / A))⇧2"
          using ‹is_real A› ‹is_real D› div_reals
          by (simp add: power2_eq_square)
        ultimately
        have "dsc ≤ ¦Re (D - A - 2 * B * cnj B / A)¦"
          by simp
        moreover
        have "Re (D - A - 2 * B * cnj B / A) ≥ 0"
        proof-
          have "Re (A⇧2 + B*cnj B) ≥ 0"
            using ‹is_real A›
            by (simp add: power2_eq_square)
          also have "Re (A⇧2 + 2*B*cnj B - A*D) ≥ Re (A⇧2 + B*cnj B)"
            using ‹Re A * Re D ≤ Re (B*cnj B)›
            using ‹is_real A› ‹is_real D›
            by simp
          finally have "Re (A⇧2 + 2*B*cnj B - A*D) ≥ 0"
            by simp
          show ?thesis
            using divide_right_mono_neg[OF ‹Re (A⇧2 + 2*B*cnj B - A*D) ≥ 0›, of "Re A"] ‹Re A < 0› ‹is_real A› ‹A ≠ 0›
            by (simp add: add_divide_distrib diff_divide_distrib) (subst Re_divide_real, auto simp add: power2_eq_square field_simps)
        qed
        ultimately
        have "dsc ≤ Re (D - A - 2 * B * cnj B / A)"
          by simp
        hence "- Re (2 * (B * cnj B) / A) ≥ Re ((A - D + cor dsc))"
          by (simp add: field_simps)
        hence "- (Re (2 * (B * cnj B)) / Re A) ≥ Re (A - D + cor dsc)"
          using ‹is_real A› ‹A ≠ 0›
          by (subst (asm) Re_divide_real, auto)
        from divide_right_mono[OF this, of "Re (2 * B * cnj B)"]
        have "- 1 / Re A ≥ Re (A - D + cor dsc) / Re (2 * B * cnj B)"
          using ‹Re A < 0› ‹B ≠ 0› ‹A ≠ 0› ‹0 < Re (2 * (B * cnj B))›
          using ‹(Re (2 * (B * cnj B)) / Re A) / Re (2 * B * cnj B) = 1 / Re A›
          by simp
        from mult_right_mono_neg[OF this, of "Re A"]
        show ?thesis
          using ‹is_real A› ‹is_real D› ‹B ≠ 0› ‹Re A < 0› ‹A ≠ 0›
          apply (subst *)
          apply (subst Re_mult_real, simp)
          apply (subst Re_divide_real, simp, simp)
          apply (simp add: field_simps)
          done
      qed
      ultimately
      show ?thesis
        by simp
    next
      case False
      hence "A = 0"
        using ‹¬ Re A > 0› ‹is_real A›
        using complex_eq_if_Re_eq by auto
      thus ?thesis
        by simp
    qed
  qed
qed


lemma chordal_circle_det_positive:
  fixes x y :: real
  assumes "x * y < 0"
  shows "x / (x - y) > 0"
proof (cases "x > 0")
  case True
  hence "y < 0"
    using ‹x * y < 0›
    by (smt mult_nonneg_nonneg)
  have "x - y > 0"
    using ‹x > 0› ‹y < 0›
    by auto
  thus ?thesis
    using ‹x > 0›
    by (metis zero_less_divide_iff)
next
  case False
  hence *: "y > 0 ∧ x < 0"
    using ‹x * y < 0›
    using mult_nonpos_nonpos[of x y]
    by (cases "x=0") force+

  have "x - y < 0"
    using *
    by auto
  thus ?thesis
    using *
    by (metis zero_less_divide_iff)
qed

lemma cor_sqrt_squared: "x ≥ 0 ⟹ (cor (sqrt x))⇧2 = cor x"
  by (simp add: power2_eq_square)

lemma chordal_circle1:
  assumes "is_real A" and "is_real D" and "Re (A * D) < 0" and "r = sqrt(Re ((4*A)/(A-D)))"
  shows "mk_circline A 0 0 D = chordal_circle ∞⇩h r"
using assms
proof (transfer, transfer)
  fix A D r
  assume *: "is_real A" "is_real D" "Re (A * D) < 0" "r = sqrt (Re ((4*A)/(A-D)))"
  hence "A ≠ 0 ∨ D ≠ 0"
    by auto
  hence "(A, 0, 0, D) ∈ hermitean_nonzero"
    using eq_cnj_iff_real[of A] eq_cnj_iff_real[of D] *
    unfolding hermitean_def
    by (simp add: mat_adj_def mat_cnj_def)
  moreover
  have "(- (cor r)⇧2, 0, 0, 4 - (cor r)⇧2) ∈ hermitean_nonzero"
    by (simp add: hermitean_def mat_adj_def mat_cnj_def power2_eq_square)
  moreover
  have "A ≠ D"
    using ‹Re (A * D) < 0› ‹is_real A› ‹is_real D›
    by auto
  have "Re ((4*A)/(A-D)) ≥ 0"
  proof-
    have "Re A / Re (A - D) ≥ 0"
      using ‹Re (A * D) < 0› ‹is_real A› ‹is_real D›
      using chordal_circle_det_positive[of "Re A" "Re D"]
      by simp
    thus ?thesis
      using ‹is_real A› ‹is_real D› ‹A ≠ D›
      by (subst Re_divide_real, auto)
  qed
  moreover
  have "- (cor (sqrt (Re (4 * A / (A - D)))))⇧2 = cor (Re (4 / (D - A))) * A"
    using ‹is_real A› ‹is_real D› ‹A ≠ D› ‹Re ((4*A)/(A-D)) ≥ 0›
    by (simp add: cor_sqrt_squared field_simps)
  moreover
  have "4 - 4 * A / (A - D) = 4 * D / (D - A)"
    using‹A ≠ D› 
    by (simp add: divide_simps split: if_split_asm) (simp add: minus_mult_right)
  hence **: "4 - (cor (sqrt (Re (4 * A / (A - D)))))⇧2 = cor (Re (4 / (D - A))) * D"
    using ‹Re ((4*A)/(A-D)) ≥ 0› ‹is_real A› ‹is_real D› ‹A ≠ D›
    by (simp add: cor_sqrt_squared field_simps)
  ultimately
  show "circline_eq_cmat (mk_circline_cmat A 0 0 D) (chordal_circle_cvec_cmat ∞⇩v r)"
    using * ‹is_real A› ‹is_real D› ‹A ≠ D› ‹r = sqrt(Re ((4*A)/(A-D)))›
    by (simp, rule_tac x="Re(4/(D-A))" in exI, auto, simp_all add: **)
qed

lemma chordal_circle2:
  assumes "is_real A" and "is_real D" and "Re (A * D) < 0" and "r = sqrt(Re ((4*D)/(D-A)))"
  shows "mk_circline A 0 0 D = chordal_circle 0⇩h r"
using assms
proof (transfer, transfer)
  fix A D r
  assume *: "is_real A" "is_real D" "Re (A * D) < 0" "r = sqrt (Re ((4*D)/(D-A)))"
  hence "A ≠ 0 ∨ D ≠ 0"
    by auto
  hence "(A, 0, 0, D) ∈ {H. hermitean H ∧ H ≠ mat_zero}"
    using eq_cnj_iff_real[of A] eq_cnj_iff_real[of D] *
    unfolding hermitean_def
    by (simp add: mat_adj_def mat_cnj_def)
  moreover
  have "(4 - (cor r)⇧2, 0, 0, - (cor r)⇧2) ∈ {H. hermitean H ∧ H ≠ mat_zero}"
    by (auto simp add: hermitean_def mat_adj_def mat_cnj_def power2_eq_square)
  moreover
  have "A ≠ D"
    using ‹Re (A * D) < 0› ‹is_real A› ‹is_real D›
    by auto
  have "Re((4*D)/(D-A)) ≥ 0"
  proof-
    have "Re D / Re (D - A) ≥ 0"
      using ‹Re (A * D) < 0› ‹is_real A› ‹is_real D›
      using chordal_circle_det_positive[of "Re D" "Re A"]
      by (simp add: field_simps)
    thus ?thesis
      using ‹is_real A› ‹is_real D› ‹A ≠ D› Re_divide_real by force
  qed
  have  "4 - 4 * D / (D - A) = 4 * A / (A - D)"
    by (simp add: divide_simps split: if_split_asm) (simp add: ‹A ≠ D› minus_mult_right)
  hence **: "4 - (cor (sqrt (Re ((4*D)/(D-A)))))⇧2 = cor (Re (4 / (A - D))) * A"
    using ‹is_real A› ‹is_real D› ‹A ≠ D› ‹Re (4 * D / (D - A)) ≥ 0›
    by (simp add: cor_sqrt_squared field_simps)
  moreover
  have "- (cor (sqrt (Re ((4*D)/(D-A)))))⇧2 = cor (Re (4 / (A - D))) * D"
    using ‹is_real A› ‹is_real D› ‹A ≠ D› ‹Re (4 * D / (D - A)) ≥ 0›
    by (simp add: cor_sqrt_squared field_simps)
  ultimately
  show "circline_eq_cmat (mk_circline_cmat A 0 0 D) (chordal_circle_cvec_cmat 0⇩v r)"
    using ‹is_real A› ‹is_real D› ‹A ≠ 0 ∨ D ≠ 0› ‹r = sqrt (Re ((4*D)/(D-A)))›
    using *
    by (simp, rule_tac x="Re (4/(A-D))" in exI, auto, simp_all add: **)
qed

lemma chordal_circle':
  assumes "B ≠ 0" and "(A, B, C, D) ∈ hermitean_nonzero" and "Re (mat_det (A, B, C, D)) ≤ 0" and
  "C * a⇧2  + (D - A) * a - B = 0" and "r = sqrt((4 - Re((-4 * a/B) * A)) / (1 + Re (a*cnj a)))"
  shows "mk_circline A B C D = chordal_circle (of_complex a) r"
using assms
proof (transfer, transfer)
  fix A B C D a :: complex and r :: real

  let ?k = "(-4) * a / B"

  assume *: "(A, B, C, D) ∈ {H. hermitean H ∧ H ≠ mat_zero}" and **: "B ≠ 0" "C * a⇧2 + (D - A) * a - B = 0" and rr: "r = sqrt ((4 - Re (?k * A)) / (1 + Re (a * cnj a)))" and det: "Re (mat_det (A, B, C, D)) ≤ 0"

  have "is_real A" "is_real D" "C = cnj B"
    using * hermitean_elems
    by auto

  from ** have a12: "let dsc = sqrt(Re ((D-A)⇧2 + 4 * (B*cnj B)))
                      in a = (A - D + cor dsc) / (2 * C) ∨ a = (A - D - cor dsc) / (2 * C)"
  proof-
    have "Re ((D-A)⇧2 + 4 * (B*cnj B)) ≥ 0"
      using ‹is_real A› ‹is_real D›
      by (subst complex_mult_cnj_cmod) (simp add: power2_eq_square)
    hence "ccsqrt ((D - A)⇧2 - 4 * C * - B) = cor (sqrt (Re ((D - A)⇧2 + 4 * (B * cnj B))))"
      using csqrt_real[of "((D - A)⇧2 + 4 * (B * cnj B))"] ‹is_real A› ‹is_real D› ‹C = cnj B›
      by (auto simp add: power2_eq_square field_simps)
    thus ?thesis
      using complex_quadratic_equation_two_roots[of C a "D - A" "-B"]
      using  ‹C * a⇧2 + (D - A) * a - B = 0› ‹B ≠ 0› ‹C = cnj B›
      by (simp add: Let_def)
  qed

  have "is_real ?k"
    using a12 ‹C = cnj B› ‹is_real A› ‹is_real D›
    by (auto simp add: Let_def)
  have "a ≠ 0"
    using **
    by auto
  hence "Re ?k ≠ 0"
    using ‹is_real (-4*a / B)› ‹B ≠ 0›
    by (metis complex.expand divide_eq_0_iff divisors_zero zero_complex.simps(1) zero_complex.simps(2) zero_neq_neg_numeral)
  moreover
  have "(-4) * a = cor (Re ?k) * B"
    using complex_of_real_Re[OF ‹is_real (-4*a/B)›] ‹B ≠ 0›
    by simp
  moreover
  have "is_real (a/B)"
    using ‹is_real ?k› is_real_mult_real[of "-4" "a / B"]
    by simp
  hence "is_real (B * cnj a)"
    using * ‹C = cnj B›
    by (metis (no_types, lifting) Im_complex_div_eq_0 complex_cnj_divide eq_cnj_iff_real hermitean_elems(3) mem_Collect_eq mult.commute)
  hence "B * cnj a = cnj B * a"
    using eq_cnj_iff_real[of "B * cnj a"]
    by simp
  hence "-4 * cnj a = cor (Re ?k) * C"
    using ‹C = cnj B›
    using complex_of_real_Re[OF ‹is_real ?k›] ‹B ≠ 0›
    by (simp, simp add: field_simps)
  moreover
  have "1 + a * cnj a ≠ 0"
    by (simp add: complex_mult_cnj_cmod)
  have "r⇧2 = (4 - Re (?k * A)) / (1 + Re (a * cnj a))"
  proof-
    have "Re (a / B * A) ≥ -1"
      using a12 chordal_circle_radius_positive[of A B C D] * ‹B ≠ 0› det
      by (auto simp add: Let_def field_simps)
    from mult_right_mono_neg[OF this, of "-4"]
    have "4 - Re (?k * A) ≥ 0"
      using Re_mult_real[of "-4" "a / B * A"]
      by (simp add: field_simps)
    moreover
    have "1 + Re (a * cnj a) > 0"
      using ‹a ≠ 0› complex_mult_cnj complex_neq_0
      by auto
    ultimately
    have "(4 - Re (?k * A)) / (1 + Re (a * cnj a)) ≥ 0"
      by (metis divide_nonneg_pos)
    thus ?thesis
      using rr
      by simp
  qed
  hence "r⇧2 = Re ((4 - ?k * A) / (1 + a * cnj a))"
    using ‹is_real ?k› ‹is_real A› ‹1 + a * cnj a ≠ 0›
    by (subst Re_divide_real, auto)
  hence "(cor r)⇧2 =  (4 - ?k * A) / (1 + a * cnj a)"
    using ‹is_real ?k› ‹is_real A› mult_reals[of ?k A] 
    by (simp add: cor_squared)
  hence "4 - (cor r)⇧2 * (a * cnj a + 1) = cor (Re ?k) * A"
    using complex_of_real_Re[OF ‹is_real (-4*a/B)›]
    using ‹1 + a * cnj a ≠ 0›
    by (simp add: field_simps)
  moreover

  have "?k = cnj ?k"
    using ‹is_real ?k›
    using eq_cnj_iff_real[of "-4*a/B"]
    by simp

  have "?k⇧2 = cor ((cmod ?k)⇧2)"
    using  cor_cmod_real[OF ‹is_real ?k›]
    unfolding power2_eq_square by force
  hence "?k⇧2 = ?k * cnj ?k"
    using complex_mult_cnj_cmod[of ?k]
    by simp
  hence ***: "a * cnj a = (cor ((Re ?k)⇧2) * B * C) / 16"
    using complex_of_real_Re[OF ‹is_real (-4*a/B)›] ‹C = cnj B› ‹is_real (-4*a/B)› ‹B ≠ 0›
    by simp
  from ** have "cor ((Re ?k)⇧2) * B * C - 4 * cor (Re ?k) * (D-A) - 16 = 0"
    using complex_of_real_Re[OF ‹is_real ?k›]
    by (simp add: power2_eq_square, simp add: field_simps, algebra)
  hence "?k * (D-A) = 4 * (cor ((Re ?k)⇧2) * B * C / 16 - 1)"
    by (subst (asm) complex_of_real_Re[OF ‹is_real ?k›]) algebra
  hence "?k * (D-A) = 4 * (a*cnj a - 1)"
    by (subst (asm)  ***[symmetric]) simp

  hence "4 * a * cnj a - (cor r)⇧2 * (a * cnj a + 1) = cor (Re ?k) * D"
    using ‹4 - (cor r)⇧2 * (a * cnj a + 1) = cor (Re ?k) * A›
    using complex_of_real_Re[OF ‹is_real (-4*a/B)›]
    by simp algebra
  ultimately
  show "circline_eq_cmat (mk_circline_cmat A B C D) (chordal_circle_cvec_cmat (of_complex_cvec a) r)"
    using * ‹a ≠ 0›
    by (simp, rule_tac x="Re (-4*a / B)" in exI, simp)
qed

end