diff --git a/thys/Budan_Fourier/BF_Misc.thy b/thys/Budan_Fourier/BF_Misc.thy
--- a/thys/Budan_Fourier/BF_Misc.thy
+++ b/thys/Budan_Fourier/BF_Misc.thy
@@ -1,1294 +1,1294 @@
 (*
     Author:     Wenda Li <wl302@cam.ac.uk / liwenda1990@hotmail.com>
 *)
 
 section "Misc results for polynomials and sign variations"
 
 theory BF_Misc imports
   "HOL-Computational_Algebra.Polynomial_Factorial"
   "HOL-Computational_Algebra.Fundamental_Theorem_Algebra"
   Sturm_Tarski.Sturm_Tarski
 begin
 
 subsection \<open>Induction on polynomial roots\<close>
 
 (*adapted from the poly_root_induct in Polynomial.thy.*)
 lemma poly_root_induct_alt [case_names 0 no_proots root]:
   fixes p :: "'a :: idom poly"
   assumes "Q 0"
   assumes "\<And>p. (\<And>a. poly p a \<noteq> 0) \<Longrightarrow> Q p"
   assumes "\<And>a p. Q p \<Longrightarrow> Q ([:-a, 1:] * p)"
   shows   "Q p"
 proof (induction "degree p" arbitrary: p rule: less_induct)
   case (less p)
   have ?case when "p=0" using \<open>Q 0\<close> that by auto
   moreover have ?case when "\<nexists>a. poly p a = 0"
     using assms(2) that by blast
   moreover have ?case when "\<exists>a. poly p a = 0" "p\<noteq>0"
   proof -
     obtain a where "poly p a =0" using \<open>\<exists>a. poly p a = 0\<close> by auto
     then obtain q where pq:"p= [:-a,1:] * q" by (meson dvdE poly_eq_0_iff_dvd)
     then have "q\<noteq>0" using \<open>p\<noteq>0\<close> by auto
     then have "degree q<degree p" unfolding pq by (subst degree_mult_eq,auto)
     then have "Q q" using less by auto
     then show ?case using assms(3) unfolding pq by auto
   qed
   ultimately show ?case by auto
 qed
 
 subsection \<open>Misc\<close>
 
 lemma lead_coeff_pderiv:
   fixes p :: "'a::{comm_semiring_1,semiring_no_zero_divisors,semiring_char_0} poly"
   shows "lead_coeff (pderiv p) = of_nat (degree p) * lead_coeff p"
   apply (auto simp:degree_pderiv coeff_pderiv)
   apply (cases "degree p")
   by (auto simp add: coeff_eq_0)
 
 lemma gcd_degree_le_min:
   assumes "p\<noteq>0" "q\<noteq>0"
   shows "degree (gcd p q) \<le> min (degree p) (degree q)"
   by (simp add: assms(1) assms(2) dvd_imp_degree_le)
 
 lemma lead_coeff_normalize_field:
   fixes p::"'a::{field,semidom_divide_unit_factor} poly"
   assumes "p\<noteq>0"
   shows "lead_coeff (normalize p) = 1"
   by (metis (no_types, lifting) assms coeff_normalize divide_self_if dvd_field_iff 
       is_unit_unit_factor leading_coeff_0_iff normalize_eq_0_iff normalize_idem)
 
 lemma smult_normalize_field_eq:
   fixes p::"'a::{field,semidom_divide_unit_factor} poly"
   shows "p = smult (lead_coeff p) (normalize p)"
 proof (rule poly_eqI)
   fix n
   have "unit_factor (lead_coeff p) = lead_coeff p"
     by (metis dvd_field_iff is_unit_unit_factor unit_factor_0)
   then show "coeff p n = coeff (smult (lead_coeff p) (normalize p)) n"
     by simp
 qed
 
 lemma lead_coeff_gcd_field:
   fixes p q::"'a::field_gcd poly"
   assumes "p\<noteq>0 \<or> q\<noteq>0"
   shows "lead_coeff (gcd p q) = 1"
   using assms by (metis gcd.normalize_idem gcd_eq_0_iff lead_coeff_normalize_field)
 
 lemma poly_gcd_0_iff:
   "poly (gcd p q) x = 0 \<longleftrightarrow> poly p x=0 \<and> poly q x=0"
   by (simp add:poly_eq_0_iff_dvd)
 
 lemma degree_eq_oneE:
   fixes p :: "'a::zero poly"
   assumes "degree p = 1"
   obtains a b where "p = [:a,b:]" "b\<noteq>0"
 proof -
   obtain a b q where p:"p=pCons a (pCons b q)"
     by (metis pCons_cases)
   with assms have "q=0" by (cases "q = 0") simp_all
   with p have "p=[:a,b:]" by auto
   moreover then have "b\<noteq>0" using assms by auto
   ultimately show ?thesis ..
 qed
 
 subsection \<open>More results about sign variations (i.e. @{term changes}\<close>
 
 lemma changes_0[simp]:"changes (0#xs) = changes xs"
   by (cases xs) auto
 
 lemma changes_Cons:"changes (x#xs) = (if filter (\<lambda>x. x\<noteq>0) xs = [] then 
                             0 
                           else if x* hd (filter (\<lambda>x. x\<noteq>0) xs) < 0 then 
                             1 + changes xs 
                           else changes xs)"
   apply (induct xs)
   by auto
 
 lemma changes_filter_eq:
   "changes (filter (\<lambda>x. x\<noteq>0) xs) = changes xs"
   apply (induct xs)
   by (auto simp add:changes_Cons)
 
 lemma changes_filter_empty:
   assumes "filter (\<lambda>x. x\<noteq>0) xs = []"
   shows "changes xs = 0" "changes (a#xs) = 0" using assms
   apply (induct xs)
   apply auto
   by (metis changes_0 neq_Nil_conv)
 
 lemma changes_append:
   assumes "xs\<noteq> [] \<and> ys\<noteq> [] \<longrightarrow> (last xs = hd ys \<and> last xs\<noteq>0)"
   shows "changes (xs@ys) = changes xs + changes ys"
   using assms
 proof (induct xs)
   case Nil
   then show ?case by simp
 next
   case (Cons a xs)
   have ?case when "xs=[]"
     using that Cons 
     apply (cases ys)
     by auto
   moreover have ?case when "ys=[]"
     using that Cons by auto
   moreover have ?case when "xs\<noteq>[]" "ys\<noteq>[]"
   proof -
     have "filter (\<lambda>x. x \<noteq> 0) xs \<noteq>[]"
       using that Cons 
       apply auto 
       by (metis (mono_tags, lifting) filter.simps(1) filter.simps(2) filter_append snoc_eq_iff_butlast)
     then have "changes (a # xs @ ys) = changes (a # xs) + changes ys"
       apply (subst (1 2) changes_Cons)
       using that Cons by auto
     then show ?thesis by auto
   qed
   ultimately show ?case by blast
 qed
 
 lemma changes_drop_dup:
   assumes "xs\<noteq> []" "ys\<noteq> [] \<longrightarrow> last xs=hd ys"
   shows "changes (xs@ys) = changes (xs@ tl ys)"
   using assms
 proof (induct xs)
   case Nil
   then show ?case by simp
 next
   case (Cons a xs)
   have ?case when "ys=[]"
     using that by simp
   moreover have ?case when "ys\<noteq>[]" "xs=[]"
     using that Cons
     apply auto
     by (metis changes.simps(3) list.exhaust_sel not_square_less_zero)
   moreover have ?case when "ys\<noteq>[]" "xs\<noteq>[]"
   proof -
     define ts ts' where "ts = filter (\<lambda>x. x \<noteq> 0) (xs @ ys)"
       and "ts' = filter (\<lambda>x. x \<noteq> 0) (xs @ tl ys)"
     have "(ts = [] \<longleftrightarrow> ts' = []) \<and> hd ts = hd ts'"
     proof (cases "filter (\<lambda>x. x \<noteq> 0) xs = []")
       case True
       then have "last xs = 0" using \<open>xs\<noteq>[]\<close> 
         by (metis (mono_tags, lifting) append_butlast_last_id append_is_Nil_conv 
             filter.simps(2) filter_append list.simps(3))
       then have "hd ys=0" using Cons(3)[rule_format, OF \<open>ys\<noteq>[]\<close>] \<open>xs\<noteq>[]\<close> by auto
       then have "filter (\<lambda>x. x \<noteq> 0) ys = filter (\<lambda>x. x \<noteq> 0) (tl ys)"
         by (metis (mono_tags, lifting) filter.simps(2) list.exhaust_sel that(1))
       then show ?thesis unfolding ts_def ts'_def by auto
     next
       case False
       then show ?thesis unfolding ts_def ts'_def by auto
     qed
     moreover have "changes (xs @ ys) = changes (xs @ tl ys)"
       apply (rule Cons(1))
       using that Cons(3) by auto
     moreover have "changes (a # xs @ ys) = (if ts = [] then 0 else if a * hd ts < 0 
             then 1 + changes (xs @ ys) else changes (xs @ ys))"
       using changes_Cons[of a "xs @ ys",folded ts_def] .
     moreover have "changes (a # xs @ tl ys) = (if ts' = [] then 0 else if a * hd ts' < 0 
             then 1 + changes (xs @ tl ys) else changes (xs @ tl ys))"
       using changes_Cons[of a "xs @ tl ys",folded ts'_def] .
     ultimately show ?thesis by auto
   qed
   ultimately show ?case by blast
 qed
 
 (*
   TODO: the following lemmas contain duplicates of some lemmas in 
           Winding_Number_Eval/Missing_Algebraic.thy
   Will resolve later.  
 *)
 
 lemma Im_poly_of_real:
   "Im (poly p (of_real x)) = poly (map_poly Im p) x"
   apply (induct p)
   by (auto simp add:map_poly_pCons)
  
 lemma Re_poly_of_real:
   "Re (poly p (of_real x)) = poly (map_poly Re p) x"
   apply (induct p)
   by (auto simp add:map_poly_pCons)
 
 subsection \<open>More about @{term map_poly} and @{term of_real}\<close>
 
 lemma of_real_poly_map_pCons[simp]:"map_poly of_real (pCons a p) = pCons (of_real a) (map_poly of_real p)" 
   by (simp add: map_poly_pCons)    
     
 lemma of_real_poly_map_plus[simp]: "map_poly of_real (p + q) = map_poly of_real p +  map_poly of_real q" 
   apply (rule poly_eqI)
   by (auto simp add: coeff_map_poly)  
  
 lemma of_real_poly_map_smult[simp]:"map_poly of_real (smult s p) = smult (of_real s) (map_poly of_real p)" 
   apply (rule poly_eqI)
   by (auto simp add: coeff_map_poly)    
 
 lemma of_real_poly_map_mult[simp]:"map_poly of_real (p*q) = map_poly of_real p * map_poly of_real q"
   by (induct p,intro poly_eqI,auto)
     
 lemma of_real_poly_map_poly:
   "of_real (poly p x) = poly (map_poly of_real p) (of_real x)" 
   by (induct p,auto)    
 
 lemma of_real_poly_map_power:"map_poly of_real (p^n) = (map_poly of_real p) ^ n"
   by (induct n,auto)
 
 (*FIXME: not duplicate*)
 lemma of_real_poly_eq_iff [simp]: "map_poly of_real p = map_poly of_real q \<longleftrightarrow> p = q"
   by (auto simp: poly_eq_iff coeff_map_poly)
 
 (*FIXME: not duplicate*)
 lemma of_real_poly_eq_0_iff [simp]: "map_poly of_real p = 0 \<longleftrightarrow> p = 0"
   by (auto simp: poly_eq_iff coeff_map_poly)
 
 subsection \<open>More about @{term order}\<close>
 
 lemma order_multiplicity_eq:
   assumes "p\<noteq>0"
   shows "order a p = multiplicity [:-a,1:] p"
   by (metis assms multiplicity_eqI order_1 order_2)
 
 lemma order_gcd:
   assumes "p\<noteq>0" "q\<noteq>0"
   shows "order x (gcd p q) = min (order x p) (order x q)"
 proof -
   have "prime [:- x, 1:]" 
     apply (auto simp add: prime_elem_linear_poly normalize_poly_def  intro!:primeI)
     by (simp add: pCons_one)
   then show ?thesis  
     using assms
     by (auto simp add:order_multiplicity_eq intro:multiplicity_gcd)
 qed
 
 lemma order_linear[simp]: "order x [:-a,1:] = (if x=a then 1 else 0)"
   by (auto simp add:order_power_n_n[where n=1,simplified] order_0I)
   
 lemma map_poly_order_of_real:
   assumes "p\<noteq>0"
   shows "order (of_real t) (map_poly of_real p) = order t p" using assms
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by simp
 next
   case (no_proots p)
   then have "order t p = 0" using order_root by blast
   moreover have "poly (map_poly of_real p) (of_real x) \<noteq>0" for x
     apply (subst of_real_poly_map_poly[symmetric])
     using no_proots order_root by simp
   then have "order (of_real t) (map_poly of_real p) = 0"
     using order_root by blast
   ultimately show ?case by auto
 next
   case (root a p)
   define a1 where "a1=[:-a,1:]"
   have [simp]:"a1\<noteq>0" "p\<noteq>0" unfolding a1_def using root(2) by auto
   have "order (of_real t) (map_poly of_real a1) = order t a1"
     unfolding a1_def by simp
   then show ?case 
     apply (fold a1_def)
     by (simp add:order_mult root)
 qed
 
 lemma order_pcompose:
   assumes "pcompose p q\<noteq>0"
   shows "order x (pcompose p q) = order x (q-[:poly q x:]) * order (poly q x) p" 
   using \<open>pcompose p q\<noteq>0\<close>
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by simp
 next
   case (no_proots p)
   have "order x (p \<circ>\<^sub>p q) = 0"
     apply (rule order_0I)
     using no_proots by (auto simp:poly_pcompose)
   moreover have "order (poly q x) p = 0"
     apply (rule order_0I)
     using no_proots by (auto simp:poly_pcompose)
   ultimately show ?case by auto
 next
   case (root a p)
   define a1 where "a1=[:-a,1:]"
   have [simp]: "a1\<noteq>0" "p\<noteq>0" "a1 \<circ>\<^sub>p q \<noteq>0" "p \<circ>\<^sub>p q \<noteq> 0" 
     subgoal using root(2) unfolding a1_def by simp
     subgoal using root(2) by auto
     using root(2) by (fold a1_def,auto simp:pcompose_mult)
   have "order x ((a1 * p) \<circ>\<^sub>p q) = order x (a1  \<circ>\<^sub>p q) + order x (p \<circ>\<^sub>p q)"
     unfolding pcompose_mult by (auto simp: order_mult)
   also have "... = order x (q-[:poly q x:]) * (order (poly q x) a1 + order (poly q x) p)"
   proof -
     have "order x (a1  \<circ>\<^sub>p q) = order x (q-[:poly q x:]) * order (poly q x) a1"
       unfolding a1_def 
       apply (auto simp: pcompose_pCons algebra_simps diff_conv_add_uminus )
       by (simp add: order_0I)
     moreover have "order x (p \<circ>\<^sub>p q) = order x (q - [:poly q x:]) * order (poly q x) p"
       apply (rule root.hyps)
       by auto
     ultimately show ?thesis by (auto simp:algebra_simps)
   qed
   also have "... =  order x (q - [:poly q x:]) * order (poly q x) (a1 * p)"
     by (auto simp:order_mult)
   finally show ?case unfolding a1_def .
 qed
 
 subsection \<open>Polynomial roots / zeros\<close>
 
 definition proots_within::"'a::comm_semiring_0 poly \<Rightarrow> 'a set \<Rightarrow> 'a set" where
   "proots_within p s = {x\<in>s. poly p x=0}"
 
 abbreviation proots::"'a::comm_semiring_0 poly \<Rightarrow> 'a set" where
   "proots p \<equiv> proots_within p UNIV"
 
 lemma proots_def: "proots p = {x. poly p x=0}" 
   unfolding proots_within_def by auto 
 
 lemma proots_within_empty[simp]:
   "proots_within p {} = {}" unfolding proots_within_def by auto
 
 lemma proots_within_0[simp]:
   "proots_within 0 s = s" unfolding proots_within_def by auto
 
 lemma proots_withinI[intro,simp]:
   "poly p x=0 \<Longrightarrow> x\<in>s \<Longrightarrow> x\<in>proots_within p s"
   unfolding proots_within_def by auto
 
 lemma proots_within_iff[simp]:
   "x\<in>proots_within p s \<longleftrightarrow> poly p x=0 \<and> x\<in>s"
   unfolding proots_within_def by auto
 
 lemma proots_within_union:
   "proots_within p A \<union> proots_within p B = proots_within p (A \<union> B)"
   unfolding proots_within_def by auto
 
 lemma proots_within_times:
   fixes s::"'a::{semiring_no_zero_divisors,comm_semiring_0} set"
   shows "proots_within (p*q) s = proots_within p s \<union> proots_within q s"
   unfolding proots_within_def by auto
 
 lemma proots_within_gcd:
   fixes s::"'a::{factorial_ring_gcd,semiring_gcd_mult_normalize} set"
   shows "proots_within (gcd p q) s= proots_within p s \<inter> proots_within q s"
   unfolding proots_within_def 
   by (auto simp add: poly_eq_0_iff_dvd) 
 
 lemma proots_within_inter:
   "NO_MATCH UNIV s \<Longrightarrow> proots_within p s = proots p \<inter> s" 
   unfolding proots_within_def by auto
 
 lemma proots_within_proots[simp]:
   "proots_within p s \<subseteq> proots p"
   unfolding proots_within_def by auto
 
 lemma finite_proots[simp]: 
   fixes p :: "'a::idom poly"
   shows "p\<noteq>0 \<Longrightarrow> finite (proots_within p s)"
   unfolding proots_within_def using poly_roots_finite by fast
 
 lemma proots_within_pCons_1_iff:
   fixes a::"'a::idom"
   shows "proots_within [:-a,1:] s = (if a\<in>s then {a} else {})"
     "proots_within [:a,-1:] s = (if a\<in>s then {a} else {})"
   by (cases "a\<in>s",auto)
 
 lemma proots_within_uminus[simp]:
   fixes p :: "'a::comm_ring poly"
   shows "proots_within (- p) s = proots_within p s"
   by auto
 
 lemma proots_within_smult:
   fixes a::"'a::{semiring_no_zero_divisors,comm_semiring_0}"
   assumes "a\<noteq>0"
   shows "proots_within (smult a p) s = proots_within p s"
   unfolding proots_within_def using assms by auto 
 
 subsection \<open>Polynomial roots counting multiplicities.\<close>
 
 (*counting the number of proots WITH MULTIPLICITIES within a set*)
 definition proots_count::"'a::idom poly \<Rightarrow> 'a set \<Rightarrow> nat" where
   "proots_count p s = (\<Sum>r\<in>proots_within p s. order r p)"  
 
 lemma proots_count_emtpy[simp]:"proots_count p {} = 0"
   unfolding proots_count_def by auto
 
 lemma proots_count_times:
   fixes s :: "'a::idom set"
   assumes "p*q\<noteq>0"
   shows "proots_count (p*q) s = proots_count p s + proots_count q s"
 proof -
   define pts where "pts=proots_within p s" 
   define qts where "qts=proots_within q s"
   have [simp]: "finite pts" "finite qts" 
     using \<open>p*q\<noteq>0\<close> unfolding pts_def qts_def by auto
   have "(\<Sum>r\<in>pts \<union> qts. order r p) =  (\<Sum>r\<in>pts. order r p)" 
   proof (rule comm_monoid_add_class.sum.mono_neutral_cong_right,simp_all)
     show "\<forall>i\<in>pts \<union> qts - pts. order i p = 0" 
       unfolding pts_def qts_def proots_within_def using order_root by fastforce
   qed
   moreover have "(\<Sum>r\<in>pts \<union> qts. order r q) = (\<Sum>r\<in>qts. order r q)" 
   proof (rule comm_monoid_add_class.sum.mono_neutral_cong_right,simp_all)
     show "\<forall>i\<in>pts \<union> qts - qts. order i q = 0" 
       unfolding pts_def qts_def proots_within_def using order_root by fastforce
   qed
   ultimately show ?thesis unfolding proots_count_def
     apply (simp add:proots_within_times order_mult[OF \<open>p*q\<noteq>0\<close>] sum.distrib)
     apply (fold pts_def qts_def)
     by auto
 qed
 
 lemma proots_count_power_n_n:
   shows "proots_count ([:- a, 1:]^n) s = (if a\<in>s \<and> n>0 then n else 0)"
 proof -
   have "proots_within ([:- a, 1:] ^ n) s= (if a\<in>s \<and> n>0 then {a} else {})"
     unfolding proots_within_def by auto
   thus ?thesis unfolding proots_count_def using order_power_n_n by auto
 qed
 
 lemma degree_proots_count:
   fixes p::"complex poly"
   shows "degree p = proots_count p UNIV"
 proof (induct "degree p" arbitrary:p )
   case 0
   then obtain c where c_def:"p=[:c:]" using degree_eq_zeroE by auto
   then show ?case unfolding proots_count_def 
     apply (cases "c=0")
     by (auto intro!:sum.infinite simp add:infinite_UNIV_char_0 order_0I)
 next
   case (Suc n)
   then have "degree p\<noteq>0" and "p\<noteq>0" by auto
   obtain z where "poly p z = 0" 
     using Fundamental_Theorem_Algebra.fundamental_theorem_of_algebra \<open>degree p\<noteq>0\<close> constant_degree[of p]
     by auto
   define onez where "onez=[:-z,1:]" 
   have [simp]: "onez\<noteq>0" "degree onez = 1" unfolding onez_def by auto
   obtain q where q_def:"p= onez * q" 
     using poly_eq_0_iff_dvd \<open>poly p z = 0\<close> dvdE unfolding onez_def by blast
   hence "q\<noteq>0" using \<open>p\<noteq>0\<close> by auto
   hence "n=degree q" using degree_mult_eq[of onez q] \<open>Suc n = degree p\<close> 
     apply (fold q_def)
     by auto
   hence "degree q = proots_count q UNIV" using Suc.hyps(1) by simp
   moreover have " Suc 0 = proots_count onez UNIV" 
     unfolding onez_def using proots_count_power_n_n[of z 1 UNIV]
     by auto
   ultimately show ?case 
     unfolding q_def using degree_mult_eq[of onez q] proots_count_times[of onez q UNIV] \<open>q\<noteq>0\<close>
     by auto
 qed
 
 lemma proots_count_smult:
   fixes a::"'a::{semiring_no_zero_divisors,idom}"
   assumes "a\<noteq>0"
   shows "proots_count (smult a p) s= proots_count p s"
 proof (cases "p=0")
   case True
   then show ?thesis by auto
 next
   case False
   then show ?thesis 
     unfolding proots_count_def
     using order_smult[OF assms] proots_within_smult[OF assms] by auto
 qed
 
 
 lemma proots_count_pCons_1_iff:
   fixes a::"'a::idom"
   shows "proots_count [:-a,1:] s = (if a\<in>s then 1 else 0)"
   unfolding proots_count_def 
   by (cases "a\<in>s",auto simp add:proots_within_pCons_1_iff order_power_n_n[of _ 1,simplified])
 
 lemma proots_count_uminus[simp]:
   "proots_count (- p) s = proots_count p s"
   unfolding proots_count_def by simp
 
 lemma card_proots_within_leq:
   assumes "p\<noteq>0"
   shows "proots_count p s \<ge> card (proots_within p s)" using assms
 proof (induct rule:poly_root_induct[of _ "\<lambda>x. x\<in>s"])
   case 0
   then show ?case unfolding proots_within_def proots_count_def by auto
 next
   case (no_roots p)
   then have "proots_within p s = {}" by auto
   then show ?case unfolding proots_count_def by auto
 next
   case (root a p)
   have "card (proots_within ([:- a, 1:] * p) s) 
       \<le> card (proots_within [:- a, 1:] s)+card (proots_within p s)" 
     unfolding proots_within_times by (auto simp add:card_Un_le)
   also have "... \<le> 1+ proots_count p s"
   proof -
     have "card (proots_within [:- a, 1:] s) \<le> 1"
     proof (cases "a\<in>s")
       case True
       then have "proots_within [:- a, 1:] s = {a}" by auto
       then show ?thesis by auto
     next
       case False
       then have "proots_within [:- a, 1:] s = {}" by auto
       then show ?thesis by auto
     qed
     moreover have "card (proots_within p s) \<le> proots_count p s"
       apply (rule root.hyps)
       using root by auto
     ultimately show ?thesis by auto
   qed
   also have "... =  proots_count ([:- a,1:] * p) s"
     apply (subst proots_count_times)
     subgoal by (metis mult_eq_0_iff pCons_eq_0_iff root.prems zero_neq_one)  
     using root by (auto simp add:proots_count_pCons_1_iff)
   finally have "card (proots_within ([:- a, 1:] * p) s) \<le> proots_count ([:- a, 1:] * p) s" .
   then show ?case 
     by (metis (no_types, opaque_lifting) add.inverse_inverse add.inverse_neutral minus_pCons 
         mult_minus_left proots_count_uminus proots_within_uminus)
 qed
 
 (*FIXME: not duplicate*)
 lemma proots_count_0_imp_empty:
   assumes "proots_count p s=0" "p\<noteq>0"
   shows "proots_within p s = {}"
 proof -
   have "card (proots_within p s) = 0"
     using card_proots_within_leq[OF \<open>p\<noteq>0\<close>,of s] \<open>proots_count p s=0\<close> by auto
   moreover have "finite (proots_within p s)" using \<open>p\<noteq>0\<close> by auto
   ultimately show ?thesis by auto
 qed
 
 lemma proots_count_leq_degree:
   assumes "p\<noteq>0"
   shows "proots_count p s\<le> degree p" using assms
 proof (induct rule:poly_root_induct[of _ "\<lambda>x. x\<in>s"])
   case 0
   then show ?case by auto
 next
   case (no_roots p)
   then have "proots_within p s = {}" by auto
   then show ?case unfolding proots_count_def by auto
 next
   case (root a p)
   have "proots_count ([:a, - 1:] * p) s = proots_count [:a, - 1:] s + proots_count p s"
     apply (subst proots_count_times)
     using root by auto
   also have "... = 1 + proots_count p s"
   proof -
     have "proots_count [:a, - 1:] s  =1"
       by (metis (no_types, lifting) add.inverse_inverse add.inverse_neutral minus_pCons 
           proots_count_pCons_1_iff proots_count_uminus root.hyps(1))
     then show ?thesis by auto
   qed
   also have "... \<le>  degree ([:a,-1:] * p)" 
     apply (subst degree_mult_eq)
     subgoal by auto
     subgoal using root by auto
     subgoal using root by (simp add: \<open>p \<noteq> 0\<close>)
     done
   finally show ?case .
 qed
 
 (*TODO: not a duplicate*)
 lemma proots_count_union_disjoint:
   assumes "A \<inter> B = {}" "p\<noteq>0"
   shows "proots_count p (A \<union> B) = proots_count p A + proots_count p B"
   unfolding proots_count_def 
   apply (subst proots_within_union[symmetric])
   apply (subst sum.union_disjoint)
   using assms by auto
 
 lemma proots_count_cong:
   assumes order_eq:"\<forall>x\<in>s. order x p = order x q" and "p\<noteq>0" and "q\<noteq>0"
   shows "proots_count p s = proots_count q s" unfolding proots_count_def
 proof (rule sum.cong)
   have "poly p x = 0 \<longleftrightarrow> poly q x = 0" when "x\<in>s" for x
     using order_eq that by (simp add: assms(2) assms(3) order_root)
   then show "proots_within p s = proots_within q s" by auto
   show "\<And>x. x \<in> proots_within q s \<Longrightarrow> order x p = order x q"
     using order_eq by auto
 qed
 
 lemma proots_count_of_real:
   assumes "p\<noteq>0"
   shows "proots_count (map_poly of_real p) ((of_real::_\<Rightarrow>'a::{real_algebra_1,idom}) ` s) 
             = proots_count p s"
 proof -
   define k where "k=(of_real::_\<Rightarrow>'a)"
   have "proots_within (map_poly of_real p) (k ` s) =k ` (proots_within p s)"
     unfolding proots_within_def k_def by (auto simp add:of_real_poly_map_poly[symmetric])
   then have "proots_count (map_poly of_real p) (k ` s) 
                 = (\<Sum>r\<in>k ` (proots_within p s). order r (map_poly of_real p))"
     unfolding proots_count_def by simp
   also have "... = sum ((\<lambda>r. order r (map_poly of_real p)) \<circ> k) (proots_within p s)"
     apply (subst sum.reindex) 
     unfolding k_def by (auto simp add: inj_on_def)
   also have "... = proots_count p s" unfolding proots_count_def
     apply (rule sum.cong)
     unfolding k_def comp_def using \<open>p\<noteq>0\<close> by (auto simp add:map_poly_order_of_real) 
   finally show ?thesis unfolding k_def .
 qed
 
 (*Is field really necessary here?*)
 lemma proots_pcompose:
   fixes p q::"'a::field poly"
   assumes "p\<noteq>0" "degree q=1"
   shows "proots_count (pcompose p q) s = proots_count p (poly q ` s)"
 proof -
   obtain a b where ab:"q=[:a,b:]" "b\<noteq>0"
     using \<open>degree q=1\<close> degree_eq_oneE by metis
   define f where "f=(\<lambda>y. (y-a)/b)"
   have f_eq:"f (poly q x) = x" "poly q (f x) = x" for x
     unfolding f_def using ab by auto
   have "proots_count (p \<circ>\<^sub>p q) s = (\<Sum>r\<in> f ` proots_within p (poly q ` s). order r (p \<circ>\<^sub>p q))" 
     unfolding proots_count_def
     apply (rule arg_cong2[where f =sum])
     apply (auto simp:poly_pcompose proots_within_def f_eq)
     by (metis (mono_tags, lifting) f_eq(1) image_eqI mem_Collect_eq)
   also have "... = (\<Sum>x\<in>proots_within p (poly q ` s). order (f x) (p \<circ>\<^sub>p q))"
     apply (subst sum.reindex)
     subgoal unfolding f_def inj_on_def using \<open>b\<noteq>0\<close> by auto
     by simp
   also have "... = (\<Sum>x\<in>proots_within p (poly q ` s). order x p)"
   proof -
     have "p \<circ>\<^sub>p q \<noteq> 0" using assms(1) assms(2) pcompose_eq_0 by force
     moreover have "order (f x) (q - [:x:]) = 1" for x 
     proof -
       have "order (f x) (q - [:x:]) = order (f x) (smult b [:-((x - a) / b),1:])"
         unfolding f_def using ab by auto
       also have "... = 1"
         apply (subst order_smult)
         using \<open>b\<noteq>0\<close> unfolding f_def by auto
       finally show ?thesis .
     qed
     ultimately have "order (f x) (p \<circ>\<^sub>p q) = order x p" for x
       apply (subst order_pcompose)
       using f_eq by auto
     then show ?thesis by auto
   qed
   also have "... =  proots_count p (poly q ` s)"
     unfolding proots_count_def by auto
   finally show ?thesis .
 qed
 
 subsection \<open>Composition of a polynomial and a rational function\<close>
 
 (*composition of a polynomial and a rational function. Maybe a more general version in the future?*)
 definition fcompose::"'a ::field poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly" where
   "fcompose p q r = fst (fold_coeffs (\<lambda>a (c,d). (d*[:a:] + q * c,r*d)) p (0,1))"
 
 lemma fcompose_0 [simp]: "fcompose 0 q r = 0"
   by (simp add: fcompose_def)
 
 lemma fcompose_const[simp]:"fcompose [:a:] q r = [:a:]"
   unfolding fcompose_def by (cases "a=0") auto
 
 lemma fcompose_pCons: 
   "fcompose (pCons a p) q1 q2 = smult a (q2^(degree (pCons a p))) + q1 * fcompose p q1 q2"
 proof (cases "p=0")
   case False
   define ff where "ff=(\<lambda>a (c, d). (d * [:a:] + q1 * c, q2 * d))"
   define fc where "fc=fold_coeffs ff p (0, 1)"
   have snd_ff:"snd fc = (if p=0 then 1 else q2^(degree p + 1))" unfolding fc_def
     apply (induct p)
     subgoal by simp
     subgoal for a p
       by (auto simp add:ff_def split:if_splits prod.splits)
     done
 
   have "fcompose (pCons a p) q1 q2 = fst (fold_coeffs ff (pCons a p) (0, 1))"
     unfolding fcompose_def ff_def by simp
   also have "... = fst (ff a fc)"
     using False unfolding fc_def by auto
   also have "... = snd fc * [:a:] + q1 * fst fc"
     unfolding ff_def by (auto split:prod.splits)
   also have "... = smult a (q2^(degree (pCons a p))) + q1 * fst fc"
     using snd_ff False by auto
   also have "... = smult a (q2^(degree (pCons a p))) + q1 * fcompose p q1 q2"
     unfolding fc_def ff_def fcompose_def by simp
   finally show ?thesis .
 qed simp
 
 lemma fcompose_uminus:
   "fcompose (-p) q r = - fcompose p q r"
   by (induct p) (auto simp:fcompose_pCons)
 
 lemma fcompose_add_less:
   assumes "degree p1 > degree p2"
   shows "fcompose (p1+p2) q1 q2 
             = fcompose p1 q1 q2 + q2^(degree p1-degree p2) * fcompose p2 q1 q2"
   using assms
 proof (induction p1 p2 rule: poly_induct2)
   case (pCons a1 p1 a2 p2)
   have ?case when "p2=0"
     using that by (simp add:fcompose_pCons smult_add_left)
   moreover have ?case when "p2\<noteq>0" "\<not> degree p2 < degree p1"
     using that pCons(2) by auto
   moreover have ?case when "p2\<noteq>0" "degree p2< degree p1"
   proof -
     define d1 d2 where "d1=degree (pCons a1 p1)" and "d2=degree (pCons a2 p2) "
     define fp1 fp2 where "fp1= fcompose p1 q1 q2" and "fp2=fcompose p2 q1 q2"
 
     have "fcompose (pCons a1 p1 + pCons a2 p2) q1 q2 
             = fcompose (pCons (a1+a2) (p1+p2)) q1 q2"
       by simp
     also have "... = smult (a1 + a2) (q2 ^ d1) + q1 * fcompose (p1 + p2) q1 q2"
     proof -
       have "degree (pCons (a1 + a2) (p1 + p2)) = d1"
         unfolding d1_def using that degree_add_eq_left by fastforce
       then show ?thesis unfolding fcompose_pCons by simp
     qed
     also have "... = smult (a1 + a2) (q2 ^ d1) + q1 * (fp1 + q2 ^ (d1 - d2) * fp2)"
     proof -
       have "degree p1 - degree p2 = d1 - d2"
         unfolding d1_def d2_def using that by simp
       then show ?thesis 
         unfolding pCons(1)[OF that(2),folded fp1_def fp2_def] by simp
     qed
     also have "... = fcompose (pCons a1 p1) q1 q2 + q2 ^ (d1 - d2) 
                         * fcompose (pCons a2 p2) q1 q2"
     proof -
       have "d1 > d2" unfolding d1_def d2_def using that by auto
       then show ?thesis
         unfolding fcompose_pCons
         apply (fold d1_def d2_def fp1_def fp2_def)
         by (simp add:algebra_simps smult_add_left power_add[symmetric])
     qed
     finally show ?thesis unfolding d1_def d2_def .
   qed
   ultimately show ?case by blast
 qed simp
 
 lemma fcompose_add_eq:
   assumes "degree p1 = degree p2"
   shows "q2^(degree p1 - degree (p1+p2)) * fcompose (p1+p2) q1 q2 
             = fcompose p1 q1 q2 + fcompose p2 q1 q2"
   using assms
 proof (induction p1 p2 rule: poly_induct2)
   case (pCons a1 p1 a2 p2)
   have ?case when "p1+p2=0"
   proof -
     have "p2=-p1" using that by algebra
     then show ?thesis by (simp add:fcompose_pCons fcompose_uminus smult_add_left)
   qed
   moreover have ?case when "p1=0"
   proof -
     have "p2=0"
       using pCons(2) that by (auto split:if_splits)
     then show ?thesis using that by simp
   qed
   moreover have ?case when "p1\<noteq>0" "p1+p2\<noteq>0"
   proof -
     define d1 d2 dp where "d1=degree (pCons a1 p1)" and "d2=degree (pCons a2 p2)"
                             and "dp = degree p1 - degree (p1+p2)"
     define fp1 fp2 where "fp1= fcompose p1 q1 q2" and "fp2=fcompose p2 q1 q2"
 
     have "q2 ^ (degree (pCons a1 p1) - degree (pCons a1 p1 + pCons a2 p2)) *
              fcompose (pCons a1 p1 + pCons a2 p2) q1 q2 
                 = q2 ^ dp * fcompose (pCons (a1+a2) (p1 +p2)) q1 q2"
       unfolding dp_def using that by auto
     also have "... = smult (a1 + a2) (q2 ^ d1) + q1 * (q2 ^ dp * fcompose (p1 + p2) q1 q2)"
     proof -
       have "degree p1 \<ge> degree (p1 + p2)" 
         by (metis degree_add_le degree_pCons_eq_if not_less_eq_eq order_refl pCons.prems zero_le)
       then show ?thesis 
         unfolding fcompose_pCons dp_def d1_def using that
         by (simp add:algebra_simps power_add[symmetric])
     qed
     also have "... =  smult (a1 + a2) (q2 ^ d1) + q1 * (fp1 + fp2)"
       apply (subst pCons(1)[folded dp_def fp1_def fp2_def])
       subgoal by (metis degree_pCons_eq_if diff_Suc_Suc diff_zero not_less_eq_eq pCons.prems zero_le)
       subgoal by simp
       done
     also have "... = fcompose (pCons a1 p1) q1 q2 + fcompose (pCons a2 p2) q1 q2"
     proof -
       have *:"d1 = degree (pCons a2 p2)"
         unfolding d1_def using pCons(2) by simp
       show ?thesis 
         unfolding fcompose_pCons
         apply (fold d1_def fp1_def fp2_def *)
         by (simp add:smult_add_left algebra_simps)
     qed
     finally show ?thesis .
   qed
   ultimately show ?case by blast
 qed simp
 
 lemma fcompose_add_const:
   "fcompose ([:a:] + p) q1 q2 = smult a (q2 ^ degree p) + fcompose p q1 q2"
   apply (cases p)
   by (auto simp add:fcompose_pCons smult_add_left)
 
 lemma fcompose_smult: "fcompose (smult a p) q1 q2 = smult a (fcompose p q1 q2)"
   by (induct p) (simp_all add:fcompose_pCons smult_add_right)
 
 lemma fcompose_mult: "fcompose (p1*p2) q1 q2 = fcompose p1 q1 q2 * fcompose p2 q1 q2"
 proof (induct p1)
   case 0
   then show ?case by simp
 next
   case (pCons a p1)
   have ?case when "p1=0 \<or> p2=0" 
     using that by (auto simp add:fcompose_smult)
   moreover have ?case when "p1\<noteq>0" "p2\<noteq>0" "a=0"
     using that by (simp add:fcompose_pCons pCons)
   moreover have ?case when "p1\<noteq>0" "p2\<noteq>0" "a\<noteq>0"
   proof -
     have "fcompose (pCons a p1 * p2) q1 q2 
             = fcompose (pCons 0 (p1 * p2) + smult a p2) q1 q2"
       by (simp add:algebra_simps)
     also have "... =  fcompose (pCons 0 (p1 * p2)) q1 q2 
                         + q2 ^ (degree p1 +1) * fcompose (smult a p2) q1 q2"
     proof -
       have "degree (pCons 0 (p1 * p2)) > degree (smult a p2)"
         using that by (simp add: degree_mult_eq)
       from fcompose_add_less[OF this,of q1 q2] that 
       show ?thesis by (simp add:degree_mult_eq)
     qed
     also have "... = fcompose (pCons a p1) q1 q2 * fcompose p2 q1 q2"
       using that by (simp add:fcompose_pCons fcompose_smult pCons algebra_simps)
     finally show ?thesis .
   qed
   ultimately show ?case by blast
 qed
 
 lemma fcompose_poly:
   assumes "poly q2 x\<noteq>0"
   shows "poly p (poly q1 x/poly q2 x) = poly (fcompose p q1 q2) x / poly (q2^(degree p)) x"
   apply (induct p)
   using assms by (simp_all add:fcompose_pCons field_simps)
 
 lemma poly_fcompose:
    assumes "poly q2 x\<noteq>0"
    shows "poly (fcompose p q1 q2) x = poly p (poly q1 x/poly q2 x) * (poly q2 x)^(degree p)"
   using fcompose_poly[OF assms] assms by (auto simp add:field_simps)
 lemma poly_fcompose_0_denominator:
   assumes "poly q2 x=0"
   shows "poly (fcompose p q1 q2) x = poly q1 x ^ degree p * lead_coeff p"
   apply (induct p)
   using assms by (auto simp add:fcompose_pCons)
 
 lemma fcompose_0_denominator:"fcompose p q1 0 = smult (lead_coeff p) (q1^degree p)"
   apply (induct p)
   by (auto simp:fcompose_pCons)
 
 lemma fcompose_nzero:
   fixes p::"'a::field poly"
   assumes "p\<noteq>0" and "q2\<noteq>0" and nconst:"\<forall>c. q1 \<noteq> smult c q2"
       and infi:"infinite (UNIV::'a set)"
   shows "fcompose p q1 q2 \<noteq> 0" using \<open>p\<noteq>0\<close>
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by simp
 next
   case (no_proots p)
   have False when "fcompose p q1 q2 = 0"
   proof -
     obtain x where "poly q2 x\<noteq>0" 
     proof -
       have "finite (proots q2)" using \<open>q2\<noteq>0\<close> by auto
       then have "\<exists>x. poly q2 x\<noteq>0" 
         by (meson UNIV_I ex_new_if_finite infi proots_withinI)
       then show ?thesis using that by auto
     qed
     define y where "y = poly q1 x / poly q2 x"
     have "poly p y = 0"
       using \<open>fcompose p q1 q2 = 0\<close> fcompose_poly[OF \<open>poly q2 x\<noteq>0\<close>,of p q1,folded y_def] 
       by simp
     then show False using no_proots(1) by auto
   qed
   then show ?case by auto
 next
   case (root a p)
   have "fcompose [:- a, 1:] q1 q2 \<noteq> 0" 
     unfolding fcompose_def using nconst[rule_format,of a]
     by simp
   moreover have "fcompose p q1 q2 \<noteq> 0" 
     using root by fastforce
   ultimately show ?case unfolding fcompose_mult by auto
 qed
 
 subsection \<open>Bijection (@{term bij_betw}) and the number of polynomial roots\<close>
 
 lemma proots_fcompose_bij_eq:
   fixes p::"'a::field poly"
   assumes bij:"bij_betw (\<lambda>x. poly q1 x/poly q2 x) A B" and "p\<noteq>0" 
       and nzero:"\<forall>x\<in>A. poly q2 x\<noteq>0"
       and max_deg: "max (degree q1) (degree q2) \<le> 1"
       and nconst:"\<forall>c. q1 \<noteq> smult c q2"
       and infi:"infinite (UNIV::'a set)"
   shows "proots_count p B = proots_count (fcompose p q1 q2) A"
   using \<open>p\<noteq>0\<close>
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by simp
 next
   case (no_proots p)
   have "proots_count p B = 0"
   proof -
     have "proots_within p B = {}"
       using no_proots by auto
     then show ?thesis unfolding proots_count_def by auto
   qed
   moreover have "proots_count (fcompose p q1 q2) A = 0"
   proof -
     have "proots_within (fcompose p q1 q2) A = {}"
       using no_proots unfolding proots_within_def
-      by (smt div_0 empty_Collect_eq fcompose_poly nzero)
+      by (smt (verit) div_0 empty_Collect_eq fcompose_poly nzero)
     then show ?thesis unfolding proots_count_def by auto
   qed
   ultimately show ?case by auto
 next
   case (root b p)
   have "proots_count ([:- b, 1:] * p) B = proots_count [:- b, 1:] B + proots_count p B"
     using proots_count_times[OF \<open>[:- b, 1:] * p \<noteq> 0\<close>] by simp
   also have "... = proots_count (fcompose [:- b, 1:] q1 q2) A 
                     + proots_count (fcompose p q1 q2) A" 
   proof -
     define g where "g=(\<lambda>x. poly q1 x/poly q2 x)"
 
     have "proots_count [:- b, 1:] B = proots_count (fcompose [:- b, 1:] q1 q2) A" 
     proof (cases "b\<in>B")
       case True
       then have "proots_count [:- b, 1:] B = 1"
         unfolding proots_count_pCons_1_iff by simp
       moreover have "proots_count (fcompose [:- b, 1:] q1 q2) A = 1"
       proof -
         obtain a where "b=g a" "a\<in>A"
           using bij[folded g_def] True 
           by (metis bij_betwE bij_betw_the_inv_into f_the_inv_into_f_bij_betw)
         define qq where "qq=q1 - smult b q2"
         have qq_0:"poly qq a=0" and qq_deg: "degree qq\<le>1" and \<open>qq\<noteq>0\<close>
           unfolding qq_def
           subgoal using \<open>b=g a\<close> nzero[rule_format,OF \<open>a\<in>A\<close>] unfolding g_def by auto
           subgoal using max_deg by (simp add: degree_diff_le)
           subgoal using nconst[rule_format,of b] by auto
           done
         have "proots_within qq A = {a}" 
         proof -
           have "a\<in>proots_within qq A" 
             using qq_0 \<open>a\<in>A\<close> by auto
           moreover have "card (proots_within qq A) = 1" 
           proof -
             have "finite (proots_within qq A)" using \<open>qq\<noteq>0\<close> by simp
             moreover have "proots_within qq A \<noteq> {}"
               using \<open>a\<in>proots_within qq A\<close> by auto
             ultimately have "card (proots_within qq A) \<noteq>0" by auto
             moreover have "card (proots_within qq A) \<le> 1" 
               by (meson \<open>qq \<noteq> 0\<close> card_proots_within_leq le_trans proots_count_leq_degree qq_deg)
             ultimately show ?thesis by auto
           qed
           ultimately show ?thesis by (metis card_1_singletonE singletonD)
         qed
         moreover have "order a qq=1" 
           by (metis One_nat_def \<open>qq \<noteq> 0\<close> le_antisym le_zero_eq not_less_eq_eq order_degree
                 order_root qq_0 qq_deg)
         ultimately show ?thesis unfolding fcompose_def proots_count_def qq_def  
           by auto
       qed
       ultimately show ?thesis by auto
     next
       case False
       then have "proots_count [:- b, 1:] B  = 0"
         unfolding proots_count_pCons_1_iff by simp
       moreover have "proots_count (fcompose [:- b, 1:] q1 q2) A = 0"
       proof -
         have "proots_within (fcompose [:- b, 1:] q1 q2) A = {}"
         proof (rule ccontr)
           assume "proots_within (fcompose [:- b, 1:] q1 q2) A \<noteq> {}"
           then obtain a where "a\<in>A" "poly q1 a = b * poly q2 a"
             unfolding fcompose_def proots_within_def by auto
           then have "b = g a"
             unfolding g_def using nzero[rule_format,OF \<open>a\<in>A\<close>] by auto
           then have "b\<in>B" using \<open>a\<in>A\<close> bij[folded g_def] using bij_betwE by blast
           then show False using False by auto
         qed  
         then show ?thesis unfolding proots_count_def by auto
       qed      
       ultimately show ?thesis by simp
     qed
     moreover have "proots_count p B = proots_count (fcompose p q1 q2) A" 
       apply (rule root.hyps)
       using mult_eq_0_iff root.prems by blast
     ultimately show ?thesis by auto
   qed
   also have "... = proots_count (fcompose ([:- b, 1:] * p) q1 q2) A"
   proof (cases "A={}")
     case False
     have "fcompose [:- b, 1:] q1 q2 \<noteq>0" 
       using nconst[rule_format,of b] unfolding fcompose_def by auto
     moreover have "fcompose p q1 q2 \<noteq> 0" 
       apply (rule fcompose_nzero[OF _ _ nconst infi])
       subgoal using \<open>[:- b, 1:] * p \<noteq> 0\<close> by auto
       subgoal using nzero False by auto
       done
     ultimately show ?thesis unfolding fcompose_mult
       apply (subst proots_count_times)
       by auto
   qed auto
   finally show ?case .
 qed
 
 lemma proots_card_fcompose_bij_eq:
   fixes p::"'a::field poly"
   assumes bij:"bij_betw (\<lambda>x. poly q1 x/poly q2 x) A B" and "p\<noteq>0" 
       and nzero:"\<forall>x\<in>A. poly q2 x\<noteq>0"
       and max_deg: "max (degree q1) (degree q2) \<le> 1"
       and nconst:"\<forall>c. q1 \<noteq> smult c q2"
       and infi:"infinite (UNIV::'a set)"
   shows "card (proots_within p B) = card (proots_within (fcompose p q1 q2) A)"
   using \<open>p\<noteq>0\<close>
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by simp
 next
   case (no_proots p)
   have "proots_within p B = {}" using no_proots by auto
   moreover have "proots_within (fcompose p q1 q2) A = {}" 
     using no_proots fcompose_poly 
-    by (smt Collect_empty_eq divide_eq_0_iff nzero proots_within_def)
+    by (smt (verit) Collect_empty_eq divide_eq_0_iff nzero proots_within_def)
   ultimately show ?case by auto
 next
   case (root b p)
   then have [simp]:"p\<noteq>0" by auto
 
   have ?case when "b\<notin>B \<or> poly p b=0"
   proof -
     have "proots_within ([:- b, 1:] * p) B = proots_within p B"
       using that by auto
     moreover have "proots_within (fcompose ([:- b, 1:] * p) q1 q2) A 
         = proots_within (fcompose p q1 q2) A"
       using that nzero unfolding fcompose_mult proots_within_times 
       apply (auto simp add: poly_fcompose)
       using bij bij_betwE by blast
     ultimately show ?thesis using root by auto
   qed
   moreover have ?case when "b\<in>B" "poly p b\<noteq>0"
   proof -
     define bb where "bb=[:- b, 1:]"
     have "card (proots_within (bb * p) B) = card {b} + card (proots_within p B)"
     proof -
       have "proots_within bb B = {b}" 
         using that unfolding bb_def by auto
       then show ?thesis unfolding proots_within_times
         apply (subst card_Un_disjoint)
         by (use that in auto)
     qed
     also have "... = 1 + card (proots_within (fcompose p q1 q2) A)"
       using root.hyps by simp
     also have "... = card (proots_within (fcompose (bb * p) q1 q2) A)"
       unfolding proots_within_times fcompose_mult
     proof (subst card_Un_disjoint) 
       obtain a where b_poly:"b=poly q1 a / poly q2 a" and "a\<in>A" 
         by (metis (no_types, lifting) \<open>b \<in> B\<close> bij bij_betwE bij_betw_the_inv_into 
             f_the_inv_into_f_bij_betw)
       define bbq pq where "bbq=fcompose bb q1 q2" and "pq=fcompose p q1 q2"
       have bbq_0:"poly bbq a=0" and bbq_deg: "degree bbq\<le>1" and "bbq\<noteq>0"
         unfolding bbq_def bb_def 
         subgoal using \<open>a \<in> A\<close> b_poly nzero poly_fcompose by fastforce
         subgoal by (metis (no_types, lifting) degree_add_le degree_pCons_eq_if degree_smult_le 
            dual_order.trans fcompose_const fcompose_pCons max.boundedE max_deg mult_cancel_left2 
            one_neq_zero one_poly_eq_simps(1) power.simps)
         subgoal by (metis \<open>a \<in> A\<close> \<open>poly (fcompose [:- b, 1:] q1 q2) a = 0\<close> fcompose_nzero infi 
                   nconst nzero one_neq_zero pCons_eq_0_iff)
         done
       show "finite (proots_within bbq A)" using \<open>bbq\<noteq>0\<close> by simp
       show "finite (proots_within pq A)" unfolding pq_def 
         by (metis \<open>a \<in> A\<close> \<open>p \<noteq> 0\<close> fcompose_nzero finite_proots infi nconst nzero poly_0 pq_def)
       have bbq_a:"proots_within bbq A = {a}"
       proof -
         have "a\<in>proots_within bbq A" 
           by (simp add: \<open>a \<in> A\<close> bbq_0)
         moreover have "card (proots_within bbq A) = 1" 
         proof -
           have "card (proots_within bbq A) \<noteq>0" 
             using \<open>a\<in>proots_within bbq A\<close> \<open>finite (proots_within bbq A)\<close>
             by auto
           moreover have "card (proots_within bbq A) \<le> 1" 
             by (meson \<open>bbq \<noteq> 0\<close> card_proots_within_leq le_trans proots_count_leq_degree bbq_deg)
           ultimately show ?thesis by auto
         qed
         ultimately show ?thesis by (metis card_1_singletonE singletonD)
       qed
       show "proots_within (bbq) A \<inter> proots_within (pq) A = {}" 
         using b_poly bbq_a fcompose_poly nzero pq_def that(2) by fastforce
       show "1 + card (proots_within pq A) = card (proots_within bbq A) + card (proots_within pq A)"
         using bbq_a by simp
     qed
     finally show ?thesis unfolding bb_def .
   qed
   ultimately show ?case by auto
 qed 
 
 lemma proots_pcompose_bij_eq:
   fixes p::"'a::idom poly"
   assumes bij:"bij_betw (\<lambda>x. poly q x) A B" and "p\<noteq>0" 
       and q_deg: "degree q = 1"
   shows "proots_count p B = proots_count (p \<circ>\<^sub>p q) A" using \<open>p\<noteq>0\<close>
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by simp
 next
   case (no_proots p)
   have "proots_count p B = 0"
   proof -
     have "proots_within p B = {}"
       using no_proots by auto
     then show ?thesis unfolding proots_count_def by auto
   qed
   moreover have "proots_count (p \<circ>\<^sub>p q) A = 0"
   proof -
     have "proots_within (p \<circ>\<^sub>p q) A = {}"
       using no_proots unfolding proots_within_def
       by (auto simp:poly_pcompose)
     then show ?thesis unfolding proots_count_def by auto
   qed
   ultimately show ?case by auto
 next
   case (root b p)
   have "proots_count ([:- b, 1:] * p) B = proots_count [:- b, 1:] B + proots_count p B"
     using proots_count_times[OF \<open>[:- b, 1:] * p \<noteq> 0\<close>] by simp
   also have "... = proots_count ([:- b, 1:] \<circ>\<^sub>p q) A + proots_count (p \<circ>\<^sub>p q) A"
   proof -
     have "proots_count [:- b, 1:] B = proots_count ([:- b, 1:] \<circ>\<^sub>p q) A"
     proof (cases "b\<in>B")
       case True
       then have "proots_count [:- b, 1:] B = 1" 
         unfolding proots_count_pCons_1_iff by simp
       moreover have "proots_count ([:- b, 1:] \<circ>\<^sub>p q) A = 1" 
       proof -
         obtain a where "b=poly q a" "a\<in>A"
           using True bij by (metis bij_betwE bij_betw_the_inv_into f_the_inv_into_f_bij_betw)
         define qq where "qq=[:- b:] + q"
         have qq_0:"poly qq a=0" and qq_deg: "degree qq\<le>1" and \<open>qq\<noteq>0\<close>
           unfolding qq_def
           subgoal using \<open>b=poly q a\<close> by auto
           subgoal using q_deg by (simp add: degree_add_le)
           subgoal using q_deg add.inverse_unique by force
           done
         have "proots_within qq A = {a}" 
         proof -
           have "a\<in>proots_within qq A" 
             using qq_0 \<open>a\<in>A\<close> by auto
           moreover have "card (proots_within qq A) = 1" 
           proof -
             have "finite (proots_within qq A)" using \<open>qq\<noteq>0\<close> by simp
             moreover have "proots_within qq A \<noteq> {}"
               using \<open>a\<in>proots_within qq A\<close> by auto
             ultimately have "card (proots_within qq A) \<noteq>0" by auto
             moreover have "card (proots_within qq A) \<le> 1" 
               by (meson \<open>qq \<noteq> 0\<close> card_proots_within_leq le_trans proots_count_leq_degree qq_deg)
             ultimately show ?thesis by auto
           qed
           ultimately show ?thesis by (metis card_1_singletonE singletonD)
         qed
         moreover have "order a qq=1" 
           by (metis One_nat_def \<open>qq \<noteq> 0\<close> le_antisym le_zero_eq not_less_eq_eq order_degree
                 order_root qq_0 qq_deg)
         ultimately show ?thesis unfolding pcompose_def proots_count_def qq_def  
           by auto
       qed
       ultimately show ?thesis by auto
     next
       case False
       then have "proots_count [:- b, 1:] B  = 0"
         unfolding proots_count_pCons_1_iff by simp
       moreover have "proots_count ([:- b, 1:] \<circ>\<^sub>p q) A = 0"
       proof -
         have "proots_within ([:- b, 1:] \<circ>\<^sub>p q) A = {}"
           unfolding pcompose_def 
           apply auto
           using False bij bij_betwE by blast
         then show ?thesis unfolding proots_count_def by auto
       qed      
       ultimately show ?thesis by simp
     qed
     moreover have "proots_count p B = proots_count (p \<circ>\<^sub>p q) A"
       apply (rule root.hyps)
       using \<open>[:- b, 1:] * p \<noteq> 0\<close> by auto
     ultimately show ?thesis by auto
   qed
   also have "... = proots_count (([:- b, 1:] * p) \<circ>\<^sub>p q) A"
     unfolding pcompose_mult
     apply (subst proots_count_times)
     subgoal by (metis (no_types, lifting) One_nat_def add.right_neutral degree_0 degree_mult_eq
       degree_pCons_eq_if degree_pcompose mult_eq_0_iff one_neq_zero one_pCons pcompose_mult
       q_deg root.prems)
     by simp
   finally show ?case .
 qed
 
 lemma proots_card_pcompose_bij_eq:
   fixes p::"'a::idom poly"
   assumes bij:"bij_betw (\<lambda>x. poly q x) A B" and "p\<noteq>0" 
       and q_deg: "degree q = 1"
   shows "card (proots_within p B) = card (proots_within (p \<circ>\<^sub>p q) A)" using \<open>p\<noteq>0\<close>
 proof (induct p rule:poly_root_induct_alt)
   case 0
   then show ?case by auto
 next
   case (no_proots p)
   have "proots_within p B = {}" using no_proots by auto
   moreover have "proots_within (p \<circ>\<^sub>p q) A = {}" using no_proots 
     by (simp add: poly_pcompose proots_within_def)
   ultimately show ?case by auto
 next
   case (root b p)
   then have [simp]:"p\<noteq>0" by auto
   have ?case when "b\<notin>B \<or> poly p b=0"
   proof -
     have "proots_within ([:- b, 1:] * p) B = proots_within p B"
       using that by auto
     moreover have "proots_within (([:- b, 1:] * p) \<circ>\<^sub>p q) A = proots_within (p \<circ>\<^sub>p q) A"
       using that unfolding pcompose_mult proots_within_times
       apply (auto simp add: poly_pcompose)
       using bij bij_betwE by blast
     ultimately show ?thesis using root.hyps[OF \<open>p\<noteq>0\<close>] by auto
   qed
   moreover have ?case when "b\<in>B" "poly p b\<noteq>0"
   proof -
     define bb where "bb=[:- b, 1:]"
     have "card (proots_within (bb * p) B) = card {b} + card (proots_within p B)"
     proof -
       have "proots_within bb B = {b}" 
         using that unfolding bb_def by auto
       then show ?thesis unfolding proots_within_times
         apply (subst card_Un_disjoint)
         by (use that in auto)
     qed
     also have "... = 1 + card (proots_within (p \<circ>\<^sub>p q) A)"
       using root.hyps by simp
     also have "... = card (proots_within ((bb * p) \<circ>\<^sub>p q) A)"
       unfolding proots_within_times pcompose_mult
     proof (subst card_Un_disjoint) 
       obtain a where "b=poly q a" "a\<in>A" 
         by (metis \<open>b \<in> B\<close> bij bij_betwE bij_betw_the_inv_into f_the_inv_into_f_bij_betw)
       define bbq pq where "bbq=bb \<circ>\<^sub>p q" and "pq=p \<circ>\<^sub>p q"
       have bbq_0:"poly bbq a=0" and bbq_deg: "degree bbq\<le>1" and "bbq\<noteq>0"
         unfolding bbq_def bb_def poly_pcompose
         subgoal using \<open>b=poly q a\<close> by auto
         subgoal using q_deg by (simp add: degree_add_le degree_pcompose)
         subgoal using pcompose_eq_0 q_deg by fastforce
         done
       show "finite (proots_within bbq A)" using \<open>bbq\<noteq>0\<close> by simp
       show "finite (proots_within pq A)" unfolding pq_def 
         by (metis \<open>p \<noteq> 0\<close> finite_proots pcompose_eq_0 q_deg zero_less_one) 
       have bbq_a:"proots_within bbq A = {a}"
       proof -
         have "a\<in>proots_within bbq A" 
           unfolding bb_def proots_within_def poly_pcompose bbq_def 
           using \<open>b=poly q a\<close> \<open>a\<in>A\<close> by simp
         moreover have "card (proots_within bbq A) = 1" 
         proof -
           have "card (proots_within bbq A) \<noteq>0" 
             using \<open>a\<in>proots_within bbq A\<close> \<open>finite (proots_within bbq A)\<close>
             by auto
           moreover have "card (proots_within bbq A) \<le> 1" 
             by (meson \<open>bbq \<noteq> 0\<close> card_proots_within_leq le_trans proots_count_leq_degree bbq_deg)
           ultimately show ?thesis by auto
         qed
         ultimately show ?thesis by (metis card_1_singletonE singletonD)
       qed
       show "proots_within (bbq) A \<inter> proots_within (pq) A = {}" 
         using bbq_a \<open>b = poly q a\<close> that(2) unfolding pq_def by (simp add:poly_pcompose)
       show "1 + card (proots_within pq A) = card (proots_within bbq A) + card (proots_within pq A)"
         using bbq_a by simp
     qed
     finally show ?thesis unfolding bb_def .
   qed
   ultimately show ?case by auto
 qed
 
 end
\ No newline at end of file
diff --git a/thys/Three_Squares/Quadratic_Forms.thy b/thys/Three_Squares/Quadratic_Forms.thy
--- a/thys/Three_Squares/Quadratic_Forms.thy
+++ b/thys/Three_Squares/Quadratic_Forms.thy
@@ -1,977 +1,978 @@
 (*  Title:      Three_Squares/Quadratic_Forms.thy
     Author:     Anton Danilkin
 *)
 
 section \<open>Properties of quadratic forms and their equivalences\<close>
 
 theory Quadratic_Forms
   imports Complex_Main Low_Dimensional_Linear_Algebra
 begin
 
 consts qf_app :: "'a \<Rightarrow> 'b \<Rightarrow> int" (infixl "$$" 65)
 
 definition qf2_app :: "mat2 \<Rightarrow> vec2 \<Rightarrow> int" where
 "qf2_app m v = <v | m $ v>"
 
 adhoc_overloading qf_app qf2_app
 
 definition qf3_app :: "mat3 \<Rightarrow> vec3 \<Rightarrow> 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 \<Rightarrow> bool"
 
 definition qf2_positive_definite :: "mat2 \<Rightarrow> bool" where
 "qf2_positive_definite m = (\<forall>v. v \<noteq> 0 \<longrightarrow> m $$ v > 0)"
 
 adhoc_overloading qf_positive_definite qf2_positive_definite
 
 definition qf3_positive_definite :: "mat3 \<Rightarrow> bool" where
 "qf3_positive_definite m = (\<forall>v. v \<noteq> 0 \<longrightarrow> 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 "\<forall>v. m $$ v \<ge> 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 "\<forall>v. m $$ v \<ge> 0"
   using assms unfolding qf3_positive_definite_def
   by (metis order_less_le order_refl qf3_app_zero)
 
 consts qf_action :: "'a \<Rightarrow> 'a \<Rightarrow> 'a" (infixl "\<cdot>" 55)
 
 definition qf2_action :: "mat2 \<Rightarrow> mat2 \<Rightarrow> mat2" where
 "qf2_action a u = u\<^sup>T * a * u"
 
 adhoc_overloading qf_action qf2_action
 
 definition qf3_action :: "mat3 \<Rightarrow> mat3 \<Rightarrow> mat3" where
 "qf3_action a u = u\<^sup>T * a * u"
 
 adhoc_overloading qf_action qf3_action
 
 lemma qf2_action_id:
   fixes a :: mat2
   shows "a \<cdot> 1 = a"
   unfolding qf2_action_def
   by simp
 
 lemma qf3_action_id:
   fixes a :: mat3
   shows "a \<cdot> 1 = a"
   unfolding qf3_action_def
   by simp
 
 lemma qf2_action_mul [simp]:
   fixes a u v :: mat2
   shows "a \<cdot> (u * v) = (a \<cdot> u) \<cdot> v"
   unfolding qf2_action_def
   by (simp add: algebra_simps)
 
 lemma qf3_action_mul [simp]:
   fixes a u v :: mat3
   shows "a \<cdot> (u * v) = (a \<cdot> u) \<cdot> v"
   unfolding qf3_action_def
   by (simp add: algebra_simps)
 
 consts qf_equiv :: "'a \<Rightarrow> 'a \<Rightarrow> bool" (infix "~" 65)
 
 definition qf2_equiv :: "mat2 \<Rightarrow> mat2 \<Rightarrow> bool" where
 "qf2_equiv a b = (\<exists>u. mat_det u = 1 \<and> a \<cdot> u = b)"
 
 adhoc_overloading qf_equiv qf2_equiv
 
 definition qf3_equiv :: "mat3 \<Rightarrow> mat3 \<Rightarrow> bool" where
 "qf3_equiv a b = (\<exists>u. mat_det u = 1 \<and> a \<cdot> u = b)"
 
 adhoc_overloading qf_equiv qf3_equiv
 
 lemma qf2_equiv_sym_impl:
   fixes a b :: mat2
   shows "a ~ b \<Longrightarrow> b ~ a"
 unfolding qf2_equiv_def qf2_action_def
 proof -
   assume "\<exists>u. mat_det u = 1 \<and> u\<^sup>T * a * u = b"
   then obtain u where "mat_det u = 1 \<and> u\<^sup>T * a * u = b" by blast
   hence "mat_det (u\<^sup>-\<^sup>1) = 1 \<and> ((u\<^sup>-\<^sup>1)\<^sup>T) * b * (u\<^sup>-\<^sup>1) = a"
     unfolding mat2_inverse_transpose[symmetric]
     using mat2_inverse_cancel_left[of "u\<^sup>T"] mat2_inverse_cancel_right
     by (auto simp add: algebra_simps)
   thus "\<exists>u. mat_det u = 1 \<and> u\<^sup>T * b * u = a" by blast
 qed
 
 lemma qf3_equiv_sym_impl:
   fixes a b :: mat3
   shows "a ~ b \<Longrightarrow> b ~ a"
 unfolding qf3_equiv_def qf3_action_def
 proof -
   assume "\<exists>u. mat_det u = 1 \<and> u\<^sup>T * a * u = b"
   then obtain u where "mat_det u = 1 \<and> u\<^sup>T * a * u = b" by blast
   hence "mat_det (u\<^sup>-\<^sup>1) = 1 \<and> ((u\<^sup>-\<^sup>1)\<^sup>T) * b * (u\<^sup>-\<^sup>1) = a"
     unfolding mat3_inverse_transpose[symmetric]
     using mat3_inverse_cancel_left[of "u\<^sup>T"] mat3_inverse_cancel_right
     by (auto simp add: algebra_simps)
   thus "\<exists>u. mat_det u = 1 \<and> u\<^sup>T * b * u = a" by blast
 qed
 
 lemma qf2_equiv_sym:
   fixes a b :: mat2
   shows "a ~ b \<longleftrightarrow> b ~ a"
   using qf2_equiv_sym_impl by blast
 
 lemma qf3_equiv_sym:
   fixes a b :: mat3
   shows "a ~ b \<longleftrightarrow> 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 \<cdot> 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 \<cdot> 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 \<longleftrightarrow> 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 \<longleftrightarrow> 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 \<longleftrightarrow> mat2_sym b"
 proof -
   obtain u where "mat_det u = 1" "u\<^sup>T * 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 "u\<^sup>T" "a\<^sup>T * u" "a * u"]
           mat2_inversable_cancel_right[of u "a\<^sup>T" a]
     by (auto simp add: algebra_simps)
 qed
 
 lemma qf3_equiv_preserves_sym:
   fixes a b :: mat3
   assumes "a ~ b"
   shows "mat3_sym a \<longleftrightarrow> mat3_sym b"
 proof -
   obtain u where "mat_det u = 1" "u\<^sup>T * 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 "u\<^sup>T" "a\<^sup>T * u" "a * u"]
           mat3_inversable_cancel_right[of u "a\<^sup>T" 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) \<subseteq> range (($$) a)"
 proof -
   obtain u where 0: "mat_det u = 1" "a \<cdot> 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) \<subseteq> range (($$) a)"
 proof -
   obtain u where 0: "mat_det u = 1" "a \<cdot> 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 \<open>Lemma 1.1 from @{cite nathanson1996}.\<close>
 
 lemma qf2_positive_definite_criterion:
   fixes a
   assumes "mat_sym a"
   shows "qf_positive_definite a \<longleftrightarrow> mat2\<^sub>1\<^sub>1 a > 0 \<and> mat_det a > 0"
 proof
   assume 0: "qf_positive_definite a"
   have "vec2 1 0 \<noteq> 0 \<longrightarrow> a $$ vec2 1 0 > 0" using 0
     unfolding qf2_positive_definite_def by blast
   hence 1: "mat2\<^sub>1\<^sub>1 a > 0"
     unfolding zero_vec2_def vec2_dot_def times_mat2_def mat2_app_def qf2_app_def
     by auto
   let ?v = "vec2 (- mat2\<^sub>1\<^sub>2 a) (mat2\<^sub>1\<^sub>1 a)"
   have "?v \<noteq> 0 \<longrightarrow> a $$ ?v > 0" using 0
     unfolding qf2_positive_definite_def by blast
   hence "mat2\<^sub>1\<^sub>1 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 "mat2\<^sub>1\<^sub>1 a > 0 \<and> mat_det a > 0" using 1 2 by blast
 next
   assume 3: "mat2\<^sub>1\<^sub>1 a > 0 \<and> mat_det a > 0"
   show "qf_positive_definite a" unfolding qf2_positive_definite_def
   proof (rule allI; rule impI)
     fix v :: vec2
     assume "v \<noteq> 0"
     hence 4: "vec2\<^sub>1 v \<noteq> 0 \<or> vec2\<^sub>2 v \<noteq> 0" unfolding zero_vec2_def
       by (metis vec2.collapse)
     let ?n = "(mat2\<^sub>1\<^sub>1 a * vec2\<^sub>1 v + mat2\<^sub>1\<^sub>2 a * vec2\<^sub>2 v)\<^sup>2 + (mat_det a) * (vec2\<^sub>2 v)\<^sup>2"
     have 5: "mat2\<^sub>1\<^sub>1 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 "\<exists>t. 2 * \<bar>k + m * t\<bar> \<le> m" (is "\<exists> t. ?P t")
 proof -
   let ?s = "k div m"
   have 0: "k - m * ?s \<ge> 0 \<and> 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) \<le> m"
     hence "?P (- ?s)" using 0 by auto
     thus ?thesis by blast
   next
     assume "\<not> (2 * (k - m * ?s) \<le> 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 \<open>Lemma 1.2 from @{cite nathanson1996}.\<close>
 
 lemma lemma_1_2:
   fixes b :: mat2
   assumes "mat_sym b"
   assumes "qf_positive_definite b"
   shows "\<exists>a. a ~ b \<and>
              2 * \<bar>mat2\<^sub>1\<^sub>2 a\<bar> \<le> mat2\<^sub>1\<^sub>1 a \<and>
              mat2\<^sub>1\<^sub>1 a \<le> (2 / sqrt 3) * sqrt (mat_det a)" (is "\<exists>a. ?P a")
 proof -
   define a\<^sub>1\<^sub>1 where "a\<^sub>1\<^sub>1 \<equiv> LEAST y. y > 0 \<and> (\<exists>x. int y = b $$ x)"
   have 0: "\<exists>y. y > 0 \<and> (\<exists>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: "a\<^sub>1\<^sub>1 > 0" "int a\<^sub>1\<^sub>1 = b $$ r"
     using a\<^sub>1\<^sub>1_def LeastI_ex[OF 0] by auto
   let ?h = "gcd (vec2\<^sub>1 r) (vec2\<^sub>2 r)"
   have "r \<noteq> 0" using assms(2) 1 by fastforce
   hence 2: "?h > 0" by (simp; metis vec2.collapse zero_vec2_def)
   let ?r' = "vec2 (vec2\<^sub>1 r div ?h) (vec2\<^sub>2 r div ?h)"
   have "?r' \<noteq> 0" using 2 unfolding zero_vec2_def
     by (simp add: algebra_simps dvd_div_eq_0_iff)
   hence "nat (b $$ ?r') > 0 \<and> (\<exists>x. int (nat (b $$ ?r')) = b $$ x)"
     using assms(2) unfolding qf2_positive_definite_def by auto
-  hence "a\<^sub>1\<^sub>1 \<le> nat (b $$ ?r')" unfolding a\<^sub>1\<^sub>1_def by (rule Least_le)
+  hence *: "a\<^sub>1\<^sub>1 \<le> nat (b $$ ?r')" unfolding a\<^sub>1\<^sub>1_def by (rule Least_le)
   hence "a\<^sub>1\<^sub>1 \<le> b $$ ?r'" using 1 by auto
-  also have "... = (b $$ r) div ?h\<^sup>2" proof -
+  also have **: "... = (b $$ r) div ?h\<^sup>2" proof -
     have "(b $$ ?r') * ?h\<^sup>2 = b $$ r" "?h\<^sup>2 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 ?h\<^sup>2" using 2 by auto
   qed
   also have "... = a\<^sub>1\<^sub>1 div ?h\<^sup>2" using 1 by auto
   finally have "a\<^sub>1\<^sub>1 \<le> a\<^sub>1\<^sub>1 div ?h\<^sup>2" .
   also have "... \<le> a\<^sub>1\<^sub>1" using 1 2
-    by (smt (z3) div_by_1 int_div_less_self of_nat_0_less_iff zero_less_power)
+    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 s\<^sub>1 s\<^sub>2 where 3: "1 = (vec2\<^sub>1 r) * s\<^sub>2 - (vec2\<^sub>2 r) * s\<^sub>1"
     by (metis bezout_int diff_minus_eq_add mult.commute mult_minus_right)
   define a'\<^sub>1\<^sub>2 where "
     a'\<^sub>1\<^sub>2 \<equiv>
       (mat2\<^sub>1\<^sub>1 b) * (vec2\<^sub>1 r) * s\<^sub>1
     + (mat2\<^sub>1\<^sub>2 b) * ((vec2\<^sub>1 r) * s\<^sub>2 + (vec2\<^sub>2 r) * s\<^sub>1)
     + (mat2\<^sub>2\<^sub>2 b) * (vec2\<^sub>2 r) * s\<^sub>2
   "
   obtain t where 4: "2 * \<bar>a'\<^sub>1\<^sub>2 + a\<^sub>1\<^sub>1 * (t :: int)\<bar> \<le> a\<^sub>1\<^sub>1"
     using 1 congruence_class_close by fastforce
   define a\<^sub>1\<^sub>2 where "a\<^sub>1\<^sub>2 \<equiv> a'\<^sub>1\<^sub>2 + a\<^sub>1\<^sub>1 * t"
   define a\<^sub>2\<^sub>2 where "a\<^sub>2\<^sub>2 \<equiv> b $$ (vec2 (s\<^sub>1 + (vec2\<^sub>1 r) * t) (s\<^sub>2 + (vec2\<^sub>2 r) * t))"
   let ?u = "
     mat2
       (vec2\<^sub>1 r) (s\<^sub>1 + (vec2\<^sub>1 r) * t)
       (vec2\<^sub>2 r) (s\<^sub>2 + (vec2\<^sub>2 r) * t)
   "
   let ?a = "
     mat2
       a\<^sub>1\<^sub>1 a\<^sub>1\<^sub>2
       a\<^sub>1\<^sub>2 a\<^sub>2\<^sub>2
   "
   have 5: "mat_det ?u = 1" unfolding mat2_det_def
     using 3 by (simp add: algebra_simps)
   have 6: "?a = b \<cdot> ?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
               a\<^sub>1\<^sub>2_def a\<^sub>1\<^sub>2_def a\<^sub>2\<^sub>2_def a'\<^sub>1\<^sub>2_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 * \<bar>mat2\<^sub>1\<^sub>2 ?a\<bar> \<le> mat2\<^sub>1\<^sub>1 ?a" using 4 unfolding a\<^sub>1\<^sub>2_def by auto
   have "a\<^sub>1\<^sub>1 \<le> int (nat a\<^sub>2\<^sub>2)"
     unfolding a\<^sub>1\<^sub>1_def a\<^sub>2\<^sub>2_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 "a\<^sub>1\<^sub>1 \<le> a\<^sub>2\<^sub>2" using 1 by linarith
   hence "4 * a\<^sub>1\<^sub>1\<^sup>2 \<le> 4 * a\<^sub>1\<^sub>1 * a\<^sub>2\<^sub>2" unfolding power2_eq_square using 1 by auto
   also have "... = 4 * (mat_det ?a + a\<^sub>1\<^sub>2\<^sup>2)"
     unfolding mat2_det_def power2_eq_square by auto
   also have "... = 4 * mat_det ?a + (2 * \<bar>a\<^sub>1\<^sub>2\<bar>)\<^sup>2"
     unfolding power2_eq_square by auto
   also have "... \<le> 4 * mat_det ?a + (int a\<^sub>1\<^sub>1)\<^sup>2"
     using 4 power2_le_iff_abs_le unfolding a\<^sub>1\<^sub>2_def by (smt (verit))
   finally have "3 * a\<^sub>1\<^sub>1\<^sup>2 \<le> 4 * (mat_det ?a)" by auto
   hence "a\<^sub>1\<^sub>1\<^sup>2 \<le> (4 / 3) * mat_det ?a" by linarith
   hence "sqrt (a\<^sub>1\<^sub>1\<^sup>2) \<le> sqrt ((4 / 3) * mat_det ?a)"
     using real_sqrt_le_mono by blast
   hence "a\<^sub>1\<^sub>1 \<le> sqrt ((4 / 3) * mat_det ?a)" by auto
   hence "a\<^sub>1\<^sub>1 \<le> (sqrt 4) / (sqrt 3) * sqrt (mat_det ?a)"
     unfolding real_sqrt_mult real_sqrt_divide by blast
   hence 9: "mat2\<^sub>1\<^sub>1 ?a \<le> (2 / sqrt 3) * sqrt (mat_det ?a)" by auto
   have "?P ?a" using 7 8 9 by blast
   thus ?thesis by blast
 qed
 
 text \<open>Theorem 1.2 from @{cite nathanson1996}.\<close>
 
 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 * \<bar>mat2\<^sub>1\<^sub>2 a\<bar> \<le> mat2\<^sub>1\<^sub>1 a"
        "mat2\<^sub>1\<^sub>1 a \<le> (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: "mat2\<^sub>1\<^sub>1 a \<ge> 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: "mat2\<^sub>1\<^sub>1 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: "mat2\<^sub>1\<^sub>1 a = 1" using 4 5 by auto
   have 7: "mat2\<^sub>1\<^sub>2 a = 0" "mat2\<^sub>2\<^sub>1 a = 0" using 0 1 6 unfolding mat2_sym_def by auto
   have 8: "mat2\<^sub>2\<^sub>2 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 \<open>Lemma 1.3 from @{cite nathanson1996}.\<close>
 
 lemma lemma_1_3:
   fixes a :: mat3
   assumes "mat_sym a"
   defines "a' \<equiv>
     mat2
       (mat3\<^sub>1\<^sub>1 a * mat3\<^sub>2\<^sub>2 a - (mat3\<^sub>1\<^sub>2 a)\<^sup>2) (mat3\<^sub>1\<^sub>1 a * mat3\<^sub>2\<^sub>3 a - mat3\<^sub>1\<^sub>2 a * mat3\<^sub>1\<^sub>3 a)
       (mat3\<^sub>1\<^sub>1 a * mat3\<^sub>2\<^sub>3 a - mat3\<^sub>1\<^sub>2 a * mat3\<^sub>1\<^sub>3 a) (mat3\<^sub>1\<^sub>1 a * mat3\<^sub>3\<^sub>3 a - (mat3\<^sub>1\<^sub>3 a)\<^sup>2)
   "
   defines "d' \<equiv>
     mat_det (
       mat2
         (mat3\<^sub>1\<^sub>1 a) (mat3\<^sub>1\<^sub>2 a)
         (mat3\<^sub>1\<^sub>2 a) (mat3\<^sub>2\<^sub>2 a)
     )
   "
   shows
     "mat_det a' = mat3\<^sub>1\<^sub>1 a * mat_det a" (is "?P")
     "\<And>x. mat3\<^sub>1\<^sub>1 a * (a $$ x) =
       (mat3\<^sub>1\<^sub>1 a * vec3\<^sub>1 x + mat3\<^sub>1\<^sub>2 a * vec3\<^sub>2 x + mat3\<^sub>1\<^sub>3 a * vec3\<^sub>3 x)\<^sup>2 +
       (a' $$ (vec2 (vec3\<^sub>2 x) (vec3\<^sub>3 x)))" (is "\<And>x. ?Q x")
     "qf_positive_definite a \<Longrightarrow> qf_positive_definite a'"
     "qf_positive_definite a \<longleftrightarrow> mat3\<^sub>1\<^sub>1 a > 0 \<and> d' > 0 \<and> 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: "\<And>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 \<Longrightarrow> mat3\<^sub>1\<^sub>1 a > 0"
   proof -
     assume 3: "qf_positive_definite a"
     show "mat3\<^sub>1\<^sub>1 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 \<Longrightarrow> qf_positive_definite a'"
   unfolding qf2_positive_definite_def
   proof
     fix v :: vec2
     assume 5: "qf_positive_definite a"
     hence 6: "mat3\<^sub>1\<^sub>1 a > 0" using 2 by blast
     have "a' $$ v \<le> 0 \<Longrightarrow> v = 0" proof -
       assume 7: "a' $$ v \<le> 0"
       obtain x\<^sub>2 x\<^sub>3 where 8: "v = vec2 x\<^sub>2 x\<^sub>3" by (rule vec2.exhaust)
       define x\<^sub>1 where "x\<^sub>1 \<equiv> - (mat3\<^sub>1\<^sub>2 a * x\<^sub>2 + mat3\<^sub>1\<^sub>3 a * x\<^sub>3)"
       have "mat3\<^sub>1\<^sub>1 a * (a $$ (vec3 x\<^sub>1 (mat3\<^sub>1\<^sub>1 a * x\<^sub>2) (mat3\<^sub>1\<^sub>1 a * x\<^sub>3))) =
         (mat3\<^sub>1\<^sub>1 a * x\<^sub>1 + mat3\<^sub>1\<^sub>2 a * mat3\<^sub>1\<^sub>1 a * x\<^sub>2 + mat3\<^sub>1\<^sub>3 a * mat3\<^sub>1\<^sub>1 a * x\<^sub>3)\<^sup>2 +
         (a' $$ (vec2 (mat3\<^sub>1\<^sub>1 a * x\<^sub>2) (mat3\<^sub>1\<^sub>1 a * x\<^sub>3)))"
         unfolding 1 by (simp add: algebra_simps)
       also have "... = a' $$ (vec2 (mat3\<^sub>1\<^sub>1 a * x\<^sub>2) (mat3\<^sub>1\<^sub>1 a * x\<^sub>3))"
         unfolding x\<^sub>1_def by (simp add: algebra_simps)
       also have "... = (mat3\<^sub>1\<^sub>1 a)\<^sup>2 * (a' $$ v)"
         unfolding 8 vec2_dot_def mat2_app_def qf2_app_def power2_eq_square
         by (simp add: algebra_simps)
       also have "... \<le> 0"
         using 7 unfolding power2_eq_square by (simp add: mult_nonneg_nonpos)
       finally have "mat3\<^sub>1\<^sub>1 a * (a $$ (vec3 x\<^sub>1 (mat3\<^sub>1\<^sub>1 a * x\<^sub>2) (mat3\<^sub>1\<^sub>1 a * x\<^sub>3))) \<le> 0"
         .
       hence "a $$ (vec3 x\<^sub>1 (mat3\<^sub>1\<^sub>1 a * x\<^sub>2) (mat3\<^sub>1\<^sub>1 a * x\<^sub>3)) \<le> 0"
         using 6 by (simp add: mult_le_0_iff)
       hence "vec3 x\<^sub>1 (mat3\<^sub>1\<^sub>1 a * x\<^sub>2) (mat3\<^sub>1\<^sub>1 a * x\<^sub>3) = 0"
         using 5 unfolding qf3_positive_definite_def by fastforce
       hence "x\<^sub>2 = 0" "x\<^sub>3 = 0" using 6 unfolding zero_vec3_def by fastforce+
       thus "v = 0" unfolding 8 zero_vec2_def by blast
     qed
     thus "v \<noteq> 0 \<longrightarrow> 0 < a' $$ v" by fastforce
   qed
   have 9: "qf_positive_definite a \<Longrightarrow> d' > 0 \<and> mat_det a > 0"
   proof -
     assume 10: "qf_positive_definite a"
     have 11: "mat3\<^sub>1\<^sub>1 a > 0" using 2 10 by blast
     have "qf_positive_definite a'" using 4 10 by blast
     hence 12: "mat2\<^sub>1\<^sub>1 a' > 0 \<and> 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 \<and> (mat_det a) > 0" using 13 14 by blast
   qed
   have 15: "mat3\<^sub>1\<^sub>1 a > 0 \<and> d' > 0 \<and> mat_det a > 0 \<Longrightarrow> qf_positive_definite a"
   proof -
     assume 16: "mat3\<^sub>1\<^sub>1 a > 0 \<and> d' > 0 \<and> mat_det a > 0"
     have 17: "mat2\<^sub>1\<^sub>1 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 \<le> 0 \<Longrightarrow> x = 0" proof -
         assume "a $$ x \<le> 0"
         hence 20: "mat3\<^sub>1\<^sub>1 a * (a $$ x) \<le> 0"
           using 16 mult_le_0_iff order_less_le by auto
         have "a' $$ (vec2 (vec3\<^sub>2 x) (vec3\<^sub>3 x)) \<le> 0"
           using 20 unfolding 1 power2_eq_square by (smt (verit) zero_le_square)
         hence 21: "a' $$ (vec2 (vec3\<^sub>2 x) (vec3\<^sub>3 x)) = 0"
           using 19 qf2_positive_definite_positive
           using nle_le by blast
         have 22: "vec3\<^sub>2 x = 0" "vec3\<^sub>3 x = 0"
           using 19 21 unfolding zero_vec2_def qf2_positive_definite_def
           by (smt (verit) vec2.inject)+
         have "mat3\<^sub>1\<^sub>1 a * vec3\<^sub>1 x + mat3\<^sub>1\<^sub>2 a * vec3\<^sub>2 x + mat3\<^sub>1\<^sub>3 a * vec3\<^sub>3 x = 0"
           using 20 21 unfolding 1 by fastforce
         hence "vec3\<^sub>1 x = 0" using 16 22 by fastforce
         thus "x = 0" using 22 unfolding zero_vec3_def by (metis vec3.collapse)
       qed
       thus "x \<noteq> 0 \<longrightarrow> 0 < a $$ x" by fastforce
     qed
   qed
   show "qf_positive_definite a \<longleftrightarrow> mat3\<^sub>1\<^sub>1 a > 0 \<and> d' > 0 \<and> mat_det a > 0"
     using 2 9 15 by blast
 qed
 
 text \<open>Lemma 1.4 from @{cite nathanson1996}.\<close>
 
 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' \<equiv>
     mat2
       (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>2\<^sub>2 b - (mat3\<^sub>1\<^sub>2 b)\<^sup>2) (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>2\<^sub>3 b - mat3\<^sub>1\<^sub>2 b * mat3\<^sub>1\<^sub>3 b)
       (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>2\<^sub>3 b - mat3\<^sub>1\<^sub>2 b * mat3\<^sub>1\<^sub>3 b) (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>3\<^sub>3 b - (mat3\<^sub>1\<^sub>3 b)\<^sup>2)
   "
   defines "a' \<equiv> b' \<cdot> v'"
   defines "v \<equiv>
     mat3
       1 r s
       0 (mat2\<^sub>1\<^sub>1 v') (mat2\<^sub>1\<^sub>2 v')
       0 (mat2\<^sub>2\<^sub>1 v') (mat2\<^sub>2\<^sub>2 v')
   "
   defines "a \<equiv> b \<cdot> v"
   shows
     "\<And>y. mat3\<^sub>1\<^sub>1 b * (b $$ y) =
       (mat3\<^sub>1\<^sub>1 b * vec3\<^sub>1 y + mat3\<^sub>1\<^sub>2 b * vec3\<^sub>2 y + mat3\<^sub>1\<^sub>3 b * vec3\<^sub>3 y)\<^sup>2 +
       (b' $$ (vec2 (vec3\<^sub>2 y) (vec3\<^sub>3 y)))" (is "\<And>y. ?P y")
     "mat3\<^sub>1\<^sub>1 a = mat3\<^sub>1\<^sub>1 b"
     "\<And>x. mat3\<^sub>1\<^sub>1 a * (a $$ x) =
       (mat3\<^sub>1\<^sub>1 a * vec3\<^sub>1 x + mat3\<^sub>1\<^sub>2 a * vec3\<^sub>2 x + mat3\<^sub>1\<^sub>3 a * vec3\<^sub>3 x)\<^sup>2 +
       (a' $$ (vec2 (vec3\<^sub>2 x) (vec3\<^sub>3 x)))" (is "\<And>x. ?Q x")
 proof -
   show "\<And>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 "mat3\<^sub>1\<^sub>1 a = mat3\<^sub>1\<^sub>1 b"
     unfolding a_def v_def times_mat3_def mat3_transpose_def qf3_action_def
     by force
   show "\<And>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 \<open>Lemma 1.5 from @{cite nathanson1996}.\<close>
 
 lemma lemma_1_5:
   fixes u\<^sub>1\<^sub>1 u\<^sub>2\<^sub>1 u\<^sub>3\<^sub>1
   assumes "Gcd {u\<^sub>1\<^sub>1, u\<^sub>2\<^sub>1, u\<^sub>3\<^sub>1} = 1"
   shows "\<exists>u. mat3\<^sub>1\<^sub>1 u = u\<^sub>1\<^sub>1 \<and> mat3\<^sub>2\<^sub>1 u = u\<^sub>2\<^sub>1 \<and> mat3\<^sub>3\<^sub>1 u = u\<^sub>3\<^sub>1 \<and> mat_det u = 1"
 proof -
   let ?a = "gcd u\<^sub>1\<^sub>1 u\<^sub>2\<^sub>1"
   show ?thesis proof cases
     assume "?a = 0"
     hence 0: "u\<^sub>1\<^sub>1 = 0" "u\<^sub>2\<^sub>1 = 0" "u\<^sub>3\<^sub>1 = 1 \<or> u\<^sub>3\<^sub>1 = -1" using assms by auto
     let ?u = "
       mat3
       0 0 (- 1)
       0 u\<^sub>3\<^sub>1 0
       u\<^sub>3\<^sub>1 0 0"
     show ?thesis
       apply (rule exI[of _ ?u])
       unfolding mat3_det_def using 0 apply auto
       done
   next
     assume 1: "?a \<noteq> 0"
     obtain u\<^sub>1\<^sub>2 u\<^sub>2\<^sub>2 where 2: "u\<^sub>1\<^sub>1 * u\<^sub>2\<^sub>2 - u\<^sub>2\<^sub>1 * u\<^sub>1\<^sub>2 = ?a" using bezout_int
       by (metis diff_minus_eq_add mult.commute mult_minus_right)
     have "gcd ?a u\<^sub>3\<^sub>1 = 1" using assms by (simp add: gcd.assoc)
     then obtain u\<^sub>3\<^sub>3 b where 3: "?a * u\<^sub>3\<^sub>3 - b * u\<^sub>3\<^sub>1 = 1" using bezout_int
       by (metis diff_minus_eq_add mult.commute mult_minus_right)
     let ?u = "
       mat3
       u\<^sub>1\<^sub>1 u\<^sub>1\<^sub>2 ((u\<^sub>1\<^sub>1 div ?a) * b)
       u\<^sub>2\<^sub>1 u\<^sub>2\<^sub>2 ((u\<^sub>2\<^sub>1 div ?a) * b)
       u\<^sub>3\<^sub>1 0 u\<^sub>3\<^sub>3"
     have "mat_det ?u =
           u\<^sub>1\<^sub>1 * u\<^sub>2\<^sub>2 * u\<^sub>3\<^sub>3
         + u\<^sub>1\<^sub>2 * ((u\<^sub>2\<^sub>1 div ?a) * b) * u\<^sub>3\<^sub>1
         - ((u\<^sub>1\<^sub>1 div ?a) * b) * u\<^sub>2\<^sub>2 * u\<^sub>3\<^sub>1
         - u\<^sub>1\<^sub>2 * u\<^sub>2\<^sub>1 * u\<^sub>3\<^sub>3"
       unfolding mat3_det_def by auto
     also have "... =
           u\<^sub>1\<^sub>1 * u\<^sub>2\<^sub>2 * u\<^sub>3\<^sub>3
         + u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a) * (b * u\<^sub>3\<^sub>1)
         - u\<^sub>2\<^sub>2 * (u\<^sub>1\<^sub>1 div ?a) * (b * u\<^sub>3\<^sub>1)
         - u\<^sub>1\<^sub>2 * u\<^sub>2\<^sub>1 * u\<^sub>3\<^sub>3"
       by auto
     also have "... =
           u\<^sub>1\<^sub>1 * u\<^sub>2\<^sub>2 * u\<^sub>3\<^sub>3
         + u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a) * (?a * u\<^sub>3\<^sub>3 - 1)
         - u\<^sub>2\<^sub>2 * (u\<^sub>1\<^sub>1 div ?a) * (?a * u\<^sub>3\<^sub>3 - 1)
         - u\<^sub>1\<^sub>2 * u\<^sub>2\<^sub>1 * u\<^sub>3\<^sub>3"
       using 3 by (simp add: algebra_simps)
     also have "... =
           u\<^sub>1\<^sub>1 * u\<^sub>2\<^sub>2 * u\<^sub>3\<^sub>3
         + u\<^sub>1\<^sub>2 * ((u\<^sub>2\<^sub>1 div ?a) * ?a) * u\<^sub>3\<^sub>3
         - u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a)
         - u\<^sub>2\<^sub>2 * ((u\<^sub>1\<^sub>1 div ?a) * ?a) * u\<^sub>3\<^sub>3
         + u\<^sub>2\<^sub>2 * (u\<^sub>1\<^sub>1 div ?a)
         - u\<^sub>1\<^sub>2 * u\<^sub>2\<^sub>1 * u\<^sub>3\<^sub>3"
       by (simp add: algebra_simps)
     also have "... =
         - u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a)
         + u\<^sub>2\<^sub>2 * (u\<^sub>1\<^sub>1 div ?a)"
       by (simp add: algebra_simps)
     also have "... =
         - u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a)
         + u\<^sub>2\<^sub>2 * u\<^sub>1\<^sub>1 div ?a"
       by (metis dvd_div_mult gcd_dvd1 mult.commute)
     also have "... =
         - u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a)
         + (?a + u\<^sub>2\<^sub>1 * u\<^sub>1\<^sub>2) div ?a"
       using 2 by (simp add: algebra_simps)
     also have "... =
         - u\<^sub>1\<^sub>2 * (u\<^sub>2\<^sub>1 div ?a)
         + 1 + (u\<^sub>2\<^sub>1 * u\<^sub>1\<^sub>2) 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 \<open>Lemma 1.6 from @{cite nathanson1996}.\<close>
 
 lemma lemma_1_6:
   fixes c :: mat3
   assumes "mat_sym c"
   assumes "qf_positive_definite c"
   shows "\<exists>a. a ~ c \<and>
              2 * (max \<bar>mat3\<^sub>1\<^sub>2 a\<bar> \<bar>mat3\<^sub>1\<^sub>3 a\<bar>) \<le> mat3\<^sub>1\<^sub>1 a \<and>
              mat3\<^sub>1\<^sub>1 a \<le> (4 / 3) * root 3 (mat_det a)"
 proof -
   define a\<^sub>1\<^sub>1 where "a\<^sub>1\<^sub>1 \<equiv> LEAST y. y > 0 \<and> (\<exists>x. int y = c $$ x)"
   have 0: "\<exists>y. y > 0 \<and> (\<exists>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: "a\<^sub>1\<^sub>1 > 0" "int a\<^sub>1\<^sub>1 = c $$ t"
     using a\<^sub>1\<^sub>1_def LeastI_ex[OF 0] by auto
   let ?h = "Gcd {vec3\<^sub>1 t, vec3\<^sub>2 t, vec3\<^sub>3 t}"
   have "t \<noteq> 0" using assms(2) 1 by fastforce
   hence 2: "?h > 0" by (simp; metis vec3.collapse zero_vec3_def)
   let ?t' = "vec3 (vec3\<^sub>1 t div ?h) (vec3\<^sub>2 t div ?h) (vec3\<^sub>3 t div ?h)"
   have "?t' \<noteq> 0" using 2 unfolding zero_vec3_def
     by (auto simp add: algebra_simps dvd_div_eq_0_iff)
   hence "nat (c $$ ?t') > 0 \<and> (\<exists>x. int (nat (c $$ ?t')) = c $$ x)"
     using assms(2) unfolding qf3_positive_definite_def by auto
   hence "a\<^sub>1\<^sub>1 \<le> nat (c $$ ?t')" unfolding a\<^sub>1\<^sub>1_def by (rule Least_le)
   hence "a\<^sub>1\<^sub>1 \<le> c $$ ?t'" using 1 by auto
   also have "... = (c $$ t) div ?h\<^sup>2" proof -
     have "?h dvd vec3\<^sub>1 t" "?h dvd vec3\<^sub>2 t" "?h dvd vec3\<^sub>3 t"
       by (meson Gcd_dvd insertCI)+
     then have "(c $$ ?t') * ?h\<^sup>2 = c $$ t" "?h\<^sup>2 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 ?h\<^sup>2" using 2 by auto
   qed
   also have "... = a\<^sub>1\<^sub>1 div ?h\<^sup>2" using 1 by auto
   finally have "a\<^sub>1\<^sub>1 \<le> a\<^sub>1\<^sub>1 div ?h\<^sup>2" .
   also have "... \<le> a\<^sub>1\<^sub>1" using 1 2
-    by (smt (z3) div_by_1 int_div_less_self of_nat_0_less_iff zero_less_power)
+    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: "mat3\<^sub>1\<^sub>1 u = vec3\<^sub>1 t" "mat3\<^sub>2\<^sub>1 u = vec3\<^sub>2 t"
                          "mat3\<^sub>3\<^sub>1 u = vec3\<^sub>3 t" "mat_det u = 1"
     using lemma_1_5 by blast
   define b where "b \<equiv> c \<cdot> 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: "a\<^sub>1\<^sub>1 = (LEAST y. y > 0 \<and> (\<exists>x. int y = b $$ x))"
     unfolding a\<^sub>1\<^sub>1_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: "a\<^sub>1\<^sub>1 = mat3\<^sub>1\<^sub>1 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' \<equiv>
     mat2
       (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>2\<^sub>2 b - (mat3\<^sub>1\<^sub>2 b)\<^sup>2) (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>2\<^sub>3 b - mat3\<^sub>1\<^sub>2 b * mat3\<^sub>1\<^sub>3 b)
       (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>2\<^sub>3 b - mat3\<^sub>1\<^sub>2 b * mat3\<^sub>1\<^sub>3 b) (mat3\<^sub>1\<^sub>1 b * mat3\<^sub>3\<^sub>3 b - (mat3\<^sub>1\<^sub>3 b)\<^sup>2)
   "
   have 8: "mat_sym b'" unfolding b'_def mat2_sym_def by auto
   have 9: "mat_det b' = mat3\<^sub>1\<^sub>1 b * mat_det b"
           "\<And>x. mat3\<^sub>1\<^sub>1 b * (b $$ x) =
             (mat3\<^sub>1\<^sub>1 b * vec3\<^sub>1 x + mat3\<^sub>1\<^sub>2 b * vec3\<^sub>2 x + mat3\<^sub>1\<^sub>3 b * vec3\<^sub>3 x)\<^sup>2 +
             (b' $$ (vec2 (vec3\<^sub>2 x) (vec3\<^sub>3 x)))"
           "qf_positive_definite b'"
     using 4 5 b'_def lemma_1_3 by blast+
   obtain a' v' where 10: "a' = b' \<cdot> v'"
                          "mat_det v' = 1"
                          "mat2\<^sub>1\<^sub>1 a' \<le> (2 / sqrt 3) * sqrt (mat3\<^sub>1\<^sub>1 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 * \<bar>(mat3\<^sub>1\<^sub>2 b) * (mat2\<^sub>1\<^sub>1 v') +
                             (mat3\<^sub>1\<^sub>3 b) * (mat2\<^sub>2\<^sub>1 v') + a\<^sub>1\<^sub>1 * (r :: int)\<bar> \<le> a\<^sub>1\<^sub>1"
                        "2 * \<bar>(mat3\<^sub>1\<^sub>2 b) * (mat2\<^sub>1\<^sub>2 v') +
                             (mat3\<^sub>1\<^sub>3 b) * (mat2\<^sub>2\<^sub>2 v') + a\<^sub>1\<^sub>1 * (s :: int)\<bar> \<le> a\<^sub>1\<^sub>1"
     using 1 congruence_class_close by fastforce
   define a\<^sub>1\<^sub>2 where "a\<^sub>1\<^sub>2 \<equiv> (mat3\<^sub>1\<^sub>2 b) * (mat2\<^sub>1\<^sub>1 v') +
                           (mat3\<^sub>1\<^sub>3 b) * (mat2\<^sub>2\<^sub>1 v') + a\<^sub>1\<^sub>1 * r"
   define a\<^sub>1\<^sub>3 where "a\<^sub>1\<^sub>3 \<equiv> (mat3\<^sub>1\<^sub>2 b) * (mat2\<^sub>1\<^sub>2 v') +
                           (mat3\<^sub>1\<^sub>3 b) * (mat2\<^sub>2\<^sub>2 v') + a\<^sub>1\<^sub>1 * s"
   define v where "v \<equiv>
     mat3
       1 r s
       0 (mat2\<^sub>1\<^sub>1 v') (mat2\<^sub>1\<^sub>2 v')
       0 (mat2\<^sub>2\<^sub>1 v') (mat2\<^sub>2\<^sub>2 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 \<equiv> b \<cdot> 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: "a\<^sub>1\<^sub>1 = (LEAST y. y > 0 \<and> (\<exists>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: "mat3\<^sub>1\<^sub>1 a = mat3\<^sub>1\<^sub>1 b"
            "\<And>x. mat3\<^sub>1\<^sub>1 a * (a $$ x) =
               (mat3\<^sub>1\<^sub>1 a * vec3\<^sub>1 x + mat3\<^sub>1\<^sub>2 a * vec3\<^sub>2 x + mat3\<^sub>1\<^sub>3 a * vec3\<^sub>3 x)\<^sup>2 +
               (a' $$ (vec2 (vec3\<^sub>2 x) (vec3\<^sub>3 x)))"
     using 4 5 10 lemma_1_4 unfolding b'_def v_def a_def by blast+
   have 16: "mat2\<^sub>1\<^sub>1 a' = mat3\<^sub>1\<^sub>1 a * mat3\<^sub>2\<^sub>2 a - (mat3\<^sub>1\<^sub>2 a)\<^sup>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 a\<^sub>2\<^sub>2 where "a\<^sub>2\<^sub>2 \<equiv> a $$ (vec3 0 1 0)"
   have 17: "a\<^sub>1\<^sub>1 = mat3\<^sub>1\<^sub>1 a" unfolding 7 15 by auto
   have 18: "a\<^sub>1\<^sub>2 = mat3\<^sub>1\<^sub>2 a" "a\<^sub>1\<^sub>3 = mat3\<^sub>1\<^sub>3 a" "a\<^sub>2\<^sub>2 = mat3\<^sub>2\<^sub>2 a" "a\<^sub>2\<^sub>2 = mat3\<^sub>2\<^sub>2 a"
     using 13(1) 17
     unfolding a_def v_def a\<^sub>1\<^sub>2_def a\<^sub>1\<^sub>3_def a\<^sub>2\<^sub>2_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 \<bar>mat3\<^sub>1\<^sub>2 a\<bar> \<bar>mat3\<^sub>1\<^sub>3 a\<bar>) \<le> mat3\<^sub>1\<^sub>1 a"
     using 11 unfolding 17[symmetric] 18 a\<^sub>1\<^sub>2_def[symmetric] a\<^sub>1\<^sub>3_def[symmetric]
     by auto
   have 21: "(2 / sqrt 3) ^ 6 = 64 / 27" unfolding power_def by auto
   have "a\<^sub>1\<^sub>1 \<le> int (nat a\<^sub>2\<^sub>2)"
     unfolding a\<^sub>1\<^sub>1_def a\<^sub>2\<^sub>2_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 "a\<^sub>1\<^sub>1 \<le> a\<^sub>2\<^sub>2" using 1 by linarith
   hence "(mat3\<^sub>1\<^sub>1 a)\<^sup>2 \<le> a\<^sub>1\<^sub>1 * a\<^sub>2\<^sub>2" unfolding power2_eq_square using 1 17 by auto
   also have "... = mat2\<^sub>1\<^sub>1 a' + a\<^sub>1\<^sub>2\<^sup>2" using 16 17 18 by auto
   also have "... \<le> (2 / sqrt 3) * sqrt (mat3\<^sub>1\<^sub>1 b * mat_det b) + a\<^sub>1\<^sub>2\<^sup>2"
     using 10 by auto
   also have "... \<le> (2 / sqrt 3) * sqrt (mat3\<^sub>1\<^sub>1 a * mat_det a) + (mat3\<^sub>1\<^sub>1 a)\<^sup>2 / 4"
     using 11 13 15 17 18 a\<^sub>1\<^sub>2_def
           power2_le_iff_abs_le[of "real_of_int (mat3\<^sub>1\<^sub>1 a)" "(mat3\<^sub>1\<^sub>2 a) * 2"]
     by auto
   finally have "(3 / 4) * (mat3\<^sub>1\<^sub>1 a)\<^sup>2 \<le> (2 / sqrt 3) * sqrt (mat3\<^sub>1\<^sub>1 a * mat_det a)"
     by (simp add: algebra_simps)
   hence "(mat3\<^sub>1\<^sub>1 a)\<^sup>2 \<le> (2 / sqrt 3) ^ 3 * sqrt (mat3\<^sub>1\<^sub>1 a * mat_det a)"
     unfolding power2_eq_square power3_eq_cube by (simp add: algebra_simps)
   hence "((mat3\<^sub>1\<^sub>1 a)\<^sup>2)\<^sup>2 \<le> ((2 / sqrt 3) ^ 3 * sqrt (mat3\<^sub>1\<^sub>1 a * mat_det a))\<^sup>2"
     using 1(1) unfolding 17[symmetric] power2_eq_square
     by (metis of_int_power power2_eq_square power_mono zero_le_square)
   hence "(mat3\<^sub>1\<^sub>1 a) ^ 4 \<le> (2 / sqrt 3) ^ 6 * (sqrt (mat3\<^sub>1\<^sub>1 a * mat_det a))\<^sup>2"
     by (simp add: power_mult_distrib)
   hence "(mat3\<^sub>1\<^sub>1 a) ^ 4 \<le> (2 / sqrt 3) ^ 6 * mat3\<^sub>1\<^sub>1 a * mat_det a"
     using 4 5 13 17 lemma_1_3 by auto
   hence "(mat3\<^sub>1\<^sub>1 a) ^ 3 \<le> (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 "(mat3\<^sub>1\<^sub>1 a) ^ 3 \<le> (64 / 27) * mat_det a" using 21 by auto
   hence "root 3 ((mat3\<^sub>1\<^sub>1 a) ^ 3) \<le> root 3 ((64 / 27) * mat_det a)" by auto
   hence "root 3 ((mat3\<^sub>1\<^sub>1 a) ^ 3) \<le> root 3 (64 / 27) * root 3 (mat_det a)"
     unfolding real_root_mult by auto
   hence "mat3\<^sub>1\<^sub>1 a \<le> root 3 (64 / 27) * root 3 (mat_det a)"
     using odd_real_root_power_cancel by auto
   hence 22: "mat3\<^sub>1\<^sub>1 a \<le> (4 / 3) * root 3 (mat_det a)"
     using real_root_divide by force
   show ?thesis using 19 20 22 by blast
 qed
 
 text \<open>Theorem 1.3 from @{cite nathanson1996}.\<close>
 
 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 \<and>
         2 * (max \<bar>mat3\<^sub>1\<^sub>2 a\<bar> \<bar>mat3\<^sub>1\<^sub>3 a\<bar>) \<le> mat3\<^sub>1\<^sub>1 a \<and>
         mat3\<^sub>1\<^sub>1 a \<le> (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: "mat3\<^sub>1\<^sub>2 a = 0" "mat3\<^sub>1\<^sub>3 a = 0" using 0 3 by auto
   have "mat3\<^sub>1\<^sub>1 a \<ge> 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: "mat3\<^sub>1\<^sub>1 a = 1" using 0 3 by auto
   define a' where "a' \<equiv>
     mat2
       (mat3\<^sub>2\<^sub>2 a) (mat3\<^sub>2\<^sub>3 a)
       (mat3\<^sub>3\<^sub>2 a) (mat3\<^sub>3\<^sub>3 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: "\<forall>v. v \<noteq> 0 \<longrightarrow> a $$ v > 0"
     show "\<forall>v. v \<noteq> 0 \<longrightarrow> a' $$ v > 0"
     proof (rule; rule)
       fix v :: vec2
       assume "v \<noteq> 0"
       hence "vec3 0 (vec2\<^sub>1 v) (vec2\<^sub>2 v) \<noteq> 0"
         unfolding zero_vec2_def zero_vec3_def by (metis vec2.collapse vec3.inject)
       hence "a $$ (vec3 0 (vec2\<^sub>1 v) (vec2\<^sub>2 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' \<cdot> u' = 1"
     using 7 8 9 qf2_det_one_equiv_canonical[of a']
     unfolding qf2_equiv_def
     by auto
   define u where "u \<equiv>
     mat3
       1 0 0
       0 (mat2\<^sub>1\<^sub>1 u') (mat2\<^sub>1\<^sub>2 u')
       0 (mat2\<^sub>2\<^sub>1 u') (mat2\<^sub>2\<^sub>2 u')
   "
   have 12: "mat_det u = 1"
     using 11 unfolding u_def mat2_det_def mat3_det_def by auto
   have 13: "a \<cdot> 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
diff --git a/thys/Three_Squares/Three_Squares.thy b/thys/Three_Squares/Three_Squares.thy
--- a/thys/Three_Squares/Three_Squares.thy
+++ b/thys/Three_Squares/Three_Squares.thy
@@ -1,1436 +1,1436 @@
 (*  Title:      Three_Squares/Three_Squares.thy
     Author:     Anton Danilkin
     Author:     Loïc Chevalier
 *)
 
 section \<open>Legendre's three squares theorem and its consequences\<close>
 
 theory Three_Squares
   imports "Dirichlet_L.Dirichlet_Theorem" Residues_Properties Quadratic_Forms
 begin
 
 subsection \<open>Legendre's three squares theorem\<close>
 
 definition quadratic_residue_alt :: "int \<Rightarrow> int \<Rightarrow> bool" where
 "quadratic_residue_alt m a = (\<exists>x y. x\<^sup>2 - a = y * m)"
 
 lemma quadratic_residue_alt_equiv: "quadratic_residue_alt = QuadRes"
   unfolding quadratic_residue_alt_def QuadRes_def
   by (metis cong_iff_dvd_diff cong_modulus_mult dvdE dvd_refl mult.commute)
 
 lemma sq_nat_abs: "(nat \<bar>v\<bar>)\<^sup>2 = nat (v\<^sup>2)"
   by (simp add: nat_power_eq[symmetric])
 
 text \<open>Lemma 1.7 from @{cite nathanson1996}.\<close>
 
 lemma three_squares_using_quadratic_residue:
   fixes n d' :: nat
   assumes "n \<ge> 2"
   assumes "d' > 0"
   assumes "QuadRes (d' * n - 1) (- d')"
   shows "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
 proof -
   define a where "a \<equiv> d' * n - 1"
   from assms(3) obtain x y where "x\<^sup>2 + int d' = y * int a"
     unfolding a_def quadratic_residue_alt_equiv[symmetric]
               quadratic_residue_alt_def
     by auto
   hence Hxy: "x\<^sup>2 + d' = y * a" by auto
   have "y \<ge> 1"
     using assms Hxy
     unfolding a_def
     by (smt (verit) bot_nat_0.not_eq_extremum mult_le_0_iff of_nat_0_le_iff
                     of_nat_le_0_iff power2_eq_square zero_le_square)
   moreover from Hxy have Hxy\<^sub>2: "d' = y * a - x\<^sup>2" by (simp add: algebra_simps)
   define M where "M \<equiv> mat3 y x 1 x a 0 1 0 n"
   moreover have Msym: "mat_sym M"
     unfolding mat3_sym_def M_def mat3.sel
     by simp
   moreover have Mdet_eq_1: "mat_det M = 1"
   proof -
     have "mat_det M = (y * a - x\<^sup>2) * n - a"
       unfolding mat3_det_def M_def mat3.sel power2_eq_square
       by (simp add: algebra_simps)
     also have "... = int d' * n - a" unfolding Hxy\<^sub>2 by simp
     also have "... = 1" unfolding a_def using assms int_ops by force
     finally show ?thesis .
   qed
   moreover have "mat_det (mat2 y x x a) > 0"
     using Hxy\<^sub>2 assms
     unfolding mat2_det_def power2_eq_square
     by simp
   ultimately have "qf_positive_definite M"
     by (auto simp add: lemma_1_3(4))
   hence "M ~ 1"
     using Msym and Mdet_eq_1
     by (simp add: qf3_det_one_equiv_canonical)
   moreover have "M $$ vec3 0 0 1 = n"
     unfolding M_def qf3_app_def mat3_app_def vec3_dot_def mat3.sel vec3.sel
     by (simp add: algebra_simps)
   hence "n \<in> range (($$) M)" by (metis rangeI)
   ultimately have "n \<in> range (($$) (1 :: mat3))"
     using qf3_equiv_preserves_range by simp
   then obtain u :: vec3 where "1 $$ u = n"
     by auto
   hence "<u | u> = n"
     unfolding qf3_app_def mat3_app_def one_mat3_def
     by simp
   hence "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. int n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
     unfolding vec3_dot_def power2_eq_square by metis
   hence "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = (nat \<bar>x\<^sub>1\<bar>)\<^sup>2 + (nat \<bar>x\<^sub>2\<bar>)\<^sup>2 + (nat \<bar>x\<^sub>3\<bar>)\<^sup>2"
     unfolding sq_nat_abs
     apply (simp add: nat_add_distrib[symmetric])
     apply (metis nat_int)
     done
   thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" by blast
 qed
 
 lemma prime_linear_combination:
   fixes a m :: nat
   assumes "m > 1"
   assumes "coprime a m"
   obtains j :: nat where "prime (a + m * j) \<and> j \<noteq> 0"
 proof -
   assume 0: "\<And>j. prime (a + m * j) \<and> j \<noteq> 0 \<Longrightarrow> thesis"
 
   have 1: "infinite {p. prime p \<and> [p = a] (mod m)}"
     using assms
     by (rule Dirichlet_Theorem.residues_nat.Dirichlet[unfolded residues_nat_def])
 
   have 2: "finite {j. prime (nat (int a + int m * j)) \<and> j \<le> 0}"
     apply (rule finite_subset[of _ "{- (int a) div (int m)..0}"])
     subgoal
       apply (rule subsetI)
       subgoal for j
       proof -
         assume 1: "j \<in> {j. prime (nat (int a + int m * j)) \<and> j \<le> 0}"
         have "int a + int m * j \<ge> 0" using 1 prime_ge_0_int by force
         hence "int m * j \<ge> - int a" by auto
         hence "j \<ge> (- int a) div int m"
           using assms 1
           by (smt (verit) unique_euclidean_semiring_class.cong_def
                           coprime_crossproduct_nat coprime_iff_invertible_int
                           coprime_int_iff int_distrib(3) int_ops(2) int_ops(7)
                           mem_Collect_eq mult_cancel_right1 zdiv_mono1
                           nonzero_mult_div_cancel_left of_nat_0_eq_iff
                           of_nat_le_0_iff prime_ge_2_int prime_nat_iff_prime)
         thus "j \<in> {- (int a) div (int m)..0}" using 1 by auto
       qed
       done
     subgoal by blast
     done
 
   have "{p. prime p \<and> [p = a] (mod m)} =
         {p. prime p \<and> (\<exists>j. int p = int a + int m * j)}"
     unfolding cong_sym_eq[of _ a]
     unfolding cong_int_iff[symmetric] cong_iff_lin
     ..
   also have "... = {p. prime p \<and> (\<exists>j. p = nat (int a + int m * j))}"
     by (metis (opaque_lifting) nat_int nat_eq_iff
                                prime_ge_0_int prime_nat_iff_prime)
   also have "... = (\<lambda>j. nat (int a + int m * j)) `
                    {j. prime (nat (int a + int m * j))}"
     by blast
   finally have "infinite ((\<lambda>j. nat (int a + int m * j)) `
                           {j. prime (nat (int a + int m * j))})"
     using 1 by metis
   hence "infinite {j. prime (nat (int a + int m * j))}"
     using finite_imageI by blast
   hence "infinite ({j. prime (nat (int a + int m * j))} -
                    {j. prime (nat (int a + int m * j)) \<and> j \<le> 0})"
     using 2 Diff_infinite_finite by blast
   hence "infinite {j. prime (nat (int a + int m * j)) \<and> j > 0}"
     by (rule back_subst[of infinite]; auto)
   hence "infinite (int ` {j. prime (nat (int a + int m * j)) \<and> j \<noteq> (0 :: nat)})"
     apply (rule back_subst[of infinite])
     unfolding image_def using zero_less_imp_eq_int apply auto
     done
   hence "infinite {j. prime (nat (int a + int m * j)) \<and> (j :: nat) \<noteq> 0}"
     using finite_imageI by blast
   hence "infinite {j. prime (a + m * j) \<and> j \<noteq> 0}"
     apply (rule back_subst[of infinite])
     apply (auto simp add: int_ops nat_plus_as_int)
     done
   thus thesis using 0 not_finite_existsD by blast
 qed
 
 text \<open>Lemma 1.8 from @{cite nathanson1996}.\<close>
 
 lemma three_squares_using_mod_four:
   fixes n :: nat
   assumes "n mod 4 = 2"
   shows "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
 proof -
   have "n > 1" using assms by auto
   have "coprime (n - 1) (4 * n)"
     by (smt (verit, del_insts) Suc_pred assms bot_nat_0.not_eq_extremum
                                coprime_commute coprime_diff_one_right_nat
                                coprime_mod_right_iff coprime_mult_left_iff
                                diff_Suc_1 mod_Suc mod_mult_self1_is_0 mult_0_right
                                numeral_2_eq_2 zero_neq_numeral)
   then obtain j where H_j:
       "prime ((n - 1) + (4 * n) * j) \<and> j \<noteq> 0"
     using prime_linear_combination[of "4 * n" "n - 1"] \<open>n > 1\<close> by auto
   have "j > 0" using H_j by blast
 
   define d' where "d' \<equiv> 4 * j + 1"
   define p where "p \<equiv> d' * n - 1"
   have "prime p"
     unfolding p_def d'_def
     using conjunct1[OF H_j] apply (rule back_subst[of prime])
     using \<open>n > 1\<close> apply (simp add: algebra_simps nat_minus_as_int nat_plus_as_int)
     done
   have "p mod 4 = 1"
     unfolding p_def
     apply (subst mod_diff_eq_nat)
     subgoal unfolding d'_def using \<open>n > 1\<close> \<open>j > 0\<close> by simp
     subgoal
       apply (subst mod_mult_eq[symmetric])
       unfolding assms d'_def apply simp
       done
     done
   have "d' * n mod 4 = 2"
     using assms p_def d'_def \<open>p mod 4 = 1\<close>
     by (metis mod_mult_cong mod_mult_self4 nat_mult_1)
   hence "d' mod 4 = 1" using assms by (simp add: d'_def)
 
   have "QuadRes p (- d')"
   proof -
     have d'_expansion: "d' = (\<Prod>q\<in>prime_factors d'. q ^ multiplicity q d')"
       using prime_factorization_nat unfolding d'_def by auto
 
     have "odd d'" unfolding d'_def by simp
     hence d'_prime_factors_odd: "q \<in> prime_factors d' \<Longrightarrow> odd q" for q
       by fastforce
 
     have d'_prime_factors_gt_2: "q \<in> prime_factors d' \<Longrightarrow> q > 2" for q
       using d'_prime_factors_odd
       by (metis even_numeral in_prime_factors_imp_prime
                 order_less_le prime_ge_2_nat)
 
     have "[p = - 1] (mod d')"
       unfolding p_def cong_iff_dvd_diff apply simp
       using \<open>n > 1\<close> apply (smt (verit) Suc_as_int Suc_pred add_gr_0 d'_def
                                        dvd_nat_abs_iff dvd_triv_left
                                        less_numeral_extra(1) mult_pos_pos
                                        of_nat_less_0_iff order_le_less_trans
                                        zero_less_one_class.zero_le_one)
       done
     hence d'_prime_factors_2_p_mod:
         "q \<in> prime_factors d' \<Longrightarrow> [p = - 1] (mod q)" for q
       by (rule cong_dvd_modulus; auto)
 
     have "d' mod 4 = (\<Prod>q\<in>prime_factors d'. q ^ multiplicity q d') mod 4"
       using d'_expansion by argo
     also have "... = (\<Prod>q\<in>prime_factors d'. (q mod 4) ^ multiplicity q d') mod 4"
       apply (subst mod_prod_eq[symmetric])
       apply (subst power_mod[symmetric])
       apply (subst mod_prod_eq)
       apply blast
       done
     also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
                       (q mod 4) ^ multiplicity q d') mod 4"
       apply (rule arg_cong[of _ _ "\<lambda>x. x mod 4"])
       apply (rule prod.mono_neutral_right)
       subgoal by blast
       subgoal by blast
       subgoal
         unfolding unique_euclidean_semiring_class.cong_def
         apply (rule ballI)
         using d'_prime_factors_odd apply simp
         apply (metis One_nat_def dvd_0_right even_even_mod_4_iff
                      even_numeral mod_exhaust_less_4)
         done
       done
     also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
                       ((int q) mod 4) ^ multiplicity q d') mod 4"
       by (simp add: int_ops)
     also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
                       ((- 1) mod 4) ^ multiplicity q d') mod 4"
       apply (rule arg_cong[of _ _ "\<lambda>x. x mod 4"])
       apply (rule prod.cong[OF refl])
       unfolding unique_euclidean_semiring_class.cong_def nat_mod_as_int apply simp
       apply (metis nat_int of_nat_mod of_nat_numeral)
       done
     also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
                       (- 1) ^ multiplicity q d') mod 4"
       apply (subst mod_prod_eq[symmetric])
       apply (subst power_mod)
       apply (subst mod_prod_eq)
       apply blast
       done
     finally have "[\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
             (- 1) ^ multiplicity q d' = 1 :: int] (mod 4)"
       using \<open>d' mod 4 = 1\<close>
       by (simp add: unique_euclidean_semiring_class.cong_def)
     hence prod_prime_factors_minus_one:
         "(\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
           (- 1) ^ multiplicity q d') = (1 :: int)"
       unfolding power_sum[symmetric]
       unfolding minus_one_power_iff unique_euclidean_semiring_class.cong_def
       by presburger
 
     have "p > 2" using \<open>prime p\<close> \<open>p mod 4 = 1\<close> nat_less_le prime_ge_2_nat by force
 
     have d'_prime_factors_Legendre:
       "q \<in> prime_factors d' \<Longrightarrow>
        Legendre p q = Legendre q p" for q
     proof -
       assume "q \<in> prime_factors d'"
       have "prime q" using \<open>q \<in> prime_factors d'\<close> by blast
       have "q > 2" using d'_prime_factors_gt_2 \<open>q \<in> prime_factors d'\<close> by blast
       show "Legendre p q = Legendre q p"
         using \<open>prime p\<close> \<open>p > 2\<close> \<open>p mod 4 = 1\<close>
               \<open>prime q\<close> \<open>q > 2\<close> Legendre_equal[of p q]
         unfolding unique_euclidean_semiring_class.cong_def
         using zmod_int[of p 4]
         by auto
     qed
 
     have "Legendre (- d') p = Legendre (- 1) p * Legendre d' p"
       using \<open>prime p\<close> Legendre_mult[of p "- 1" d'] by auto
     also have "... = Legendre d' p"
       using \<open>prime p\<close> \<open>p > 2\<close> \<open>p mod 4 = 1\<close> Legendre_minus_one[of p]
       unfolding unique_euclidean_semiring_class.cong_def nat_mod_as_int
       by auto
     also have "... = (\<Prod>q\<in>prime_factors d'. Legendre (q ^ multiplicity q d') p)"
       apply (subst d'_expansion)
       using \<open>prime p\<close> \<open>p > 2\<close> Legendre_prod[of p] apply auto
       done
     also have "... = (\<Prod>q\<in>prime_factors d'. (Legendre q p) ^ multiplicity q d')"
       using \<open>prime p\<close> \<open>p > 2\<close> Legendre_power by auto
     also have "... = (\<Prod>q\<in>prime_factors d'. (Legendre p q) ^ multiplicity q d')"
       using d'_prime_factors_Legendre by auto
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       (Legendre (- 1) q) ^ multiplicity q d')"
       apply (rule prod.cong[OF refl])
       using d'_prime_factors_2_p_mod apply (metis Legendre_cong)
       done
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       (if [q = 1] (mod 4) then 1 else - 1) ^ multiplicity q d')"
       apply (rule prod.cong[OF refl])
       subgoal for q
         using Legendre_minus_one_alt[of q] d'_prime_factors_gt_2[of q]
         by (smt (verit) cong_int_iff in_prime_factors_iff int_eq_iff_numeral
                         less_imp_of_nat_less numeral_Bit0 numeral_One of_nat_1
                         prime_nat_int_transfer)
       done
     also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
                       (if [q = 1] (mod 4) then 1 else - 1) ^ multiplicity q d')"
       apply (rule prod.mono_neutral_right)
       subgoal by blast
       subgoal by blast
       subgoal
         unfolding unique_euclidean_semiring_class.cong_def
         apply (rule ballI)
         using d'_prime_factors_odd apply simp
         apply (metis One_nat_def dvd_0_right even_even_mod_4_iff
                      even_numeral mod_exhaust_less_4)
         done
       done
     also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 4)}.
                       (- 1) ^ multiplicity q d')"
       by (rule prod.cong[OF refl];
           simp add: unique_euclidean_semiring_class.cong_def)
     also have "... = 1" using prod_prime_factors_minus_one .
     finally show "QuadRes p (- d')"
       unfolding Legendre_def
       by (metis one_neq_neg_one one_neq_zero)
   qed
   thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
     using \<open>n > 1\<close> three_squares_using_quadratic_residue[of n d']
     unfolding d'_def p_def
     by auto
 qed
 
 lemma three_mod_eight_power_iff:
   fixes n :: nat
   shows "(3 :: int) ^ n mod 8 = (if even n then 1 else 3)"
 proof (induction n)
   case 0
   thus ?case by auto
 next
   case (Suc n)
   thus ?case
     apply (cases "even n")
     subgoal
       using mod_mult_left_eq[of 3 8 "3 ^ n"] apply simp
       apply presburger
       done
     subgoal
       using mod_mult_left_eq[of 3 8 "3 * 3 ^ n"] apply simp
       apply presburger
       done
     done
 qed
 
 text \<open>Lemma 1.9 from @{cite nathanson1996}.\<close>
 
 lemma three_squares_using_mod_eight:
   fixes n :: nat
   assumes "n mod 8 \<in> {1, 3, 5}"
   shows "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
 proof cases
   assume "n = 1"
   hence "n = 1\<^sup>2 + 0\<^sup>2 + 0\<^sup>2" unfolding power2_eq_square by auto
   thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" by blast
 next
   assume "n \<noteq> 1"
   hence "n > 1" using assms by auto
 
   have H_n:
     "(n mod 8 = 1 \<Longrightarrow> P) \<Longrightarrow>
      (n mod 8 = 3 \<Longrightarrow> P) \<Longrightarrow>
      (n mod 8 = 5 \<Longrightarrow> P) \<Longrightarrow> P" for P
     using assms by auto
 
   define c :: nat where "c \<equiv> if n mod 8 = 3 then 1 else 3"
   have "c * n \<ge> 1" unfolding c_def using \<open>n > 1\<close> by auto
 
   obtain k where H_k: "2 * k = c * n - 1"
     using H_n
     unfolding c_def
     by (smt (verit, ccfv_threshold) dvd_mod even_mult_iff even_numeral
                                     odd_numeral odd_one odd_two_times_div_two_nat)
   have k_mod_4: "k mod 4 = (if n mod 8 = 5 then 3 else 1)" (is "k mod 4 = ?v")
   proof -
     have "c * n mod 8 = (if n mod 8 = 5 then 7 else 3)"
     using H_n
     proof cases
       assume "n mod 8 = 1"
       have "3 * n mod 8 = 3"
         using \<open>n mod 8 = 1\<close> mod_mult_right_eq[of 3 n 8]
         by auto
       thus ?thesis unfolding c_def using \<open>n mod 8 = 1\<close> by auto
     next
       assume "n mod 8 = 3"
       thus ?thesis unfolding c_def by auto
     next
       assume "n mod 8 = 5"
       have "3 * n mod 8 = 7"
         using \<open>n mod 8 = 5\<close> mod_mult_right_eq[of 3 n 8]
         by auto
       thus ?thesis unfolding c_def using \<open>n mod 8 = 5\<close> by auto
     qed
     hence "2 * k mod 8 = (if n mod 8 = 5 then 6 else 2)"
       unfolding H_k using \<open>c * n \<ge> 1\<close> mod_diff_eq_nat by simp
     hence "2 * (k mod 4) = 2 * ?v" by (simp add: mult_mod_right)
     thus ?thesis by simp
   qed
 
   have "coprime k (4 * n)"
     using k_mod_4 H_k \<open>c * n \<ge> 1\<close>
     by (metis One_nat_def coprime_Suc_left_nat coprime_commute
               coprime_diff_one_right_nat coprime_mod_left_iff
               coprime_mult_right_iff mult_2_right numeral_2_eq_2 numeral_3_eq_3
               numeral_Bit0 order_less_le_trans zero_less_one zero_neq_numeral)
   then obtain j where H_j:
       "prime (k + (4 * n) * j) \<and> j \<noteq> 0"
     using prime_linear_combination[of k "n - 1"] \<open>n > 1\<close>
     by (metis One_nat_def Suc_leI bot_nat_0.not_eq_extremum mult_is_0
               nat_1_eq_mult_iff nat_less_le prime_linear_combination
               zero_neq_numeral)
   have "j > 0" using H_j by blast
 
   define p where "p \<equiv> k + (4 * n) * j"
   have "prime p"
     unfolding p_def
     using conjunct1[OF H_j] apply (rule back_subst[of prime])
     apply (simp add: int_ops nat_plus_as_int)
     done
   have "[p = k] (mod 4 * n)"
     unfolding p_def unique_euclidean_semiring_class.cong_def by auto
 
   have "p > 2"
     using \<open>prime p\<close> \<open>[p = k] (mod 4 * n)\<close> \<open>coprime k (4 * n)\<close>
     by (metis cong_dvd_iff cong_dvd_modulus_nat coprime_common_divisor_nat
               dvd_mult2 even_numeral le_neq_implies_less odd_one prime_ge_2_nat)
 
   define d' where "d' \<equiv> 8 * j + c"
   have "d' > 1" unfolding d'_def using \<open>j > 0\<close> by simp
   have H_2_p: "2 * p = d' * n - 1"
     unfolding p_def d'_def
     using \<open>c * n \<ge> 1\<close> H_k
     by (smt (verit, del_insts) Nat.add_diff_assoc add.commute
                                add_mult_distrib mult.commute mult_2 numeral_Bit0)
 
   have "QuadRes p (- d')"
   proof -
     have d'_expansion: "d' = (\<Prod>q\<in>prime_factors d'. q ^ multiplicity q d')"
       using \<open>j > 0\<close> prime_factorization_nat unfolding d'_def by auto
 
     have "odd d'" unfolding c_def d'_def by simp
     hence d'_prime_factors_odd: "q \<in> prime_factors d' \<Longrightarrow> odd q" for q
       by fastforce
 
     have d'_prime_factors_gt_2: "q \<in> prime_factors d' \<Longrightarrow> q > 2" for q
       using d'_prime_factors_odd
       by (metis even_numeral in_prime_factors_imp_prime
                 order_less_le prime_ge_2_nat)
 
     have "[2 * p = - 1] (mod d')"
       using \<open>n > 1\<close> \<open>d' > 1\<close>
       unfolding H_2_p cong_iff_dvd_diff
       by (simp add: int_ops less_1_mult order_less_imp_not_less)
     hence d'_prime_factors_2_p_mod:
         "q \<in> prime_factors d' \<Longrightarrow> [2 * p = - 1] (mod q)" for q
       by (rule cong_dvd_modulus; auto)
 
     have "q \<in> prime_factors d' \<Longrightarrow> coprime (2 * int p) q" for q
       using d'_prime_factors_2_p_mod
       by (metis cong_imp_coprime cong_sym coprime_1_left
                 coprime_minus_left_iff mult_2 of_nat_add)
     hence d'_prime_factors_coprime:
         "q \<in> prime_factors d' \<Longrightarrow> coprime (int p) q" for q
       using d'_expansion by auto
 
     have Legendre_using_quadratic_reciprocity:
       "Legendre (- d') p =
        (\<Prod>q\<in>prime_factors d'. (Legendre p q) ^ multiplicity q d')"
     proof cases
       assume "n mod 8 \<in> {1, 3}"
 
       have "k mod 4 = 1" using \<open>n mod 8 \<in> {1, 3}\<close> k_mod_4 by auto
       hence "p mod 4 = 1"
         using \<open>[p = k] (mod 4 * n)\<close>
         by (metis unique_euclidean_semiring_class.cong_def cong_modulus_mult_nat)
       hence "[int p = 1] (mod 4)"
         by (metis cong_mod_left cong_refl int_ops(2) int_ops(3) of_nat_mod)
 
       have d'_prime_factors_Legendre:
         "q \<in> prime_factors d' \<Longrightarrow>
          Legendre p q = Legendre q p" for q
       proof -
         assume "q \<in> prime_factors d'"
         have "prime q" using \<open>q \<in> prime_factors d'\<close> by blast
         have "q > 2" using d'_prime_factors_gt_2 \<open>q \<in> prime_factors d'\<close> by blast
         show "Legendre p q = Legendre q p"
           using \<open>prime p\<close> \<open>p > 2\<close> \<open>[int p = 1] (mod 4)\<close>
                 \<open>prime q\<close> \<open>q > 2\<close> Legendre_equal[of p q]
           by auto
       qed
 
       have "Legendre (- d') p = Legendre (- 1) p * Legendre d' p"
         using \<open>prime p\<close> Legendre_mult[of p "- 1" d'] by auto
       also have "... = Legendre d' p"
         using \<open>prime p\<close> \<open>p > 2\<close> \<open>[int p = 1] (mod 4)\<close> Legendre_minus_one
         by auto
       also have "... = (\<Prod>q\<in>prime_factors d'. Legendre (q ^ multiplicity q d') p)"
         apply (subst d'_expansion)
         using \<open>prime p\<close> \<open>p > 2\<close> Legendre_prod[of p] apply auto
         done
       also have "... = (\<Prod>q\<in>prime_factors d'. (Legendre q p) ^ multiplicity q d')"
         using \<open>prime p\<close> \<open>p > 2\<close> Legendre_power by auto
       also have "... = (\<Prod>q\<in>prime_factors d'. (Legendre p q) ^ multiplicity q d')"
         using d'_prime_factors_Legendre by auto
       finally show ?thesis .
     next
       assume "n mod 8 \<notin> {1, 3}"
       hence "n mod 8 = 5" using assms by auto
 
       have "[p = 3] (mod 4)"
         using \<open>n mod 8 = 5\<close> k_mod_4 \<open>[p = k] (mod 4 * n)\<close>
         by (metis cong_mod_right cong_modulus_mult_nat)
       have "[d' = 3] (mod 8)"
         using \<open>n mod 8 = 5\<close>
         unfolding d'_def c_def cong_iff_dvd_diff
         by (simp add: unique_euclidean_semiring_class.cong_def)
       have "[d' = - 1] (mod 4)"
         using \<open>[d' = 3] (mod 8)\<close> cong_modulus_mult[of d' 3 4 2]
         unfolding unique_euclidean_semiring_class.cong_def nat_mod_as_int
         by auto
 
       have d'_prime_factors_cases:
         "q \<in> prime_factors d' \<Longrightarrow>
          multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)" for q
       proof -
         assume "q \<in> prime_factors d'"
         consider "[q = 0] (mod 4)"
                | "[q = 1] (mod 4)"
                | "[q = 2] (mod 4)"
                | "[q = 3] (mod 4)"
           unfolding unique_euclidean_semiring_class.cong_def by (simp; linarith)
         thus "multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)"
         proof cases
           assume "[q = 0] (mod 4)"
           hence False
             using d'_prime_factors_odd \<open>q \<in> prime_factors d'\<close>
             by (meson cong_0_iff dvd_trans even_numeral)
           thus "multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)"
             by blast
         next
           assume "[q = 1] (mod 4)"
           thus "multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)"
             by blast
         next
           assume "[q = 2] (mod 4)"
           hence "q = 2"
             using d'_prime_factors_odd \<open>q \<in> prime_factors d'\<close>
             by (metis unique_euclidean_semiring_class.cong_def
                       dvd_mod_iff even_numeral)
           thus "multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)"
             by (simp add: \<open>odd d'\<close> not_dvd_imp_multiplicity_0)
         next
           assume "[q = 3] (mod 4)"
           thus "multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)"
             by blast
         qed
       qed
 
       have "d' = (\<Prod>q\<in>{q\<in>prime_factors d'. True}. q ^ multiplicity q d')"
         using d'_expansion by auto
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'.
                         multiplicity q d' = 0 \<or> [q = 1] (mod 4) \<or> [q = 3] (mod 4)}.
                         q ^ multiplicity q d')"
         using d'_prime_factors_cases by meson
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. multiplicity q d' = 0} \<union>
                         {q\<in>prime_factors d'. [q = 1] (mod 4) \<or>
                         [q = 3] (mod 4)}. q ^ multiplicity q d')"
         by (rule prod.cong; blast)
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4) \<or>
                         [q = 3] (mod 4)}. q ^ multiplicity q d')"
         by (rule prod.mono_neutral_left[symmetric]; auto)
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)} \<union>
                             {q\<in>prime_factors d'. [q = 3] (mod 4)}.
                         q ^ multiplicity q d')"
         by (rule prod.cong; blast)
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                         q ^ multiplicity q d') *
                        (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                         q ^ multiplicity q d')"
         by (rule prod.union_disjoint;
             auto simp add: unique_euclidean_semiring_class.cong_def)
       finally have d'_expansion_mod_4:
         "d' = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                q ^ multiplicity q d') *
               (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                q ^ multiplicity q d')" .
 
       have "int d' mod 4 = int ((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                                  q ^ multiplicity q d') *
                                 (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                                  q ^ multiplicity q d') mod 4)"
         using d'_expansion_mod_4
         by presburger
       also have "... = ((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                          ((q mod 4) ^ multiplicity q d') mod 4) mod 4) *
                        ((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                          ((q mod 4) ^ multiplicity q d') mod 4) mod 4) mod 4"
         unfolding mod_mult_eq mod_prod_eq power_mod ..
       also have "... = int (((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                               (1 ^ multiplicity q d') mod 4) mod 4) *
                             ((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                               (3 ^ multiplicity q d') mod 4) mod 4) mod 4)"
         unfolding unique_euclidean_semiring_class.cong_def by auto
       also have "... = ((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                          (((int 3) mod 4) ^ multiplicity q d') mod 4) mod 4) mod 4"
         by (simp add: int_ops)
       also have "... = ((\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                          (((- 1) mod 4) ^ multiplicity q d') mod 4) mod 4) mod 4"
         by auto
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                         ((- 1) ^ multiplicity q d')) mod 4"
         unfolding power_mod mod_prod_eq mod_mod_trivial ..
       finally have "[d' = \<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                           ((- 1) ^ multiplicity q d')] (mod 4)"
         unfolding unique_euclidean_semiring_class.cong_def .
       hence "[\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
               ((- 1) ^ multiplicity q d') = - 1 :: int] (mod 4)"
         using \<open>[d' = - 1] (mod 4)\<close>
         unfolding unique_euclidean_semiring_class.cong_def
         by argo
       hence prod_d'_prime_factors_q_3_mod_4_minus_one:
           "(\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
             ((- 1) ^ multiplicity q d')) = (- 1 :: int)"
         unfolding power_sum[symmetric]
         unfolding minus_one_power_iff unique_euclidean_semiring_class.cong_def
         by auto
 
       have d'_prime_factors_q_1_mod_4_Legendre:
         "q \<in> prime_factors d' \<Longrightarrow>
          [q = 1] (mod 4) \<Longrightarrow>
          Legendre p q = Legendre q p" for q
       proof -
         assume "q \<in> prime_factors d'"
         assume "[q = 1] (mod 4)"
         have "prime q" using \<open>q \<in> prime_factors d'\<close> by blast
         have "q > 2" using d'_prime_factors_gt_2 \<open>q \<in> prime_factors d'\<close> by blast
         show "Legendre p q = Legendre q p"
           using \<open>prime p\<close> \<open>p > 2\<close> \<open>[p = 3] (mod 4)\<close> \<open>[q = 1] (mod 4)\<close>
                 \<open>prime q\<close> \<open>q > 2\<close> Legendre_equal[of p q]
           unfolding unique_euclidean_semiring_class.cong_def
           using zmod_int[of q 4]
           by auto
       qed
 
       have d'_prime_factors_q_3_mod_4_Legendre:
         "q \<in> prime_factors d' \<Longrightarrow>
          [q = 3] (mod 4) \<Longrightarrow>
          Legendre p q = - Legendre q p" for q
       proof -
         assume "q \<in> prime_factors d'"
         assume "[q = 3] (mod 4)"
         have "prime q" using \<open>q \<in> prime_factors d'\<close> by blast
         have "q > 2" using d'_prime_factors_gt_2 \<open>q \<in> prime_factors d'\<close> by blast
         have "p \<noteq> q"
           using d'_prime_factors_coprime[of q] \<open>q \<in> prime_factors d'\<close> \<open>prime p\<close>
           by auto
         show "Legendre p q = - Legendre q p"
           using \<open>prime p\<close> \<open>p > 2\<close> \<open>[p = 3] (mod 4)\<close> \<open>[q = 3] (mod 4)\<close>
                 \<open>prime q\<close> \<open>q > 2\<close> \<open>p \<noteq> q\<close> Legendre_opposite[of p q]
           unfolding unique_euclidean_semiring_class.cong_def
           using zmod_int[of p 4] zmod_int[of q 4]
           by fastforce
       qed
 
       have d'_prime_factors_q_0_2_mod_4:
           "q \<in> prime_factors d' \<Longrightarrow>
            ([q = 0] (mod 4) \<or> [q = 2] (mod 4)) \<Longrightarrow>
            Legendre p q = 1" for q
         unfolding unique_euclidean_semiring_class.cong_def
         using d'_prime_factors_odd mod_mod_cancel[of 2 4 q]
         by fastforce
 
       have "Legendre (- d') p = Legendre (- 1) p * Legendre d' p"
         using \<open>prime p\<close> Legendre_mult[of p "- 1" d'] by auto
       also have "... = - Legendre d' p"
         using \<open>prime p\<close> \<open>p > 2\<close> \<open>[p = 3] (mod 4)\<close> Legendre_minus_one[of p]
         unfolding unique_euclidean_semiring_class.cong_def nat_mod_as_int
         by (auto simp add: cong_0_iff Legendre_def)
       also have "... = - (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                           (Legendre q p) ^ multiplicity q d') *
                          (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                           (Legendre q p) ^ multiplicity q d')"
         apply (subst d'_expansion_mod_4)
         using \<open>prime p\<close> \<open>p > 2\<close> Legendre_mult[of p]
               Legendre_prod[of p "\<lambda>q. q ^ multiplicity q d'"] Legendre_power[of p]
         apply simp
         done
       also have "... = - (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                           (Legendre p q) ^ multiplicity q d') *
                          (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                           (- Legendre p q) ^ multiplicity q d')"
         using d'_prime_factors_q_1_mod_4_Legendre
               d'_prime_factors_q_3_mod_4_Legendre
         by auto
       also have "... = - (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                           (Legendre p q) ^ multiplicity q d') *
                          (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                           (Legendre p q * (- 1)) ^ multiplicity q d')"
         by auto
       also have "... = - (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                           (Legendre p q) ^ multiplicity q d') *
                          (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                           (Legendre p q) ^ multiplicity q d') *
                          (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                           (- 1) ^ multiplicity q d')"
         unfolding power_mult_distrib prod.distrib by auto
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d')"
         unfolding prod_d'_prime_factors_q_3_mod_4_minus_one by auto
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 0] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 1] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 2] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 3] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d')"
         using d'_prime_factors_q_0_2_mod_4 by auto
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 0] (mod 4)} \<union>
                             {q\<in>prime_factors d'. [q = 1] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 2] (mod 4)} \<union>
                             {q\<in>prime_factors d'. [q = 3] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d')"
         using prod.union_disjoint[of "{q\<in>prime_factors d'. [q = 0] (mod 4)}"
                                      "{q\<in>prime_factors d'. [q = 1] (mod 4)}"
                                      "\<lambda>q. (Legendre p (int q)) ^
                                           multiplicity q d'"]
               prod.union_disjoint[of "{q\<in>prime_factors d'. [q = 2] (mod 4)}"
                                      "{q\<in>prime_factors d'. [q = 3] (mod 4)}"
                                      "\<lambda>q. (Legendre p (int q)) ^
                                           multiplicity q d'"]
         by (fastforce simp add: unique_euclidean_semiring_class.cong_def)
       also have "... = (\<Prod>q\<in>({q\<in>prime_factors d'. [q = 0] (mod 4)} \<union>
                              {q\<in>prime_factors d'. [q = 1] (mod 4)}) \<union>
                             ({q\<in>prime_factors d'. [q = 2] (mod 4)} \<union>
                              {q\<in>prime_factors d'. [q = 3] (mod 4)}).
                         (Legendre p q) ^ multiplicity q d')"
         by (rule prod.union_disjoint[symmetric];
             auto simp add: unique_euclidean_semiring_class.cong_def)
       also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'.
                              [q = 0] (mod 4) \<or> [q = 1] (mod 4) \<or>
                              [q = 2] (mod 4) \<or> [q = 3] (mod 4)}.
                         (Legendre p q) ^ multiplicity q d')"
         by (rule prod.cong; auto)
       also have "... = (\<Prod>q\<in>prime_factors d'. (Legendre p q) ^ multiplicity q d')"
         by (rule prod.cong;
             auto simp add: unique_euclidean_semiring_class.cong_def)
       finally show ?thesis .
     qed
 
     have "q \<in> prime_factors d' \<Longrightarrow> Legendre 4 q = 1" for q
     proof -
       assume "q \<in> prime_factors d'"
       have "q dvd 4 \<Longrightarrow> q \<le> 4" by (simp add: dvd_imp_le)
       hence "q dvd 4 \<Longrightarrow> q \<in> {0, 1, 2, 3, 4}" by auto
       hence "q dvd 4 \<Longrightarrow> q \<in> {1, 2, 4}" by auto
       hence "\<not> q dvd 4" using \<open>q \<in> prime_factors d'\<close> d'_prime_factors_odd[of q]
         by (metis empty_iff even_numeral in_prime_factors_imp_prime
                   insert_iff not_prime_1)
       hence "\<not> int q dvd 4" by presburger
       thus "Legendre 4 q = 1"
         unfolding Legendre_def QuadRes_def cong_0_iff power2_eq_square
         by (metis cong_refl mult_2 numeral_Bit0)
     qed
     hence "Legendre (- d') p = (\<Prod>q\<in>prime_factors d'.
            (Legendre (2 * 2) q * Legendre p q) ^ multiplicity q d')"
       using Legendre_using_quadratic_reciprocity by auto
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       (Legendre 2 q * Legendre (2 * p) q) ^ multiplicity q d')"
       apply (rule prod.cong[OF refl])
       using d'_prime_factors_gt_2 Legendre_mult in_prime_factors_imp_prime
       by (metis int_ops(7) of_nat_numeral prime_nat_int_transfer mult.assoc)
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       (Legendre 2 q * Legendre (- 1) q) ^ multiplicity q d')"
       apply (rule prod.cong[OF refl])
       using d'_prime_factors_2_p_mod Legendre_cong
       unfolding unique_euclidean_semiring_class.cong_def
       apply metis
       done
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       ((if [q = 1] (mod 8) \<or> [q = 7] (mod 8) then 1 else - 1) *
                        (if [q = 1] (mod 4) then 1 else - 1)) ^ multiplicity q d')"
       apply (rule prod.cong[OF refl])
       subgoal for q
         apply (rule arg_cong2[of _ _ _ _ "\<lambda>a b. (a * b) ^ multiplicity q d'"])
         subgoal
           using Legendre_two_alt[of q] d'_prime_factors_gt_2[of q]
           unfolding unique_euclidean_semiring_class.cong_def nat_mod_as_int
           by force
         subgoal
           using Legendre_minus_one_alt[of q] d'_prime_factors_gt_2[of q]
           unfolding unique_euclidean_semiring_class.cong_def nat_mod_as_int
           by force
         done
       done
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       ((if [q = 5] (mod 8) \<or> [q = 7] (mod 8) then - 1 else 1)) ^
                        multiplicity q d')"
       apply (rule prod.cong)
       subgoal by blast
       subgoal for q
         apply (rule arg_cong[of _ _ "\<lambda>a. a ^ multiplicity q d'"])
         unfolding unique_euclidean_semiring_class.cong_def apply (simp; presburger)
         done
       done
     also have "... = (\<Prod>q\<in>prime_factors d'.
                       (if [q = 5] (mod 8) \<or> [q = 7] (mod 8)
                        then (- 1) ^ multiplicity q d' else 1))"
       by (rule prod.cong; auto)
     also have "... = (\<Prod>q\<in>{q\<in>prime_factors d'. [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                       (- 1) ^ multiplicity q d')"
       using prod.inter_filter[symmetric] by fast
     also have "... = (- 1) ^ (\<Sum>q\<in>{q\<in>prime_factors d'.
                                    [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                               multiplicity q d')"
       by (simp add: power_sum)
     finally have Legendre_using_sum:
       "Legendre (- d') p =
          (- 1) ^ (\<Sum>q\<in>{q\<in>prime_factors d'. [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                   multiplicity q d')" .
 
     have "[\<Sum>q\<in>{q\<in>prime_factors d'. [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
            multiplicity q d' = 0] (mod 2)"
     proof -
       have "d' = (\<Prod>q\<in>{q \<in> prime_factors d'.
                        [q = 1] (mod 8) \<or> [q = 3] (mod 8) \<or>
                        [q = 5] (mod 8) \<or> [q = 7] (mod 8)}. q ^ multiplicity q d')"
         apply (subst d'_expansion)
         apply (rule prod.cong)
         subgoal
           apply (rule Set.set_eqI)
           subgoal for q
             apply (rule iffI)
             subgoal
               using d'_prime_factors_odd[of q]
               unfolding unique_euclidean_semiring_class.cong_def
               apply simp
               apply presburger
               done
             subgoal by blast
             done
           done
         subgoal by blast
         done
       also have "... = (\<Prod>q\<in>({q \<in> prime_factors d'. [q = 1] (mod 8)} \<union>
                              {q \<in> prime_factors d'. [q = 3] (mod 8)}) \<union>
                             ({q \<in> prime_factors d'. [q = 5] (mod 8)} \<union>
                              {q \<in> prime_factors d'. [q = 7] (mod 8)}).
                         q ^ multiplicity q d')"
         by (rule prod.cong; auto)
       also have "... = (\<Prod>q\<in>({q \<in> prime_factors d'. [q = 1] (mod 8)} \<union>
                              {q \<in> prime_factors d'. [q = 3] (mod 8)}).
                         q ^ multiplicity q d') *
                        (\<Prod>q\<in>({q \<in> prime_factors d'. [q = 5] (mod 8)} \<union>
                              {q \<in> prime_factors d'. [q = 7] (mod 8)}).
                         q ^ multiplicity q d')"
         by (rule prod.union_disjoint;
             force simp add: unique_euclidean_semiring_class.cong_def)
       also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                         q ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                         q ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                         q ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                         q ^ multiplicity q d')"
         using prod.union_disjoint[of "{q\<in>prime_factors d'. [q = 1] (mod 8)}"
                                      "{q\<in>prime_factors d'. [q = 3] (mod 8)}"
                                      "\<lambda>q. q ^ multiplicity q d'"]
               prod.union_disjoint[of "{q\<in>prime_factors d'. [q = 5] (mod 8)}"
                                      "{q\<in>prime_factors d'. [q = 7] (mod 8)}"
                                      "\<lambda>q. q ^ multiplicity q d'"]
         by (force simp add: unique_euclidean_semiring_class.cong_def)
       finally have "int (d' mod 8) = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                                       q ^ multiplicity q d') *
                                      (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                                       q ^ multiplicity q d') *
                                      (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                                       q ^ multiplicity q d') *
                                      (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                                       q ^ multiplicity q d') mod 8"
         by auto
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          q ^ multiplicity q d') mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          q ^ multiplicity q d') mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          q ^ multiplicity q d') mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          q ^ multiplicity q d') mod 8) mod 8"
         by (metis (no_types, lifting) mod_mult_eq mod_mod_trivial)
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          (q ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          (q ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          (q ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          (q ^ multiplicity q d') mod 8) mod 8) mod 8"
         unfolding mod_prod_eq ..
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          ((q mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          ((q mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          ((q mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          ((q mod 8) ^ multiplicity q d') mod 8) mod 8) mod 8"
         unfolding power_mod ..
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          (((int q) mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          (((int q) mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          (((int q) mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          (((int q) mod 8) ^ multiplicity q d') mod 8) mod 8) mod 8"
         by (simp add: int_ops)
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          ((1 mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          ((3 mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          (((- 3) mod 8) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          (((- 1) mod 8) ^ multiplicity q d') mod 8) mod 8) mod 8"
         unfolding cong_int_iff[symmetric] unique_euclidean_semiring_class.cong_def
         by auto
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          (1 ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          (3 ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          ((- 3) ^ multiplicity q d') mod 8) mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          ((- 1) ^ multiplicity q d') mod 8) mod 8) mod 8"
         unfolding power_mod ..
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                          1 ^ multiplicity q d') mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          3 ^ multiplicity q d') mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          (- 3) ^ multiplicity q d') mod 8) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          (- 1) ^ multiplicity q d') mod 8) mod 8"
         unfolding mod_prod_eq ..
       also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 1] (mod 8)}.
                         1 ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                         3 ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                         (- 3) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                         (- 1) ^ multiplicity q d') mod 8"
         by (metis (no_types, lifting) mod_mult_eq mod_mod_trivial)
       also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                         3 ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                         (- 3) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                         (- 1) ^ multiplicity q d') mod 8"
         by auto
       also have "... = (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                         3 ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                         3 ^ multiplicity q d' * (- 1) ^ multiplicity q d') *
                        (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                         (- 1) ^ multiplicity q d') mod 8"
         unfolding power_mult_distrib[symmetric] by auto
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)}.
                          3 ^ multiplicity q d') *
                         (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          3 ^ multiplicity q d')) *
                        ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          (- 1) ^ multiplicity q d') *
                         (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          (- 1) ^ multiplicity q d')) mod 8"
         unfolding prod.distrib by (simp add: algebra_simps)
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'. [q = 3] (mod 8)} \<union>
                              {q \<in> prime_factors d'. [q = 5] (mod 8)}.
                          3 ^ multiplicity q d') *
                         (\<Prod>q\<in>{q \<in> prime_factors d'. [q = 5] (mod 8)} \<union>
                              {q \<in> prime_factors d'. [q = 7] (mod 8)}.
                          (- 1) ^ multiplicity q d')) mod 8"
         apply (subst
           prod.union_disjoint[of "{q\<in>prime_factors d'. [q = 5] (mod 8)}"
                                  "{q\<in>prime_factors d'. [q = 7] (mod 8)}"
                                  "\<lambda>q. (- 1) ^ multiplicity q d'"]
         )
         apply ((force simp add: unique_euclidean_semiring_class.cong_def)+)[3]
         apply (subst
           prod.union_disjoint[of "{q\<in>prime_factors d'. [q = 3] (mod 8)}"
                                  "{q\<in>prime_factors d'. [q = 5] (mod 8)}"
                                  "\<lambda>q. 3 ^ multiplicity q d'"]
         )
         apply ((force simp add: unique_euclidean_semiring_class.cong_def)+)[3]
         apply blast
         done
       also have "... = ((\<Prod>q\<in>{q \<in> prime_factors d'.
                               [q = 3] (mod 8) \<or> [q = 5] (mod 8)}.
                          3 ^ multiplicity q d') *
                         (\<Prod>q\<in>{q \<in> prime_factors d'.
                               [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                          (- 1) ^ multiplicity q d')) mod 8"
         by (rule arg_cong2[of _ _ _ _ "\<lambda>A B. ((\<Prod>q\<in>A. _ q) * (\<Prod>q\<in>B. _ q)) mod 8"];
             auto)
       also have "... = (((\<Prod>q\<in>{q \<in> prime_factors d'.
                                [q = 3] (mod 8) \<or> [q = 5] (mod 8)}.
                           3 ^ multiplicity q d') mod 8) *
                         ((\<Prod>q\<in>{q \<in> prime_factors d'.
                                [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                           (- 1) ^ multiplicity q d') mod 8)) mod 8"
         unfolding mod_mult_eq ..
       also have "... = ((3 ^ (\<Sum>q\<in>{q \<in> prime_factors d'.
                                    [q = 3] (mod 8) \<or> [q = 5] (mod 8)}.
                               multiplicity q d') mod 8) *
                         ((- 1) ^ (\<Sum>q\<in>{q \<in> prime_factors d'.
                                        [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                                   multiplicity q d') mod 8)) mod 8"
         unfolding power_sum ..
       also have "... =
           int (case ((\<Sum>q\<in>{q \<in> prime_factors d'.
                            [q = 3] (mod 8) \<or> [q = 5] (mod 8)}.
                       multiplicity q d') mod 2,
                      (\<Sum>q\<in>{q \<in> prime_factors d'.
                            [q = 5] (mod 8) \<or> [q = 7] (mod 8)}.
                       multiplicity q d') mod 2) of
                  (0    , 0    ) \<Rightarrow> 1
                | (0    , Suc 0) \<Rightarrow> 7
                | (Suc 0, 0    ) \<Rightarrow> 3
                | (Suc 0, Suc 0) \<Rightarrow> 5)" (is "_ = int ?d'_mod_8")
         unfolding three_mod_eight_power_iff minus_one_power_iff
         by (simp; simp add: odd_iff_mod_2_eq_one)
       finally have d'_mod_8: "d' mod 8 = ?d'_mod_8"by linarith
 
       have "[d' = 1] (mod 8) \<or> [d' = 3] (mod 8)"
         unfolding d'_def c_def unique_euclidean_semiring_class.cong_def
         using assms
         by auto
       hence "?d'_mod_8 = 1 \<or> ?d'_mod_8 = 3"
         unfolding unique_euclidean_semiring_class.cong_def d'_mod_8 by auto
       thus ?thesis
         unfolding unique_euclidean_semiring_class.cong_def
-        by (smt (z3) Collect_cong One_nat_def cong_exp_iff_simps(11)
+        by (smt (verit) Collect_cong One_nat_def cong_exp_iff_simps(11)
                      even_mod_2_iff even_numeral nat.case(2) numeral_eq_iff
                      numerals(1) old.nat.simps(4) parity_cases prod.simps(2)
                      semiring_norm(84))
     qed
     hence "Legendre (- d') p = 1"
       using Legendre_using_sum
       unfolding minus_one_power_iff cong_0_iff
       by argo
     thus "QuadRes p (- d')"
       unfolding Legendre_def
       by (metis one_neq_neg_one one_neq_zero)
   qed
 
   from \<open>QuadRes p (- d')\<close> obtain x\<^sub>0 y where "x\<^sub>0\<^sup>2 - (- d') = y * (int p)"
     unfolding quadratic_residue_alt_equiv[symmetric]
               quadratic_residue_alt_def
     by auto
   hence "(int p) dvd (x\<^sub>0\<^sup>2 - - d')" by simp
 
   define x :: int where "x \<equiv> if odd x\<^sub>0 then x\<^sub>0 else (x\<^sub>0 + p)"
 
   have "even (4 * int n * j)" by simp
   moreover have "odd k" using \<open>coprime k (4 * n)\<close> by auto
   ultimately have "odd (int p)" unfolding p_def by simp
 
   have "odd x" unfolding x_def using \<open>odd (int p)\<close> by auto
 
   have "QuadRes (2 * p) (- d')"
   unfolding quadratic_residue_alt_equiv[symmetric]
             quadratic_residue_alt_def
   proof -
     have "2 dvd (x\<^sup>2 - - d')" unfolding d'_def c_def using \<open>odd x\<close> by auto
     moreover from \<open>(int p) dvd (x\<^sub>0\<^sup>2 - - d')\<close>
     have "(int p) dvd (x\<^sup>2 - - d')"
       unfolding x_def power2_eq_square
       apply (simp add: algebra_simps)
       unfolding add.assoc[symmetric, of d' "x\<^sub>0 * x\<^sub>0"]
       apply auto
       done
     moreover have "coprime 2 (int p)" using \<open>odd (int p)\<close> by auto
     ultimately have "(int (2 * p)) dvd (x\<^sup>2 - - d')" by (simp add: divides_mult)
     hence "(x\<^sup>2 - - d') mod (int (2 * p)) = 0" by simp
     hence "\<exists>y. x\<^sup>2 - - d' = int (2 * p) * y" by auto
     hence "\<exists>y. x\<^sup>2 - - d' = y * int (2 * p)" by (simp add: algebra_simps)
     thus "\<exists>x y. x\<^sup>2 - - d' = y * int (2 * p)" by (rule exI[where ?x=x])
   qed
 
   have "n \<ge> 2" using \<open>1 < n\<close> by auto
   moreover have "0 < nat d'" unfolding d'_def using \<open>j > 0\<close> by simp
   moreover have "QuadRes (int (nat d' * n - 1)) (- d')"
     using \<open>d' > 1\<close> H_2_p \<open>QuadRes (2 * p) (- d')\<close> by (simp add: int_ops)
   ultimately show "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
     using three_squares_using_quadratic_residue[of n "nat d'"]
     by auto
 qed
 
 lemma power_two_mod_eight:
   fixes n :: nat
   shows "n\<^sup>2 mod 8 \<in> {0, 1, 4}"
 proof -
   have 0: "n\<^sup>2 mod 8 = (n mod 8)\<^sup>2 mod 8"
     unfolding power2_eq_square by (simp add: mod_mult_eq)
   have "n mod 8 \<in> {0, 1, 2, 3, 4, 5, 6, 7}" by auto
   hence "(n mod 8)\<^sup>2 mod 8 \<in> {0, 1, 4}"
     unfolding power2_eq_square by auto
   thus "n\<^sup>2 mod 8 \<in> {0, 1, 4}" unfolding 0 .
 qed
 
 lemma power_two_mod_four:
   fixes n :: nat
   shows "n\<^sup>2 mod 4 \<in> {0, 1}"
   using power_two_mod_eight[of n] mod_mod_cancel[of 4 8 "n\<^sup>2"] by auto
 
 text \<open>Theorem 1.4 from @{cite nathanson1996}.\<close>
 
 theorem three_squares_iff:
   fixes n :: nat
   shows "(\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2) \<longleftrightarrow> (\<nexists>a k. n = 4 ^ a * (8 * k + 7))"
 proof
   assume "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
   then obtain x\<^sub>1 x\<^sub>2 x\<^sub>3 where 0: "n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" by blast
   show "\<nexists>a k. n = 4 ^ a * (8 * k + 7)"
   proof
     assume "\<exists>a k. n = 4 ^ a * (8 * k + 7)"
     then obtain a k where 1: "n = 4 ^ a * (8 * k + 7)" by blast
     from 0 1 show False
     proof (induction a arbitrary: n x\<^sub>1 x\<^sub>2 x\<^sub>3 rule: nat.induct)
       fix n x\<^sub>1 x\<^sub>2 x\<^sub>3 :: nat
       assume 2: "n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
       assume "n = 4 ^ 0 * (8 * k + 7)"
       hence 3: "n mod 8 = 7" unfolding 1 by auto
       have "(x\<^sub>1\<^sup>2 mod 8 + x\<^sub>2\<^sup>2 mod 8 + x\<^sub>3\<^sup>2 mod 8) mod 8 = 7"
         unfolding 2 3[symmetric]
         by (meson mod_add_cong mod_mod_trivial)
       thus False
         using power_two_mod_eight[of x\<^sub>1]
               power_two_mod_eight[of x\<^sub>2]
               power_two_mod_eight[of x\<^sub>3]
         by auto
     next
       fix a' n x\<^sub>1 x\<^sub>2 x\<^sub>3 :: nat
       assume 4: "\<And>n' x\<^sub>1' x\<^sub>2' x\<^sub>3' :: nat.
                  n' = x\<^sub>1'\<^sup>2 + x\<^sub>2'\<^sup>2 + x\<^sub>3'\<^sup>2 \<Longrightarrow> n' = 4 ^ a' * (8 * k + 7) \<Longrightarrow> False"
       assume 5: "n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
       assume "n = 4 ^ Suc a' * (8 * k + 7)"
       hence "n = 4 * (4 ^ a' * (8 * k + 7))" (is "n = 4 * ?m") by auto
       hence 6: "4 * ?m = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" unfolding 5 by auto
       have "(x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2) mod 4 = 0" using 6 by presburger
       hence "((x\<^sub>1\<^sup>2 mod 4) + (x\<^sub>2\<^sup>2 mod 4) + (x\<^sub>3\<^sup>2 mod 4)) mod 4 = 0"
         by presburger
       hence "x\<^sub>1\<^sup>2 mod 4 = 0 \<and> x\<^sub>2\<^sup>2 mod 4 = 0 \<and> x\<^sub>3\<^sup>2 mod 4 = 0"
         using power_two_mod_four[of x\<^sub>1]
               power_two_mod_four[of x\<^sub>2]
               power_two_mod_four[of x\<^sub>3]
         by (auto; presburger)
       hence "even x\<^sub>1 \<and> even x\<^sub>2 \<and> even x\<^sub>3"
         by (metis dvd_0_right even_even_mod_4_iff even_power)
       hence "4 * ?m = 4 * ((x\<^sub>1 div 2)\<^sup>2 + (x\<^sub>2 div 2)\<^sup>2 + (x\<^sub>3 div 2)\<^sup>2)"
         unfolding 6 by fastforce
       hence "?m = (x\<^sub>1 div 2)\<^sup>2 + (x\<^sub>2 div 2)\<^sup>2 + (x\<^sub>3 div 2)\<^sup>2" by auto
       thus False using 4 by blast
     qed
   qed
 next
   assume 7: "\<nexists>a k. n = 4 ^ a * (8 * k + 7)"
   show "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
   proof cases
     assume "n = 0"
     thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" by auto
   next
     assume 8: "n \<noteq> 0"
     have "n > 0 \<Longrightarrow> \<exists>a m. n = 4 ^ a * m \<and> \<not> 4 dvd m"
     proof (induction n rule: less_induct)
       fix n :: nat
       assume 9: "\<And>n'. n' < n \<Longrightarrow> n' > 0 \<Longrightarrow> \<exists>a' m'. n' = 4 ^ a' * m' \<and> \<not> 4 dvd m'"
       assume 10: "n > 0"
       show "\<exists>a m. n = 4 ^ a * m \<and> \<not> 4 dvd m"
       proof cases
         assume 11: "4 dvd n"
         have "n div 4 < n" "n div 4 > 0" using 10 11 by auto
         then obtain a' m' where 12: "n div 4 = 4 ^ a' * m' \<and> \<not> 4 dvd m'"
           using 9 by blast
         have "n = 4 ^ (Suc a') * m' \<and> \<not> 4 dvd m'"
           using 11 12 by auto
         thus "\<exists>a m. n = 4 ^ a * m \<and> \<not> 4 dvd m" by blast
       next
         assume "\<not> 4 dvd n"
         hence "n = 4 ^ 0 * n \<and> \<not> 4 dvd n" by auto
         thus "\<exists>a m. n = 4 ^ a * m \<and> \<not> 4 dvd m" by blast
       qed
     qed
     then obtain a m where 13: "n = 4 ^ a * m" "\<not> 4 dvd m" using 8 by auto
     have 14: "m mod 8 \<noteq> 7"
     proof
       assume "m mod 8 = 7"
       then obtain k where "m = 8 * k + 7" by (metis div_mod_decomp mult.commute)
       hence "n = 4 ^ a * (8 * k + 7)" unfolding 13 by blast
       thus False using 7 by blast
     qed
     have "m mod 4 = 2 \<or> m mod 8 \<in> {1, 3, 5, 7}"
       using 13 by (simp; presburger)
     hence "m mod 4 = 2 \<or> m mod 8 \<in> {1, 3, 5}"
       using 14 by blast
     hence "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. m = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
       using three_squares_using_mod_four three_squares_using_mod_eight
       by blast
     hence "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = (2 ^ a * x\<^sub>1)\<^sup>2 + (2 ^ a * x\<^sub>2)\<^sup>2 + (2 ^ a * x\<^sub>3)\<^sup>2"
       unfolding 13 power2_eq_square
       unfolding mult.assoc[of "2 ^ a"]
       unfolding mult.commute[of "2 ^ a"]
       unfolding mult.assoc
       unfolding power_add[symmetric]
       unfolding mult_2[symmetric]
       unfolding power_mult
       unfolding mult.assoc[symmetric]
       unfolding add_mult_distrib[symmetric]
       unfolding mult.commute[of "4 ^ a"]
       by simp
     thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" by blast
   qed
 qed
 
 text \<open>Theorem 1.5 from @{cite nathanson1996}.\<close>
 
 theorem odd_three_squares_using_mod_eight:
   fixes n :: nat
   assumes "n mod 8 = 3"
   shows "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. odd x\<^sub>1 \<and> odd x\<^sub>2 \<and> odd x\<^sub>3 \<and> n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
 proof -
   obtain x\<^sub>1 x\<^sub>2 x\<^sub>3 where 0: "n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
     using assms three_squares_using_mod_eight by blast
   have "(x\<^sub>1\<^sup>2 mod 8 + x\<^sub>2\<^sup>2 mod 8 + x\<^sub>3\<^sup>2 mod 8) mod 8 = 3"
     unfolding 0 assms[symmetric]
     by (meson mod_add_cong mod_mod_trivial)
   hence "x\<^sub>1\<^sup>2 mod 8 = 1 \<and> x\<^sub>2\<^sup>2 mod 8 = 1 \<and> x\<^sub>3\<^sup>2 mod 8 = 1"
     using power_two_mod_eight[of x\<^sub>1]
           power_two_mod_eight[of x\<^sub>2]
           power_two_mod_eight[of x\<^sub>3]
     by auto
   hence "odd x\<^sub>1 \<and> odd x\<^sub>2 \<and> odd x\<^sub>3"
     by (metis dvd_mod even_numeral even_power odd_one pos2)
   hence "odd x\<^sub>1 \<and> odd x\<^sub>2 \<and> odd x\<^sub>3 \<and> n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" using 0 by blast
   thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. odd x\<^sub>1 \<and> odd x\<^sub>2 \<and> odd x\<^sub>3 \<and> n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" by blast
 qed
 
 subsection \<open>Consequences\<close>
 
 lemma four_decomposition:
   fixes n :: nat
   shows "\<exists>x y z. n = x\<^sup>2 + y\<^sup>2 + z\<^sup>2 + z"
 proof -
   have "(4 * n + 1) mod 8 \<in> {1, 3, 5}" by (simp; presburger)
   then obtain x y z where 0: "4 * n + 1 = x\<^sup>2 + y\<^sup>2 + z\<^sup>2"
     using three_squares_using_mod_eight by blast
   hence 1: "1 = (x\<^sup>2 + y\<^sup>2 + z\<^sup>2) mod 4"
     by (metis Suc_0_mod_numeral(2) Suc_eq_plus1 mod_add_left_eq
               mod_mult_self1_is_0)
   show ?thesis
   proof cases
     assume "even x"
     then obtain x' where H_x: "x = 2 * x'" by blast
     show ?thesis
     proof cases
       assume "even y"
       then obtain y' where H_y: "y = 2 * y'" by blast
       show ?thesis
       proof cases
         assume "even z"
         then obtain z' where H_z: "z = 2 * z'" by blast
         show ?thesis using 1 unfolding H_x H_y H_z by auto
       next
         assume "odd z"
         then obtain z' where H_z: "z = 2 * z' + 1" using oddE by blast
         have "n = x'\<^sup>2 + y'\<^sup>2 + z'\<^sup>2 + z'"
           using 0 unfolding H_x H_y H_z power2_eq_square by auto
         thus ?thesis by blast
       qed
     next
       assume "odd y"
       then obtain y' where H_y: "y = 2 * y' + 1" using oddE by blast
       show ?thesis
       proof cases
         assume "even z"
         then obtain z' where H_z: "z = 2 * z'" by blast
         have "n = x'\<^sup>2 + z'\<^sup>2 + y'\<^sup>2 + y'"
           using 0 unfolding H_x H_y H_z power2_eq_square by auto
         thus ?thesis by blast
       next
         assume "odd z"
         then obtain z' where H_z: "z = 2 * z' + 1" using oddE by blast
         show ?thesis
           using 1
           unfolding H_x H_y H_z power2_eq_square
           by (metis dvd_mod even_add even_mult_iff even_numeral odd_one)
       qed
     qed
   next
     assume "odd x"
     then obtain x' where H_x: "x = 2 * x' + 1" using oddE by blast
     show ?thesis
     proof cases
       assume "even y"
       then obtain y' where H_y: "y = 2 * y'" by blast
       show ?thesis
       proof cases
         assume "even z"
         then obtain z' where H_z: "z = 2 * z'" by blast
         have "n = y'\<^sup>2 + z'\<^sup>2 + x'\<^sup>2 + x'"
           using 0 unfolding H_x H_y H_z power2_eq_square by auto
         thus ?thesis by blast
       next
         assume "odd z"
         then obtain z' where H_z: "z = 2 * z' + 1" using oddE by blast
         show ?thesis
           using 1
           unfolding H_x H_y H_z power2_eq_square
           by (metis dvd_mod even_add even_mult_iff even_numeral odd_one)
       qed
     next
       assume "odd y"
       then obtain y' where H_y: "y = 2 * y' + 1" using oddE by blast
       show ?thesis
       proof cases
         assume "even z"
         then obtain z' where H_z: "z = 2 * z'" by blast
         show ?thesis
           using 1
           unfolding H_x H_y H_z power2_eq_square
           by (metis dvd_mod even_add even_mult_iff even_numeral odd_one)
       next
         assume "odd z"
         then obtain z' where H_z: "z = 2 * z' + 1" using oddE by blast
         show ?thesis
           using 1
           unfolding H_x H_y H_z power2_eq_square
           by (simp add: mod_add_eq[symmetric])
       qed
     qed
   qed
 qed
 
 theorem four_decomposition_int:
   fixes n :: int
   shows "(\<exists>x y z. n = x\<^sup>2 + y\<^sup>2 + z\<^sup>2 + z) \<longleftrightarrow> n \<ge> 0"
 proof
   assume "\<exists>x y z. n = x\<^sup>2 + y\<^sup>2 + z\<^sup>2 + z"
   then obtain x y z where 0: "n = x\<^sup>2 + y\<^sup>2 + z\<^sup>2 + z" by blast
   show "n \<ge> 0"
     unfolding 0 power2_eq_square
     by (smt (verit) div_pos_neg_trivial mult_le_0_iff
                     nonzero_mult_div_cancel_right sum_squares_ge_zero)
 next
   assume "n \<ge> 0"
   thus "\<exists>x y z. n = x\<^sup>2 + y\<^sup>2 + z\<^sup>2 + z"
     using four_decomposition[of "nat n"]
     by (smt (verit) int_eq_iff int_plus of_nat_power)
 qed
 
 theorem four_squares:
   fixes n :: nat
   shows "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3 x\<^sub>4. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2 + x\<^sub>4\<^sup>2"
 proof cases
   assume "\<exists>a k. n = 4 ^ a * (8 * k + 7)"
   then obtain a k where "n = 4 ^ a * (8 * k + 7)" by blast
   hence 0: "n = 4 ^ a * (8 * k + 6) + (2 ^ a)\<^sup>2"
     by (metis add_mult_distrib left_add_mult_distrib mult.commute mult_numeral_1
         numeral_Bit0 numeral_plus_numeral power2_eq_square power_mult_distrib
         semiring_norm(5))
   have "\<nexists>a' k'. 4 ^ a * (8 * k + 6) = 4 ^ a' * (8 * k' + 7)"
   proof
     assume "\<exists>a' k'. 4 ^ a * (8 * k + 6) = 4 ^ a' * (8 * k' + 7)"
     then obtain a' k' where 1: "4 ^ a * (8 * k + 6) = 4 ^ a' * (8 * k' + 7)"
       by blast
     show False
     proof (cases rule: linorder_cases[of a a'])
       assume "a < a'"
       hence 2: "a' = a + (a' - a)" "a' - a > 0" by auto
       have 3: "4 ^ a * (8 * k + 6) = 4 ^ a * 4 ^ (a' - a) * (8 * k' + 7)"
         using 1 2 by (metis power_add)
       have "2 = (8 * k + 6) mod 4" by presburger
       also have "... = (4 ^ (a' - a) * (8 * k' + 7)) mod 4" using 3 by auto
       also have "... = 0" using 2 by auto
       finally show False by auto
     next
       assume "a = a'"
       hence "8 * k + 6 = 8 * k' + 7" using 1 by auto
       thus False by presburger
     next
       assume "a > a'"
       hence 4: "a = a' + (a - a')" "a - a' > 0" by auto
       have 5: "4 ^ a' * 4 ^ (a - a') * (8 * k + 6) = 4 ^ a' * (8 * k' + 7)"
         using 1 4 by (metis power_add)
       have "0 = (4 ^ (a - a') * (8 * k + 6)) mod 4" using 4 by auto
       also have "... = (8 * k' + 7) mod 4" using 5 by auto
       also have "... = 3" by presburger
       finally show False by auto
     qed
   qed
   then obtain x\<^sub>1 x\<^sub>2 x\<^sub>3 where "4 ^ a * (8 * k + 6) = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2"
     using three_squares_iff by blast
   thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3 x\<^sub>4. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2 + x\<^sub>4\<^sup>2" unfolding 0 by auto
 next
   assume "\<nexists>a k. n = 4 ^ a * (8 * k + 7)"
   hence "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2" using three_squares_iff by blast
   thus "\<exists>x\<^sub>1 x\<^sub>2 x\<^sub>3 x\<^sub>4. n = x\<^sub>1\<^sup>2 + x\<^sub>2\<^sup>2 + x\<^sub>3\<^sup>2 + x\<^sub>4\<^sup>2" by (metis zero_power2 add_0)
 qed
 
 end