Theory FourSquares

(*  Title:      FourSquares.thy
Author:     Roelof Oosterhuis, 2007  Rijksuniversiteit Groningen
*)

section ‹Lagrange's four-square theorem›

theory FourSquares
imports "HOL-Number_Theory.Number_Theory"
begin

context

fixes sum4sq_nat :: "nat ⇒ nat ⇒ nat ⇒ nat ⇒ nat"
defines "sum4sq_nat a b c d ≡ a^2+b^2+c^2+d^2"

fixes is_sum4sq_nat :: "nat ⇒ bool"
defines "is_sum4sq_nat n ≡ (∃ a b c d. n = sum4sq_nat a b c d)"

begin

(* TODO: move this lemma to the distribution (?). (It is used also in QuadForm and TwoSquares) *)
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"
have h: "a - k * n = a mod n" by (simp add: div_mult_mod_eq algebra_simps k_def)
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

text ‹Shows that all nonnegative integers can be written as the sum of four squares. The proof consists of the following steps:
\begin{itemize}\item For every prime $p=2n+1$ the two sets of residue classes $$\{x^2 \bmod p\mid 0\le x\le n\} ~{\rm and}~ \{-1-y^2 \bmod p \mid 0 \le y \le n\}$$ both contain $n+1$ different elements and therefore they must have at least one element in common.
\item Hence there exist $x,y$ such that $x^2+y^2+1^2+0^2$ is a multiple of $p$.
\item The next step is to show, by an infinite descent, that $p$ itself can be written as the sum of four squares.
\item Finally, using the multiplicity of this form, the same holds for all positive numbers.
\end{itemize}›

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

private definition
is_sum4sq_int :: "int ⇒ bool" where
"is_sum4sq_int n ⟷ (∃ a b c d. n = sum4sq_int(a,b,c,d))"

private lemma mult_sum4sq_int: "sum4sq_int(a,b,c,d) * sum4sq_int(p,q,r,s) =
sum4sq_int(a*p+b*q+c*r+d*s, a*q-b*p-c*s+d*r,
a*r+b*s-c*p-d*q, a*s-b*r+c*q-d*p)"
by (unfold sum4sq_int_def, simp add: eval_nat_numeral field_simps)

private lemma sum4sq_int_nat_eq: "sum4sq_nat a b c d = sum4sq_int (a, b, c, d)"
unfolding sum4sq_nat_def sum4sq_int_def by simp

private lemma is_sum4sq_int_nat_eq: "is_sum4sq_nat n = is_sum4sq_int (int n)"
unfolding is_sum4sq_nat_def is_sum4sq_int_def
proof
assume "∃a b c d. n = sum4sq_nat a b c d"
thus "∃a b c d. int n = sum4sq_int (a, b, c, d)" using sum4sq_int_nat_eq by force
next
assume "∃a b c d. int n = sum4sq_int (a, b, c, d)"
then obtain a b c d where "int n = sum4sq_int (a, b, c, d)" by blast
hence "int n = sum4sq_int (int (nat ¦a¦), int (nat ¦b¦), int (nat ¦c¦), int (nat ¦d¦))"
unfolding sum4sq_int_def by simp
hence "int n = int (sum4sq_nat (nat ¦a¦) (nat ¦b¦) (nat ¦c¦) (nat ¦d¦))"
using sum4sq_int_nat_eq by presburger
thus "∃a b c d. n = sum4sq_nat a b c d" by auto
qed

private lemma is_mult_sum4sq_int: "is_sum4sq_int x ⟹ is_sum4sq_int y ⟹ is_sum4sq_int (x*y)"
by (unfold is_sum4sq_int_def, auto simp only: mult_sum4sq_int, blast)

private lemma is_mult_sum4sq_nat: "is_sum4sq_nat x ⟹ is_sum4sq_nat y ⟹ is_sum4sq_nat (x*y)"
using is_mult_sum4sq_int is_sum4sq_int_nat_eq by simp

private lemma mult_oddprime_is_sum4sq: "⟦ prime (nat p); odd p ⟧ ⟹
∃ t. 0<t ∧ t<p ∧ is_sum4sq_int (p*t)"
proof -
assume p1: "prime (nat p)"
then have p0: "p > 1" and "prime p"
assume p2: "odd p"
then obtain n where n: "p = 2*n+1" using oddE by blast
with p1 have n0: "n > 0" by (auto simp add: prime_nat_iff)
let ?C = "{0 ..< p}"
let ?D = "{0 .. n}"
let ?f = "%x. x^2 mod p"
let ?g = "%x. (-1-x^2) mod p"
let ?A = "?f  ?D"
let ?B = "?g  ?D"
have finC: "finite ?C" by simp
have finD: "finite ?D" by simp
from p0 have AsubC: "?A ⊆ ?C" and BsubC: "?B ⊆ ?C"
by auto
with finC have finA: "finite ?A" and finB: "finite ?B"
from AsubC BsubC have AunBsubC: "?A ∪ ?B ⊆ ?C" by (rule Un_least)
from p0 have cardC: "card ?C = nat p" using card_atLeastZeroLessThan_int by blast
from n0 have cardD: "card ?D = 1+ nat n" by simp
have cardA: "card ?A = card ?D"
proof -
have "inj_on ?f ?D"
proof (unfold inj_on_def, auto)
fix x fix y
assume x0: "0 ≤ x" and xn: "x ≤ n" and y0: "0 ≤ y" and yn: " y ≤ n"
and xyp: "x^2 mod p = y^2 mod p"
with p0 have "[x^2 = y^2] (mod p)" using cong_def by blast
hence "p dvd x^2-y^2" using cong_iff_dvd_diff by blast
hence "p dvd (x+y)*(x-y)" by (simp add: power2_eq_square algebra_simps)
hence "p dvd x+y ∨ p dvd x-y" using ‹prime p› p0
by (auto dest: prime_dvd_multD)
moreover
{ assume "p dvd x+y"
moreover from xn yn n have "x+y < p" by auto
ultimately have "¬ x+y > 0" by (auto simp add: zdvd_not_zless)
with x0 y0 have "x = y" by auto } ― ‹both are zero›
moreover
{ assume ass: "p dvd x-y"
have "x = y"
proof (rule ccontr, case_tac "x-y ≥ 0")
assume "x-y ≥ 0" and "x ≠ y" hence "x-y > 0" by auto
with ass have "¬ x-y < p" by (auto simp add: zdvd_not_zless)
with xn y0 n p0 show "False" by auto
next
assume "¬ 0 ≤ x-y" hence "y-x > 0" by auto
moreover from x0 yn n p0 have "y-x < p" by auto
ultimately have "¬ p dvd y-x" by (auto simp add: zdvd_not_zless)
moreover from ass have "p dvd -(x-y)" by (simp only: dvd_minus_iff)
ultimately show "False" by auto
qed }
ultimately show "x=y" by auto
qed
with finD show ?thesis by (simp only: inj_on_iff_eq_card)
qed
have cardB: "card ?B = card ?D"
proof -
have "inj_on ?g ?D"
proof (unfold inj_on_def, auto)
fix x fix y
assume x0: "0 ≤ x" and xn: "x ≤ n" and y0: "0 ≤ y" and yn: " y ≤ n"
and xyp: "(-1-x^2) mod p = (-1-y^2) mod p"
with p0 have "[-1-y^2 = -1-x^2] (mod p)" by (simp only: cong_def)
hence "p dvd (-1-y^2) - (-1-x^2)" by (simp only: cong_iff_dvd_diff)
moreover have "-1-y^2 - (-1-x^2) = x^2 - y^2" by arith
ultimately have "p dvd x^2-y^2" by simp
hence "p dvd (x+y)*(x-y)" by (simp add: power2_eq_square algebra_simps)
with p1 have "p dvd x+y ∨ p dvd x-y" using ‹prime p› p0
by (auto dest: prime_dvd_multD)
moreover
{ assume "p dvd x+y"
moreover from xn yn n have "x+y < p" by auto
ultimately have "¬ x+y > 0" by (auto simp add: zdvd_not_zless)
with x0 y0 have "x = y" by auto } ― ‹both are zero›
moreover
{ assume ass: "p dvd x-y"
have "x = y"
proof (rule ccontr, case_tac "x-y ≥ 0")
assume "x-y ≥ 0" and "x ≠ y" hence "x-y > 0" by auto
with ass have "¬ x-y < p" by (auto simp add: zdvd_not_zless)
with xn y0 n p0 show "False" by auto
next
assume "¬ 0 ≤ x-y" hence "y-x > 0" by auto
moreover from x0 yn n p0 have "y-x < p" by auto
ultimately have "¬ p dvd y-x" by (auto simp add: zdvd_not_zless)
moreover from ass have "p dvd -(x-y)" by (simp only: dvd_minus_iff)
ultimately show "False" by auto
qed }
ultimately show "x=y" by auto
qed
with finD show ?thesis by (simp only: inj_on_iff_eq_card)
qed
have "?A ∩ ?B ≠ {}"
proof (rule ccontr, auto)
assume ABdisj: "?A ∩ ?B = {}"
from cardA cardB cardD have "2 + 2*(nat n) = card ?A + card ?B" by auto
also with finA finB ABdisj have "… = card (?A ∪ ?B)"
by (simp only: card_Un_disjoint)
also with finC AunBsubC have "… ≤ card ?C" by (simp only: card_mono)
also with cardC have "… = nat p" by simp
finally have "2 + 2*(nat n) ≤ nat p" by simp
with n show "False" by arith
qed
then obtain z where "z ∈ ?A ∧ z ∈ ?B" by auto
then obtain x y where xy: "x ∈ ?D ∧ y ∈ ?D ∧ z = x^2 mod p ∧ z = (-1-y^2) mod p" by blast
with p0 have "[x^2=-1-y^2](mod p)" by (simp add: cong_def)
hence "p dvd x^2-(-1-y^2)" by (simp only: cong_iff_dvd_diff)
moreover have "x^2-(-1-y^2)=x^2+y^2+1" by arith
ultimately have "p dvd sum4sq_int(x,y,1,0)" by (auto simp add: sum4sq_int_def)
then obtain t where t: "p * t = sum4sq_int(x,y,1,0)" by (auto simp only: dvd_def eq_refl)
hence "is_sum4sq_int (p*t)" by (unfold is_sum4sq_int_def, auto)
moreover have "t > 0 ∧ t < p"
proof
have "x^2 ≥ 0 ∧ y^2 ≥ 0" by simp
hence "x^2+y^2+1 > 0" by arith
with t have "p*t > 0" by (unfold sum4sq_int_def, auto)
moreover
{ assume "t < 0" with p0 have "p*t < p*0" by (simp only: zmult_zless_mono2)
hence "p*t < 0" by simp }
moreover
{ assume "t = 0" hence "p*t = 0" by simp }
ultimately have "¬ t < 0 ∧ t ≠ 0" by auto
thus "t > 0" by simp
from xy have "x^2 ≤ n^2 ∧ y^2 ≤ n^2" by (auto simp add: power_mono)
hence "x^2+y^2+1 ≤ 2*n^2 + 1" by auto
with t have contr: "p*t ≤ 2*n^2+1" by (simp add: sum4sq_int_def)
moreover
{ assume "t > n+1"
with p0 have "p*(n+1) < p*t" by (simp only: zmult_zless_mono2)
with n have "p*t > (2*n+1)*n + (2*n+1)*1" by (simp only: distrib_left)
hence "p*t > 2*n*n + n + 2*n + 1" by (simp only: distrib_right mult_1_left)
with n0 have "p*t > 2*n^2 + 1" by (simp add: power2_eq_square) }
ultimately have "¬ t > n+1" by auto
with n0 n show "t < p" by auto
qed
ultimately show ?thesis by blast
qed

private lemma zprime_is_sum4sq: "prime (nat p) ⟹ is_sum4sq_int p"
proof (cases)
assume p2: "p=2"
hence "p = sum4sq_int(1,1,0,0)" by (auto simp add: sum4sq_int_def)
thus ?thesis by (auto simp add: is_sum4sq_int_def)
next
assume "¬ p =2" and prp: "prime (nat p)"
hence "¬ nat p = 2" by simp
with prp have "2 < nat p" using prime_nat_iff by force
moreover with prp have "odd (nat p)" using prime_odd_nat[of "nat p"] by blast
ultimately have "odd p" by (simp add: even_nat_iff)
with prp have "∃ t. 0<t ∧ t<p ∧ is_sum4sq_int (p*t)" by (rule mult_oddprime_is_sum4sq)
then obtain a b c d t where pt_sol: "0<t ∧ t<p ∧ p*t = sum4sq_int(a,b,c,d)"
by (unfold is_sum4sq_int_def, blast)
hence Qt: "0<t ∧ t<p ∧ (∃ a1 a2 a3 a4. p*t = sum4sq_int(a1,a2,a3,a4))"
(is "?Q t") by blast
have "?Q 1"
proof (rule ccontr)
assume nQ1: "¬ ?Q 1"
have "¬ ?Q t"
proof (induct t rule: infinite_descent0_measure[where V="λx. (nat x)- 1"], clarify)
fix x a b c d
assume "nat x - 1 = 0" and "x > 0" and s: "p*x=sum4sq_int(a,b,c,d)" and "x < p"
moreover hence "x = 1" by arith
ultimately have "?Q 1" by auto
with nQ1 show False by auto
next
fix x
assume "0 < nat x - 1" and "¬ ¬ ?Q x"
then obtain a1 a2 a3 a4 where ass: "1<x ∧ x<p ∧ p*x = sum4sq_int(a1,a2,a3,a4)"
by auto
have "∃ y. nat y - 1 < nat x - 1 ∧ ?Q y"
proof (cases)
assume evx: "even x"
hence "even (x*p)" by simp
with ass have ev1234: "even (a1^2+a2^2 + a3^2+a4^2)"
by (auto simp add: sum4sq_int_def ac_simps)
have "∃ b1 b2 b3 b4. p*x=sum4sq_int(b1,b2,b3,b4) ∧ even (b1+b2) ∧ even (b3+b4)"
proof (cases)
assume ev12: "even (a1^2+a2^2)"
with ev1234 ass show ?thesis by auto
next
assume "¬ even (a1^2+a2^2)"
hence odd12: "odd (a1^2+a2^2)" by simp
with ev1234 have odd34: "odd (a3^2+a4^2)" by auto
show ?thesis
proof (cases)
assume ev1: "even (a1^2)"
with odd12 have odd2: "odd (a2^2)" by simp
show ?thesis
proof (cases)
assume "even (a3^2)"
moreover from ass have "p*x = sum4sq_int(a1,a3,a2,a4)" by (auto simp add: sum4sq_int_def)
ultimately show ?thesis using odd2 odd34 ev1 by auto
next
assume "¬ even (a3^2)"
moreover from ass have "p*x = sum4sq_int(a1,a4,a2,a3)" by (auto simp add: sum4sq_int_def)
ultimately show ?thesis using odd34 odd2 ev1 by auto
qed
next
assume odd1: "¬ even (a1^2)"
with odd12 have ev2: "even (a2^2)" by simp
show ?thesis
proof (cases)
assume "even (a3^2)"
moreover from ass have "sum4sq_int(a1,a4,a2,a3)=p*x" by (auto simp add: sum4sq_int_def)
ultimately show ?thesis using odd34 odd1 ev2 by force
next
assume "¬ even (a3^2)"
moreover from ass have "sum4sq_int(a1,a3,a2,a4)=p*x" by (auto simp add: sum4sq_int_def)
ultimately show ?thesis using odd34 odd1 ev2 by force
qed
qed
qed
then obtain b1 b2 b3 b4
where b: "p*x = sum4sq_int(b1,b2,b3,b4) ∧ even (b1+b2) ∧ even (b3+b4)" by auto
then obtain c1 c3 where c13: "b1+b2 = 2*c1 ∧ b3+b4 = 2*c3"
using evenE[of "b1+b2"] evenE[of "b3+b4"] by meson
from b have "even (b1-b2) ∧ even (b3-b4)" by simp
then obtain c2 c4 where c24: "b1-b2 = 2*c2 ∧ b3-b4 = 2*c4"
using evenE[of "b1-b2"] evenE[of "b3-b4"] by meson
from evx obtain y where y: "x = 2*y" using evenE by blast
hence "4*(p*y) = 2*(p*x)" by (simp add: ac_simps)
also from b have "… = 2*b1^2 + 2*b2^2 + 2*b3^2 + 2*b4^2"
by (auto simp: sum4sq_int_def)
also have "… = (b1 + b2)^2 + (b1 - b2)^2 + (b3 + b4)^2 + (b3 - b4)^2"
by (auto simp add: power2_eq_square algebra_simps)
also with c13 c24 have "… = 4*(c1^2 + c2^2 + c3^2 + c4^2)"
finally have "p * y = sum4sq_int(c1,c2,c3,c4)" by (auto simp add: sum4sq_int_def)
moreover from y ass have "0 < y ∧ y < p ∧ (nat y) - 1 < (nat x) - 1" by arith
ultimately show ?thesis by blast
next
assume xodd: "¬ even x"
with ass have "∃ c1 c2 c3 c4. 2*¦a1-c1*x¦≤x ∧ 2*¦a2-c2*x¦≤x ∧ 2*¦a3-c3*x¦≤x ∧ 2*¦a4-c4*x¦≤x"
then obtain b1 c1 b2 c2 b3 c3 b4 c4 where
bc_def: "b1 = a1-c1*x ∧ b2 = a2-c2*x ∧ b3 = a3-c3*x ∧ b4 = a4-c4*x"
and "2*¦b1¦≤x ∧ 2*¦b2¦≤x ∧ 2*¦b3¦≤x ∧ 2*¦b4¦≤x"
by blast
moreover have "2*¦b1¦≠x ∧ 2*¦b2¦≠x ∧ 2*¦b3¦≠x ∧ 2*¦b4¦≠x" using xodd by fastforce
ultimately have bc_abs: "2*¦b1¦<x ∧ 2*¦b2¦<x ∧ 2*¦b3¦<x ∧ 2*¦b4¦<x" by auto
let ?B = "b1^2 + b2^2 + b3^2 + b4^2"
let ?C = "c1^2 + c2^2 + c3^2 + c4^2"
have "x dvd ?B"
proof
from bc_def ass have
"?B = p*x - 2*(a1*c1+a2*c2+a3*c3+a4*c4)*x + ?C*x^2"
unfolding sum4sq_int_def by (auto simp add: power2_eq_square algebra_simps)
thus "?B = x*(p - 2*(a1*c1+a2*c2+a3*c3+a4*c4) + ?C*x)"
by (auto simp add: ac_simps power2_eq_square
distrib_left right_diff_distrib)
qed
then obtain y where y: "?B = x * y" by (auto simp add: dvd_def)
let ?A1 = "a1*b1 + a2*b2 + a3*b3 + a4*b4"
let ?A2 = "a1*b2 - a2*b1 - a3*b4 + a4*b3"
let ?A3 = "a1*b3 + a2*b4 - a3*b1 - a4*b2"
let ?A4 = "a1*b4 - a2*b3 + a3*b2 - a4*b1"
let ?A = "sum4sq_int(?A1,?A2,?A3,?A4)"
have "x dvd ?A1 ∧ x dvd ?A2 ∧ x dvd ?A3 ∧ x dvd ?A4"
proof (safe)
from bc_def have
"?A1 = (b1+c1*x)*b1 + (b2+c2*x)*b2 + (b3+c3*x)*b3 + (b4+c4*x)*b4"
by simp
also with y have "… = x*(y + c1*b1 + c2*b2 + c3*b3 + c4*b4)"
by (auto simp add: distrib_left power2_eq_square ac_simps)
finally show "x dvd ?A1" by auto
from bc_def have
"?A2 = (b1+c1*x)*b2 - (b2+c2*x)*b1 - (b3+c3*x)*b4 + (b4+c4*x)*b3"
by simp
also have "… = x*(c1*b2 - c2*b1 - c3*b4 + c4*b3)"
by (auto simp add: distrib_left right_diff_distrib ac_simps)
finally show "x dvd ?A2" by auto
from bc_def have
"?A3 = (b1+c1*x)*b3 + (b2+c2*x)*b4 - (b3+c3*x)*b1 - (b4+c4*x)*b2"
by simp
also have "… = x*(c1*b3 + c2*b4 - c3*b1 - c4*b2)"
by (auto simp add: distrib_left right_diff_distrib ac_simps)
finally show "x dvd ?A3" by auto
from bc_def have
"?A4 = (b1+c1*x)*b4 - (b2+c2*x)*b3 + (b3+c3*x)*b2 - (b4+c4*x)*b1"
by simp
also have "… = x*(c1*b4 - c2*b3 + c3*b2 - c4*b1)"
by (auto simp add: distrib_left right_diff_distrib ac_simps)
finally show "x dvd ?A4" by auto
qed
then obtain d1 d2 d3 d4 where d:
"?A1=x*d1 ∧ ?A2=x*d2 ∧ ?A3=x*d3 ∧ ?A4=x*d4"
let ?D = "sum4sq_int(d1,d2,d3,d4)"
from d have "x^2*?D = ?A"
by (auto simp only: sum4sq_int_def power_mult_distrib distrib_left)
also have "… = sum4sq_int(a1,a2,a3,a4)*sum4sq_int(b1,b2,b3,b4)"
by (simp only: mult_sum4sq_int)
also with y ass have "… = (p*x)*(x*y)" by (auto simp add: sum4sq_int_def)
also have "… = x^2*(p*y)" by (simp only: power2_eq_square ac_simps)
finally have "x^2*(?D - p*y) = 0" by (auto simp add: right_diff_distrib)
with ass have "p*y = ?D" by auto
moreover have y_l_x: "y < x"
proof -
have "4*b1^2 = (2*¦b1¦)^2 ∧ 4*b2^2 = (2*¦b2¦)^2 ∧
4*b3^2 = (2*¦b3¦)^2 ∧ 4*b4^2 = (2*¦b4¦)^2" by simp
with bc_abs have "4*b1^2<x^2 ∧ 4*b2^2<x^2 ∧ 4*b3^2<x^2 ∧ 4*b4^2<x^2"
using power_strict_mono [of "2*¦b¦" x 2 for b]
by auto
hence "?B < x^2" by auto
with y have "x*(x-y) > 0"
by (auto simp add: power2_eq_square right_diff_distrib)
moreover from ass have "x > 0" by simp
ultimately show ?thesis using zero_less_mult_pos by fastforce
qed
moreover have "y > 0"
proof -
have b2pos: "b1^2 ≥ 0 ∧ b2^2 ≥ 0 ∧ b3^2 ≥ 0 ∧ b4^2 ≥ 0" by simp
hence "?B = 0 ∨ ?B > 0" by arith
moreover
{ assume "?B = 0"
moreover from b2pos have
"?B-b1^2 ≥ 0 ∧ ?B-b2^2 ≥ 0 ∧ ?B-b3^2 ≥ 0 ∧ ?B-b4^2 ≥ 0" by arith
ultimately have "b1^2 ≤ 0 ∧ b2^2 ≤ 0 ∧ b3^2 ≤ 0 ∧ b4^2 ≤ 0" by auto
with b2pos have  "b1^2 = 0 ∧ b2^2 = 0 ∧ b3^2 = 0 ∧ b4^2 = 0" by arith
hence "b1 = 0 ∧ b2 = 0 ∧ b3 = 0 ∧ b4 = 0" by auto
with bc_def have "x dvd a1 ∧ x dvd a2 ∧ x dvd a3 ∧ x dvd a4"
by auto
hence "x^2 dvd a1^2 ∧ x^2 dvd a2^2 ∧ x^2 dvd a3^2 ∧ x^2 dvd a4^2" by simp
hence "x^2 dvd a1^2+a2^2+a3^2+a4^2" by (simp only: dvd_add)
with ass have "x^2 dvd p*x" by (auto simp only: sum4sq_int_def)
hence "x*x dvd x*p" by (simp only: power2_eq_square ac_simps)
with ass have "nat x dvd nat p"
moreover from ass prp have "x ≥ 0 ∧ x ≠ 1 ∧ x ≠ p ∧ prime (nat p)" by simp
ultimately have "False" unfolding prime_nat_iff by auto }
moreover
{ assume "?B > 0"
with y have "x*y > 0" by simp
moreover from ass have "x > 0" by simp
ultimately have ?thesis using zero_less_mult_pos by blast }
ultimately show ?thesis by auto
qed
moreover with y_l_x have "(nat y) - 1 < (nat x) - 1" by arith
moreover from y_l_x ass have "y < p" by auto
ultimately show ?thesis by blast
qed
thus "∃ y. nat y - 1 < nat x - 1 ∧ ¬ ¬ ?Q y" by blast
qed
with Qt show "False" by simp
qed
thus "is_sum4sq_int p" by (auto simp add: is_sum4sq_int_def)
qed

private lemma prime_is_sum4sq: "prime p ⟹ is_sum4sq_nat p"
using zprime_is_sum4sq is_sum4sq_int_nat_eq by simp

theorem sum_of_four_squares: "is_sum4sq_nat n"
proof (induction n rule: nat_less_induct)
case (1 n)
show ?case
proof (cases "n>1")
case False
hence "n = 0 ∨ n = 1" by auto
moreover have "0 = sum4sq_nat 0 0 0 0" "1 = sum4sq_nat 1 0 0 0" unfolding sum4sq_nat_def by auto
ultimately show ?thesis unfolding is_sum4sq_nat_def by blast
next
case True
then obtain p m where dec: "prime p ∧ n = p * m" using prime_factor_nat[of n]
by (auto elim: dvdE)
moreover hence "m<n" using n_less_m_mult_n[of m p] prime_gt_Suc_0_nat[of p] True by linarith
ultimately have "is_sum4sq_nat m" "is_sum4sq_nat p" using 1 prime_is_sum4sq by blast+
thus ?thesis using dec is_mult_sum4sq_nat by blast
qed
qed

end

end