module linop; % Linear operator package.
% Author: Anthony C. Hearn.
% Copyright (c) 1987 The RAND Corporation. All rights reserved.
fluid '(!*intstr);
symbolic procedure linear u;
for each x in u do
if not idp x then typerr(x,'operator) else flag(list x,'linear);
rlistat '(linear);
symbolic procedure formlnr u;
begin scalar x,y,z;
x := car u;
if null cdr u or null cddr u
then rerror(alg,29,list("Linear operator",
x,"called with too few arguments"));
y := cadr u;
z := !*a2k caddr u . cdddr u;
return if y = 1 then u
else if not depends(y,car z)
then list('times,y,x . 1 . z)
else if atom y then u
else if car y eq 'plus
then 'plus . for each j in cdr y collect formlnr(x . j. z)
else if car y eq 'minus
then list('minus,formlnr(x . cadr y . z))
else if car y eq 'difference
then list('difference,formlnr(x . cadr y . z),
formlnr(x . caddr y . z))
else if car y eq 'times then formlntms(x,cdr y,z,u)
else if car y eq 'quotient then formlnquot(x,cdr y,z,u)
else if car y eq 'recip
then formlnrecip(x,carx(cdr y,'recip),z,u)
else if y := expt!-separate(y,car z)
then list('times,car y,x . cdr y . z)
else u
end;
symbolic procedure formseparate(u,v);
%separates U into two parts, and returns a dotted pair of them: those
%which are not commutative and do not depend on V, and the remainder;
begin scalar w,x,y;
for each z in u do
if not noncomp z and not depends(z,v) then x := z . x
else if (w := expt!-separate(z,v))
then <<x := car w . x; y := cdr w . y>>
else y := z . y;
return reversip!* x . reversip!* y
end;
symbolic procedure expt!-separate(u,v);
%determines if U is an expression in EXPT that can be separated into
%two parts, one that does not depend on V and one that does,
%except if there is no non-dependent part, NIL is returned;
if not eqcar(u,'expt) or depends(cadr u,v)
or not eqcar(caddr u,'plus)
then nil
else expt!-separate1(cdaddr u,cadr u,v);
symbolic procedure expt!-separate1(u,v,w);
begin scalar x;
x := formseparate(u,w);
return if null car x then nil
else list('expt,v,replus car x) .
if null cdr x then 1 else list('expt,v,replus cdr x)
end;
symbolic procedure formlntms(u,v,w,x);
%U is a linear operator, V its first argument with TIMES removed,
%W the rest of the arguments and X the whole expression.
%Value is the transformed expression;
begin scalar y;
y := formseparate(v,car w);
return if null car y then x
else 'times . aconc!*(car y,
if null cddr y then formlnr(u . cadr y . w)
else u . ('times . cdr y) . w)
end;
symbolic procedure formlnquot(fn,quotargs,rest,whole);
%FN is a linear operator, QUOTARGS its first argument with QUOTIENT
%removed, REST the remaining arguments, WHOLE the whole expression.
%Value is the transformed expression;
begin scalar x;
return if not depends(cadr quotargs,car rest)
then list('quotient,formlnr(fn . car quotargs . rest),
cadr quotargs)
else if not depends(car quotargs,car rest)
and car quotargs neq 1
then list('times,car quotargs,
formlnr(fn . list('recip,cadr quotargs) . rest))
else if eqcar(car quotargs,'plus)
then 'plus . for each j in cdar quotargs
collect formlnr(fn . ('quotient . j . cdr quotargs)
. rest)
else if eqcar(car quotargs,'minus)
then list('minus,formlnr(fn .
('quotient . cadar quotargs . cdr quotargs)
. rest))
else if eqcar(car quotargs,'times)
and car(x := formseparate(cdar quotargs,car rest))
then 'times . aconc!*(car x,
formlnr(fn . list('quotient,mktimes cdr x,
cadr quotargs) . rest))
else if eqcar(cadr quotargs,'times)
and car(x := formseparate(cdadr quotargs,car rest))
then list('times,list('recip,mktimes car x),
formlnr(fn . list('quotient,car quotargs,mktimes cdr x)
. rest))
else if x := expt!-separate(car quotargs,car rest)
then list('times,car x,formlnr(fn . list('quotient,cdr x,cadr
quotargs) . rest))
else if x := expt!-separate(cadr quotargs,car rest)
then list('times,list('recip,car x),
formlnr(fn . list('quotient,car quotargs,cdr x)
. rest))
else if (x := reval!* cadr quotargs) neq cadr quotargs
then formlnquot(fn,list(car quotargs,x),rest,whole)
else whole
end;
symbolic procedure formlnrecip(fn,reciparg,rest,whole);
% FN is a linear operator, RECIPARG the RECIP argument, REST the
% remaining arguments, WHOLE the whole expression. Value is the
% transformed expression.
begin scalar x;
return if not depends(reciparg,car rest)
then list('quotient,fn . 1 . rest,reciparg)
else if eqcar(reciparg,'minus)
then list('minus,formlnr(fn . ('recip . cdr reciparg) . rest))
else if eqcar(reciparg,'times)
and car(x := formseparate(cdr reciparg,car rest))
then list('times,list('recip,mktimes car x),
formlnr(fn . list('recip,mktimes cdr x)
. rest))
else if x := expt!-separate(reciparg,car rest)
then list('times,list('recip,car x),
formlnr(fn . list('recip,cdr x)
. rest))
else if (x := reval!* reciparg) neq reciparg
then formlnrecip(fn,x,rest,whole)
else whole
end;
symbolic procedure mktimes u;
if null cdr u then car u else 'times . u;
symbolic procedure reval!* u;
%like REVAL, except INTSTR is always ON;
begin scalar !*intstr;
!*intstr := t;
return reval u
end;
endmodule;
end;