Remove LoadPath "/usr/local/lib/coq/theories/Reals".

Require Import ssreflect ssrfun eqtype ssrbool ssrnat seq choice finset div
  fintype tuple finfun bigop prime binomial ssralg morphism fingroup perm
  automorphism quotient action finalg zmodp matrix mxalgebra vector zmodp.

Set Implicit Arguments.
Unset Strict Implicit.
Unset Printing Implicit Defensive.

Import GroupScope.
Import GRing.Theory.
Open Local Scope ring_scope.

Open Scope fun_scope.

Section comp_det_def.

Variable R : comRingType.

Fixpoint expand sign (mat : seq (seq R)) det acc : R :=
  if mat is a :: tl then
    sign * head 0 a * det (rev acc ++ (map behead tl)) +
    expand (- sign) tl det (behead a :: acc)
  else
    0.

Fixpoint sdet n (mat : seq (seq R)) : R :=
  if n is p.+1 then expand 1 mat (sdet p) [::] else 1.

Definition sseq2mat n m (mat: seq (seq R)) :=
  \matrix_(i < n, j < m) nth 0 (nth [::] mat i) j.

Definition mat2sseq n m (mat : 'M[R]_(n,m)) : seq (seq R) :=
map (fun i => map (fun j => mat i j) (enum 'I_m)) (enum 'I_n).

Lemma can_mat2sseq : forall n m mat, sseq2mat n m (mat2sseq mat) = mat.
Proof.
move =>  n m mat; apply/matrixP => x y.
rewrite /sseq2mat /mat2sseq mxE.
rewrite (nth_map x); last by rewrite size_enum_ord ltn_ord.
rewrite (nth_map y); last by rewrite size_enum_ord ltn_ord.
by rewrite !nth_ord_enum.
Qed.

Lemma castmx_det :
  forall n p (m : 'M[R]_n) (h : n = p), \det m = \det (castmx (h, h) m).
Proof.
 by move => n p m h; move: (h) m; rewrite h => h' m; rewrite castmx_id.
Qed.

Lemma col'0_sseq2mat_behead :
  forall l p mat, col' ord0 (sseq2mat l p.+1 mat) =
       sseq2mat l p (map behead mat).
Proof.
move => l p mat; rewrite /col' /sseq2mat; apply/matrixP => i j; rewrite !mxE /=.
case w: (i < size mat).
 rewrite (nth_map ([::]:seq R)) // /behead.
 by  case: (nth _ _ i); case j => [ [| ?] ?].
by rewrite !nth_default // ?size_map leqNgt negbT.
Qed.

Lemma col'0_sseq2mat :
  forall l1 p n mat (h : p.+1 = n), 
    col' (cast_ord h ord0) (sseq2mat l1 n mat) =
    sseq2mat l1 n.-1 (map behead mat).
Proof.
by move => l1 p n mat h; case h; rewrite cast_ord_id col'0_sseq2mat_behead.
Qed.

Lemma cast_sseq2mat :
  forall n1 n2 n3 n4 (q : n1 = n2) (h : n3 = n4) m,
    castmx (q, h) (sseq2mat n1 n3 m) = sseq2mat n2 n4 m.
Proof. 
by move => n1 n2 n3 n4 q h; move: (q) (h); rewrite q h => q' h' m;
   rewrite castmx_id.
Qed.

Lemma sseq2mat_cat :
  forall k l p m1 m2, size m1 = k -> 
    sseq2mat (k + l) p (m1 ++ m2) = col_mx (sseq2mat k p m1) (sseq2mat l p m2).
Proof.
move => k l p m1 m2 s1; apply/matrixP => i j; rewrite /sseq2mat /col_mx !mxE.
rewrite nth_cat s1; case: splitP => [i' -> | j' ->]; rewrite mxE //.
by rewrite -{2}[k]addn0 subn_add2l subn0. 
Qed.

Lemma zip_rcons : forall T1 T2 (l1:seq T1)(l2:seq T2) a b,
  size l1 = size l2 -> zip (rcons l1 a) (rcons l2 b) = rcons (zip l1 l2) (a,b).
Proof.
move => T1 T2; elim => [ | c l1 Hl1]; first by case.
by move => [ | d l2] a b //; move => [q] /=; rewrite Hl1.
Qed.

Lemma rev_zip : forall T1 T2 (l1:seq T1)(l2:seq T2),
   size l1 = size l2 -> rev (zip l1 l2) = zip (rev l1) (rev l2).
Proof.
move => T1 T2; elim => [|a l1 Hl1] [|b l2] // [sl].
by rewrite !rev_cons zip_rcons ?size_rev // Hl1.
Qed.

Lemma row'0_sseq2mat : forall m1 r m2 k l p, size m1 = k ->
   row' (rshift k 0) (sseq2mat (k + l.+1) p (m1 ++ (cons r m2))) =
   sseq2mat (k + l.+1).-1 p (m1 ++ m2).
Proof.
move => m1 r m2 k l p sm1; apply/matrixP => i j; rewrite !mxE /= /bump addn0.
case q: (k <= i).
 by rewrite !nth_cat add1n sm1 !ltnNge leqW //= leq_subS q.
by rewrite !nth_cat add0n ltnNge sm1 q.
Qed.

Definition czip (v:seq R) (l: seq (seq R)) :=
  map (fun p => let: (x,y) := p in x::y) (zip v l).

Lemma expand_correct :
  forall mat p (k l:nat) (f : seq (seq R) -> R)
   s acc v (h :(p.+1 = k + l)%N),
   size acc = k -> size mat <= l -> size v = k ->
    (forall mat', size mat' <= p -> f mat' = \det (sseq2mat p p mat')) ->
    let r := sseq2mat (k+l) p.+1 (rev (czip v acc) ++ mat) in
    expand s mat f acc =
      s * (-1)^+ k * \sum_(i < l) r (rshift k i) 0 *
      cofactor (castmx (refl_equal (k + l)%N, h) r)
         (rshift k i) (cast_ord h ord0).
Proof.
elim => [ | r mat' Ih] p k l f s acc v h sa sm sv fc.
 set ref := sseq2mat _ _ _; rewrite /= big1 ?mulr0 // => i _.
 suffices -> :  ref (rshift k i) 0 = 0 by rewrite mul0r.
 rewrite /ref sseq2mat_cat; first by rewrite col_mxEd mxE !nth_default.
 by rewrite size_rev /czip size_map size_zip sa sv minnn. 
move: h sm; case ld : l => [ | l'] // h /= sm; set ref := sseq2mat _ _ _.
have qp : p.+1 = (k.+1 + l')%N by rewrite h addnS.
have qp' : p = (k + l'.+1).-1 by rewrite addnS -addSn -qp.
have sa' : size ((behead r)::acc) = k.+1 by rewrite /= sa.
have sv' : size ((head 0 r)::v) = k.+1 by rewrite /= sv.
have sz : size (rev (czip v acc)) = k.
 by rewrite size_rev /czip size_map size_zip sv sa minnn.
move: (Ih p _ l' f (-s) _ _ qp sa' sm  sv' fc) => -> {Ih}.
rewrite big_ord_recl {1}/ref.
have -> : sseq2mat (k + l'.+1) p.+1 (rev (czip v acc) ++ (r::mat'))
              (rshift k ord0) 0 = head 0 r.
 by rewrite /sseq2mat mxE /= nth_cat addn0 sz ltnn subnn /=.
have -> : cofactor (castmx (Logic.eq_refl (k + l'.+1)%N, h) ref)
      (rshift k ord0) (cast_ord h ord0) =
        (-1)^+k * f (rev acc ++ map behead mat').
 rewrite /cofactor /ref /= !addn0; congr (_ * _); rewrite fc; last first.
  by rewrite size_cat size_rev size_map sa qp' addnS /= leq_add2l -ltnS.
 rewrite cast_sseq2mat col'0_sseq2mat map_cat map_rev /czip -map_comp.
 have unzs : {in (zip v acc), funcomp tt behead
                       (fun p0 : R * seq R => let '(x, y) := p0 in cons x y) =1
                 @snd  R _} by move => [a b] _ /=.
 rewrite (eq_in_map unzs) -/unzip2 unzip2_zip; last by rewrite sa sv leqnn.
 rewrite /= row'0_sseq2mat; last by rewrite size_rev.
 by rewrite -(cast_sseq2mat (esym qp') (esym qp')) -castmx_det.
rewrite mulr_addr; congr (_ + _).
 rewrite !mulrA; congr (_ * _); rewrite -!mulrA; congr (_ * _).
 by rewrite  mulrC -mulrA -exprn_addr -signr_odd addnn odd_double expr0 mulr1.
rewrite -mulN1r exprS -(mulrA (-1)) !(mulrC (-1)) !mulrA
 -(mulrA _ (-1) (-1)) mulrNN !mulr1 (mulrC s); congr (_ * _).
apply eq_big => // i _;rewrite !mxE /ref.
have -> : rshift k (lift ord0 i) = (k + i.+1)%N :> nat by [].
have -> : rshift k.+1 i = (k.+1 + i)%N :> nat by [].
congr (_ * _).
 by rewrite !nth_cat !ltnNge !size_rev /czip !size_map !size_zip sa' sv'
   sa sv !minnn leq_addr addnS addKn (leqW (leq_addr i k)) -(addnS k i) addKn.
congr (_ * _); first by rewrite /= !addnS.
rewrite !cast_sseq2mat !col'0_sseq2mat.
have h' : (k + l'.+1).-1 = (k.+1 + l').-1 by rewrite addnS.
rewrite (castmx_det _ h'); apply: f_equal.
apply/matrixP => x y; rewrite /row' castmxE !mxE /czip -!map_rev.
have sc : size (cons _ v) = size (cons _ acc) by move => *; rewrite /= sv.
have sav : size v = size acc by rewrite sv.
rewrite rev_zip // !rev_cons zip_rcons ?map_rcons ?size_rev //.
by rewrite cat_rcons !map_cat !map_cons rev_zip //= /bump leq0n addnS.
Qed.

Lemma sdet_correct :
  forall n mat, size mat <= n -> sdet n mat = \det (sseq2mat n n mat).
Proof.
elim => [ | n Ih] mat sm /=.
 by rewrite /sdet /sseq2mat (flatmx0 (\matrix_(i,j) _)) -(flatmx0 1%:M) det1.
rewrite (expand_correct _ (erefl _ : n.+1 = (0 + n.+1)%N)
           (erefl (size nil)) sm (erefl (size nil)) Ih).
rewrite !castmx_id (expand_det_col _ 0) expr0 !mul1r; apply: eq_bigr => i _.
by rewrite cast_ord_id (_ : rshift 0 i = i) //; apply: val_inj.
Qed.

Lemma sdet_compute :
  forall n (m:'M[R]_n), \det m = sdet n (mat2sseq m).
Proof.
move => n m; rewrite sdet_correct ?(can_mat2sseq m) //.
 by rewrite /mat2sseq size_map size_enum_ord.
Qed.

End comp_det_def.

Require Import ZArith.

(* The next part is about adding a commutative ring structure to Z. *)

Lemma Zeq_axiom : Equality.axiom Zeq_bool.
rewrite /Equality.axiom => x y; case h: (Zeq_bool x y).
  by apply: ReflectT; apply/Zeq_is_eq_bool. 
by apply ReflectF; apply/Zeq_is_eq_bool; rewrite h.
Qed.

Definition Zpickle (x : Z) : nat :=
  match x with
    Z0 => 0%nat
  | Zpos p => ((nat_of_P p).*2 - 1)%nat
  | Zneg p => ((nat_of_P p).*2)%nat
  end.

Definition Zunpickle (x : nat) : Z :=
  if x == 0%nat then Z0 else
  if odd x then Zpos (P_of_succ_nat (half x)) else
   Zneg (P_of_succ_nat (half (x - 1))).

Lemma SZ_pickleK : cancel Zpickle Zunpickle.
move => [ | p | p]; first done; rewrite /Zpickle /Zunpickle;
 move/ltP: (lt_O_nat_of_P p); case pk : (nat_of_P p) => [ | k] _ //.
 have -> : (k.+1.*2 - 1 == 0)%nat = false by rewrite doubleS subSS subn0.
 have -> : odd (k.+1.*2 - 1) by rewrite odd_sub // odd_double.
 have -> : (k.+1.*2 -1 = 1 + k.*2)%nat by rewrite doubleS subSS subn0 addSn add0n.
 rewrite (half_bit_double _ true); apply f_equal.
 by apply: nat_of_P_inj; rewrite nat_of_P_o_P_of_succ_nat_eq_succ.
rewrite (_ : (k.+1.*2 == 0)%nat = false) // (_ : (k.+1.*2 - 1 = 1 + k.*2)%nat) //.
rewrite odd_double (half_bit_double _ true); apply f_equal.
by apply: nat_of_P_inj; rewrite nat_of_P_o_P_of_succ_nat_eq_succ.
Qed.

Lemma n10Z : ~~Zeq_bool 1 0.
Proof. by []. Qed.

Definition SZ_eqMixin :=  Equality.Mixin Zeq_axiom.
Canonical Structure SZ_eqType : eqType := Equality.pack SZ_eqMixin.
Definition SZ_countable_mixin := CanCountMixin SZ_pickleK.
Definition SZ_choiceMixin := CountChoiceMixin SZ_countable_mixin.
Canonical Structure SZ_choiceType := Eval hnf in ChoiceType Z SZ_choiceMixin.
Canonical Structure SZ_countType := Eval hnf in CountType Z SZ_countable_mixin.
Definition SZ_ZmodMixin :=
  ZmodMixin Zplus_assoc Zplus_comm Zplus_0_l Zplus_opp_l.
Canonical Structure SZ_ZmodType := ZmodType Z SZ_ZmodMixin.
Definition SZ_comRingMixin : GRing.Ring.mixin_of SZ_ZmodType := 
  ComRingMixin Zmult_assoc Zmult_comm Zmult_1_l Zmult_plus_distr_l n10Z.
Canonical Structure SZ_Ring := Eval hnf in RingType Z SZ_comRingMixin.
Canonical Structure SZ_comRing := Eval hnf in ComRingType Z Zmult_comm.

Eval vm_compute in sdet 2 ([:: [:: 1; 2]; [::3; 4]])%Z.

Eval vm_compute in sdet 4
 ([::[::  1;  2;  3;  4];
     [::  5;  6;  7;  8];
     [::  9; 10; 28;  0];
     [::  3; 2;   2; -1]])%Z.

(* The following computation takes 8 seconds. *)
Time Eval vm_compute in sdet 9
 ([::[::  1;  2;  3;  4; -2;  3; 7;  8; 1];
     [::  5;  6;  7;  8; 14; 21; 5; 12; 2];
     [::  9; 10; 28;  0;  9;  7; 6; 11; 3];
     [::  3; 2;   2; -1;  1;  2; 3;  4; 4];
     [::  1;  2;  3;  4; -2;  3; 9;  8; 5];
     [::  5;  6;  7;  8;  3;  13;5; 12; 6];
     [::  9; 10; 29;  0;  5; 25; 6; 11; 7];
     [::  3; 2;   2; -1;  0;  0; 0;  0; 8];
     [::  3; 2;   2; -1;  0;  0; 3;  0; 9]])%Z.