Library smith

Require Import ssreflect ssrfun ssrbool eqtype ssrnat div seq path.
Require Import ssralg ssrint ssrnum fintype.
Require Import dvdring matrix mxalgebra bigop zmodp perm mxstructure.

Section smith_def.

Local Open Scope ring_scope.

Variable R : euclidRingType.

Local Notation "''M_' ( m , n )" := 'M[R]_(m, n) : type_scope.
Local Notation "''M_' ( n )" := 'M[R]_(n, n) : type_scope.
Local Notation "''M_' n" := 'M[R]_(n, n) : type_scope.

Definition combine_mx (a b c d : R) (m : nat) (k : 'I_m) :=
  let k' := rshift 1 k in
  let d := \row_j (a *+ (j == 0) + d *+ (j == k') +
            ((j != 0) && (j != k'))%:R) in
  diag_mx d + c *: delta_mx k' 0 + b *: delta_mx 0 k'.

Import GRing.Theory.

Lemma combine_mxE (a b c d : R) (m n : nat) (M : 'M_(1 + m,1 + n)) (k : 'I_m) :
  let k' := rshift 1 k in
  let r0 := a *: row 0 M + b *: row k' M in
  let rk := c *: row 0 M + d *: row k' M in
  combine_mx a b c d m k *m M =
    \matrix_i (r0 *+ (i == 0) + rk *+ (i == k') +
               row i M *+ ((i != 0) && (i != k'))).

Lemma combine_mx_inv (a b c d : R) m k :
  a * d - b * c = 1 ->
  combine_mx a b c d m k *m combine_mx d (-b) (-c) a m k= 1%:M.

Definition bz_mx (a b : R) (m : nat) (k : 'I_m) : 'M_(1+m) :=
  let:(_,u,v,a1,b1) := egcdr a b in combine_mx u v (-b1) a1 m k.

Lemma bz_mx_00 {m n : nat} {M : 'M_(1 + m,1 + n)} {k : 'I_m} :
 (bz_mx (M 0 0) (M (rshift 1 k) 0) m k *m M) 0 0 %=
  gcdr (M 0 0) (M (rshift 1 k) 0).

Lemma sdvd_bz_mx (m n : nat) (M : 'M_(1 + m,1 + n)) (k : 'I_m) :
 ~~ (M 0 0 %| M (rshift 1 k) 0) ->
 (bz_mx (M 0 0) (M (rshift 1 k) 0) m k *m M) 0 0 %<| M 0 0.

Lemma unit_bz_mx m n (M : 'M_(1 + m,1 + n)) k :
  (bz_mx (M 0 0) (M (rshift 1 k) 0) m k) \in unitmx.

Local Open Scope ring_scope.

Fixpoint improve_pivot k {m n} : 'M_(1 + m, 1 + n) ->
  'M_(1 + m) * 'M_(1 + m, 1 + n) * 'M_(1 + n) :=
  match k with
  | 0 => fun A => (1%:M,A,1%:M)
  | p.+1 => fun A =>
    let a := A 0 0 in
      if [pick i | ~~(a %| A (rshift 1 i) 0)] is Some i then
        let L := bz_mx (A 0 0) (A (rshift 1 i) 0) m i in
        let A := L *m A in
        let:(L',A,R') := improve_pivot p A in
        (L' *m L, A, R')
      else let u := dlsubmx A in let v := ursubmx A in
      let u' := map_mx (fun x => 1 - odflt 0 (x %/? a)) u in
      let L := block_mx 1 0 u' 1%:M in
      let A := block_mx a%:M v (const_mx a) (drsubmx A + u' *m v) in
      if [pick ij | ~~(a %| A ij.1 (rshift 1 ij.2))] is Some ij then
        let A := xrow 0 ij.1 A in
        let R := (bz_mx (A 0 0) (A 0 (rshift 1 ij.2)) n ij.2)^T in
        let A := A *m R in
        let:(L',A,R') := improve_pivot p A in
        (L' *m xrow 0 ij.1 L, A, R *m R')
      else (L, A, 1%:M)
  end.

Section improve_pivot_correct.

Variables m n : nat.

CoInductive improve_pivot_spec M :
  'M_(1 + m) * 'M_(1 + m,1 + n) * 'M_(1 + n) -> Type :=
  ImprovePivotSpec L A R of L *m M *m R = (A : 'M_(1 + m, 1 + n))
  & (forall i j, A 0 0 %| A i j)
  & (forall i, A i 0 = A 0 0)
  & A 0 0 %| M 0 0
  & L \in unitmx & R \in unitmx
  : improve_pivot_spec M (L,A,R).

Lemma improve_pivotP : forall k (M : 'M_(1 + m,1 + n)),
  (enorm (M 0%R 0%R) <= k)%N -> M 0 0 != 0 ->
    improve_pivot_spec M (improve_pivot k M).

End improve_pivot_correct.

Fixpoint Smith {m n} : 'M_(m,n) -> 'M_(m) * seq R * 'M_(n) :=
  match m, n return 'M_(m, n) -> 'M_(m) * seq R * 'M_(n) with
  | _.+1, _.+1 => fun A : 'M_(1 + _, 1 + _) =>
      if [pick ij | A ij.1 ij.2 != 0] is Some (i, j) then
      let a := A i j in let A := xrow i 0 (xcol j 0 A) in
      (* this is where Euclidean norm eases termination argument *)
      let:(L,A,R) := improve_pivot (enorm a) A in
      let a := A 0 0 in
      let u := dlsubmx A in let v := ursubmx A in
      let v' := map_mx (fun x => odflt 0 (x %/? a)) v in
      let A := (drsubmx A) - (const_mx 1 *m v) in
      let:(L', d, R') := Smith (map_mx (fun x => odflt 0 (x %/? a)) A) in
      (lift0_mx L' *m block_mx 1 0 (-const_mx 1) 1%:M *m (xcol i 0 L),
       a :: [seq x * a | x <- d],
       (xrow j 0 R) *m block_mx 1 (-v') 0 1%:M *m lift0_mx R')
    else (1%:M, [::], 1%:M)
  | _, _ => fun A => (1%:M, [::], 1%:M)
  end.

CoInductive Smith_spec {m n} M : 'M[R]_m * seq R * 'M[R]_n -> Type :=
  SmithSpec L0 d R0 of L0 *m M *m R0 = diag_mx_seq m n d
  & sorted (@dvdr R) d
  & L0 \in unitmx & R0 \in unitmx : Smith_spec M (L0, d, R0).

Lemma SmithP : forall (m n : nat) (M : 'M_(m,n)), Smith_spec M (Smith M).

End smith_def.