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.
Require Import ssralg ssrint ssrnum fintype.
Require Import dvdring matrix mxalgebra bigop zmodp perm mxstructure.
Formalization of Anders' Haskell algorithm for general
Smith Normal form using ssreflect matrices.
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'))).
Try to rewrite without case analysis
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.
a * d - b * c = 1 ->
combine_mx a b c d m k *m combine_mx d (-b) (-c) a m k= 1%:M.
If the following Bezout identity holds: u * a1 + v * b1 = 1,
bz_mx a b n k represents the following matrix (dots are zeros):
(kth column)
/ u .... v ..... \
| . 1 .......... |
| ....1......... |
| -b1 .. a1 .... | (kth row)
| ..........1... |
\ .............1 /
(determinant is +/-1)
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.
Two-steps approach:
1) Make the first column look like
___
/ g \
| g |
| . |
| . |
| g |
\___/
2)
For any i j s.t. ~~ g %| M i j, xrow 0 i M, bezout step on the first row
and back to 1)
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.
This page has been generated by coqdoc