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.