module matsm; % Simplification of matrices.
% Author: Anthony C. Hearn.
% Copyright (c) 1998 Anthony C. Hearn. All rights reserved.
% This module is rife with essential references to RPLAC-based
% functions.
symbolic procedure matsm!*(u,v);
% Matrix expression simplification function.
matsm!*1 matsm u;
% symbolic procedure matsm!*1 u;
% begin scalar sub2;
% sub2 := !*sub2; % Since we need value for each element.
% u := 'mat . for each j in u collect
% for each k in j
% collect <<!*sub2 := sub2; !*q2a subs2 k>>;
% !*sub2 := nil; % Since all substitutions done.
% return u
% end;
symbolic procedure matsm!*1 u;
begin
% We use subs2!* to make sure each element simplified fully.
u := 'mat . for each j in u collect
for each k in j collect !*q2a subs2!* k;
!*sub2 := nil; % Since all substitutions done.
return u
end;
symbolic procedure matsm u;
begin scalar x,y;
for each j in nssimp(u,'matrix) do
<<y := multsm(car j,matsm1 cdr j);
x := if null x then y else addm(x,y)>>;
return x
end;
symbolic procedure matsm1 u;
%returns matrix canonical form for matrix symbol product U;
begin scalar x,y,z; integer n;
a: if null u then return z
else if eqcar(car u,'!*div) then go to d
else if atom car u then go to er
else if caar u eq 'mat then go to c1
else <<x := lispapply(caar u,cdar u); %
if eqcar(x,'mat) then x := matsm x>>; %
b: z := if null z then x
else if null cdr z and null cdar z then multsm(caar z,x)
else multm(x,z);
c: u := cdr u;
go to a;
c1: if not lchk cdar u then rerror(matrix,3,"Matrix mismatch");
x := for each j in cdar u collect
for each k in j collect xsimp k;
go to b;
d: y := matsm cadar u;
if (n := length car y) neq length y
then rerror(matrix,4,"Non square matrix")
else if (z and n neq length z)
then rerror(matrix,5,"Matrix mismatch")
else if cddar u then go to h
else if null cdr y and null cdar y then go to e;
x := subfg!*;
subfg!* := nil;
if null z then z := apply1(get('mat,'inversefn),y)
else if null(x := get('mat,'lnrsolvefn))
then z := multm(apply1(get('mat,'inversefn),y),z)
else z := apply2(get('mat,'lnrsolvefn),y,z);
subfg!* := x;
% Make sure there are no power substitutions.
z := for each j in z collect for each k in j collect
<<!*sub2 := t; subs2 k>>;
go to c;
e: if null caaar y then rerror(matrix,6,"Zero divisor");
y := revpr caar y;
z := if null z then list list y else multsm(y,z);
go to c;
h: if null z then z := generateident n;
go to c;
er: rerror(matrix,7,list("Matrix",car u,"not set"))
end;
symbolic procedure lchk u;
begin integer n;
if null u or atom car u then return nil;
n := length car u;
repeat u := cdr u
until null u or atom car u or length car u neq n;
return null u
end;
symbolic procedure addm(u,v);
% Returns sum of two matrix canonical forms U and V.
% Returns U + 0 as U. Patch by Francis Wright.
if v = '(((nil . 1))) then u else % FJW.
for each j in addm1(u,v,function cons)
collect addm1(car j,cdr j,function addsq);
symbolic procedure addm1(u,v,w);
if null u and null v then nil
else if null u or null v then rerror(matrix,8,"Matrix mismatch")
else apply2(w,car u,car v) . addm1(cdr u,cdr v,w);
symbolic procedure tp u; tp1 matsm u;
put('tp,'rtypefn,'getrtypecar);
symbolic procedure tp1 u;
%returns transpose of the matrix canonical form U;
%U is destroyed in the process;
begin scalar v,w,x,y,z;
v := w := list nil;
while car u do
<<x := u;
y := z := list nil;
while x do
<<z := cdr rplacd(z,list caar x);
x := cdr rplaca(x,cdar x)>>;
w := cdr rplacd(w,list cdr y)>>;
return cdr v
end;
symbolic procedure scalprod(u,v);
%returns scalar product of two lists (vectors) U and V;
if null u and null v then nil ./ 1
else if null u or null v then rerror(matrix,9,"Matrix mismatch")
else addsq(multsq(car u,car v),scalprod(cdr u,cdr v));
symbolic procedure multm(u,v);
%returns matrix product of two matrix canonical forms U and V;
(for each y in u
collect for each k in x collect subs2 scalprod(y,k))
where x = tp1 v;
symbolic procedure multsm(u,v);
%returns product of standard quotient U and matrix standard form V;
if u = (1 ./ 1) then v
else for each j in v collect for each k in j collect multsq(u,k);
% Explicit substitution code for matrices.
symbolic procedure matsub(u,v);
'mat . for each x in cdr v collect
for each y in x collect subeval1(u,y);
put('matrix,'subfn,'matsub);
endmodule;
end;