Library Coqtail.Arith.Lagrange_four_square

Require Import ZArith Omega Znumtheory.
Require Import Natsets MyNat Ztools Zeqm Znumfacts.

Definition of the predicate : being the sum of four squares

Definition foursquare (n : Z) : Type :=
{a : _ & { b: _ & { c : _ & { d |
n = a × a + b × b + c × c + d × d
}}}}.

Euler's identity and compatibility of the 4-square property wrt. multiplication

Lemma euler's_identity : a b c d w x y z : Z,
let s z := z × z in
(a × a + b × b + c × c + d × d) × (w × w + x × x + y × y + z × z) =
s (a×w + b×x + c×y + d×z) +
s (a×x - b×w - c×z + d×y) +
s (a×y + b×z - c×w - d×x) +
s (a×z - b×y + c×x - d×w).
Proof.
intros; unfold s.
ring.
Qed.

Lemma mult_foursquare_compat : n m : Z,
foursquare nfoursquare mfoursquare (n × m)%Z.
Proof.
intros n m (a, (b, (c, (d, Hn)))) (x, (y, (z, (t, Hm)))); subst.
rewrite euler's_identity.
repeat eexists.
Defined.

Lemma Zmod_sqrt_eq_compat : p i j, prime p
(0 i → 0 j → 2 × i < p → 2 × j < p
i × i j × j [p]i = j)%Z.
Proof.
intros p i j Pp Pi Pj Bi Bj ES.
destruct (Z_eq_dec ((i + j) mod p) 0) as [ E | E ].
cut (i + j = 0); [ omega | ].
cut (i + j < p); [ intros H | omega ].
rewrite <- E; symmetry; apply Zmod_small.
omega.

destruct (modular_inverse _ _ Pp E) as (a, Ha).

cut (i j [p]).
intros M; red in M.
rewrite (Zmod_small i p), (Zmod_small j p) in M; omega.

apply eqm_minus_0.
replace (i - j) with ((i - j) × 1) by ring.
rewrite <- Ha.
apply eqm_minus_0 in ES.
replace 0 with (0 × a) by ring.
rewrite <- ES.
apply eq_eqm; ring.
Qed.

Lemma prime_odd : p, 2 pprime pp mod 2 = 1.
Proof.
intros p N2 Pp.
rem (p mod 2) r Er.
pose proof Z_mod_lt p 2 as help_omega.
cut (r 0).
omega.

clear help_omega; intros Rn; subst.
apply Zmod_divide in Rn; [ | omega ].
refine (N2 (prime_div_prime _ _ _ _ _)); auto.
apply prime_2.
Qed.

Lemma odd_bound_1 : p i, p mod 2 = 1 → i < Zsucc (p / 2) → 2 × i < p.
Proof.
intros p i Op Bi.
unfold Zsucc in Bi.
apply Zle_lt_trans with (2 × (p / 2)).
omega.

rewrite Z_mult_div_mod; [ | auto with × ].
rewrite Op; auto with ×.
Qed.

modsym x m is the unique y congruent to x such that -m/2≤x<m/2

Definition modsym x m := (x + m / 2) mod m - m / 2.

Lemma modsym_bounds : x m, 0 < m- m 2 × modsym x m < m.
Proof.
intros x m Pm; unfold modsym.
pose proof Z_mod_lt (x + m / 2) m.
pose proof Z_mod_lt m 2.
split; ring_simplify; rewrite Z_mult_div_mod; omega.
Qed.

Lemma modsym_mod_compat : x m, (modsym x m) mod m = x mod m.
Proof.
intros x m; unfold modsym.
rewrite Zminus_mod_idemp_l.
f_equal; ring.
Qed.

Lemma mod_modsym_compat : x m, modsym (x mod m) m = modsym x m.
Proof.
intros x m; unfold modsym.
rewrite Zplus_mod_idemp_l; auto.
Qed.

Lemma modsym_mod_diff : x m, 0 < m{ k | modsym x m = x mod m + m × k }.
Proof.
intros x m.
((modsym x m - x mod m) / m).
rewrite Z_mult_div_mod; [ | auto with × ].
rewrite Zminus_mod, modsym_mod_compat.
rewrite Zminus_mod_idemp_l.
do 2 rewrite Zminus_mod_idemp_r.
rewrite Zminus_diag, Zmod_0_l.
ring.
Qed.

Lemma modsym_eqm : x m, modsym x m x [ m ].
Proof.
intros x m.
apply modsym_mod_compat.
Qed.

Lemma prime_div_false : a p, prime p(a | p) → 1 < a < pFalse.
Proof.
intros a p Pp D Bp.
destruct (prime_divisors p Pp a D) as [|[|[|]]]; omega.
Qed.

Lemma Zbounding_square : x m, 0 < m-m x mx ^ 2 m ^ 2.
Proof.
intros x m Pm Bx.
simpl; unfold Zpower_pos; simpl.
rewrite Zmult_assoc.
rewrite <- Zabs_square.
apply Zle_trans with (Zabs x × m).
rewrite <- Zmult_assoc.
apply Zmult_le_compat_l; zify; omega.

rewrite Zmult_comm.
apply Zmult_le_compat_l; zify; omega.
Qed.

All prime numbers divides some 1+l²+m² (plus convenient conditions on l,m)

Lemma prime_dividing_sum_of_two_squares_plus_one : p,
prime p → 3 p
{l : _ & {m : _ & { k |
p × k = 1 + l × l + m × m
2 × m < p 2 × l < p 0 < k 0 l 0 m (0 < l 0 < m)}}}.
Proof.
intros p Pp Op.

pose (np := Zabs_nat p).

assert (p_odd : p mod 2 = 1) by (apply prime_odd; auto || omega).

pose (s := fun x : Zx × x).
assert (s_pos : x, 0 s x).
intros; unfold s; rewrite <- Zabs_square.
apply Zmult_le_0_compat; auto with ×.

assert (mod_pos : x, 0 x mod p).
intros; apply Z_mod_lt; destruct Pp; omega.

assert (hp_pos : 0 p / 2) by (apply Z_div_pos; omega).

assert (autobound : i, (i < S (Zabs_nat (p / 2)))%nat → 2 × Z_of_nat i < p).
intros i Li.
apply inj_lt in Li.
rewrite inj_S, inj_Zabs_nat, Zabs_eq in Li; auto.
apply odd_bound_1; auto.

pose (FL := fun lZabs_nat (s (Z_of_nat l) mod p)).
pose (FM := fun mZabs_nat ((-1 - s (Z_of_nat m)) mod p)).
assert (IFL : injective (S (Zabs_nat (p / 2))) FL).
intros i j Li Lj E.
unfold FL in E.
apply Zabs_nat_inj in E; eauto.
apply inj_eq_rev.
apply (Zmod_sqrt_eq_compat p (Z_of_nat i) (Z_of_nat j)); auto with ×.

assert (IFM : injective (S (Zabs_nat (p / 2))) FM).
intros i j Li Lj E.
unfold FM in E.
apply inj_eq_rev.
apply (Zmod_sqrt_eq_compat p (Z_of_nat i) (Z_of_nat j)); auto with ×.

apply Zabs_nat_inj in E; eauto.
symmetry; apply eqm_minus_0.
apply eqm_minus_0 in E.
rewrite <- E.
apply eq_eqm; unfold s; ring.

assert (BFL : bounded (S (Zabs_nat (p / 2))) (Zabs_nat p) FL).
intros i _.
unfold FL.
apply inj_lt_rev; do 2 rewrite inj_Zabs_nat.
repeat rewrite Zabs_eq; auto with ×.
apply Z_mod_lt; auto with ×.

assert (BFM : bounded (S (Zabs_nat (p / 2))) (Zabs_nat p) FM).
intros i _.
unfold FM.
apply inj_lt_rev; do 2 rewrite inj_Zabs_nat.
repeat rewrite Zabs_eq; auto with ×.
apply Z_mod_lt; auto with ×.

assert (CL := count_image_injective (S (Zabs_nat (p / 2))) np FL IFL BFL).
assert (CM := count_image_injective (S (Zabs_nat (p / 2))) np FM IFM BFM).

remember (image FL (S (Zabs_nat (p / 2)))) as L.
remember (image FM (S (Zabs_nat (p / 2)))) as M.

destruct (count_drawers L M np) as (i, (Bi, (Li, Mi))).
rewrite CL, CM.
apply inj_lt_rev.
rewrite inj_plus, inj_S, inj_Zabs_nat.
rewrite Zabs_eq; auto.
unfold np; rewrite inj_Zabs_nat, Zabs_eq; auto with ×.
unfold Zsucc.
ring_simplify.
pose proof Z_mult_div_bounds p 2.
omega.

rewrite HeqL in Li; rewrite HeqM in Mi.
destruct (image_true _ _ _ Li) as (l, (Bl, Hl)).
destruct (image_true _ _ _ Mi) as (m, (Bm, Hm)).
(Z_of_nat l); (Z_of_nat m).
unfold FL in Hl; unfold FM in Hm.
clear HeqL CL HeqM CM Li Mi IFL IFM BFL BFM FL FM.
subst.
apply inj_eq in Hm; do 2 rewrite inj_Zabs_nat in Hm.
do 2 rewrite Zabs_eq in Hm; auto.
symmetry in Hm; apply eqm_minus_0 in Hm.
rewrite divide_eqm in Hm; notzero.
pose proof s_pos (Z_of_nat l).
pose proof s_pos (Z_of_nat m).
destruct (Zdivide_inf _ _ Hm) as (k, Hk).
k.
assert (0 < k × p) by omega.
assert (0 < k) by (apply Zmult_lt_0_reg_r with p; assumption || omega).
assert (3 k × p).
replace 3 with (1 × 3) by reflexivity.
apply Zmult_le_compat; omega.
repeat split; auto with ×.
rewrite Zmult_comm, <- Hk.
unfold s; ring.

rem (Z_of_nat l) lz Elz.
rem (Z_of_nat m) mz Emz.
cut (0 lz 0 mz); [ omega | ].
cut (0 s lz 0 s mz).
cut ( a, 0 s a → 0 a).         intros Hyp [?|?]; [left|right]; apply Hyp; auto.
clear; intros [] H; tauto || zify; auto with ×.
omega.
Defined.

Building bounds to prove things like x1²+x2²+x3²+x4²=m² /\ -m/2≤xi<m/2 => xi=-m/2

Lemma egality_case_sum_of_four : a b c d M,
a Mb Mc Md M
a + b + c + d = 4 × M((a = M) (b = M)) ((c = M) (d = M)).
Proof.
intros.
omega.
Qed.

Lemma square_bound_equality_case : a M,
-M a < MM × M a × aa = - M.
Proof.
intros a M Bounds Bsqr.
destruct (Z_le_dec 0 M); [ | omega ].
cut (Zabs M Zabs a).
intro; zify; omega.

rewrite <- (Zabs_square M), <- (Zabs_square a) in Bsqr.
apply sqrt_le_compat; auto with ×.
Qed.

Lemma square_bound : x m, -m x mx × x m × m.
Proof.
intros x m Bx.
assert (Pm : 0 m) by omega.
rewrite <- Zabs_square.
apply Zle_trans with (Zabs x × m).
apply Zmult_le_compat_l; auto with ×.
zify; omega.

apply Zmult_le_compat_r; auto with ×.
zify; omega.
Qed.

Lemma square_bound_opp : x m, m x -mx × x m × m.
Proof.
intros x m Bx.
rewrite <- (Zmult_opp_opp m m).
apply square_bound; auto with ×.
Qed.

Lemma egality_case_sum_of_four_squares : a b c d M,
-M a < M-M b < M-M c < M-M d < M
a × a + b × b + c × c + d × d = 4 × (M × M)
((a = -M) (b = -M)) ((c = -M) (d = -M)).
Proof.
intros a b c d m Ba Bb Bc Bd E.
assert (P : x, -m x < m-m x m) by (intros; omega).
pose proof square_bound _ _ (P _ Ba).
pose proof square_bound _ _ (P _ Bb).
pose proof square_bound _ _ (P _ Bc).
pose proof square_bound _ _ (P _ Bd).
assert (Pm : 0 < m) by omega.
cut ((Zabs a = m Zabs b = m) (Zabs c = m Zabs d = m)).
intros; zify; omega.

split; split; (rewrite <- (Zabs_eq m); [ | auto with × ]);
apply sqrt_eq_compat_abs;
omega.
Qed.

If mp is the sum of four squares, then so is np for a smaller n

Lemma foursquare_prime_factor_decreasing :
p, prime p m, (1 < m m < p)%Z
foursquare (m × p) →
sigT (fun n ⇒ ((0 < n n < m)%Z × foursquare (n × p))%type).
Proof.
intros p Pp m (LBm, UBm) FSmp.
assert (help0 : m 0); auto with ×.
assert (help1 : 0 < m); auto with ×.
destruct FSmp as (x1, (x2, (x3, (x4, Hx)))).
pose (y1 := modsym x1 m).
pose (y2 := modsym x2 m).
pose (y3 := modsym x3 m).
pose (y4 := modsym x4 m).
assert (eqm1 : y1 x1 [m]) by apply modsym_eqm.
assert (eqm2 : y2 x2 [m]) by apply modsym_eqm.
assert (eqm3 : y3 x3 [m]) by apply modsym_eqm.
assert (eqm4 : y4 x4 [m]) by apply modsym_eqm.

assert (Dm : (m | y1 × y1 + y2 × y2 + y3 × y3 + y4 × y4)).
apply divide_eqm; auto.
rewrite eqm1, eqm2, eqm3, eqm4.
eapply multiple_eqm; eauto.

assert (Dmp : ¬ (m | p)).
intros D.
eapply prime_div_false; eauto.

destruct (Zdivide_inf _ _ Dm) as (r, Hr).
r.
split.
assert (Pr : 0 r).
apply Zmult_le_0_reg_r with m; auto with ×.
rewrite <- Hr.
pose proof Zle_0_square y1.
pose proof Zle_0_square y2.
pose proof Zle_0_square y3.
pose proof Zle_0_square y4.
omega.

assert (Nr : 0 r).
intros Zr; subst.

ring_simplify (0 × m) in Hr.
assert (y0 : (y1 = 0 y2 = 0) (y3 = 0 y4 = 0)).
pose proof Zle_0_square y1.
pose proof Zeq_0_square y1.
pose proof Zle_0_square y2.
pose proof Zeq_0_square y2.
pose proof Zle_0_square y3.
pose proof Zeq_0_square y3.
pose proof Zle_0_square y4.
pose proof Zeq_0_square y4.
omega.

apply Dmp.
apply Zmult_divide_compat_rev_l with m; auto.
rewrite Hx.
repeat apply Zdivide_plus_r;
apply Zdivide_square, divide_eqm; auto.
rewrite <- eqm1; apply eq_eqm; tauto.
rewrite <- eqm2; apply eq_eqm; tauto.
rewrite <- eqm3; apply eq_eqm; tauto.
rewrite <- eqm4; apply eq_eqm; tauto.

assert (Lrm : r m).
apply Zmult_le_reg_r with m; auto with ×.
rewrite <- Hr.
apply Zmult_le_reg_r with 4; auto with ×.
ring_simplify.
clear -help1 LBm.
assert (B : x, 4 × (modsym x m) ^ 2 m ^2).
clear -help1; intros x.
pose proof modsym_bounds x m help1 as Bms.
set (y := modsym x m); fold y in Bms; clearbody y.
replace (4 × y ^ 2) with ((2 × y) ^ 2) by ring.
apply Zbounding_square; auto.
now auto with ×.

pose proof B x1 as Ey1; fold y1 in Ey1.
pose proof B x2 as Ey2; fold y2 in Ey2.
pose proof B x3 as Ey3; fold y3 in Ey3.
pose proof B x4 as Ey4; fold y4 in Ey4.
omega.

assert (Nmr : r m).
intros Emr; subst.
assert (Eyi : (2 × - y1 = m 2 × - y2 = m) (2 × - y3 = m 2 × - y4 = m)).
assert (Eyi : (2 × y1 = - m 2 × y2 = - m) (2 × y3 = - m 2 × y4 = - m)).
apply egality_case_sum_of_four_squares;
try (apply modsym_bounds; notzero).
rewrite <- Hr; ring.
omega.

destruct Eyi as ((Ey1, Ey2), (Ey3, Ey4)).
apply Dmp.
apply Zmult_divide_compat_rev_l with m; auto.
rewrite Hx.
apply Zmult_divide_compat_rev_l with 4; notzero.
rewrite <- divide_eqm; notzero.
ring_simplify.
assert (REW : x, 4 × x ^ 2 = (2 × x) × (2 × x)) by (intros; ring).
repeat rewrite REW; clear REW.
transitivity (m × m + m × m + m × m + m × m).
cut ( x y, 2 × -y = my x [m]
(2 × x) × (2 × x) m × m [4 × (m × m)]).
intros HYP.
repeat apply Zplus_eqm;
(eapply HYP; [ | eassumption ]; eassumption).

intros x y Ey Mxy.
apply eqm_square_half; notzero.
rewrite <- Ey at 2.
apply eqm_mult_compat_l.
rewrite <- Mxy.
replace (- y) with (y + 2 × - y) by ring.
rewrite Ey.
rewrite eqm_minus_0.
apply divide_eqm; notzero.
(- 1); ring.

apply divide_eqm; notzero.
1; ring.

omega.

assert (Erpm : (r × p) × (m × m) = (r × m) × (m × p)) by ring.
rewrite Hx, <-Hr in Erpm.
rewrite euler's_identity in Erpm.
rem (y1 × x1 + y2 × x2 + y3 × x3 + y4 × x4) t1 Et1.
rem (y1 × x2 - y2 × x1 - y3 × x4 + y4 × x3) t2 Et2.
rem (y1 × x3 + y2 × x4 - y3 × x1 - y4 × x2) t3 Et3.
rem (y1 × x4 - y2 × x3 + y3 × x2 - y4 × x1) t4 Et4.

assert (D1 : (m | t1)).
rewrite <- divide_eqm; auto; rewrite Et1.
unfold y1, y2, y3, y4; repeat rewrite modsym_eqm.
eapply multiple_eqm; eauto.

assert (D2 : (m | t2)).
rewrite <- divide_eqm; auto; rewrite Et2.
unfold y1, y2, y3, y4; repeat rewrite modsym_eqm.
try ring .
red; f_equal; ring.

assert (D3 : (m | t3)).
rewrite <- divide_eqm; auto; rewrite Et3.
unfold y1, y2, y3, y4; repeat rewrite modsym_eqm.
red; f_equal; ring.

assert (D4 : (m | t4)).
rewrite <- divide_eqm; auto; rewrite Et4.
unfold y1, y2, y3, y4; repeat rewrite modsym_eqm.
red; f_equal; ring.

(t1 / m); (t2 / m); (t3 / m); (t4 / m).
apply Zmult_reg_r with (m × m).
clear -help0; destruct m; auto; intro H; inversion H.
rewrite Erpm.
ring_simplify.
cut ( a, (m | a)(a / m) ^ 2 × m ^ 2 = a ^ 2).
intros Q.
do 4 (rewrite Q; auto).

clear -help0.
intros a D.
destruct (Zdivide_inf m a D) as (k, E); subst.
rewrite Z_div_mult_full; auto.
ring.
Defined.

Induction scheme to use the previous lemma
Definition lt_wf_rect :=
fun p (P : natType) F
well_founded_induction_type
(well_founded_ltof nat (fun mm)) P F p.

Using the previous lemma, one can decrease mp into np to eventually get p

Lemma foursquare_prime : p, prime pfoursquare p.
Proof.
intros p Pp.

destruct (Z_eq_dec p 2) as [E | E].
do 2 1%Z; do 2 0%Z; auto.

assert (3 p) as Op by
(pose proof (prime_ge_2 p Pp); zify; omega); clear E.

pose (fs_mult := fun p m
prod (foursquare (Z_of_nat (S m) × p))
(Z_of_nat (S m) < p)).
assert (sigT (fs_mult p)) as Hm.
destruct (prime_dividing_sum_of_two_squares_plus_one p Pp Op) as
(l, (m, (k, (Ep, (Bm, (Bl, (Pk, Plm))))))).
assert (tech1 : Z_of_nat (S (Zabs_nat (k - 1))) = k).
rewrite inj_S, inj_Zabs_nat, Zabs_eq; [ | omega ]; auto with ×.
(Zabs_nat (k - 1)); split.
0%Z; 1%Z; l; m.
transitivity (p × k).
rewrite tech1; ring.
rewrite Ep; auto.

rewrite tech1; clear tech1.
apply Zmult_lt_reg_r with (4 × p); [ auto with × | ].
replace (p × (4 × p)) with (2 × (p × p) + 2 × (p × p)) by ring.
replace (k × (4 × p)) with (4 × (p × k)) by ring.
rewrite Ep.
assert (tech2 : a b, 0 < a → 2 × a < b → 4 × a × a + 4 < b × b).
clear; intros a b Pa LT.
assert (LE : 2 × a + 1 b) by omega.
assert (Pda : 0 2 × a + 1) by omega.
assert (LE2 := Zmult_le_compat _ _ _ _ LE LE Pda Pda).
eapply Zlt_le_trans; [ | apply LE2 ].
ring_simplify.
omega.

assert (tech3 : a b, 0 a → 2 × a < b → 4 × a × a < b × b).
clear; intros a b Pa LT.
assert (LE : 2 × a + 1 b) by omega.
assert (Pda : 0 2 × a + 1) by omega.
assert (LE2 := Zmult_le_compat _ _ _ _ LE LE Pda Pda).
eapply Zlt_le_trans; [ | apply LE2 ].
ring_simplify.
omega.

assert (p2_pos : 0 < p × p).
transitivity (1 × p).
omega.
apply Zmult_lt_compat_r; omega.

rem (l × l) ll Ell.
rem (m × m) mm Emm.
rem (p × p) pp Epp.
destruct Plm as (NNl, (NNm, [Pl | Pm])).
specialize (tech2 _ p Pl Bl).
specialize (tech3 _ p NNm Bm).
rewrite <- Zmult_assoc, <-Ell, <-Emm, <-Epp in ×.
clear -tech2 tech3 p2_pos.
ring_simplify.
omega.

specialize (tech2 _ p Pm Bm).
specialize (tech3 _ p NNl Bl).
rewrite <- Zmult_assoc, <-Ell, <-Emm, <-Epp in ×.
clear -tech2 tech3 p2_pos.
ring_simplify.
omega.

clear Op.
destruct Hm as (m, Hm).
generalize p Pp m Hm; clear Hm m Pp p.
intros p Pp m.
refine (lt_wf_rect m (fun nfs_mult p nfoursquare p) _); clear m.
intros [|m] IH (FSm, UBm).
destruct p; assumption.

assert (LBm : 1 < Z_of_nat (S (S m))) by (zify; omega).
destruct (foursquare_prime_factor_decreasing p Pp _ (conj LBm UBm) FSm) as (n, ((LBn, UBn), FSn)).
apply IH with (Zabs_nat (n - 1)).
unfold ltof.
zify; omega.

unfold fs_mult.
rewrite inj_S, inj_Zabs_nat, Zabs_eq; auto with *; unfold Zsucc.
replace (n - 1 + 1) with n by ring.
split; auto; omega.
Defined.

Trivial application of the previous lemma and euler's identity

Theorem lagrange_4_square_theorem : n, 0 nfoursquare n.
Proof.
intros n; eapply Z_prime_rect; auto.
repeat ( 0); auto.

1; repeat ( 0); auto.

apply foursquare_prime.

apply mult_foursquare_compat.
Defined.

Definition lagrange_fun (n : Z) : (Z × Z) × (Z × Z) :=
let (a, ha) := lagrange_4_square_theorem (Zabs n) (Zabs_pos n) in
let (b, hb) := ha in
let (c, hc) := hb in
let (d, _ ) := hc in
((a, b), (c, d)).

Lemma lagrange_fun_spec (n : Z) :
(let (ab, cd) := lagrange_fun n in let (a, b) := ab in let (c, d) := cd in
Zabs n = a × a + b × b + c × c + d × d).
Proof.
unfold lagrange_fun.
destruct (lagrange_4_square_theorem (Zabs n) (Zabs_pos n)) as (a, (b, (c, (d, pr)))).
exact pr.
Qed.

Extraction "Arith/Lagrange_four_square.ml" lagrange_fun.