# Theory TwoSquares

```(*  Title:      TwoSquares.thy
Authors:    Roelof Oosterhuis, 2007  Rijksuniversiteit Groningen
Jaime Mendizabal Roche, 2016
*)

theory TwoSquares
imports
"HOL-Number_Theory.Number_Theory"
begin

context

fixes sum2sq_nat :: "nat ⇒ nat ⇒ nat"
defines "sum2sq_nat a b ≡ a^2+b^2"

fixes is_sum2sq_nat :: "nat ⇒ bool"
defines "is_sum2sq_nat n ≡ (∃ a b. n = sum2sq_nat a b)"

begin

(* TODO: move this lemma to the distribution (?). (It is used also in QuadForm and FourSquares) *)
private lemma best_division_abs: "(n::int) > 0 ⟹ ∃ k. 2 * ¦a - k*n¦ ≤ n"
proof -
assume a: "n > 0"
define k where "k = a div n"
hence h: "a - k * n = a mod n" by (simp add: mod_div_mult_eq algebra_simps)
thus ?thesis
proof (cases "2 * (a mod n) ≤ n")
case True
hence "2 * ¦a - k*n¦ ≤ n" using h pos_mod_sign a by auto
thus ?thesis by blast
next
case False
hence "2 * (n - a mod n) ≤ n" by auto
have "a - (k+1)*n = a mod n - n" using h by (simp add: algebra_simps)
hence "2 * ¦a - (k+1)*n¦ ≤ n" using h pos_mod_bound[of n a] a False by fastforce
thus ?thesis by blast
qed
qed

private definition
sum2sq_int :: "int × int ⇒ int" where
"sum2sq_int = (λ(a,b). a^2+b^2)"

private definition
is_sum2sq_int :: "int ⇒ bool" where
"is_sum2sq_int n ⟷ (∃ a b. n = sum2sq_int(a,b))"

private lemma sum2sq_int_nat_eq: "sum2sq_nat a b = sum2sq_int (a, b)"
unfolding sum2sq_nat_def sum2sq_int_def by simp

private lemma is_sum2sq_int_nat_eq: "is_sum2sq_nat n = is_sum2sq_int (int n)"
unfolding is_sum2sq_nat_def is_sum2sq_int_def
proof
assume "∃a b. n = sum2sq_nat a b"
thus "∃a b. int n = sum2sq_int (a, b)" using sum2sq_int_nat_eq by force
next
assume "∃a b. int n = sum2sq_int (a, b)"
then obtain a b where "int n = sum2sq_int (a, b)" by blast
hence "int n = sum2sq_int (int (nat ¦a¦), int (nat ¦b¦))" unfolding sum2sq_int_def by simp
hence "int n = int (sum2sq_nat (nat ¦a¦) (nat ¦b¦))" using sum2sq_int_nat_eq by presburger
thus "∃a b. n = sum2sq_nat a b" by auto
qed

private lemma product_two_squares_aux: "sum2sq_int(a, b) * sum2sq_int(c, d) = sum2sq_int(a*c - b*d, a*d + b*c)"
unfolding power2_eq_square sum2sq_int_def by (simp add: algebra_simps)

private lemma product_two_squares_int: "is_sum2sq_int m ⟹ is_sum2sq_int n ⟹ is_sum2sq_int (m*n)"
by (unfold is_sum2sq_int_def, auto simp only: product_two_squares_aux, blast)

private lemma product_two_squares_nat: "is_sum2sq_nat m ⟹ is_sum2sq_nat n ⟹ is_sum2sq_nat (m*n)"
using product_two_squares_int is_sum2sq_int_nat_eq by simp

private lemma sots1_aux:
assumes "prime (4*k+3)"
assumes "odd (multiplicity (4*k+3) n)"
shows "¬ is_sum2sq_nat n"
proof
assume "is_sum2sq_nat n"
then obtain a b where h1: "n = a^2 + b^2"  unfolding is_sum2sq_nat_def sum2sq_nat_def by blast
have ab_nz: "a ≠ 0 ∨ b ≠ 0" by (rule ccontr) (insert assms, auto simp: h1)
let ?p = "4*k+3"
let ?g = "gcd a b"
have h2: "?g ≠ 0" using assms(2) h1 odd_pos by fastforce
then obtain a' b' where h3: "a = a' * ?g" "b = b' * ?g" "coprime a' b'"
using gcd_coprime_exists by blast
have "?g^2 dvd n" using dvd_add h1 by auto
then obtain m where h4: "m * ?g^2 = n" using dvd_div_mult_self by blast
also have "… = (a' * ?g)^2 + (b' * ?g)^2" unfolding h1 using h3 by presburger
also have "… = ?g^2 * a'^2 + ?g^2 * b'^2" unfolding power2_eq_square by simp
finally have "?g^2 * m = ?g^2 * (a'^2 + b'^2)" by (simp add: distrib_left mult.commute)
hence h5: "m = a'^2 + b'^2" using h2 by auto
let ?mul = "multiplicity ?p ?g"
have "multiplicity ?p (?g^2) = ?mul + ?mul"
unfolding power2_eq_square using h2 assms
by (subst prime_elem_multiplicity_mult_distrib) simp_all
hence "even (multiplicity ?p (?g^2))" by simp
moreover have "m ≠ 0" using assms(2) h4 odd_pos by fastforce
ultimately have "odd (multiplicity ?p m)"
using assms ab_nz by (simp_all add: h4 [symmetric] prime_elem_multiplicity_mult_distrib)
hence "?p dvd m" using not_dvd_imp_multiplicity_0 by force
hence h6: "?p dvd a'^2+b'^2" using h5 by auto
{
assume "?p dvd a'^2"
moreover hence "?p dvd b'^2" using h6 dvd_add_right_iff by blast
ultimately have "?p dvd a'" "?p dvd b'" using assms(1) prime_dvd_power_nat by blast+
hence "False"
using assms(1) h3(3) coprime_common_divisor_nat[of a' b' ?p] not_prime_1 by linarith
}
hence "¬ (?p dvd a'^2)" ..
hence h7: "¬ (?p dvd a')" using assms(1)
hence "coprime ?p a'"
using assms(1) by (simp add: prime_imp_coprime)
thm prime_imp_coprime_nat
moreover have "a' ≠ 0" using h7 dvd_0_right[of ?p] by meson
ultimately obtain ainv aux where "a' * ainv = ?p * aux + 1"
using bezout_nat[of "a'" ?p]
by (auto simp: ac_simps)
hence "[a' * ainv = 1] (mod ?p)" using cong_to_1'_nat by auto
from cong_mult [OF this this] have h11: "[1 = ainv^2 * a'^2] (mod ?p)"
unfolding power2_eq_square by (simp add: algebra_simps cong_sym)
let ?bdiva = "ainv * b'"
have "[ainv^2 * (a'^2 + b'^2) = 0] (mod ?p)"
using h6 cong_dvd_modulus_nat cong_mult_self_right by blast
from cong_add [OF h11 this] have "[1 + ainv^2 * b'^2 = 0] (mod ?p)"
by fastforce
hence h8: "[?bdiva^2 + 1 = 0] (mod ?p)" by (simp add: power_mult_distrib)
{
assume "?p dvd ?bdiva"
hence "?p dvd (?bdiva^2)" by (simp add: assms(1) prime_dvd_power_nat_iff)
hence "[?bdiva^2 = 0] (mod ?p)" using cong_altdef_nat by auto
hence "[?bdiva^2 +1 = 1] (mod ?p)" using cong_add_rcancel_0_nat by blast
from this h8 have "[0 = 1] (mod ?p)" using cong_sym cong_trans by blast
hence "?p dvd 1" using cong_0_1_nat by auto
hence "False" using assms(1) by simp
}
hence "¬ (?p dvd ?bdiva)" ..
hence h9: "[?bdiva^(?p-1) = 1] (mod ?p)"
using assms(1) fermat_theorem [of ?p ?bdiva] by simp
have h10: "?p≥3" by simp
have h11: "[?bdiva^(4*k+2) = 1] (mod ?p)" using h9 by auto
have "[(?bdiva^2 + 1)^2 = 0] (mod ?p)" using h8 cong_pow [of "?bdiva^2 + 1" 0 ?p 2] by auto
moreover have "?bdiva ^ 4 = (?bdiva ^ 2) ^ 2" by auto
hence "(?bdiva^2 + 1)^2 = ?bdiva^4 + ?bdiva^2 + ?bdiva^2 + 1"
by (auto simp: algebra_simps power2_eq_square)
ultimately have "[?bdiva^4 + ?bdiva^2 + ?bdiva^2 + 1 = 0] (mod ?p)" by simp
moreover have "[?bdiva^4 + ?bdiva^2 + (?bdiva^2 + 1) = ?bdiva^4 + ?bdiva^2 + 0] (mod ?p)"
ultimately have "[?bdiva^4 + ?bdiva^2 = 0] (mod ?p)" by (simp add: cong_def)
hence "[?bdiva^4 + ?bdiva^2 + 1 = 0 + 1] (mod ?p)" using cong_add_rcancel_nat by blast
moreover have "[?bdiva^4 + (?bdiva^2 + 1) = ?bdiva^4 + 0] (mod ?p)"
ultimately have "[?bdiva^4 = 1] (mod ?p)" by (simp add: cong_def)
hence "[(?bdiva^4)^k = 1^k] (mod ?p)" using cong_pow by blast
hence h12: "[?bdiva^(4*k) = 1] (mod ?p)" by (simp add: power_mult)
hence h13: "[?bdiva^(4*k)*(?bdiva^2 + 1) = 1*(?bdiva^2 + 1)] (mod ?p)"
using cong_scalar_right by blast
have "?bdiva^(4*k)*(?bdiva^2 + 1) = ?bdiva^(4*k+2)+?bdiva^(4*k)"
hence "[?bdiva^(4*k+2)+?bdiva^(4*k) = ?bdiva^2 + 1] (mod ?p)"
using h13 unfolding nat_mult_1 by presburger
moreover have "[?bdiva^(4*k+2) + ?bdiva^(4*k) = 1 + 1] (mod ?p)"
using h11 h12 cong_add by blast
ultimately have "[?bdiva^2 + 1 = 2] (mod ?p)"
hence "[0 = 2] (mod ?p)" using h8 by (simp add: cong_def)
then have "?p dvd 2" by (auto dest: cong_dvd_iff)
then show False
by (auto dest: dvd_imp_le)
qed

private lemma sots1: assumes "is_sum2sq_nat n"
shows "⋀ k. prime (4*k+3) ⟶ even (multiplicity (4*k+3) n)"
using sots1_aux assms by blast

private lemma aux_lemma: assumes "[(a::nat) = b] (mod c)" "b < c"
shows "∃ k. a = c*k + b"
proof -
have "a mod c = b" using assms by (simp add: cong_def mod_if)
hence "b ≤ a" using assms by auto
thus ?thesis using cong_le_nat assms(1) by auto
qed

private lemma Legendre_1mod4: "prime (4*k+1::nat) ⟹ (Legendre (-1) (4*k+1)) = 1"
proof -
let ?p = "4*k+1"
let ?L = "Legendre (-1) ?p"
assume p: "prime ?p"
from p have "k ≠ 0" by (intro notI) simp_all
hence p2: "?p > 2" by simp
with p have "[?L = (-1)^((?p - 1) div 2)] (mod ?p)"
by (rule euler_criterion)
hence "[?L = (-1)^(2 * nat k)] (mod ?p)" by auto
hence "[?L = 1] (mod ?p)" unfolding power_mult by simp
hence "?p dvd 1-?L"
using cong_iff_dvd_diff dvd_minus_iff[of ?p "?L-1"] by auto
moreover have "?L=1 ∨ ?L=0 ∨ ?L=-1" by (simp add: Legendre_def)
ultimately have "?L = 1 ∨ ?p dvd 1 ∨ ?p dvd (2::int)" by auto
moreover
{ assume "?p dvd 1 ∨ ?p dvd (2::int)"
with p2 have False by (auto simp add: zdvd_not_zless) }
ultimately show ?thesis by auto
qed

private lemma qf1_prime_exists: "prime (4*k+1) ⟹ is_sum2sq_nat (4*k+1)"
proof -
let ?p = "4*k+1"
assume p: "prime ?p"
hence "Legendre (-1) ?p = 1" by (rule Legendre_1mod4)
moreover
{ assume "¬ QuadRes ?p (-1)"
hence "Legendre (-1) ?p ≠ 1" by (unfold Legendre_def, auto) }
ultimately have "QuadRes ?p (-1)" by auto
then obtain s1 where s1: "[s1^2 = -1] (mod ?p)" by (auto simp add: QuadRes_def)
hence s1': "[s1^2 + 1 = 0] (mod ?p)" by (simp add: cong_iff_dvd_diff)
define s2 where "s2 = nat ¦s1¦"
hence "int (s2^2 + 1) = s1^2 + 1" by auto
with s1' have "[int (s2^2 + 1) = 0] (mod ?p)" by presburger
hence s2: "[s2^2 + 1 = 0] (mod ?p)"
using cong_int_iff by fastforce
from p have p0: "?p > 0" by simp
then obtain s where s0p: "0 ≤ s ∧ s < ?p ∧ [s2 = s] (mod ?p)"
using cong_less_unique_nat[of ?p] by fastforce
then have "[s^2 = s2^2] (mod ?p)"
with s2 have s: "[s^2 + 1 = 0] (mod ?p)"
hence "?p dvd s^2 + 1" using cong_altdef_nat by auto
then obtain t where t: "s^2 + 1 = ?p*t" by (auto simp add: dvd_def)
hence "?p*t = sum2sq_nat s 1" by (simp add: sum2sq_nat_def)
hence qf1pt: "is_sum2sq_nat (?p*t)" by (auto simp add: is_sum2sq_nat_def)
have t_l_p: "t < ?p"
proof (rule ccontr)
assume "¬ t < ?p"
hence "t > ?p - 1" by simp
with p0 have "?p*(?p - 1) < ?p*t" by (simp only: mult_less_mono2)
also with t have "… = s^2 + 1" by simp
also have "… ≤ ?p*(?p - 1) - ?p + 2"
proof -
from s0p have "s ≤ ?p - 1" by (auto simp add: less_le)
with s0p have "s^2 ≤ (?p - 1)^2" by (simp only: power_mono)
also have "… = ?p*(?p - 1) - 1*(?p - 1)" by (simp only: power2_eq_square diff_mult_distrib)
finally show ?thesis by auto
qed
finally have "?p < 2" by simp
with p show False by (unfold prime_def, auto)
qed
have tpos: "t ≥ 1"
proof (rule ccontr)
assume "¬ t ≥ 1"
hence "t < 1" by auto
moreover
{ assume "t = 0" with t have "s^2 + 1 = 0" by simp }
moreover
{ assume "t < 0"
with p0 have "?p*t < ?p*0" by (simp only: zmult_zless_mono2)
with t have "s^2 + 1 < 0" by auto }
moreover have "s^2 ≥ 0" by (simp only: zero_le_power2)
ultimately show False by (auto simp add: less_le)
qed
moreover
{ assume t1: "t > 0"
then obtain tn where tn: "tn = t - 1" by auto
have "is_sum2sq_nat (?p*(1+ 0))" (is "?Q 0")
― ‹So, \$Q~n =\$ there exist \$x,y\$ such that \$x^2+y^2 =(p*(1+ int(n)))\$›
proof (rule ccontr)
assume nQ1: "¬ ?Q 0"
have "(1 + tn) < ?p ⟹ ¬ ?Q tn"
proof (induct tn rule: infinite_descent0)
case 0
from nQ1 show "1+ 0 < ?p ⟹ ¬ ?Q 0" by simp
next
case (smaller n)
hence n0: "n > 0" and IH: "1+ n < ?p ∧ ?Q n" by auto
then obtain x y where "x^2 + y^2 = int (?p*(1+ n))"
using is_sum2sq_int_nat_eq by (unfold is_sum2sq_int_def sum2sq_int_def, auto)
hence xy: "x^2 + y^2 = (int ?p)*(int (1+ n))" unfolding of_nat_mult by presburger
let ?n1 = "int (1 + n)"
from n0 have n1pos: "?n1 > 0" by simp
then obtain r v where rv: "v = x-r*?n1 ∧ 2*¦v¦ ≤ ?n1"
by (frule_tac n="?n1" in best_division_abs, auto)
from n1pos obtain s w where sw: "w = y-s*?n1 ∧ 2*¦w¦ ≤ ?n1"
by (frule_tac n="?n1" in best_division_abs, auto)
let ?C = "v^2 + w^2"
have "?n1 dvd ?C"
proof
from rv sw have "?C = (x-r*?n1)^2 + (y-s*?n1)^2" by simp
also have "… = x^2 + y^2 - 2*x*(r*?n1) - 2*y*(s*?n1) + (r*?n1)^2 + (s*?n1)^2"
unfolding power2_eq_square by (simp add: algebra_simps)
also with xy have "… = ?n1*?p - ?n1*(2*x*r) - ?n1*(2*y*s) + ?n1^2*r^2 + ?n1^2*s^2"
by (simp only: ac_simps power_mult_distrib)
finally show "?C = ?n1*(?p - 2*x*r - 2*y*s + ?n1*(r^2 + s^2))"
by (simp only: power_mult_distrib distrib_left ac_simps
left_diff_distrib right_diff_distrib power2_eq_square)
qed
then obtain m1 where m1: "?C = ?n1*m1" by (auto simp add: dvd_def)
have mn: "m1 < ?n1"
proof (rule ccontr)
assume "¬ m1 < ?n1" hence "?n1-m1 ≤ 0" by simp
hence "4*?n1 - 4*m1 ≤ 0" by simp
with n1pos have "2*?n1 - 4*m1 < 0" by simp
with n1pos have "?n1*(2*?n1 - 4*m1) < ?n1*0" by (simp only: zmult_zless_mono2)
hence contr: "?n1*(2*?n1- 4*m1) < 0" by simp
have hlp: "2*¦v¦ ≥ 0 ∧ 2*¦w¦ ≥ 0" by simp
from m1 have "4*?n1*m1 = 4*v^2 + 4*w^2" by arith
also have "… = (2*¦v¦)^2 + (2*¦w¦)^2"
also from rv hlp have "… ≤ ?n1^2 + (2*¦w¦)^2"
using power_mono [of "2*¦b¦" "1 + int n" 2 for b] by auto
also from sw hlp have "… ≤ ?n1^2 + ?n1^2"
using power_mono [of "2*¦b¦" "1 + int n" 2 for b] by auto
finally have "?n1*m1*4 ≤ ?n1*?n1*2" by (simp add: power2_eq_square ac_simps)
hence "?n1*(2*?n1- 4*m1) ≥ 0" by (simp only: right_diff_distrib ac_simps)
with contr show False by auto
qed
have "?p*m1 = (r*v + s*w + m1)^2 + (r*w - s*v)^2"
proof -
from m1 xy have "(?p*?n1)*?C = (x^2+y^2)*(v^2+w^2)" by simp
also have "… = (x*v + y*w)^2 + (x*w - y*v)^2"
also with rv sw have "… = ((r*?n1+v)*v + (s*?n1+w)*w)^2 + ((r*?n1+v)*w - (s*?n1+w)*v)^2"
by simp
also have "… = (?n1*(r*v) + ?n1*(s*w) + (v^2+w^2))^2 + (?n1*(r*w) - ?n1*(s*v))^2"
also from m1 have "… = (?n1*(r*v) + ?n1*(s*w) + ?n1*m1)^2 + (?n1*(r*w) - ?n1*(s*v))^2"
by simp
finally have "(?p*?n1)*?C = ?n1^2*(r*v + s*w + m1)^2 + ?n1^2*(r*w - s*v)^2"
with m1 have "?n1^2*(?p*m1) = ?n1^2*((r*v + s*w + m1)^2 + (r*w - s*v)^2)"
by (simp only: ac_simps power2_eq_square, simp add: distrib_left)
hence "?n1^2*(?p*m1 - (r*v+s*w+m1)^2 - (r*w-s*v)^2) = 0"
by (auto simp add: distrib_left right_diff_distrib)
moreover from n1pos have "?n1^2 ≠ 0" by (simp add: power2_eq_square)
ultimately show ?thesis by simp
qed
hence qf1pm1: "is_sum2sq_int ((int ?p)*m1)"
by (unfold is_sum2sq_int_def sum2sq_int_def, auto)
have m1pos: "m1 > 0"
proof -
{ assume "v^2 + w^2 = 0"
hence "v = 0 ∧ w = 0" using sum_power2_eq_zero_iff by blast
with rv sw have "?n1 dvd x ∧ ?n1 dvd y" by (unfold dvd_def, auto)
hence "?n1^2 dvd x^2 ∧ ?n1^2 dvd y^2" by simp
hence "?n1^2 dvd x^2 + y^2" by (simp only: dvd_add)
with xy have "?n1*?n1 dvd ?n1*?p" by (simp only: power2_eq_square ac_simps)
moreover from n1pos have "?n1 ≠ 0" by simp
ultimately have "?n1 dvd ?p" by (rule zdvd_mult_cancel)
with n1pos have "?n1 ≥ 0 ∧ ?n1 dvd ?p" by simp
with p have "?n1 = 1 ∨ ?n1 = ?p" unfolding prime_nat_iff by presburger
with IH have "?Q 0" by auto
with nQ1 have False by simp }
moreover
{ assume "v^2 + 1*w^2 ≠ 0"
moreover have "v^2 + w^2 ≥ 0" by simp
ultimately have vwpos: "v^2 + w^2 > 0" by arith
with m1 have "m1 ≠ 0" by auto
moreover have "m1 ≥ 0"
proof (rule ccontr)
assume "¬ m1 ≥ 0"
hence "m1 < 0" by simp
with n1pos have "?n1*m1 < ?n1*0" by (simp only: zmult_zless_mono2)
with m1 vwpos show False by simp
qed
ultimately have ?thesis by auto }
ultimately show ?thesis by auto
qed
hence "1 + int((nat m1) - 1) = m1" by arith
with qf1pm1 have Qm1: "?Q ((nat m1) - 1)"
using is_sum2sq_int_nat_eq by (simp add: algebra_simps)
then obtain mm where tmp: "mm = (nat m1) - 1 ∧ ?Q mm" by auto
moreover have "mm<n" using tmp mn m1pos by arith
moreover with IH have "1 + int mm < ?p" by auto
ultimately show ?case by auto
qed
hence "¬ is_sum2sq_nat (?p*t)" using tn tpos t_l_p by auto
with qf1pt show False by simp
qed
hence ?thesis by auto }
ultimately show ?thesis by (auto simp add: less_le)
qed

private lemma fermat_two_squares: assumes "prime p" "(¬ [p = 3] (mod 4))"
shows "is_sum2sq_nat p"
proof (cases "p=2")
case True
have "(2::nat)=1^2+1^2" using power2_eq_square by simp
thus ?thesis unfolding is_sum2sq_nat_def sum2sq_nat_def using True by fast
next
case False
hence "p > 2" using assms(1) unfolding prime_nat_iff by auto
hence h1: "odd p" using assms(1) prime_odd_nat by simp
hence h2: "¬ [p = 0] (mod 4)" unfolding cong_def by fastforce
have h3: "¬ [p = 2] (mod 4)" using h1 cong_dvd_iff [of p 2 2] cong_dvd_modulus_nat by auto
obtain x where h4: "[p = x] (mod 4) ∧ x<4" by (meson cong_less_unique_nat zero_less_numeral)
from h1 h2 h3 h4 assms have "x≠0 ∧ x≠2 ∧ x≠3 ∧ x<4" by meson
hence "x=1" by linarith
from this h4 have "[p = 1] (mod 4)" by simp
then obtain k where "p = 4*k+1" using aux_lemma by fastforce
thus ?thesis using qf1_prime_exists assms by blast
qed

private lemma sots2: assumes "⋀ k. prime (4*k+3) ⟶ even (multiplicity (4*k+3) n)"
shows "is_sum2sq_nat n" using assms
proof (induction n rule: nat_less_induct)
case (1 n)
thus ?case
proof (cases "n>1")
case f: False
thus ?thesis
proof (cases "n=1")
case True
have "(1::nat)=0^2+1^2" by (simp add: power2_eq_square)
thus ?thesis using True unfolding is_sum2sq_nat_def sum2sq_nat_def by blast
next
case False
hence "n=0" using f by simp
moreover have "(0::nat)=0^2+0^2" by (simp add: power2_eq_square)
ultimately show ?thesis unfolding is_sum2sq_nat_def sum2sq_nat_def by blast
qed
next
case True
then obtain p m where h1: "prime p ∧ n = p * m" using prime_divisor_exists[of n]
by (auto elim: dvdE)
with True have m_nz: "m ≠ 0" by (intro notI) auto
from h1 have h2: "m<n" using n_less_m_mult_n[of m p] prime_gt_Suc_0_nat[of p] True by linarith
{
assume a1: "[p = 3] (mod 4)"
then obtain kp where "p = 4*kp+3" using aux_lemma by fastforce
hence "even (multiplicity p n)" using "1.prems" h1 by auto
moreover have "multiplicity p n ≠ 0" using h1 True m_nz
by (subst multiplicity_eq_zero_iff) (auto simp: prime_gt_0_nat)
ultimately have h3: "multiplicity p n ≥ 2" by presburger
have "p dvd m"
proof (rule ccontr)
assume a2: "¬ p dvd m"
hence "multiplicity p m = 0" by (rule not_dvd_imp_multiplicity_0)
moreover from h1 have "multiplicity p p = 1" by (intro multiplicity_prime) auto
moreover have "m > 0" using h1 True by (cases "m = 0") simp_all
ultimately have "multiplicity p n = 1" using h1
using prime_elem_multiplicity_mult_distrib [of p p m] m_nz prime_gt_0_nat
by auto
thus "False" using h3 by simp
qed
then obtain m' where h4: "m = p * m'" using dvdE by blast
with h1 have h5: "n = p^2 * m'" by (simp add: power2_eq_square)
have h6: "m'<n"
using dual_order.strict_trans h1 h2 h4 nat_mult_less_cancel1 prime_gt_0_nat[of p] by blast
have "⋀ kq. prime (4*kq + 3) ⟹ even (multiplicity (4*kq + 3) m')"
proof -
fix kq::nat
let ?q = "4*kq + 3"
assume a2: "prime ?q"
{
assume p: "p=?q"
hence h7: "multiplicity ?q (p^2) = 2" using h1
by (auto intro!: multiplicity_prime_power)
have "even (multiplicity ?q n)" using 1(2)[of kq] a2 by blast
also note h5
also from p h1 h4 m_nz
have "multiplicity (4 * kq + 3) (p^2 * m') =
Suc (Suc (multiplicity (4 * kq + 3) m'))"
by (subst prime_elem_multiplicity_mult_distrib) auto
finally have "even (multiplicity ?q m')" by simp
}
moreover {
assume "p≠?q"
from a2 h4 m_nz have "multiplicity ?q n =
multiplicity (4 * kq + 3) (p⇧2) + multiplicity (4 * kq + 3) m'"
unfolding h5 by (subst prime_elem_multiplicity_mult_distrib) simp_all
also from ‹p ≠ ?q› a2 h1 have "multiplicity ?q (p^2) = 0"
by (intro multiplicity_distinct_prime_power) simp_all
finally have "multiplicity ?q n = multiplicity ?q m'" by simp
moreover have "even (multiplicity ?q n)" using 1(2)[of kq] a2 by blast
ultimately have "even (multiplicity ?q m')" by simp
}
ultimately show "even (multiplicity ?q m')" by blast
qed
hence "is_sum2sq_nat m'" by (simp add: 1 h6)
moreover have "p^2 = p^2 + 0^2" by simp
hence "is_sum2sq_nat (p^2)" unfolding is_sum2sq_nat_def sum2sq_nat_def by blast
ultimately have ?thesis using product_two_squares_nat h5 by blast
} moreover
{
assume a1: "¬ [p = 3] (mod 4)"
have "⋀ kq. prime (4*kq+3) ⟹ even (multiplicity (4*kq+3) m)"
proof -
fix kq
let ?q = "4*(kq::nat) + 3"
assume a2: "prime ?q"
{ assume "p = ?q"
then have False using a1 cong_add_rcancel_0_nat [of "4 * kq" 3 4]
}
hence "p≠?q" ..

have "n = p * m" using h1 by simp
also from h1 a2 m_nz have "multiplicity ?q … =
multiplicity (4 * kq + 3) p + multiplicity (4 * kq + 3) m"
by (subst prime_elem_multiplicity_mult_distrib) (simp_all add: prime_gt_0_nat)
also from ‹p ≠ ?q› a2 h1 have "multiplicity ?q p = 0"
by (intro prime_multiplicity_other) simp_all
finally have "multiplicity ?q n = multiplicity ?q m" by simp
moreover have "even (multiplicity ?q n)" using 1(2)[of kq] a2 by blast
ultimately show "even (multiplicity ?q m)" by simp
qed
hence "is_sum2sq_nat m" by (simp add: 1 h2)
moreover have "is_sum2sq_nat p" using fermat_two_squares a1 h1 by blast
ultimately have ?thesis using product_two_squares_nat h1 by blast
} ultimately
show ?thesis by blast
qed
qed

theorem sum_of_two_squares:
"is_sum2sq_nat n ⟷ (∀ k. prime (4*k+3) ⟶ even (multiplicity (4*k+3) n))"
using sots1[of n] sots2[of n] by blast

private lemma k_mod_eq: "(∀p::nat. prime p ∧ [p = 3] (mod 4) ⟶ P p) = (∀k. prime (4*k+3) ⟶ P (4*k+3))"
proof
assume a1: "∀p. prime p ∧ [p = 3] (mod 4) ⟶ P p"
{
fix k :: nat
assume "prime (4 * k + 3)"
moreover hence "[4*k+3 = 3] (mod 4)"
ultimately have "P (4 * k + 3)" using a1 by blast
}
thus "∀k. prime (4 * k + 3) ⟶ P (4 * k + 3)" by blast
next
assume a1: "∀k. prime (4 * k + 3) ⟶ P (4 * k + 3)"
{
fix p :: nat
assume "prime p" "[p = 3] (mod 4)"
moreover with aux_lemma obtain k where "p = 4*k+3" by fastforce
ultimately have "P p" using a1 by blast
}
thus "∀p. prime p ∧ [p = 3] (mod 4) ⟶ P p" by blast
qed

theorem sum_of_two_squares':
"is_sum2sq_nat n ⟷ (∀ p. prime p ∧ [p = 3] (mod 4) ⟶ even (multiplicity p n))"
using sum_of_two_squares k_mod_eq by presburger

theorem sum_of_two_squares_prime: assumes "prime p"
shows "is_sum2sq_nat p = [p≠3] (mod 4)"
proof (cases "[p=3] (mod 4)")
case True
have "odd (multiplicity p p)" using assms by simp
hence "¬ (is_sum2sq_nat p)" using assms True sum_of_two_squares' by blast
with True show ?thesis by simp