Theory Quaternions

section‹Theory of Quaternions›
text‹This theory is inspired by the HOL Light development of quaternions, but follows its own route. 
Quaternions are developed coinductively, as in the existing formalisation of the complex numbers. 
Quaternions are quickly shown to belong to the type classes of real normed division algebras 
and real inner product spaces. And therefore they inherit a great body of facts involving algebraic 
laws, limits, continuity, etc., which must be proved explicitly in the HOL Light version. 
The development concludes with the geometric interpretation of the product of imaginary quaternions.›
theory "Quaternions"
  imports
    "HOL-Analysis.Multivariate_Analysis"
begin

subsection‹Basic definitions›

text‹As with the complex numbers, coinduction is convenient›

codatatype quat = Quat (Re: real) (Im1: real) (Im2: real) (Im3: real)

lemma quat_eqI [intro?]: "⟦Re x = Re y; Im1 x = Im1 y; Im2 x = Im2 y; Im3 x = Im3 y⟧ ⟹ x = y"
  by (rule quat.expand) simp

lemma quat_eq_iff: "x = y ⟷ Re x = Re y ∧ Im1 x = Im1 y ∧ Im2 x = Im2 y ∧ Im3 x = Im3 y"
  by (auto intro: quat.expand)

context
begin
no_notation Complex.imaginary_unit ("𝗂")

primcorec quat_ii :: quat  ("𝗂")
  where "Re 𝗂 = 0" | "Im1 𝗂 = 1" | "Im2 𝗂 = 0" | "Im3 𝗂 = 0"

primcorec quat_jj :: quat  ("𝗃")
  where "Re 𝗃 = 0" | "Im1 𝗃 = 0" | "Im2 𝗃 = 1" | "Im3 𝗃 = 0"

primcorec quat_kk :: quat  ("𝗄")
  where "Re 𝗄 = 0" | "Im1 𝗄 = 0" | "Im2 𝗄 = 0" | "Im3 𝗄 = 1"

end

bundle quaternion_syntax begin
notation quat_ii ("𝗂")
no_notation Complex.imaginary_unit ("𝗂")
end

bundle no_quaternion_syntax begin
no_notation quat_ii ("𝗂")
notation Complex.imaginary_unit ("𝗂")
end

unbundle quaternion_syntax

subsection ‹Addition and Subtraction: An Abelian Group›

instantiation quat :: ab_group_add
begin

primcorec zero_quat
  where "Re 0 = 0" |"Im1 0 = 0" | "Im2 0 = 0" | "Im3 0 = 0"

primcorec plus_quat
  where
    "Re (x + y) = Re x + Re y"
  | "Im1 (x + y) = Im1 x + Im1 y"
  | "Im2 (x + y) = Im2 x + Im2 y"
  | "Im3 (x + y) = Im3 x + Im3 y"

primcorec uminus_quat
  where
    "Re (- x) = - Re x"
  | "Im1 (- x) = - Im1 x"
  | "Im2 (- x) = - Im2 x"
  | "Im3 (- x) = - Im3 x"

primcorec minus_quat
  where
    "Re (x - y) = Re x - Re y"
  | "Im1 (x - y) = Im1 x - Im1 y"
  | "Im2 (x - y) = Im2 x - Im2 y"
  | "Im3 (x - y) = Im3 x - Im3 y"

instance
  by standard (simp_all add: quat_eq_iff)

end


subsection ‹A Vector Space›

instantiation quat :: real_vector

begin

primcorec scaleR_quat
  where
    "Re (scaleR r x) = r * Re x"
  | "Im1 (scaleR r x) = r * Im1 x"
  | "Im2 (scaleR r x) = r * Im2 x"
  | "Im3 (scaleR r x) = r * Im3 x"

instance
  by standard (auto simp: quat_eq_iff distrib_left distrib_right  scaleR_add_right)

end

instantiation quat :: real_algebra_1

begin

primcorec one_quat
  where "Re 1 = 1" | "Im1 1 = 0" | "Im2 1 = 0" | "Im3 1 = 0"

primcorec times_quat
  where
    "Re (x * y) = Re x * Re y - Im1 x * Im1 y - Im2 x * Im2 y - Im3 x * Im3 y"
  | "Im1 (x * y) = Re x * Im1 y + Im1 x *  Re y + Im2 x * Im3 y - Im3 x * Im2 y"
  | "Im2 (x * y) = Re x * Im2 y - Im1 x * Im3 y + Im2 x *  Re y + Im3 x * Im1 y"
  | "Im3 (x * y) = Re x * Im3 y + Im1 x * Im2 y - Im2 x * Im1 y + Im3 x *  Re y"

instance
  by standard (auto simp: quat_eq_iff distrib_left distrib_right Rings.right_diff_distrib Rings.left_diff_distrib)

end


subsection ‹Multiplication and Division: A Real Division Algebra›

instantiation quat :: real_div_algebra
begin

primcorec inverse_quat
  where
    "Re (inverse x) = Re x / ((Re x)⇧2 + (Im1 x)⇧2 + (Im2 x)⇧2 + (Im3 x)⇧2)"
  | "Im1 (inverse x) = - (Im1 x) / ((Re x)⇧2 + (Im1 x)⇧2 + (Im2 x)⇧2 + (Im3 x)⇧2)"
  | "Im2 (inverse x) = - (Im2 x) / ((Re x)⇧2 + (Im1 x)⇧2 + (Im2 x)⇧2 + (Im3 x)⇧2)"
  | "Im3 (inverse x) = - (Im3 x) / ((Re x)⇧2 + (Im1 x)⇧2 + (Im2 x)⇧2 + (Im3 x)⇧2)"

definition "x div y = x * inverse y" for x y :: quat

instance
proof
  show "⋀x::quat. x ≠ 0 ⟹ inverse x * x = 1"
    by (auto simp: quat_eq_iff add_nonneg_eq_0_iff
        power2_eq_square add_divide_distrib [symmetric] diff_divide_distrib [symmetric])
  show "⋀x::quat. x ≠ 0 ⟹ x * inverse x = 1"
    by (auto simp: quat_eq_iff add_nonneg_eq_0_iff power2_eq_square add_divide_distrib [symmetric])
  show "⋀x y::quat. x div y = x * inverse y"
    by (simp add: divide_quat_def)
  show "inverse 0 = (0::quat)"
    by (auto simp: quat_eq_iff)
qed

end


subsection ‹Multiplication and Division: A Real Normed Division Algebra›

fun quat_proj
  where
    "quat_proj x 0 = Re x"
  | "quat_proj x (Suc 0) = Im1 x"
  | "quat_proj x (Suc (Suc 0)) = Im2 x"
  | "quat_proj x (Suc (Suc (Suc 0))) = Im3 x"

lemma quat_proj_add:
  assumes "i ≤ 3"
  shows "quat_proj (x+y) i = quat_proj x i + quat_proj y i"
proof -
  consider "i = 0" | "i = 1" | "i = 2" | "i = 3"
    using assms by linarith
  then show ?thesis
    by cases (auto simp: numeral_2_eq_2 numeral_3_eq_3)
qed

instantiation quat :: real_normed_div_algebra
begin

definition "norm z = sqrt ((Re z)⇧2 + (Im1 z)⇧2 + (Im2 z)⇧2 + (Im3 z)⇧2)"

definition "sgn x = x /⇩R norm x" for x :: quat

definition "dist x y = norm (x - y)" for x y :: quat

definition [code del]:
  "(uniformity :: (quat × quat) filter) = (INF e∈{0 <..}. principal {(x, y). dist x y < e})"

definition [code del]:
  "open (U :: quat set) ⟷ (∀x∈U. eventually (λ(x', y). x' = x ⟶ y ∈ U) uniformity)"

lemma norm_eq_L2: "norm z = L2_set (quat_proj z) {..3}"
  by (simp add: norm_quat_def L2_set_def numeral_3_eq_3)

instance
proof
  fix r :: real and x y :: quat and S :: "quat set"
  show "(norm x = 0) ⟷ (x = 0)"
    by (simp add: norm_quat_def quat_eq_iff add_nonneg_eq_0_iff)
  have eq: "L2_set (quat_proj (x + y)) {..3} = L2_set (λi. quat_proj x i + quat_proj y i) {..3}"
    by (rule L2_set_cong) (auto simp: quat_proj_add)
  show "norm (x + y) ≤ norm x + norm y"
    by (simp add: norm_eq_L2 eq L2_set_triangle_ineq)
  show "norm (scaleR r x) = ¦r¦ * norm x"
    by (simp add: norm_quat_def quat_eq_iff power_mult_distrib distrib_left [symmetric] real_sqrt_mult)
  show "norm (x * y) = norm x * norm y"
    by (simp add: norm_quat_def quat_eq_iff real_sqrt_mult [symmetric]
        power2_eq_square algebra_simps)
qed (rule sgn_quat_def dist_quat_def open_quat_def uniformity_quat_def)+

end


instantiation quat :: real_inner
begin

definition inner_quat_def:
  "inner x y = Re x * Re y + Im1 x * Im1 y + Im2 x * Im2 y + Im3 x * Im3 y"

instance
proof
  fix x y z :: quat and r :: real
  show "inner x y = inner y x"
    unfolding inner_quat_def by (simp add: mult.commute)
  show "inner (x + y) z = inner x z + inner y z"
    unfolding inner_quat_def by (simp add: distrib_right)
  show "inner (scaleR r x) y = r * inner x y"
    unfolding inner_quat_def by (simp add: distrib_left)
  show "0 ≤ inner x x"
    unfolding inner_quat_def by simp
  show "inner x x = 0 ⟷ x = 0"
    unfolding inner_quat_def by (simp add: add_nonneg_eq_0_iff quat_eq_iff)
  show "norm x = sqrt (inner x x)"
    unfolding inner_quat_def norm_quat_def
    by (simp add: power2_eq_square)
qed

end

lemma quat_inner_1 [simp]: "inner 1 x = Re x"
  unfolding inner_quat_def by simp

lemma quat_inner_1_right [simp]: "inner x 1 = Re x"
  unfolding inner_quat_def by simp

lemma quat_inner_i_left [simp]: "inner 𝗂 x = Im1 x"
  unfolding inner_quat_def by simp

lemma quat_inner_i_right [simp]: "inner x 𝗂 = Im1 x"
  unfolding inner_quat_def by simp

lemma quat_inner_j_left [simp]: "inner 𝗃 x = Im2 x"
  unfolding inner_quat_def by simp

lemma quat_inner_j_right [simp]: "inner x 𝗃 = Im2 x"
  unfolding inner_quat_def by simp

lemma quat_inner_k_left [simp]: "inner 𝗄 x = Im3 x"
  unfolding inner_quat_def by simp

lemma quat_inner_k_right [simp]: "inner x 𝗄 = Im3 x"
  unfolding inner_quat_def by simp

abbreviation quat_of_real :: "real ⇒ quat"
  where "quat_of_real ≡ of_real"

lemma Re_quat_of_real [simp]: "Re(quat_of_real a) = a"
  by (simp add: of_real_def)

lemma Im1_quat_of_real [simp]: "Im1(quat_of_real a) = 0"
  by (simp add: of_real_def)

lemma Im2_quat_of_real [simp]: "Im2(quat_of_real a) = 0"
  by (simp add: of_real_def)

lemma Im3_quat_of_real [simp]: "Im3(quat_of_real a) = 0"
  by (simp add: of_real_def)

lemma quat_eq_0_iff: "q = 0 ⟷ (Re q)⇧2 + (Im1 q)⇧2 + (Im2 q)⇧2 + (Im3 q)⇧2 = 0"
proof
  assume "(quat.Re q)⇧2 + (Im1 q)⇧2 + (Im2 q)⇧2 + (Im3 q)⇧2 = 0"
  then have "∀qa. qa - q = qa"
    by (simp add: add_nonneg_eq_0_iff minus_quat.ctr)
  then show "q = 0"
    by simp
qed auto

lemma quat_of_real_times_commute: "quat_of_real r * q = q * of_real r"
  by (simp add: of_real_def)

lemma quat_of_real_times_left_commute: "quat_of_real r * (p * q) = p * (of_real r * q)"
  by (simp add: of_real_def)

lemma quat_norm_units [simp]: "norm quat_ii = 1" "norm (𝗃::quat) = 1" "norm (𝗄::quat) = 1"
  by (auto simp: norm_quat_def)

lemma ii_nz [simp]: "quat_ii ≠ 0"
  using quat_ii.simps(2) by fastforce

lemma jj_nz [simp]: "𝗃 ≠ 0"
  using quat_jj.sel(3) by fastforce

lemma kk_nz [simp]: "𝗄 ≠ 0"
  using quat_kk.sel(4) by force

text‹An "expansion" theorem into the traditional notation›

lemma quat_unfold:
   "q = of_real(Re q) + 𝗂 * of_real(Im1 q) + 𝗃 * of_real(Im2 q) + 𝗄 * of_real(Im3 q)"
  by (simp add: quat_eq_iff)

lemma quat_trad: "Quat x y z w = of_real x + 𝗂 * of_real y + 𝗃 * of_real z + 𝗄 * of_real w"
  by (simp add: quat_eq_iff)

lemma of_real_eq_Quat: "of_real a = Quat a 0 0 0"
  by (simp add: quat_trad)

lemma ii_squared [simp]: "quat_ii⇧2 = -1"
  by (simp add: power2_eq_square quat.expand)

lemma jj_squared [simp]: "𝗃^2 = -1"
  by (simp add: power2_eq_square quat.expand)

lemma kk_squared [simp]: "𝗄^2 = -1"
  by (simp add: power2_eq_square quat.expand)

lemma inverse_ii [simp]: "inverse quat_ii = -quat_ii"
  by (simp add: power2_eq_square quat.expand)

lemma inverse_jj [simp]: "inverse 𝗃 = -𝗃"
  by (simp add: power2_eq_square quat.expand)

lemma inverse_kk [simp]: "inverse 𝗄 = -𝗄"
  by (simp add: power2_eq_square quat.expand)

lemma inverse_mult: "inverse (p * q) = inverse q * inverse p" for p::quat
  by (metis inverse_zero mult_not_zero nonzero_inverse_mult_distrib)

lemma quat_of_real_inverse_collapse [simp]:
  assumes "c ≠ 0"
  shows "quat_of_real c * quat_of_real (inverse c) = 1" "quat_of_real (inverse c) * quat_of_real c = 1"
  using assms by auto

subsection‹Conjugate of a quaternion›

primcorec cnj :: "quat ⇒ quat"
  where
    "Re (cnj z) = Re z"
  | "Im1 (cnj z) = - Im1 z"
  | "Im2 (cnj z) = - Im2 z"
  | "Im3 (cnj z) = - Im3 z"


lemma cnj_cancel_iff [simp]: "cnj x = cnj y ⟷ x = y"
proof
  show "cnj x = cnj y ⟹ x = y"
    by (simp add: quat_eq_iff)
qed auto

lemma cnj_cnj [simp]:
   "cnj(cnj q) = q"
  by (simp add: quat_eq_iff)

lemma cnj_of_real [simp]: "cnj(quat_of_real x) = quat_of_real x"
  by (simp add: quat_eqI)

lemma cnj_zero [simp]: "cnj 0 = 0"
  by (simp add: quat_eq_iff)

lemma cnj_zero_iff [iff]: "cnj z = 0 ⟷ z = 0"
  by (simp add: quat_eq_iff)

lemma cnj_one [simp]: "cnj 1 = 1"
  by (simp add: quat_eq_iff)

lemma cnj_one_iff [simp]: "cnj z = 1 ⟷ z = 1"
  by (simp add: quat_eq_iff)

lemma quat_norm_cnj [simp]: "norm(cnj q) = norm q"
  by (simp add: norm_quat_def)

lemma cnj_add [simp]: "cnj (x + y) = cnj x + cnj y"
  by (simp add: quat_eq_iff)

lemma cnj_sum [simp]: "cnj (sum f S) = (∑x∈S. cnj (f x))"
  by (induct S rule: infinite_finite_induct) auto

lemma cnj_diff [simp]: "cnj (x - y) = cnj x - cnj y"
  by (simp add: quat_eq_iff)

lemma cnj_minus [simp]: "cnj (- x) = - cnj x"
  by (simp add: quat_eq_iff)

lemma cnj_mult [simp]: "cnj (x * y) = cnj y * cnj x"
  by (simp add: quat_eq_iff)

lemma cnj_inverse [simp]: "cnj (inverse x) = inverse (cnj x)"
  by (simp add: quat_eq_iff)

lemma cnj_divide [simp]: "cnj (x / y) = inverse (cnj y) * cnj x"
  by (simp add: divide_quat_def)

lemma cnj_power [simp]: "cnj (x^n) = cnj x ^ n"
  by (induct n) (auto simp: power_commutes)

lemma cnj_of_nat [simp]: "cnj (of_nat n) = of_nat n"
  by (metis cnj_of_real of_real_of_nat_eq)

lemma cnj_of_int [simp]: "cnj (of_int z) = of_int z"
  by (metis cnj_of_real of_real_of_int_eq)

lemma cnj_numeral [simp]: "cnj (numeral w) = numeral w"
  by (metis of_nat_numeral cnj_of_nat)

lemma cnj_neg_numeral [simp]: "cnj (- numeral w) = - numeral w"
  by simp

lemma cnj_scaleR [simp]: "cnj (scaleR r x) = scaleR r (cnj x)"
  by (simp add: quat_eq_iff)

lemma cnj_units [simp]: "cnj quat_ii = -𝗂" "cnj 𝗃 = -𝗃" "cnj 𝗄 = -𝗄"
  by (simp_all add: quat_eq_iff)

lemma cnj_eq_of_real: "cnj q = quat_of_real x ⟷ q = quat_of_real x"
proof
  show "cnj q = quat_of_real x ⟹ q = quat_of_real x"
    by (metis cnj_of_real cnj_cnj)
qed auto

lemma quat_add_cnj: "q + cnj q = quat_of_real(2 * Re q)" "cnj q + q = quat_of_real(2 * Re q)"
  by simp_all (simp_all add: mult.commute mult_2 plus_quat.code)

lemma quat_divide_numeral:
  fixes x::quat shows "x / numeral w = x /⇩R numeral w"
  unfolding divide_quat_def
  by (metis mult.right_neutral mult_scaleR_right of_real_def of_real_inverse of_real_numeral)

lemma Re_divide_numeral [simp]: "Re (x / numeral w) = Re x / numeral w"
  by (metis divide_inverse_commute quat_divide_numeral scaleR_quat.simps(1))

lemma Im1_divide_numeral [simp]: "Im1 (x / numeral w) = Im1 x / numeral w"
  unfolding quat_divide_numeral by simp

lemma Im2_divide_numeral [simp]: "Im2 (x / numeral w) = Im2 x / numeral w"
  unfolding quat_divide_numeral by simp

lemma Im3_divide_numeral [simp]: "Im3 (x / numeral w) = Im3 x / numeral w"
  unfolding quat_divide_numeral by simp


lemma of_real_quat_re_cnj: "quat_of_real(Re q) = inverse(quat_of_real 2) * (q + cnj q)"
  by (simp add: quat_eq_iff)

lemma quat_mult_cnj_commute: "cnj q * q = q * cnj q"
  by (simp add: quat_eq_iff)

lemma quat_norm_pow_2: "quat_of_real(norm q) ^ 2 = q * cnj q"
  by (simp add: quat_eq_iff norm_quat_def flip: of_real_power) (auto simp: power2_eq_square)

lemma quat_norm_pow_2_alt: "quat_of_real(norm q) ^ 2 = cnj q * q"
  by (simp add: quat_mult_cnj_commute quat_norm_pow_2)

lemma quat_inverse_cnj: "inverse q = inverse (quat_of_real((norm q)⇧2)) * cnj q"
  by (simp add: quat_eq_iff norm_quat_def numeral_Bit0 flip: of_real_power)

lemma quat_inverse_eq_cnj: "norm q = 1 ⟹ inverse q = cnj q"
  by (metis inverse_1 mult_cancel_left norm_eq_zero norm_one cnj_one quat_norm_pow_2 right_inverse)

subsection‹Linearity and continuity of the components›

lemma bounded_linear_Re: "bounded_linear Re"
  and bounded_linear_Im1: "bounded_linear Im1"
  and bounded_linear_Im2: "bounded_linear Im2"
  and bounded_linear_Im3: "bounded_linear Im3"
  by (simp_all add: bounded_linear_intro [where K=1] norm_quat_def real_le_rsqrt add.assoc)

lemmas Cauchy_Re = bounded_linear.Cauchy [OF bounded_linear_Re]
lemmas Cauchy_Im1 = bounded_linear.Cauchy [OF bounded_linear_Im1]
lemmas Cauchy_Im2 = bounded_linear.Cauchy [OF bounded_linear_Im2]
lemmas Cauchy_Im3 = bounded_linear.Cauchy [OF bounded_linear_Im3]
lemmas tendsto_Re [tendsto_intros] = bounded_linear.tendsto [OF bounded_linear_Re]
lemmas tendsto_Im1 [tendsto_intros] = bounded_linear.tendsto [OF bounded_linear_Im1]
lemmas tendsto_Im2 [tendsto_intros] = bounded_linear.tendsto [OF bounded_linear_Im2]
lemmas tendsto_Im3 [tendsto_intros] = bounded_linear.tendsto [OF bounded_linear_Im3]
lemmas isCont_Re [simp] = bounded_linear.isCont [OF bounded_linear_Re]
lemmas isCont_Im1 [simp] = bounded_linear.isCont [OF bounded_linear_Im1]
lemmas isCont_Im2 [simp] = bounded_linear.isCont [OF bounded_linear_Im2]
lemmas isCont_Im3 [simp] = bounded_linear.isCont [OF bounded_linear_Im3]
lemmas continuous_Re [simp] = bounded_linear.continuous [OF bounded_linear_Re]
lemmas continuous_Im1 [simp] = bounded_linear.continuous [OF bounded_linear_Im1]
lemmas continuous_Im2 [simp] = bounded_linear.continuous [OF bounded_linear_Im2]
lemmas continuous_Im3 [simp] = bounded_linear.continuous [OF bounded_linear_Im3]
lemmas continuous_on_Re [continuous_intros] = bounded_linear.continuous_on[OF bounded_linear_Re]
lemmas continuous_on_Im1 [continuous_intros] = bounded_linear.continuous_on[OF bounded_linear_Im1]
lemmas continuous_on_Im2 [continuous_intros] = bounded_linear.continuous_on[OF bounded_linear_Im2]
lemmas continuous_on_Im3 [continuous_intros] = bounded_linear.continuous_on[OF bounded_linear_Im3]
lemmas has_derivative_Re [derivative_intros] = bounded_linear.has_derivative[OF bounded_linear_Re]
lemmas has_derivative_Im1 [derivative_intros] = bounded_linear.has_derivative[OF bounded_linear_Im1]
lemmas has_derivative_Im2 [derivative_intros] = bounded_linear.has_derivative[OF bounded_linear_Im2]
lemmas has_derivative_Im3 [derivative_intros] = bounded_linear.has_derivative[OF bounded_linear_Im3]
lemmas sums_Re = bounded_linear.sums [OF bounded_linear_Re]
lemmas sums_Im1 = bounded_linear.sums [OF bounded_linear_Im1]
lemmas sums_Im2 = bounded_linear.sums [OF bounded_linear_Im2]
lemmas sums_Im3 = bounded_linear.sums [OF bounded_linear_Im3]


subsection‹Quaternionic-specific theorems about sums›

lemma Re_sum [simp]: "Re(sum f S) = sum (λx.  Re(f x)) S" for f :: "'a ⇒ quat"
  by (induct S rule: infinite_finite_induct) auto

lemma Im1_sum [simp]: "Im1(sum f S) = sum (λx. Im1(f x)) S"
  by (induct S rule: infinite_finite_induct) auto

lemma Im2_sum [simp]: "Im2(sum f S) = sum (λx. Im2(f x)) S"
  by (induct S rule: infinite_finite_induct) auto

lemma Im3_sum [simp]: "Im3(sum f S) = sum (λx. Im3(f x)) S"
  by (induct S rule: infinite_finite_induct) auto

lemma in_Reals_iff_Re: "q ∈ Reals ⟷ quat_of_real(Re q) = q"
  by (metis Re_quat_of_real Reals_cases Reals_of_real)

lemma in_Reals_iff_cnj: "q ∈ Reals ⟷ cnj q = q"
  by (simp add: in_Reals_iff_Re quat_eq_iff)

lemma real_norm: "q ∈ Reals ⟹ norm q = abs(Re q)"
  by (metis in_Reals_iff_Re norm_of_real)


lemma norm_power2: "(norm q)⇧2 = Re(cnj q * q)"
  by (metis Re_quat_of_real of_real_power quat_mult_cnj_commute quat_norm_pow_2)

lemma quat_norm_imaginary: "Re x = 0 ⟹ x⇧2 = - (quat_of_real (norm x))⇧2"
  unfolding quat_norm_pow_2
  by (cases x) (auto simp: power2_eq_square cnj.ctr times_quat.ctr uminus_quat.ctr)


subsection‹Bound results for real and imaginary components of limits›

lemma Re_tendsto_upperbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. quat.Re (f x) ≤ b; net ≠ bot⟧ ⟹ Re l ≤ b"
  by (blast intro: tendsto_upperbound [OF tendsto_Re])

lemma Im1_tendsto_upperbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. Im1 (f x) ≤ b; net ≠ bot⟧ ⟹ Im1 l ≤ b"
  by (blast intro: tendsto_upperbound [OF tendsto_Im1])

lemma Im2_tendsto_upperbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. Im2 (f x) ≤ b; net ≠ bot⟧ ⟹ Im2 l ≤ b"
  by (blast intro: tendsto_upperbound [OF tendsto_Im2])

lemma Im3_tendsto_upperbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. Im3 (f x) ≤ b; net ≠ bot⟧ ⟹ Im3 l ≤ b"
  by (blast intro: tendsto_upperbound [OF tendsto_Im3])

lemma Re_tendsto_lowerbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. b ≤ quat.Re (f x); net ≠ bot⟧ ⟹ b ≤ Re l"
  by (blast intro: tendsto_lowerbound [OF tendsto_Re])

lemma Im1_tendsto_lowerbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. b ≤ Im1 (f x); net ≠ bot⟧ ⟹ b ≤ Im1 l"
  by (blast intro: tendsto_lowerbound [OF tendsto_Im1])

lemma Im2_tendsto_lowerbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. b ≤ Im2 (f x); net ≠ bot⟧ ⟹ b ≤ Im2 l"
  by (blast intro: tendsto_lowerbound [OF tendsto_Im2])

lemma Im3_tendsto_lowerbound:
   "⟦(f ⤏ l) net; ∀⇩F x in net. b ≤ Im3 (f x); net ≠ bot⟧ ⟹ b ≤ Im3 l"
  by (blast intro: tendsto_lowerbound [OF tendsto_Im3])

lemma of_real_continuous_iff: "continuous net (λx. quat_of_real(f x)) ⟷ continuous net f"
  by (simp add: continuous_def tendsto_iff)

lemma of_real_continuous_on_iff:
   "continuous_on S (λx. quat_of_real(f x)) ⟷ continuous_on S f"
  using continuous_on_Re continuous_on_of_real by fastforce

subsection‹Quaternions for describing 3D isometries›

subsubsection‹The ‹HIm› operator›

definition HIm :: "quat ⇒ real^3" where
  "HIm q ≡ vector[Im1 q, Im2 q, Im3 q]"

lemma HIm_Quat: "HIm (Quat w x y z) = vector[x,y,z]"
  by (simp add: HIm_def)

lemma him_eq: "HIm p = HIm q ⟷ Im1 p = Im1 q ∧ Im2 p = Im2 q ∧ Im3 p = Im3 q"
  by (metis HIm_def vector_3)

lemma him_of_real [simp]: "HIm(of_real a) = 0"
  by (simp add: of_real_eq_Quat HIm_Quat vec_eq_iff vector_def)

lemma him_0 [simp]: "HIm 0 = 0"
  by (metis him_of_real of_real_0)

lemma him_1 [simp]: "HIm 1 = 0"
  by (metis him_of_real of_real_1)

lemma him_cnj: "HIm(cnj q) = - HIm q"
  by (simp add: HIm_def vec_eq_iff vector_def)

lemma him_mult_left [simp]: "HIm(of_real a * q) = a *⇩R HIm q"
  by (simp add: HIm_def vec_eq_iff vector_def)

lemma him_mult_right [simp]: "HIm(q * of_real a) = a *⇩R HIm q"
  by (simp add: HIm_def vec_eq_iff vector_def)

lemma him_add [simp]: "HIm(p + q) = HIm p + HIm q"
  by (simp add: HIm_def vec_eq_iff vector_def)

lemma him_minus [simp]: "HIm(-q) = - HIm q"
  by (simp add: HIm_def vec_eq_iff vector_def)

lemma him_diff [simp]: "HIm(p - q) = HIm p - HIm q"
  by (simp add: HIm_def vec_eq_iff vector_def)

lemma him_sum [simp]: "HIm (sum f S) = (∑x∈S. HIm (f x))"
  by (induct S rule: infinite_finite_induct) auto

lemma linear_him: "linear HIm"
  by (metis him_add him_mult_right linearI mult.right_neutral mult_scaleR_right of_real_def)

subsubsection‹The ‹Hv› operator›

definition Hv :: "real^3 ⇒ quat" where
  "Hv v ≡ Quat 0 (v$1) (v$2) (v$3)"

lemma Re_Hv [simp]: "Re(Hv v) = 0"
  by (simp add: Hv_def)

lemma Im1_Hv [simp]: "Im1(Hv v) = v$1"
  by (simp add: Hv_def)

lemma Im2_Hv [simp]: "Im2(Hv v) = v$2"
  by (simp add: Hv_def)

lemma Im3_Hv [simp]: "Im3(Hv v) = v$3"
  by (simp add: Hv_def)

lemma hv_vec: "Hv(vec r) = Quat 0 r r r"
  by (simp add: Hv_def)

lemma hv_eq_zero [simp]: "Hv v = 0 ⟷ v = 0"
  by (simp add: quat_eq_iff vec_eq_iff) (metis exhaust_3)

lemma hv_zero [simp]: "Hv 0 = 0"
  by simp

lemma hv_vector [simp]: "Hv(vector[x,y,z]) = Quat 0 x y z"
  by (simp add: Hv_def)

lemma hv_basis: "Hv(axis 1 1) = 𝗂" "Hv(axis 2 1) = 𝗃" "Hv(axis 3 1) = 𝗄"
  by (auto simp: Hv_def axis_def quat_ii.code quat_jj.code quat_kk.code)

lemma hv_add [simp]: "Hv(x + y) = Hv x + Hv y"
  by (simp add: Hv_def quat_eq_iff)

lemma hv_minus [simp]: "Hv(-x) = -Hv x"
  by (simp add: Hv_def quat_eq_iff)

lemma hv_diff [simp]: "Hv(x - y) = Hv x - Hv y"
  by (simp add: Hv_def quat_eq_iff)

lemma hv_cmult [simp]: "Hv(a *⇩R x) = of_real a * Hv x"
  by (simp add: Hv_def quat_eq_iff)

lemma hv_sum [simp]: "Hv (sum f S) = (∑x∈S. Hv (f x))"
  by (induct S rule: infinite_finite_induct) auto

lemma hv_inj: "Hv x = Hv y ⟷ x = y"
  by (simp add: Hv_def quat_eq_iff vec_eq_iff) (metis (full_types) exhaust_3)

lemma linear_hv: "linear Hv"
  by unfold_locales (auto simp: of_real_def)

lemma him_hv [simp]: "HIm(Hv x) = x"
  using HIm_def hv_inj quat_eq_iff by fastforce

lemma cnj_hv [simp]: "cnj(Hv v) = -Hv v"
  using Hv_def cnj.code hv_minus by auto

lemma hv_him: "Hv(HIm q) = Quat 0 (Im1 q) (Im2 q) (Im3 q)"
  by (simp add: HIm_def)

lemma hv_him_eq: "Hv(HIm q) = q ⟷ Re q = 0"
  by (simp add: hv_him quat_eq_iff)

lemma dot_hv [simp]: "Hv u ∙ Hv v = u ∙ v"
  by (simp add: Hv_def inner_quat_def inner_vec_def sum_3)

lemma norm_hv [simp]: "norm (Hv v) = norm v"
  by (simp add: norm_eq)

subsection‹Geometric interpretation of the product of imaginary quaternions›

context includes cross3_syntax
begin

lemma mult_hv_eq_cross_dot: "Hv x * Hv y = Hv(x × y) - of_real (x ∙ y)"
  by (simp add: quat_eq_iff cross3_simps)

text‹Representing orthogonal transformations as conjugation or congruence with a quaternion›

lemma orthogonal_transformation_quat_congruence:
  assumes "norm q = 1"
  shows "orthogonal_transformation (λx. HIm(cnj q * Hv x * q))"
proof -
  have nq: "(quat.Re q)⇧2 + (Im1 q)⇧2 + (Im2 q)⇧2 + (Im3 q)⇧2 = 1"
    using assms norm_quat_def by auto
  have "Vector_Spaces.linear (*⇩R) (*⇩R) (λx. HIm (cnj q * Hv x * q))"
  proof
    show "⋀r b. HIm (cnj q * Hv (r *⇩R b) * q) = r *⇩R HIm (cnj q * Hv b * q)"
      by (metis him_mult_left hv_cmult mult_scaleR_left mult_scaleR_right scaleR_conv_of_real)
  qed (simp add: distrib_left distrib_right)
  moreover have "HIm (cnj q * Hv v * q) ∙ HIm (cnj q * Hv w * q) = ((quat.Re q)⇧2 + (Im1 q)⇧2 + (Im2 q)⇧2 + (Im3 q)⇧2)⇧2 * (v ∙ w)" for v w
    by (simp add: HIm_def inner_vec_def sum_3 power2_eq_square algebra_simps) (*SLOW, like 15s*)
  ultimately show ?thesis
    by (simp add: orthogonal_transformation_def linear_def nq)
qed

lemma orthogonal_transformation_quat_conjugation:
  assumes "q ≠ 0"
  shows "orthogonal_transformation (λx. HIm(inverse q * Hv x * q))"
proof -
  obtain c p where eq: "q = of_real c * p" and 1: "norm p = 1"
  proof
    show "q = quat_of_real (norm q) * (inverse (of_real (norm q)) * q)"
      by (metis assms mult.assoc mult.left_neutral norm_eq_zero of_real_eq_0_iff right_inverse)
    show "norm (inverse (quat_of_real (norm q)) * q) = 1"
      using assms by (simp add: norm_mult norm_inverse)
  qed
  have "c ≠ 0"
    using assms eq by auto
  then have "HIm (cnj p * Hv x * p) = HIm (inverse (quat_of_real c * p) * Hv x * (quat_of_real c * p))" for x
    apply (simp add: inverse_mult mult.assoc flip: quat_inverse_eq_cnj [OF 1] of_real_inverse)
    using quat_of_real_times_commute quat_of_real_times_left_commute quat_of_real_inverse_collapse
    by (simp add: of_real_def)
  then show ?thesis
    using orthogonal_transformation_quat_congruence [OF 1]
    by (simp add: eq)
qed

unbundle no_quaternion_syntax

end

end