(* :Title: Decomposing Matrices into Elementary Matrices, Permutation Matrices and a Diagonal Matrix.*) (* :Author: Frank Zizza Willamette University zizza@willamette.edu *) (* :Summary: Decomposes matrices with entries in a Euclidean ring as a product of elementary matrices, permutation matrices and a diagonal matrix. Allowed Euclidean rings are the integers, the Gaussian integers, and the polynomials in one variable over the fields of rational numbers, real numbers or complex numbers. *) (* :Context: ElementaryDecompositions` *) (* :Package Version: 1.0 *) (* :Copyright: Copyright 1992 *) (* :History: Created by Frank Zizza at MSRI, November 1992 *) (* :Keywords: Euclidean ring, elementary matrices, permutation matrices, invariant factors, elementary divisors. *) (* :Source: Jacobson: "Elementary Algebra I", W. H. Freeman & Company. *) (* :Mathematica Version: 2.1*) (* :Limitations: None known. *) BeginPackage["ElementaryDecompositions`"] (*Usage*) ElementaryDecomposition::usage = "ElementaryDecomposition[M] yields the factorization of the matrix M into elementary matrices, permutation matrices and a diagonal matrix. The matrix M must have entries in the ring of integers, or the ring of Gaussian integers. The result is a list {{e[i, j, a], p[i, j] ...}, {{{d[1],0,..}, {0, d[2], 0, ..}}}, {e[i, j, a], p[i, j] ...}}. e[i, j, x] is the identity matrix except in position (i, j) where is has the value x. p[i, j] is the identity matrix with the i-th and j-th columns interchanged. The middle term in the result is a diagonal matrix {{d[1],0,..}, {0, d[2], 0, ..}}. The orginal matrix M is the product of all the matrices in the list in the order in which they are printed. \n \nElementaryDecomposition[M, t] yields the factorization for matrices with entries in the ring of polynomials in one variable t over a field. The field is determined by the coefficients of the polynomials in the input matrix. \n \nElementaryDecomposition[M, Padded -> True] yields the same factorization as ElementaryDecomposition[M] with the elementary matrices and permutation matrices expanded as usual Mathematica matrices." NormalDecomposition::usage = "NormalDecomposition[M] yields a factorization of the r x s matrix M into {P, D, Q} where P is r x r and invertable, D is r x s and diagonal and Q is s x s and invertable. The original matrix M equals P . D . Q. The entries of M must all be integers or Gaussian integers. NormalDecomposition[M, t] works polynomials in a single variable over a field." InvariantFactors::usage = "InvariantFactors[M] gives the diagonal elements of the normal form of the matrix M under the equivalence relation A ~ B iff A == P . B . Q where P and Q are invertable matrices. The entries of M must all be integers or Gaussian integers.InvariantFactors[M, t] works for polynomials in a single variable over a field. It is a standard theorem that all such matrices have a diagonal normal form. The result is a list of elements {d[1], d[2], ... , d[k], 0, .. 0} of the appropriate ring with the divisibility condition d[i] | d[j] for all 1 <= i <= j <= k." ElementaryMatrix::usage = "ElementaryMatrix[i, j, x, n] is a short form for n x n matrix which equals the identity matrix except in position (i, j) where is has the value x. In formatted output, the dimension n is suppressed and i and j are subscripted." PermutationMatrix::usage = "PermutationMatrix[i, j, n] is a short form for n x n matrix that equals the identity matrix with the i-th and j-th columns interchanged. In formatted output, the dimension n is suppressed and i and j are subscripted." Pad::usage = "Pad[elementary or permutation matrix] expands ElementaryMatrix[i, j, x, n] and PermutationMatrix[i, j, n] as actual matrices, i.e. a list of lists." Padded::usage = "Padded is an option for ElementaryDecomposition which specfies whether the elementary matrices and the permutation matrices in the factorization should be given in compressed form as e[i,j,a] and p[i,j] or as actual matrices. Padded -> True gives the actual matrices. The default is Padded -> False." CheckIt::usage = "CheckIt[factorization] multiplies the list of matrices obtained from ElementaryDecomposition." (*Options*) Options[ElementaryDecomposition] = {Padded -> False}; (*Messages*) ElementaryDecomposition::notint = "Matrix does not have entries in the integers or the Gaussian integers. For polynomial entries, append the variable as in ElementaryDecomposition[M, t]." ElementaryDecomposition::notpoly = "Matrix does not have entries in the ring of polynomials in one variable." Begin["`Private`"] SetAttributes[{LeftFactors, RightFactors}, {HoldAll}] MakeFactorList[A_] := Factors[LeftFactors[], A, RightFactors[]] (*Matrix Operations*) FindMinPosition[A_?MatrixQ] := Module[{B, min}, B = $Delta[ A ]; min = Min[ Flatten[B] ]; Position[ B, min, 2, 1]] SubMatrix[A_?MatrixQ, n_Integer] := Map[Drop[#, n-1]&, Drop[A, n-1]] FindMinPosition[A_?MatrixQ, n_Integer] := FindMinPosition[ SubMatrix[A, n] ] + (n-1) SwitchRows[{i_, i_}, fac_Factors] := fac SwitchColumns[{i_, i_}, fac_Factors] := fac SwitchRows[{i_, j_}, A_?MatrixQ] := A[[ Range[ Dimensions[A][[1]] ] /. {i -> j, j -> i} ]] SwitchRows[{i_, j_}, Factors[llist_, A_, rlist_]] := Module[{m}, m = Dimensions[A][[1]]; Factors[Append[llist,PermutationMatrix[i,j,m]],SwitchRows[{i,j}, A], rlist] ] SwitchColumns[{i_, i_}, A_?MatrixQ] := A SwitchColumns[{i_, j_}, A_?MatrixQ] := Transpose[ SwitchRows[{i, j}, Transpose[A] ] ] SwitchColumns[{i_, j_}, Factors[llist_, A_, rlist_]] := Module[{n}, n = Dimensions[A][[2]]; Factors[llist, SwitchColumns[{i, j}, A], Prepend[rlist, PermutationMatrix[i,j,n]] ] ] Reposition[fac:Factors[llist_, A_, rlist_], n_Integer] := Module[{i, j}, {{i, j}} = FindMinPosition[A, n]; SwitchRows[ {n, i}, SwitchColumns[ {n, j}, fac] ] ] RowOperations[fac:Factors[llist_, A_, rlist_], n_Integer] := Module[{ann = A[[n, n]], B = A, m, an, i, q, nllist = llist}, If[ann === 0, fac, m = Dimensions[A][[1]]; For[i=n+1; an = A[[n]], i<=m, i++, If[(q = $Quotient[ A[[i,n]], ann]) =!= 0, AppendTo[nllist, ElementaryMatrix[i, n, q, m]]; B[[i]] = Chop[Simplify[B[[i]] - q*an]]]]; Factors[nllist, B, rlist] ] ] ColumnOperations[fac:Factors[llist_, A_, rlist_], n_Integer] := Module[{ann = A[[n, n]], B = Transpose[A], m, an, i, q, nrlist = rlist}, If[ann === 0, fac, m = Dimensions[A][[2]]; For[i=n+1; an = B[[n]], i<=m, i++, If[(q = $Quotient[B[[i,n]], ann]) =!= 0, PrependTo[nrlist, ElementaryMatrix[n, i, q, m]]; B[[i]] = Chop[Simplify[B[[i]] - q*an]]]]; Factors[llist, Transpose[B], nrlist] ] ] RCOperations[fac_Factors, n_Integer] := ColumnOperations[ RowOperations[ Reposition[fac, n], n], n] OneStep[fac_Factors, n_Integer] := FixedPoint[RCOperations[#, n]&, fac] DivisibleQ[A_?MatrixQ] := Module[{a}, If[ (a = Simplify[A[[1,1]]]) =!= 0, Union[ $Mod[ Flatten[A], a ]] === {0}, True] ] AddRow[{i_, j_}, Factors[llist_, A_, rlist_]] := Module[{m}, m = Dimensions[A][[1]]; Factors[Append[llist, ElementaryMatrix[i, j, -1, m]], ReplacePart[A, A[[i]]+A[[j]], i], rlist]] DivisibilityReposition[fac:Factors[llist_, A_, rlist_], n_] := Module[{B, i, u}, B = SubMatrix[A, n]; If[ DivisibleQ[B], fac, {u} = Select[ Flatten[B], $Mod[#, B[[1,1]] ] =!= 0&, 1]; i = Position[ B, u, 2, 1][[1,1]]; AddRow[{n, n+i-1}, fac] ]] DivisibleOneStep[fac_Factors, n_] := FixedPoint[DivisibilityReposition[OneStep[#1, n], n] &, fac] Diagonalize[(A_)?MatrixQ] := Fold[DivisibleOneStep, MakeFactorList[A], Range[Min[Dimensions[A]]]] (*User Interface*) GaussianIntegerQ[z_] := MatchQ[z, Complex[_Integer, _Integer] | _Integer]; ElementaryDecomposition[matrix_?MatrixQ, opts___Rule] := Module[{ans}, Which[ MatrixQ[matrix, IntegerQ], $Delta = Function[n, If[n=!=0, Abs[n], Infinity], {Listable}]; $Quotient = Function[{n, m}, Quotient[n,m]]; $Mod = Function[{n, m}, Mod[n,m], {Listable}];, MatrixQ[matrix, GaussianIntegerQ], $Delta = Function[n, If[n=!=0, Abs[n]^2, Infinity], {Listable}]; $Quotient = Function[{n, m}, Round[n / m]]; $Mod = Function[{n, m}, n - Round[n / m] * m, {Listable}];, True, Message[ElementaryDecomposition::notint]; Abort[]]; ans = Diagonalize[matrix]; ans = Apply[List, MapAt[List, ans, 2], {0, 1}]; padded = Padded /. {opts} /. Options[ElementaryDecomposition]; If[padded === False, ans, Pad[ans]]] ElementaryDecomposition[matrix_?MatrixQ, t_, opts___Rule] := Module[{ans}, If[Not[MatrixQ[matrix, PolynomialQ[#, t]&]], Message[ElementaryDecomposition::notpoly]; Abort[]]; $Delta = Function[p, If[p=!=0, Exponent[p, t], Infinity], {Listable}]; $Quotient = Function[{p, q}, PolynomialQuotient[p, q, t]]; $Mod = Function[{p,q}, PolynomialRemainder[p, q, t], {Listable}]; ans = Diagonalize[matrix]; ans = Apply[List, MapAt[List, ans, 2], {0, 1}]; padded = Padded /. {opts} /. Options[ElementaryDecomposition]; If[padded === False, ans, Pad[ans]]] NormalDecomposition[matrix_?MatrixQ] := Module[{factors}, factors = ElementaryDecomposition[matrix, Padded -> True]; { Dot @@ factors[[1]], factors[[2, 1]], Dot @@ factors[[3]] } ] NormalDecomposition[matrix_?MatrixQ, t_] := Module[{factors}, factors = ElementaryDecomposition[matrix, t, Padded -> True]; { Dot @@ factors[[1]], factors[[2, 1]], Dot @@ factors[[3]] } ] InvariantFactors[matrix_?MatrixQ] := Module[{mat}, mat = ElementaryDecomposition[matrix][[2, 1]]; Table[mat[[i,i]], {i, 1, Min[Dimensions[matrix]]}]] InvariantFactors[matrix_?MatrixQ, t_] := Module[{mat}, mat = ElementaryDecomposition[matrix, t][[2, 1]]; Table[mat[[i,i]], {i, 1, Min[Dimensions[matrix]]}]] (*Elementary Matrix Padding*) SetAttributes[Pad, {Listable}]; Pad[PermutationMatrix[i_, j_, n_]] := IdentityMatrix[n][[ Range[n] /. {i -> j, j -> i} ]] Pad[ElementaryMatrix[i_, j_, a_, n_]] := ReplacePart[IdentityMatrix[n], a, {i, j}] Pad[n_] := n (*Testing*) CheckIt[factors_] := With[{fac = Pad[factors]}, (Dot @@ fac[[1]]) . fac[[2,1]] . (Dot @@ fac[[3]])] (*Formats*) Format[PermutationMatrix[i_,j_,n_]] := Subscripted["p"[i,j]] Format[ElementaryMatrix[i_,j_,a_, n_]] := Subscripted["e"[i, j, a], {1, 2}] End[] (*`private`*) Protect[ElementaryDecomposition, NormalDecomposition, InvariantFactors] EndPackage[] (*ElementaryDecompositions`*) (**********************************************************) (*:Extras: MinimumPolynomial[A_?MatrixQ, t_Symbol] := Module[{m, p, c}, m = Length[A]; p = Last[ InvariantFactors[t IdentityMatrix[m] - A, t] ]; p/Last[CoefficientList[p, t]]] /; SameQ[Dimensions[A]] *)