Theory Quadratic_Forms

(*  Title:      Three_Squares/Quadratic_Forms.thy
    Author:     Anton Danilkin
*)

section ‹Properties of quadratic forms and their equivalences›

theory Quadratic_Forms
  imports Complex_Main Low_Dimensional_Linear_Algebra
begin

consts qf_app :: "'a ⇒ 'b ⇒ int" (infixl ‹$$› 65)

definition qf2_app :: "mat2 ⇒ vec2 ⇒ int" where
"qf2_app m v = <v | m $ v>"

adhoc_overloading qf_app qf2_app

definition qf3_app :: "mat3 ⇒ vec3 ⇒ int" where
"qf3_app m v = <v | m $ v>"

adhoc_overloading qf_app qf3_app

lemma qf2_app_zero [simp]:
  fixes m :: mat2
  shows "m $$ 0 = 0"
  unfolding qf2_app_def by auto

lemma qf3_app_zero [simp]:
  fixes m :: mat3
  shows "m $$ 0 = 0"
  unfolding qf3_app_def by auto

consts qf_positive_definite :: "'a ⇒ bool"

definition qf2_positive_definite :: "mat2 ⇒ bool" where
"qf2_positive_definite m = (∀v. v ≠ 0 ⟶ m $$ v > 0)"

adhoc_overloading qf_positive_definite qf2_positive_definite

definition qf3_positive_definite :: "mat3 ⇒ bool" where
"qf3_positive_definite m = (∀v. v ≠ 0 ⟶ m $$ v > 0)"

adhoc_overloading qf_positive_definite qf3_positive_definite

lemma qf2_positive_definite_positive:
  fixes m :: mat2
  assumes "qf_positive_definite m"
  shows "∀v. m $$ v ≥ 0"
  using assms unfolding qf2_positive_definite_def
  by (metis order_less_le order_refl qf2_app_zero)

lemma qf3_positive_definite_positive:
  fixes m :: mat3
  assumes "qf_positive_definite m"
  shows "∀v. m $$ v ≥ 0"
  using assms unfolding qf3_positive_definite_def
  by (metis order_less_le order_refl qf3_app_zero)

consts qf_action :: "'a ⇒ 'a ⇒ 'a" (infixl ‹⋅› 55)

definition qf2_action :: "mat2 ⇒ mat2 ⇒ mat2" where
"qf2_action a u = uT * a * u"

adhoc_overloading qf_action qf2_action

definition qf3_action :: "mat3 ⇒ mat3 ⇒ mat3" where
"qf3_action a u = uT * a * u"

adhoc_overloading qf_action qf3_action

lemma qf2_action_id:
  fixes a :: mat2
  shows "a ⋅ 1 = a"
  unfolding qf2_action_def
  by simp

lemma qf3_action_id:
  fixes a :: mat3
  shows "a ⋅ 1 = a"
  unfolding qf3_action_def
  by simp

lemma qf2_action_mul [simp]:
  fixes a u v :: mat2
  shows "a ⋅ (u * v) = (a ⋅ u) ⋅ v"
  unfolding qf2_action_def
  by (simp add: algebra_simps)

lemma qf3_action_mul [simp]:
  fixes a u v :: mat3
  shows "a ⋅ (u * v) = (a ⋅ u) ⋅ v"
  unfolding qf3_action_def
  by (simp add: algebra_simps)

consts qf_equiv :: "'a ⇒ 'a ⇒ bool" (infix ‹~› 65)

definition qf2_equiv :: "mat2 ⇒ mat2 ⇒ bool" where
"qf2_equiv a b = (∃u. mat_det u = 1 ∧ a ⋅ u = b)"

adhoc_overloading qf_equiv qf2_equiv

definition qf3_equiv :: "mat3 ⇒ mat3 ⇒ bool" where
"qf3_equiv a b = (∃u. mat_det u = 1 ∧ a ⋅ u = b)"

adhoc_overloading qf_equiv qf3_equiv

lemma qf2_equiv_sym_impl:
  fixes a b :: mat2
  shows "a ~ b ⟹ b ~ a"
unfolding qf2_equiv_def qf2_action_def
proof -
  assume "∃u. mat_det u = 1 ∧ uT * a * u = b"
  then obtain u where "mat_det u = 1 ∧ uT * a * u = b" by blast
  hence "mat_det (u-1) = 1 ∧ ((u-1)T) * b * (u-1) = a"
    unfolding mat2_inverse_transpose[symmetric]
    using mat2_inverse_cancel_left[of "uT"] mat2_inverse_cancel_right
    by (auto simp add: algebra_simps)
  thus "∃u. mat_det u = 1 ∧ uT * b * u = a" by blast
qed

lemma qf3_equiv_sym_impl:
  fixes a b :: mat3
  shows "a ~ b ⟹ b ~ a"
unfolding qf3_equiv_def qf3_action_def
proof -
  assume "∃u. mat_det u = 1 ∧ uT * a * u = b"
  then obtain u where "mat_det u = 1 ∧ uT * a * u = b" by blast
  hence "mat_det (u-1) = 1 ∧ ((u-1)T) * b * (u-1) = a"
    unfolding mat3_inverse_transpose[symmetric]
    using mat3_inverse_cancel_left[of "uT"] mat3_inverse_cancel_right
    by (auto simp add: algebra_simps)
  thus "∃u. mat_det u = 1 ∧ uT * b * u = a" by blast
qed

lemma qf2_equiv_sym:
  fixes a b :: mat2
  shows "a ~ b ⟷ b ~ a"
  using qf2_equiv_sym_impl by blast

lemma qf3_equiv_sym:
  fixes a b :: mat3
  shows "a ~ b ⟷ b ~ a"
  using qf3_equiv_sym_impl by blast

lemma qf2_equiv_trans:
  fixes a b c :: mat2
  assumes "a ~ b"
  assumes "b ~ c"
  shows "a ~ c"
  using assms by (metis mult_1 mat2_mul_det qf2_action_mul qf2_equiv_def)

lemma qf3_equiv_trans:
  fixes a b c :: mat3
  assumes "a ~ b"
  assumes "b ~ c"
  shows "a ~ c"
  using assms by (metis mult_1 mat3_mul_det qf3_action_mul qf3_equiv_def)

lemma qf2_action_app [simp]:
  fixes a u :: mat2
  fixes v :: vec2
  shows "(a ⋅ u) $$ v = a $$ (u $ v)"
  unfolding qf2_action_def qf2_app_def
  using vec2_dot_transpose_right by auto

lemma qf3_action_app [simp]:
  fixes a u :: mat3
  fixes v :: vec3
  shows "(a ⋅ u) $$ v = a $$ (u $ v)"
  unfolding qf3_action_def qf3_app_def
  using vec3_dot_transpose_right by auto

lemma qf2_equiv_preserves_positive_definite:
  fixes a b :: mat2
  assumes "a ~ b"
  shows "qf_positive_definite a ⟷ qf_positive_definite b"
  unfolding qf2_positive_definite_def
  by (metis assms mat2_special_preserves_zero qf2_action_app
            qf2_equiv_def qf2_equiv_sym)

lemma qf3_equiv_preserves_positive_definite:
  fixes a b :: mat3
  assumes "a ~ b"
  shows "qf_positive_definite a ⟷ qf_positive_definite b"
  unfolding qf3_positive_definite_def
  by (metis assms mat3_special_preserves_zero qf3_action_app
            qf3_equiv_def qf3_equiv_sym)

lemma qf2_equiv_preserves_sym:
  fixes a b :: mat2
  assumes "a ~ b"
  shows "mat2_sym a ⟷ mat2_sym b"
proof -
  obtain u where "mat_det u = 1" "uT * a * u = b"
    using assms unfolding qf2_action_def qf2_equiv_def by auto
  thus ?thesis
    unfolding mat2_sym_criterion
    using mat2_inversable_cancel_left[of "uT" "aT * u" "a * u"]
          mat2_inversable_cancel_right[of u "aT" a]
    by (auto simp add: algebra_simps)
qed

lemma qf3_equiv_preserves_sym:
  fixes a b :: mat3
  assumes "a ~ b"
  shows "mat3_sym a ⟷ mat3_sym b"
proof -
  obtain u where "mat_det u = 1" "uT * a * u = b"
    using assms unfolding qf3_action_def qf3_equiv_def by auto
  thus ?thesis
    unfolding mat3_sym_criterion
    using mat3_inversable_cancel_left[of "uT" "aT * u" "a * u"]
          mat3_inversable_cancel_right[of u "aT" a]
    by (auto simp add: algebra_simps)
qed

lemma qf2_equiv_preserves_det:
  fixes a b :: mat2
  assumes "a ~ b"
  shows "mat_det a = mat_det b"
  using assms unfolding qf2_action_def qf2_equiv_def by auto

lemma qf3_equiv_preserves_det:
  fixes a b :: mat3
  assumes "a ~ b"
  shows "mat_det a = mat_det b"
  using assms unfolding qf3_action_def qf3_equiv_def by auto

lemma qf2_equiv_preserves_range_subset:
  fixes a b :: mat2
  assumes "a ~ b"
  shows "range (($$) b) ⊆ range (($$) a)"
proof -
  obtain u where 0: "mat_det u = 1" "a ⋅ u = b"
    using assms unfolding qf2_equiv_def by auto
  show ?thesis unfolding 0[symmetric] image_def by auto
qed

lemma qf3_equiv_preserves_range_subset:
  fixes a b :: mat3
  assumes "a ~ b"
  shows "range (($$) b) ⊆ range (($$) a)"
proof -
  obtain u where 0: "mat_det u = 1" "a ⋅ u = b"
    using assms unfolding qf3_equiv_def by auto
  show ?thesis unfolding 0[symmetric] image_def by auto
qed

lemma qf2_equiv_preserves_range:
  fixes a b :: mat2
  assumes "a ~ b"
  shows "range (($$) a) = range (($$) b)"
  using assms qf2_equiv_sym qf2_equiv_preserves_range_subset by blast

lemma qf3_equiv_preserves_range:
  fixes a b :: mat3
  assumes "a ~ b"
  shows "range (($$) a) = range (($$) b)"
  using assms qf3_equiv_sym qf3_equiv_preserves_range_subset by blast

text ‹Lemma 1.1 from @{cite nathanson1996}.›

lemma qf2_positive_definite_criterion:
  fixes a
  assumes "mat_sym a"
  shows "qf_positive_definite a ⟷ mat211 a > 0 ∧ mat_det a > 0"
proof
  assume 0: "qf_positive_definite a"
  have "vec2 1 0 ≠ 0 ⟶ a $$ vec2 1 0 > 0" using 0
    unfolding qf2_positive_definite_def by blast
  hence 1: "mat211 a > 0"
    unfolding zero_vec2_def vec2_dot_def times_mat2_def mat2_app_def qf2_app_def
    by auto
  let ?v = "vec2 (- mat212 a) (mat211 a)"
  have "?v ≠ 0 ⟶ a $$ ?v > 0" using 0
    unfolding qf2_positive_definite_def by blast
  hence "mat211 a * mat_det a > 0" using assms 1
    unfolding zero_vec2_def vec2_dot_def times_mat2_def mat2_app_def
              mat2_det_def mat2_sym_def qf2_app_def
    by (simp add: mult.commute)
  hence 2: "mat_det a > 0" using 1 zero_less_mult_pos by blast
  show "mat211 a > 0 ∧ mat_det a > 0" using 1 2 by blast
next
  assume 3: "mat211 a > 0 ∧ mat_det a > 0"
  show "qf_positive_definite a" unfolding qf2_positive_definite_def
  proof (rule allI; rule impI)
    fix v :: vec2
    assume "v ≠ 0"
    hence 4: "vec21 v ≠ 0 ∨ vec22 v ≠ 0" unfolding zero_vec2_def
      by (metis vec2.collapse)
    let ?n = "(mat211 a * vec21 v + mat212 a * vec22 v)2 + (mat_det a) * (vec22 v)2"
    have 5: "mat211 a * (a $$ v) = ?n"
      using assms
      unfolding vec2_dot_def times_mat2_def mat2_app_def mat2_det_def
                mat2_sym_def qf2_app_def power2_eq_square
      by (simp add: algebra_simps)
    have "?n > 0" using 3 4
      by (metis add.commute add_0 add_pos_nonneg mult_eq_0_iff
          mult_pos_pos power_zero_numeral zero_le_power2 zero_less_power2)
    thus "a $$ v > 0" using 3 5 zero_less_mult_pos by metis
  qed
qed

lemma congruence_class_close:
  fixes k m :: int
  assumes "m > 0"
  shows "∃t. 2 * ¦k + m * t¦ ≤ m" (is "∃ t. ?P t")
proof -
  let ?s = "k div m"
  have 0: "k - m * ?s ≥ 0 ∧ k - m * ?s < m" using assms
    by (metis pos_mod_sign pos_mod_bound add.commute add_diff_cancel_right'
              div_mod_decomp_int mult.commute)
  show ?thesis proof cases
    assume "2 * (k - m * ?s) ≤ m"
    hence "?P (- ?s)" using 0 by auto
    thus ?thesis by blast
  next
    assume "¬ (2 * (k - m * ?s) ≤ m)"
    hence "2 * (k - m * ?s) > m" by auto
    hence "?P (- (?s + 1))" using 0 by (simp add: algebra_simps)
    thus ?thesis by blast
  qed
qed

text ‹Lemma 1.2 from @{cite nathanson1996}.›

lemma lemma_1_2:
  fixes b :: mat2
  assumes "mat_sym b"
  assumes "qf_positive_definite b"
  shows "∃a. a ~ b ∧
             2 * ¦mat212 a¦ ≤ mat211 a ∧
             mat211 a ≤ (2 / sqrt 3) * sqrt (mat_det a)" (is "∃a. ?P a")
proof -
  define a11 where "a11 ≡ LEAST y. y > 0 ∧ (∃x. int y = b $$ x)"
  have 0: "∃y. y > 0 ∧ (∃x. int y = b $$ x)"
    apply (rule exI[of _ "nat (b $$ (vec2 1 1))"])
    using assms(2) unfolding qf2_positive_definite_def zero_vec2_def
    apply (metis nat_0_le order_less_le vec2.inject zero_less_nat_eq zero_neq_one)
    done
  obtain r where 1: "a11 > 0" "int a11 = b $$ r"
    using a11_def LeastI_ex[OF 0] by auto
  let ?h = "gcd (vec21 r) (vec22 r)"
  have "r ≠ 0" using assms(2) 1 by fastforce
  hence 2: "?h > 0" by (simp; metis vec2.collapse zero_vec2_def)
  let ?r' = "vec2 (vec21 r div ?h) (vec22 r div ?h)"
  have "?r' ≠ 0" using 2 unfolding zero_vec2_def
    by (simp add: algebra_simps dvd_div_eq_0_iff)
  hence "nat (b $$ ?r') > 0 ∧ (∃x. int (nat (b $$ ?r')) = b $$ x)"
    using assms(2) unfolding qf2_positive_definite_def by auto
  hence *: "a11 ≤ nat (b $$ ?r')" unfolding a11_def by (rule Least_le)
  hence "a11 ≤ b $$ ?r'" using 1 by auto
  also have **: "... = (b $$ r) div ?h2" proof -
    have "(b $$ ?r') * ?h2 = b $$ r" "?h2 dvd b $$ r"
      unfolding qf2_app_def vec2_dot_def mat2_app_def power2_eq_square
      using 1
      by (auto simp add: algebra_simps mult_dvd_mono)
    thus "b $$ ?r' = (b $$ r) div ?h2" using 2 by auto
  qed
  also have "... = a11 div ?h2" using 1 by auto
  finally have "a11 ≤ a11 div ?h2" .
  also have "... ≤ a11" using 1 2
    by (metis * ** assms(2) div_le_dividend int_eq_iff nat_div_distrib order_refl
        qf2_positive_definite_positive verit_la_disequality)
  finally have "?h = 1" using 1 2
    by (smt (verit) int_div_less_self of_nat_0_less_iff power2_eq_square
                    zero_less_power zmult_eq_1_iff)
  then obtain s1 s2 where 3: "1 = (vec21 r) * s2 - (vec22 r) * s1"
    by (metis bezout_int diff_minus_eq_add mult.commute mult_minus_right)
  define a'12 where "
    a'12 ≡
      (mat211 b) * (vec21 r) * s1
    + (mat212 b) * ((vec21 r) * s2 + (vec22 r) * s1)
    + (mat222 b) * (vec22 r) * s2
  "
  obtain t where 4: "2 * ¦a'12 + a11 * (t :: int)¦ ≤ a11"
    using 1 congruence_class_close by fastforce
  define a12 where "a12 ≡ a'12 + a11 * t"
  define a22 where "a22 ≡ b $$ (vec2 (s1 + (vec21 r) * t) (s2 + (vec22 r) * t))"
  let ?u = "
    mat2
      (vec21 r) (s1 + (vec21 r) * t)
      (vec22 r) (s2 + (vec22 r) * t)
  "
  let ?a = "
    mat2
      a11 a12
      a12 a22
  "
  have 5: "mat_det ?u = 1" unfolding mat2_det_def
    using 3 by (simp add: algebra_simps)
  have 6: "?a = b ⋅ ?u" using assms(1) 1
    unfolding qf2_action_def mat2_transpose_def qf2_app_def vec2_dot_def
              mat2_app_def times_mat2_def mat2_sym_def
              a12_def a12_def a22_def a'12_def
    by (simp add: algebra_simps)
  have "b ~ ?a" unfolding qf2_equiv_def using 5 6 by auto
  hence 7: "?a ~ b" using qf2_equiv_sym by blast
  have 8: "2 * ¦mat212 ?a¦ ≤ mat211 ?a" using 4 unfolding a12_def by auto
  have "a11 ≤ int (nat a22)"
    unfolding a11_def a22_def
    apply (rule rev_iffD1[OF _ nat_int_comparison(3)])
    apply (rule wellorder_Least_lemma(2))
    using assms(2) 5 unfolding qf2_positive_definite_def zero_vec2_def mat2_det_def
    apply (metis add.right_neutral diff_add_cancel mat2.sel
                 mult_zero_left mult_zero_right nat_0_le order_less_le vec2.inject
                 zero_less_nat_eq zero_neq_one)
    done
  hence "a11 ≤ a22" using 1 by linarith
  hence "4 * a112 ≤ 4 * a11 * a22" unfolding power2_eq_square using 1 by auto
  also have "... = 4 * (mat_det ?a + a122)"
    unfolding mat2_det_def power2_eq_square by auto
  also have "... = 4 * mat_det ?a + (2 * ¦a12¦)2"
    unfolding power2_eq_square by auto
  also have "... ≤ 4 * mat_det ?a + (int a11)2"
    using 4 power2_le_iff_abs_le unfolding a12_def by (smt (verit))
  finally have "3 * a112 ≤ 4 * (mat_det ?a)" by auto
  hence "a112 ≤ (4 / 3) * mat_det ?a" by linarith
  hence "sqrt (a112) ≤ sqrt ((4 / 3) * mat_det ?a)"
    using real_sqrt_le_mono by blast
  hence "a11 ≤ sqrt ((4 / 3) * mat_det ?a)" by auto
  hence "a11 ≤ (sqrt 4) / (sqrt 3) * sqrt (mat_det ?a)"
    unfolding real_sqrt_mult real_sqrt_divide by blast
  hence 9: "mat211 ?a ≤ (2 / sqrt 3) * sqrt (mat_det ?a)" by auto
  have "?P ?a" using 7 8 9 by blast
  thus ?thesis by blast
qed

text ‹Theorem 1.2 from @{cite nathanson1996}.›

theorem qf2_det_one_equiv_canonical:
  fixes f :: mat2
  assumes "mat_sym f"
  assumes "qf_positive_definite f"
  assumes "mat_det f = 1"
  shows "f ~ 1"
proof -
  obtain a where
    0: "f ~ a"
       "2 * ¦mat212 a¦ ≤ mat211 a"
       "mat211 a ≤ (2 / sqrt 3) * sqrt (mat_det a)"
    using assms lemma_1_2[of f] qf2_equiv_sym by auto
  have 1: "mat_sym a"
    using assms 0 qf2_equiv_preserves_sym by auto
  have 2: "qf_positive_definite a"
    using assms 0 qf2_equiv_preserves_positive_definite by auto
  have 3: "mat_det a = 1" using assms 0 qf2_equiv_preserves_det by auto
  have 4: "mat211 a ≥ 1"
    apply (rule allE[OF 2[unfolded qf2_positive_definite_def], of "vec2 1 0"])
    unfolding zero_vec2_def vec2_dot_def mat2_app_def qf2_app_def
              qf2_positive_definite_def
    apply auto
    done
  have 5: "mat211 a < 2" using 0 unfolding 3
    by (smt (verit, best) divide_le_eq int_less_real_le mult_cancel_left1
                          of_int_1 real_sqrt_four real_sqrt_le_iff
                          real_sqrt_mult_self real_sqrt_one)
  have 6: "mat211 a = 1" using 4 5 by auto
  have 7: "mat212 a = 0" "mat221 a = 0" using 0 1 6 unfolding mat2_sym_def by auto
  have 8: "mat222 a = 1" using 3 6 7 unfolding mat2_det_def by auto
  have "a = 1" using 6 7 8 unfolding one_mat2_def using mat2.collapse by metis
  thus ?thesis using 0 by metis
qed

text ‹Lemma 1.3 from @{cite nathanson1996}.›

lemma lemma_1_3:
  fixes a :: mat3
  assumes "mat_sym a"
  defines "a' ≡
    mat2
      (mat311 a * mat322 a - (mat312 a)2) (mat311 a * mat323 a - mat312 a * mat313 a)
      (mat311 a * mat323 a - mat312 a * mat313 a) (mat311 a * mat333 a - (mat313 a)2)
  "
  defines "d' ≡
    mat_det (
      mat2
        (mat311 a) (mat312 a)
        (mat312 a) (mat322 a)
    )
  "
  shows
    "mat_det a' = mat311 a * mat_det a" (is "?P")
    "⋀x. mat311 a * (a $$ x) =
      (mat311 a * vec31 x + mat312 a * vec32 x + mat313 a * vec33 x)2 +
      (a' $$ (vec2 (vec32 x) (vec33 x)))" (is "⋀x. ?Q x")
    "qf_positive_definite a ⟹ qf_positive_definite a'"
    "qf_positive_definite a ⟷ mat311 a > 0 ∧ d' > 0 ∧ mat_det a > 0"
proof -
  show 0: ?P using assms
    unfolding a'_def mat2_det_def mat3_det_def mat3_sym_def power2_eq_square
    by (simp add: algebra_simps)
  show 1: "⋀x. ?Q x" using assms
    unfolding a'_def vec2_dot_def vec3_dot_def mat2_app_def mat3_app_def
              mat3_sym_def qf2_app_def qf3_app_def power2_eq_square
    by (simp add: algebra_simps)
  have 2: "qf_positive_definite a ⟹ mat311 a > 0"
  proof -
    assume 3: "qf_positive_definite a"
    show "mat311 a > 0"
      using allE[OF iffD1[OF qf3_positive_definite_def 3], of "vec3 1 0 0"]
      unfolding zero_vec3_def vec3_dot_def mat3_app_def qf3_app_def
      by auto
  qed
  show 4: "qf_positive_definite a ⟹ qf_positive_definite a'"
  unfolding qf2_positive_definite_def
  proof
    fix v :: vec2
    assume 5: "qf_positive_definite a"
    hence 6: "mat311 a > 0" using 2 by blast
    have "a' $$ v ≤ 0 ⟹ v = 0" proof -
      assume 7: "a' $$ v ≤ 0"
      obtain x2 x3 where 8: "v = vec2 x2 x3" by (rule vec2.exhaust)
      define x1 where "x1 ≡ - (mat312 a * x2 + mat313 a * x3)"
      have "mat311 a * (a $$ (vec3 x1 (mat311 a * x2) (mat311 a * x3))) =
        (mat311 a * x1 + mat312 a * mat311 a * x2 + mat313 a * mat311 a * x3)2 +
        (a' $$ (vec2 (mat311 a * x2) (mat311 a * x3)))"
        unfolding 1 by (simp add: algebra_simps)
      also have "... = a' $$ (vec2 (mat311 a * x2) (mat311 a * x3))"
        unfolding x1_def by (simp add: algebra_simps)
      also have "... = (mat311 a)2 * (a' $$ v)"
        unfolding 8 vec2_dot_def mat2_app_def qf2_app_def power2_eq_square
        by (simp add: algebra_simps)
      also have "... ≤ 0"
        using 7 unfolding power2_eq_square by (simp add: mult_nonneg_nonpos)
      finally have "mat311 a * (a $$ (vec3 x1 (mat311 a * x2) (mat311 a * x3))) ≤ 0"
        .
      hence "a $$ (vec3 x1 (mat311 a * x2) (mat311 a * x3)) ≤ 0"
        using 6 by (simp add: mult_le_0_iff)
      hence "vec3 x1 (mat311 a * x2) (mat311 a * x3) = 0"
        using 5 unfolding qf3_positive_definite_def by fastforce
      hence "x2 = 0" "x3 = 0" using 6 unfolding zero_vec3_def by fastforce+
      thus "v = 0" unfolding 8 zero_vec2_def by blast
    qed
    thus "v ≠ 0 ⟶ 0 < a' $$ v" by fastforce
  qed
  have 9: "qf_positive_definite a ⟹ d' > 0 ∧ mat_det a > 0"
  proof -
    assume 10: "qf_positive_definite a"
    have 11: "mat311 a > 0" using 2 10 by blast
    have "qf_positive_definite a'" using 4 10 by blast
    hence 12: "mat211 a' > 0 ∧ mat_det a' > 0"
      using qf2_positive_definite_criterion
      unfolding a'_def mat2_sym_def by fastforce
    have 13: "d' > 0"
      using 12 unfolding a'_def d'_def mat2_det_def power2_eq_square by fastforce
    have 14: "mat_det a > 0"
      using 11 12 unfolding 0 by (simp add: zero_less_mult_iff)
    show "d' > 0 ∧ (mat_det a) > 0" using 13 14 by blast
  qed
  have 15: "mat311 a > 0 ∧ d' > 0 ∧ mat_det a > 0 ⟹ qf_positive_definite a"
  proof -
    assume 16: "mat311 a > 0 ∧ d' > 0 ∧ mat_det a > 0"
    have 17: "mat211 a' > 0"
      using 16 unfolding a'_def d'_def mat2_det_def power2_eq_square
      by (simp add: algebra_simps)
    have 18: "mat_det a' > 0" using 16 unfolding 0 by fastforce
    have 19: "qf_positive_definite a'"
      using qf2_positive_definite_criterion 17 18
      unfolding a'_def mat2_sym_def by fastforce
    show "qf_positive_definite a"
    unfolding qf3_positive_definite_def
    proof
      fix x :: vec3
      have "a $$ x ≤ 0 ⟹ x = 0" proof -
        assume "a $$ x ≤ 0"
        hence 20: "mat311 a * (a $$ x) ≤ 0"
          using 16 mult_le_0_iff order_less_le by auto
        have "a' $$ (vec2 (vec32 x) (vec33 x)) ≤ 0"
          using 20 unfolding 1 power2_eq_square by (smt (verit) zero_le_square)
        hence 21: "a' $$ (vec2 (vec32 x) (vec33 x)) = 0"
          using 19 qf2_positive_definite_positive
          using nle_le by blast
        have 22: "vec32 x = 0" "vec33 x = 0"
          using 19 21 unfolding zero_vec2_def qf2_positive_definite_def
          by (smt (verit) vec2.inject)+
        have "mat311 a * vec31 x + mat312 a * vec32 x + mat313 a * vec33 x = 0"
          using 20 21 unfolding 1 by fastforce
        hence "vec31 x = 0" using 16 22 by fastforce
        thus "x = 0" using 22 unfolding zero_vec3_def by (metis vec3.collapse)
      qed
      thus "x ≠ 0 ⟶ 0 < a $$ x" by fastforce
    qed
  qed
  show "qf_positive_definite a ⟷ mat311 a > 0 ∧ d' > 0 ∧ mat_det a > 0"
    using 2 9 15 by blast
qed

text ‹Lemma 1.4 from @{cite nathanson1996}.›

lemma lemma_1_4:
  fixes b :: mat3
  fixes v' :: mat2
  fixes r s :: int
  assumes "mat_sym b"
  assumes "qf_positive_definite b"
  assumes "mat_det v' = 1"
  defines "b' ≡
    mat2
      (mat311 b * mat322 b - (mat312 b)2) (mat311 b * mat323 b - mat312 b * mat313 b)
      (mat311 b * mat323 b - mat312 b * mat313 b) (mat311 b * mat333 b - (mat313 b)2)
  "
  defines "a' ≡ b' ⋅ v'"
  defines "v ≡
    mat3
      1 r s
      0 (mat211 v') (mat212 v')
      0 (mat221 v') (mat222 v')
  "
  defines "a ≡ b ⋅ v"
  shows
    "⋀y. mat311 b * (b $$ y) =
      (mat311 b * vec31 y + mat312 b * vec32 y + mat313 b * vec33 y)2 +
      (b' $$ (vec2 (vec32 y) (vec33 y)))" (is "⋀y. ?P y")
    "mat311 a = mat311 b"
    "⋀x. mat311 a * (a $$ x) =
      (mat311 a * vec31 x + mat312 a * vec32 x + mat313 a * vec33 x)2 +
      (a' $$ (vec2 (vec32 x) (vec33 x)))" (is "⋀x. ?Q x")
proof -
  show "⋀y. ?P y" using assms
    unfolding b'_def vec2_dot_def vec3_dot_def mat2_app_def mat3_app_def
              mat3_sym_def qf2_app_def qf3_app_def power2_eq_square
    by (simp add: algebra_simps)
  show "mat311 a = mat311 b"
    unfolding a_def v_def times_mat3_def mat3_transpose_def qf3_action_def
    by force
  show "⋀y. ?Q y" using assms
    by (simp add: algebra_simps power2_eq_square
                  a_def v_def a'_def b'_def vec2_dot_def vec3_dot_def
                  times_mat2_def times_mat3_def mat2_app_def mat3_app_def
                  mat2_transpose_def mat3_transpose_def mat3_sym_def
                  qf2_app_def qf3_app_def qf2_action_def qf3_action_def)
qed

text ‹Lemma 1.5 from @{cite nathanson1996}.›

lemma lemma_1_5:
  fixes u11 u21 u31
  assumes "Gcd {u11, u21, u31} = 1"
  shows "∃u. mat311 u = u11 ∧ mat321 u = u21 ∧ mat331 u = u31 ∧ mat_det u = 1"
proof -
  let ?a = "gcd u11 u21"
  show ?thesis proof cases
    assume "?a = 0"
    hence 0: "u11 = 0" "u21 = 0" "u31 = 1 ∨ u31 = -1" using assms by auto
    let ?u = "
      mat3
      0 0 (- 1)
      0 u31 0
      u31 0 0"
    show ?thesis
      apply (rule exI[of _ ?u])
      unfolding mat3_det_def using 0 apply auto
      done
  next
    assume 1: "?a ≠ 0"
    obtain u12 u22 where 2: "u11 * u22 - u21 * u12 = ?a" using bezout_int
      by (metis diff_minus_eq_add mult.commute mult_minus_right)
    have "gcd ?a u31 = 1" using assms by (simp add: gcd.assoc)
    then obtain u33 b where 3: "?a * u33 - b * u31 = 1" using bezout_int
      by (metis diff_minus_eq_add mult.commute mult_minus_right)
    let ?u = "
      mat3
      u11 u12 ((u11 div ?a) * b)
      u21 u22 ((u21 div ?a) * b)
      u31 0 u33"
    have "mat_det ?u =
          u11 * u22 * u33
        + u12 * ((u21 div ?a) * b) * u31
        - ((u11 div ?a) * b) * u22 * u31
        - u12 * u21 * u33"
      unfolding mat3_det_def by auto
    also have "... =
          u11 * u22 * u33
        + u12 * (u21 div ?a) * (b * u31)
        - u22 * (u11 div ?a) * (b * u31)
        - u12 * u21 * u33"
      by auto
    also have "... =
          u11 * u22 * u33
        + u12 * (u21 div ?a) * (?a * u33 - 1)
        - u22 * (u11 div ?a) * (?a * u33 - 1)
        - u12 * u21 * u33"
      using 3 by (simp add: algebra_simps)
    also have "... =
          u11 * u22 * u33
        + u12 * ((u21 div ?a) * ?a) * u33
        - u12 * (u21 div ?a)
        - u22 * ((u11 div ?a) * ?a) * u33
        + u22 * (u11 div ?a)
        - u12 * u21 * u33"
      by (simp add: algebra_simps)
    also have "... =
        - u12 * (u21 div ?a)
        + u22 * (u11 div ?a)"
      by (simp add: algebra_simps)
    also have "... =
        - u12 * (u21 div ?a)
        + u22 * u11 div ?a"
      by (metis dvd_div_mult gcd_dvd1 mult.commute)
    also have "... =
        - u12 * (u21 div ?a)
        + (?a + u21 * u12) div ?a"
      using 2 by (simp add: algebra_simps)
    also have "... =
        - u12 * (u21 div ?a)
        + 1 + (u21 * u12) div ?a"
      using 1 by auto
    also have "... = 1" by (simp add: algebra_simps div_mult1_eq)
    finally have 4: "mat_det ?u = 1" .
    show ?thesis
      apply (rule exI[of _ ?u])
      using 4 apply auto
      done
  qed
qed

text ‹Lemma 1.6 from @{cite nathanson1996}.›

lemma lemma_1_6:
  fixes c :: mat3
  assumes "mat_sym c"
  assumes "qf_positive_definite c"
  shows "∃a. a ~ c ∧
             2 * (max ¦mat312 a¦ ¦mat313 a¦) ≤ mat311 a ∧
             mat311 a ≤ (4 / 3) * root 3 (mat_det a)"
proof -
  define a11 where "a11 ≡ LEAST y. y > 0 ∧ (∃x. int y = c $$ x)"
  have 0: "∃y. y > 0 ∧ (∃x. int y = c $$ x)"
    apply (rule exI[of _ "nat (c $$ (vec3 1 1 1))"])
    using assms(2) unfolding qf3_positive_definite_def zero_vec3_def
    apply (metis nat_0_le order_less_le vec3.inject zero_less_nat_eq zero_neq_one)
    done
  obtain t where 1: "a11 > 0" "int a11 = c $$ t"
    using a11_def LeastI_ex[OF 0] by auto
  let ?h = "Gcd {vec31 t, vec32 t, vec33 t}"
  have "t ≠ 0" using assms(2) 1 by fastforce
  hence 2: "?h > 0" by (simp; metis vec3.collapse zero_vec3_def)
  let ?t' = "vec3 (vec31 t div ?h) (vec32 t div ?h) (vec33 t div ?h)"
  have "?t' ≠ 0" using 2 unfolding zero_vec3_def
    by (auto simp add: algebra_simps dvd_div_eq_0_iff)
  hence "nat (c $$ ?t') > 0 ∧ (∃x. int (nat (c $$ ?t')) = c $$ x)"
    using assms(2) unfolding qf3_positive_definite_def by auto
  hence "a11 ≤ nat (c $$ ?t')" unfolding a11_def by (rule Least_le)
  hence "a11 ≤ c $$ ?t'" using 1 by auto
  also have "... = (c $$ t) div ?h2" proof -
    have "?h dvd vec31 t" "?h dvd vec32 t" "?h dvd vec33 t"
      by (meson Gcd_dvd insertCI)+
    then have "(c $$ ?t') * ?h2 = c $$ t" "?h2 dvd c $$ t"
      unfolding qf3_app_def vec3_dot_def mat3_app_def power2_eq_square
      using 1
      by (auto simp add: algebra_simps mult_dvd_mono)
    thus "c $$ ?t' = (c $$ t) div ?h2" using 2 by auto
  qed
  also have "... = a11 div ?h2" using 1 by auto
  finally have "a11 ≤ a11 div ?h2" .
  also have "... ≤ a11" using 1 2
    by (metis div_le_dividend nat_div_distrib nat_int nat_le_iff of_nat_0_le_iff)
  finally have "?h = 1" using 1 2
    by (smt (verit) int_div_less_self of_nat_0_less_iff power2_eq_square
                    zero_less_power zmult_eq_1_iff)
  then obtain u where 3: "mat311 u = vec31 t" "mat321 u = vec32 t"
                         "mat331 u = vec33 t" "mat_det u = 1"
    using lemma_1_5 by blast
  define b where "b ≡ c ⋅ u"
  have 4: "mat_sym b"
    using 3 assms(1) qf3_equiv_preserves_sym
    unfolding b_def qf3_equiv_def
    by auto
  have 5: "qf_positive_definite b"
    using 3 assms(2) qf3_equiv_preserves_positive_definite
    unfolding b_def qf3_equiv_def
    by auto
  have 6: "a11 = (LEAST y. y > 0 ∧ (∃x. int y = b $$ x))"
    unfolding a11_def apply (rule arg_cong[of _ _ Least])
    using 3 qf3_equiv_preserves_range[of c b]
    unfolding b_def image_def qf3_equiv_def
    apply fast
    done
  have 7: "a11 = mat311 b"
    using 1 3
    by (simp add: algebra_simps
                  b_def times_mat3_def vec3_dot_def mat3_app_def
                  mat3_transpose_def qf3_app_def qf3_action_def)
  define b' where "b' ≡
    mat2
      (mat311 b * mat322 b - (mat312 b)2) (mat311 b * mat323 b - mat312 b * mat313 b)
      (mat311 b * mat323 b - mat312 b * mat313 b) (mat311 b * mat333 b - (mat313 b)2)
  "
  have 8: "mat_sym b'" unfolding b'_def mat2_sym_def by auto
  have 9: "mat_det b' = mat311 b * mat_det b"
          "⋀x. mat311 b * (b $$ x) =
            (mat311 b * vec31 x + mat312 b * vec32 x + mat313 b * vec33 x)2 +
            (b' $$ (vec2 (vec32 x) (vec33 x)))"
          "qf_positive_definite b'"
    using 4 5 b'_def lemma_1_3 by blast+
  obtain a' v' where 10: "a' = b' ⋅ v'"
                         "mat_det v' = 1"
                         "mat211 a' ≤ (2 / sqrt 3) * sqrt (mat311 b * mat_det b)"
    using 8 9 qf2_equiv_sym qf2_equiv_preserves_det lemma_1_2[of b']
    unfolding qf2_equiv_def by metis
  obtain r s where 11: "2 * ¦(mat312 b) * (mat211 v') +
                            (mat313 b) * (mat221 v') + a11 * (r :: int)¦ ≤ a11"
                       "2 * ¦(mat312 b) * (mat212 v') +
                            (mat313 b) * (mat222 v') + a11 * (s :: int)¦ ≤ a11"
    using 1 congruence_class_close by fastforce
  define a12 where "a12 ≡ (mat312 b) * (mat211 v') +
                          (mat313 b) * (mat221 v') + a11 * r"
  define a13 where "a13 ≡ (mat312 b) * (mat212 v') +
                          (mat313 b) * (mat222 v') + a11 * s"
  define v where "v ≡
    mat3
      1 r s
      0 (mat211 v') (mat212 v')
      0 (mat221 v') (mat222 v')
  "
  have 12: "mat_det v = 1"
    using 10 unfolding v_def mat2_det_def mat3_det_def by (simp add: algebra_simps)
  define a where "a ≡ b ⋅ v"
  have 13: "mat_det a = mat_det b"
    using 12 qf3_equiv_preserves_det
    unfolding a_def qf3_equiv_def
    by metis
  have 14: "a11 = (LEAST y. y > 0 ∧ (∃x. int y = a $$ x))"
    unfolding 6 apply (rule arg_cong[of _ _ Least])
    using 12 qf3_equiv_preserves_range[of b a]
    unfolding a_def image_def qf3_equiv_def
    apply fast
    done
  have 15: "mat311 a = mat311 b"
           "⋀x. mat311 a * (a $$ x) =
              (mat311 a * vec31 x + mat312 a * vec32 x + mat313 a * vec33 x)2 +
              (a' $$ (vec2 (vec32 x) (vec33 x)))"
    using 4 5 10 lemma_1_4 unfolding b'_def v_def a_def by blast+
  have 16: "mat211 a' = mat311 a * mat322 a - (mat312 a)2"
    using 15(2)[of "vec3 0 1 0"]
    unfolding vec3_dot_def vec2_dot_def mat3_app_def mat2_app_def
              qf3_app_def qf2_app_def
    by simp
  define a22 where "a22 ≡ a $$ (vec3 0 1 0)"
  have 17: "a11 = mat311 a" unfolding 7 15 by auto
  have 18: "a12 = mat312 a" "a13 = mat313 a" "a22 = mat322 a" "a22 = mat322 a"
    using 13(1) 17
    unfolding a_def v_def a12_def a13_def a22_def
    by (auto simp add: algebra_simps
                       times_mat3_def vec3_dot_def mat3_app_def
                       mat3_transpose_def qf3_app_def qf3_action_def)
  have 19: "a ~ c"
    unfolding qf3_equiv_sym[of a c]
    unfolding a_def b_def qf3_equiv_def qf3_action_mul[symmetric]
    using 3 12 by (metis mat3_mul_det mult_1)
  have 20: "2 * (max ¦mat312 a¦ ¦mat313 a¦) ≤ mat311 a"
    using 11 unfolding 17[symmetric] 18 a12_def[symmetric] a13_def[symmetric]
    by auto
  have 21: "(2 / sqrt 3) ^ 6 = 64 / 27" unfolding power_def by auto
  have "a11 ≤ int (nat a22)"
    unfolding a11_def a22_def
    apply (rule rev_iffD1[OF _ nat_int_comparison(3)])
    apply (rule wellorder_Least_lemma(2))
    using assms(2) 5 12
    apply (metis a_def b_def nat_0_le nat_int of_nat_0 qf3_action_app
                 qf3_equiv_def qf3_equiv_preserves_positive_definite
                 qf3_positive_definite_def qf3_positive_definite_positive
                 vec3.sel(2) zero_neq_one zero_vec3_def zless_nat_conj)
    done
  hence "a11 ≤ a22" using 1 by linarith
  hence "(mat311 a)2 ≤ a11 * a22" unfolding power2_eq_square using 1 17 by auto
  also have "... = mat211 a' + a122" using 16 17 18 by auto
  also have "... ≤ (2 / sqrt 3) * sqrt (mat311 b * mat_det b) + a122"
    using 10 by auto
  also have "... ≤ (2 / sqrt 3) * sqrt (mat311 a * mat_det a) + (mat311 a)2 / 4"
    using 11 13 15 17 18 a12_def
          power2_le_iff_abs_le[of "real_of_int (mat311 a)" "(mat312 a) * 2"]
    by auto
  finally have "(3 / 4) * (mat311 a)2 ≤ (2 / sqrt 3) * sqrt (mat311 a * mat_det a)"
    by (simp add: algebra_simps)
  hence "(mat311 a)2 ≤ (2 / sqrt 3) ^ 3 * sqrt (mat311 a * mat_det a)"
    unfolding power2_eq_square power3_eq_cube by (simp add: algebra_simps)
  hence "((mat311 a)2)2 ≤ ((2 / sqrt 3) ^ 3 * sqrt (mat311 a * mat_det a))2"
    using 1(1) unfolding 17[symmetric] power2_eq_square
    by (metis of_int_power power2_eq_square power_mono zero_le_square)
  hence "(mat311 a) ^ 4 ≤ (2 / sqrt 3) ^ 6 * (sqrt (mat311 a * mat_det a))2"
    by (simp add: power_mult_distrib)
  hence "(mat311 a) ^ 4 ≤ (2 / sqrt 3) ^ 6 * mat311 a * mat_det a"
    using 4 5 13 17 lemma_1_3 by auto
  hence "(mat311 a) ^ 3 ≤ (2 / sqrt 3) ^ 6 * mat_det a"
    using 1(1) unfolding 17[symmetric]
    unfolding power_def apply (simp add: algebra_simps)
    apply (metis mult_left_le_imp_le of_nat_0_less_iff times_divide_eq_right)
    done
  hence "(mat311 a) ^ 3 ≤ (64 / 27) * mat_det a" using 21 by auto
  hence "root 3 ((mat311 a) ^ 3) ≤ root 3 ((64 / 27) * mat_det a)" by auto
  hence "root 3 ((mat311 a) ^ 3) ≤ root 3 (64 / 27) * root 3 (mat_det a)"
    unfolding real_root_mult by auto
  hence "mat311 a ≤ root 3 (64 / 27) * root 3 (mat_det a)"
    using odd_real_root_power_cancel by auto
  hence 22: "mat311 a ≤ (4 / 3) * root 3 (mat_det a)"
    using real_root_divide by force
  show ?thesis using 19 20 22 by blast
qed

text ‹Theorem 1.3 from @{cite nathanson1996}.›

theorem qf3_det_one_equiv_canonical:
  fixes f :: mat3
  assumes "mat_sym f"
  assumes "qf_positive_definite f"
  assumes "mat_det f = 1"
  shows "f ~ 1"
proof -
  obtain a where
    0: "f ~ a ∧
        2 * (max ¦mat312 a¦ ¦mat313 a¦) ≤ mat311 a ∧
        mat311 a ≤ (4 / 3) * root 3 (mat_det a)"
    using assms lemma_1_6[of f] qf3_equiv_sym by auto
  have 1: "mat_sym a"
    using assms 0 qf3_equiv_preserves_sym by auto
  have 2: "qf_positive_definite a"
    using assms 0 qf3_equiv_preserves_positive_definite by auto
  have 3: "mat_det a = 1" using assms 0 qf3_equiv_preserves_det by auto
  have 4: "mat312 a = 0" "mat313 a = 0" using 0 3 by auto
  have "mat311 a ≥ 1"
    apply (rule allE[OF 2[unfolded qf3_positive_definite_def], of "vec3 1 0 0"])
    unfolding zero_vec3_def vec3_dot_def mat3_app_def qf3_app_def
              qf3_positive_definite_def
    apply auto
    done
  hence 6: "mat311 a = 1" using 0 3 by auto
  define a' where "a' ≡
    mat2
      (mat322 a) (mat323 a)
      (mat332 a) (mat333 a)
  "
  have 7: "mat_det a' = 1"
    using 3 4 6 unfolding a'_def mat2_det_def mat3_det_def by auto
  have 8: "mat_sym a'" using 1 unfolding a'_def mat2_sym_def mat3_sym_def by auto
  have 9: "qf_positive_definite a'"
  using 2 unfolding qf2_positive_definite_def qf3_positive_definite_def
  proof -
    assume 10: "∀v. v ≠ 0 ⟶ a $$ v > 0"
    show "∀v. v ≠ 0 ⟶ a' $$ v > 0"
    proof (rule; rule)
      fix v :: vec2
      assume "v ≠ 0"
      hence "vec3 0 (vec21 v) (vec22 v) ≠ 0"
        unfolding zero_vec2_def zero_vec3_def by (metis vec2.collapse vec3.inject)
      hence "a $$ (vec3 0 (vec21 v) (vec22 v)) > 0" using 10 by auto
      thus "a' $$ v > 0"
        unfolding a'_def vec2_dot_def vec3_dot_def
                  mat2_app_def mat3_app_def qf2_app_def qf3_app_def
                  qf2_positive_definite_def qf3_positive_definite_def
        by auto
    qed
  qed
  obtain u' where 11: "mat_det u' = 1" "a' ⋅ u' = 1"
    using 7 8 9 qf2_det_one_equiv_canonical[of a']
    unfolding qf2_equiv_def
    by auto
  define u where "u ≡
    mat3
      1 0 0
      0 (mat211 u') (mat212 u')
      0 (mat221 u') (mat222 u')
  "
  have 12: "mat_det u = 1"
    using 11 unfolding u_def mat2_det_def mat3_det_def by auto
  have 13: "a ⋅ u = 1"
    using 1 4 6 11
    by (simp add: algebra_simps
                  a'_def u_def one_mat2_def one_mat3_def
                  times_mat2_def times_mat3_def
                  mat2_transpose_def mat3_transpose_def
                  qf2_action_def qf3_action_def mat3_sym_def)
  have "a ~ 1" using 12 13 unfolding qf3_equiv_def by blast
  thus ?thesis using 0 qf3_equiv_trans by blast
qed

end