proposition: Lemma int_Smith_normal_form m n (M : 'M[int]_(m, n)) : {L : 'M[int]_m & L \in unitmx & {R : 'M[int]_n & R \in unitmx & {d : seq int | sorted dvdz d & M = L *m (\matrix_(i, j) (d`_i *+ (i == j :> nat))) *m R}}}.
proof: Proof. move: {2}_.+1 (ltnSn (m + n)) => mn. elim: mn => // mn IHmn in m n M *; rewrite ltnS => le_mn. have [[i j] nzMij | no_ij] := pickP (fun k => M k.1 k.2 != 0); last first. do 2![exists 1%:M; first exact: unitmx1]; exists nil => //=. apply/matrixP=> i j; apply/eqP; rewrite mulmx1 mul1mx mxE nth_nil mul0rn. exact: negbFE (no_ij (i, j)). do [case: m i => [[]//|m] i; case: n j => [[]//|n] j /=] in M nzMij le_mn *. wlog Dj: j M nzMij / j = 0; last rewrite {j}Dj in nzMij. case/(_ 0 (xcol j 0 M)); rewrite ?mxE ?tpermR // => L uL [R uR [d dvD dM]]. exists L => //; exists (xcol j 0 R); last exists d => //=. by rewrite xcolE unitmx_mul uR unitmx_perm. by rewrite xcolE !mulmxA -dM xcolE -mulmxA -perm_mxM tperm2 perm_mx1 mulmx1. move Da: (M i 0) nzMij => a nz_a. have [A leA] := ubnP `|a|; elim: A => // A IHa in a leA m n M i Da nz_a le_mn *. wlog [j a'Mij]: m n M i Da le_mn / {j | ~~ (a %| M i j)%Z}; last first. have nz_j: j != 0 by apply: contraNneq a'Mij => ->; rewrite Da. case: n => [[[]//]|n] in j le_mn nz_j M a'Mij Da *. wlog{nz_j} Dj: j M a'Mij Da / j = 1; last rewrite {j}Dj in a'Mij. case/(_ 1 (xcol j 1 M)); rewrite ?mxE ?tpermR ?tpermD //. move=> L uL [R uR [d dvD dM]]; exists L => //. exists (xcol j 1 R); first by rewrite xcolE unitmx_mul uR unitmx_perm. exists d; rewrite //= xcolE !mulmxA -dM xcolE -mulmxA -perm_mxM tperm2. by rewrite perm_mx1 mulmx1. have [u [v]] := Bezoutz a (M i 1); set b := gcdz _ _ => Db. have{leA} ltA: (`|b| < A)%N. rewrite -ltnS (leq_trans _ leA) // ltnS ltn_neqAle andbC. rewrite dvdn_leq ?absz_gt0 ? dvdn_gcdl //=. by rewrite (contraNneq _ a'Mij) ?dvdzE // => <-; apply: dvdn_gcdr. pose t2 := [fun j : 'I_2 => [tuple _; _]`_j : int]; pose a1 := M i 1. pose Uul := \matrix_(k, j) t2 (t2 u (- (a1 %/ b)%Z) j) (t2 v (a %/ b)%Z j) k. pose U : 'M_(2 + n) := block_mx Uul 0 0 1%:M; pose M1 := M *m U. have{nz_a} nz_b: b != 0 by rewrite gcdz_eq0 (negPf nz_a). have uU: U \in unitmx. rewrite unitmxE det_ublock det1 (expand_det_col _ 0) big_ord_recl big_ord1. do 2!rewrite /cofactor [row' _ _]mx11_scalar !mxE det_scalar1 /=. rewrite mulr1 mul1r mulN1r opprK -[_ + _](mulzK _ nz_b) mulrDl. by rewrite -!mulrA !divzK ?dvdz_gcdl ?dvdz_gcdr // Db divzz nz_b unitr1. have{} Db: M1 i 0 = b. rewrite /M1 -(lshift0 n 1) [U]block_mxEh mul_mx_row row_mxEl. rewrite -[M](@hsubmxK _ _ 2) (@mul_row_col _ _ 2) mulmx0 addr0 !mxE /=. rewrite big_ord_recl big_ord1 !mxE /= [lshift _ _]((_ =P 0) _) // Da. by rewrite [lshift _ _]((_ =P 1) _) // mulrC -(mulrC v). have [L uL [R uR [d dvD dM1]]] := IHa b ltA _ _ M1 i Db nz_b le_mn. exists L => //; exists (R *m invmx U); last exists d => //. by rewrite unitmx_mul uR unitmx_inv. by rewrite mulmxA -dM1 mulmxK. move=> {A leA}IHa; wlog Di: i M Da / i = 0; last rewrite {i}Di in Da. case/(_ 0 (xrow i 0 M)); rewrite ?mxE ?tpermR // => L uL [R uR [d dvD dM]]. exists (xrow i 0 L); first by rewrite xrowE unitmx_mul unitmx_perm. exists R => //; exists d; rewrite //= xrowE -!mulmxA (mulmxA L) -dM xrowE. by rewrite mulmxA -perm_mxM tperm2 perm_mx1 mul1mx. without loss /forallP a_dvM0: / [forall j, a %| M 0%R j]%Z. case: (altP forallP) => [_ IH|/forallPn/sigW/IHa IH _]; exact: IH. without loss{Da a_dvM0} Da: M / forall j, M 0 j = a. pose Uur := col' 0 (\row_j (1 - (M 0%R j %/ a)%Z)). pose U : 'M_(1 + n) := block_mx 1 Uur 0 1%:M; pose M1 := M *m U. have uU: U \in unitmx by rewrite unitmxE det_ublock !det1 mulr1. case/(_ (M *m U)) => [j | L uL [R uR [d dvD dM]]]. rewrite -(lshift0 m 0) -[M](@submxK _ 1 _ 1) (@mulmx_block _ 1 m 1). rewrite (@col_mxEu _ 1) !mulmx1 mulmx0 addr0 [ulsubmx _]mx11_scalar. rewrite mul_scalar_mx !mxE !lshift0 Da. case: splitP => [j0 _ | j1 Dj]; rewrite ?ord1 !mxE // lshift0 rshift1. by rewrite mulrBr mulr1 mulrC divzK ?subrK. exists L => //; exists (R * U^-1); first by rewrite unitmx_mul uR unitmx_inv. by exists d; rewrite //= mulmxA -dM mulmxK. without loss{IHa} /forallP/(_ (_, _))/= a_dvM: / [forall k, a %| M k.1 k.2]%Z. case: (altP forallP) => [_|/forallPn/sigW [[i j] /= a'Mij] _]; first exact. have [|||L uL [R uR [d dvD dM]]] := IHa _ _ M^T j; rewrite ?mxE 1?addnC //. by exists i; rewrite mxE. exists R^T; last exists L^T; rewrite ?unitmx_tr //; exists d => //. rewrite -[M]trmxK dM !trmx_mul mulmxA; congr (_ *m _ *m _). by apply/matrixP=> i1 j1 /[!mxE]; case: eqVneq => // ->. without loss{nz_a a_dvM} a1: M a Da / a = 1. pose M1 := map_mx (divz^~ a) M; case/(_ M1 1)=> // [k|L uL [R uR [d dvD dM]]]. by rewrite !mxE Da divzz nz_a. exists L => //; exists R => //; exists [seq a * x | x <- d]. case: d dvD {dM} => //= x d; elim: d x => //= y d IHd x /andP[dv_xy /IHd]. by rewrite [dvdz _ _]dvdz_mul2l ?[_ \in _]dv_xy. have ->: M = a *: M1 by apply/matrixP=> i j; rewrite !mxE mulrC divzK ?a_dvM. rewrite dM scalemxAl scalemxAr; congr (_ *m _ *m _). apply/matrixP=> i j; rewrite !mxE mulrnAr; congr (_ *+ _). have [lt_i_d | le_d_i] := ltnP i (size d); first by rewrite (nth_map 0). by rewrite !nth_default ?size_map ?mulr0. rewrite {a}a1 -[m.+1]/(1 + m)%N -[n.+1]/(1 + n)%N in M Da *. pose Mu := ursubmx M; pose Ml := dlsubmx M. have{} Da: ulsubmx M = 1 by rewrite [_ M]mx11_scalar !mxE !lshift0 Da. pose M1 := - (Ml *m Mu) + drsubmx M. have [|L uL [R uR [d dvD dM1]]] := IHmn m n M1; first by rewrite -addnS ltnW. exists (block_mx 1 0 Ml L). by rewrite unitmxE det_lblock det_scalar1 mul1r. exists (block_mx 1 Mu 0 R). by rewrite unitmxE det_ublock det_scalar1 mul1r. exists (1 :: d); set D1 := \matrix_(i, j) _ in dM1. by rewrite /= path_min_sorted //; apply/allP => g _; apply: dvd1n. rewrite [D in _ *m D *m _](_ : _ = block_mx 1 0 0 D1); last first. by apply/matrixP=> i j; do 3?[rewrite ?mxE ?ord1 //=; case: splitP => ? ->]. rewrite !mulmx_block !(mul0mx, mulmx0, addr0) !mulmx1 add0r mul1mx -Da -dM1. by rewrite addNKr submxK. Qed.
imports: HB(structures) mathcomp(ssreflect ssrbool ssrfun eqtype ssrnat seq path) mathcomp(div choice fintype tuple finfun bigop prime order) mathcomp(ssralg poly ssrnum ssrint archimedean rat matrix) mathcomp(polydiv perm zmodp mxalgebra vector)