Here's a smaller & slightly more efficient version of the matrix_factor function.
/* Helper function for matrix_factor(m).
Returns 4 items:
overall gcd,
list of row gcds,
reduced matrix,
list of column gcds.
matrix_factor1(m) assumes that m is factored, and has NO DENOMINATORS. */
matrix_factor1(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:apply('ezgcd,list_matrix_entries(col(m,i))),
colgcdi:coli[1], coli:rest(coli),
push(colgcdi,colgcds),
push(coli,cols)),
m:apply('matrix,cols),
m:transpose(m),
colgcds:apply('ezgcd,colgcds),
tgcd:colgcds[1], colgcds:rest(colgcds),
rowgcds:[],rows:[],
for i:matrix_size(m)[1] step -1 thru 1 do
(rowi:apply('ezgcd,list_matrix_entries(row(m,i))),
rowgcdi:rowi[1], rowi:rest(rowi),
push(rowgcdi,rowgcds),
push(rowi,rows)),
[tgcd,
rowgcds,
apply('matrix,rows),
colgcds]);
/* matrix_factor(m).
Returns 4 items:
G overall gcd,
P row gcds in a diagonal matrix,
M reduced matrix,
Q column gcds in a diagonal matrix.
m = G*P.M.Q */
matrix_factor(m):=
block([r],
m:factor(m),
r:map("/",nfactors:matrix_factor1(matrixmap('num,m)),
dfactors:matrix_factor1(matrixmap('denom,m))),
[r[1],apply('diag_matrix,r[2]),r[3],apply('diag_matrix,r[4])]);
At 06:41 PM 5/6/2013, Henry Baker wrote:
>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)