I'm sure this matrix_factor code can be improved, but so far it seems to work well for me.
/* Maxima's built-in gcd doesn't work so well for fractions, so we gcd the numerator & denominator separately. */
list_gcd(l):=
block([entries,numgcd,denomgcd],
entries:factor(l),
numgcd:apply('ezgcd,map('num,entries))[1],
denomgcd:apply('ezgcd,map('denom,entries))[1],
factor(numgcd/denomgcd));
/* Incase anyone just wants the overall gcd. */
matrix_gcd(m):=list_gcd(list_matrix_entries(m));
/* Return list of 4 items:
overall gcd of all the elements,
diagonal matrix of row gcd's,
reduced matrix,
diagonal matrix of column gcd's.
The input matrix can be recovered by multiplying these 4 factors together. */
matrix_factor(m):=
block([i,colgcds,cols,coli,colgcdi,rowgcds,rows,rowi,rowgcdi,tgcd],
colgcds:[],cols:[],
for i:matrix_size(m)[2] step -1 thru 1 do
(coli:list_matrix_entries(col(m,i)),
colgcdi:list_gcd(coli),
coli:factor(coli/colgcdi),
push(colgcdi,colgcds),
push(coli,cols)),
m:apply('matrix,cols),
m:transpose(m),
tgcd:list_gcd(colgcds),
colgcds:factor(colgcds/tgcd),
rowgcds:[],rows:[],
for i:matrix_size(m)[1] step -1 thru 1 do
(rowi:list_matrix_entries(row(m,i)),
rowgcdi:list_gcd(rowi),
rowi:factor(rowi/rowgcdi),
push(rowgcdi,rowgcds),
push(rowi,rows)),
[tgcd,
apply('diag_matrix,rowgcds),
apply('matrix,rows),
apply('diag_matrix,colgcds)]);
At 11:58 AM 5/6/2013, Henry Baker wrote:
>If anyone is interested in this sort of capability, I'll post the source code, which isn't very long.
>
>This is not a complete factorization, but only pulls out the total gcd, and the row & column gcds.
>
>Maxima 5.28.0-2 http://maxima.sourceforge.net
>using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
>Distributed under the GNU Public License. See the file COPYING.
>Dedicated to the memory of William Schelter.
>The function bug_report() provides bug reporting information.
>(%i1) load("f:\\cg.mac");
>(%o0) f:\cg.mac
>(%i1) M:matrix([a,b],[c,d]);
> [ a b ]
>(%o1) [ ]
> [ c d ]
>(%i2) P:diag_matrix(u,v);
> [ u 0 ]
>(%o2) [ ]
> [ 0 v ]
>(%i3) Q:diag_matrix(x,y);
> [ x 0 ]
>(%o3) [ ]
> [ 0 y ]
>(%i4) n*P.M.Q;
> [ a n u x b n u y ]
>(%o4) [ ]
> [ c n v x d n v y ]
>(%i5) matrix_factor(%);
> [ u 0 ] [ a b ] [ x 0 ]
>(%o5) [n, [ ], [ ], [ ]]
> [ 0 v ] [ c d ] [ 0 y ]
>(%i6)