%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% P H Y S O P %
% %
% A Package for Operator Calculus %
% in Physics %
% %
% Author: Mathias Warns %
% Physics Institute %
% University of Bonn %
% Nussallee 12 %
% D-5300 BONN 1 (F.R.G.) %
% <warns@unicc.bitnet> %
% %
% Version: 1.5 06 Jan 1992 %
% %
% Designed for: REDUCE version 3.4 %
% Tested on : - Intel 386/486 AT compatible computers %
% PSL implementation of REDUCE 3.4 %
% - IBM 3084/9000-620 MVS/XA %
% PSL implementation of REDUCE 3.4 %
% %
% CAUTION: (i) The NONCOM2 package is needed to run this package %
% (ii) This package cannot be used simultaneously with %
% packages modifying the standard GETRTYPE procedure %
% %
% Copyright (c) Mathias Warns 1990 - 1992 %
% %
% Permission is granted to any individual or institution to %
% use, copy or re-distribute this software as long as it is %
% not sold for profit, provided that this copyright notice %
% is retained and the file is not altered. %
% %
% *** Revision history since issue of Version 0.99 *** %
% %
% - sloppy use of CAR on atoms corrected in various procedures %
% - MUL and TSTACK added in PHYSOPTIMES %
% - Bug in CLEARPHYSOP corrected %
% - ordering procedures recoded for greater efficiency %
% - handling of PROG expressions included via %
% procedure PHYSOPPROG %
% - procedures PHYSOPTIMES and MULTOPOP!* modified %
% - extended error handling inclued via REDERR2 %
% - PHYSOPTYPELET recoded %
% - PHYSOPCONTRACT modified for new pattern natcher %
% - EQ changed to = in MULTF and MULTFNC %
% - PHYSOPCOMMUTE/PHYSOPANTICOMMUTE and COMM2 corrected %
% - Handling of SUB and output printing adapted to 3.4 %
% %
% 1.1 130791 mw %
% - Modifications for greater efficiency in procedures ORDOP, %
% ISANINDEX and ISAVARINDEX %
% - PHYSOP2SQ slightly modified for greater efficiency %
% - Procedure COLLECTPHYSTYPE added %
% - handling of inverse and adjoint operators modified %
% procedures INV and ADJ2 modified %
% procedures INVP and ADJP recoded %
% - procedures GETPHYSTYPE!*SQ and GETPHYSTYPESF added for greater %
% efficiency in type checking of !*SQ expressions %
% - procedure GETPHYTYPE modified accordingly %
% - SIMP!* changed to SUBS2 in procedure PHYSOPSUBS %
% - Bug in EXPTEXPAND and PHYSOPEXPT corrected %
% - PHYSOPORDCHK and PHYSOPSIMP slightly enhanced %
% - PHYSOPTYPELET enhanced (COND treatment) %
% - phystypefn for PLUS and DIFFERENCE changed to GETPHYSTYPEALL %
% - GETPHYSTYPEALL added %
% - GETPHYSTYPETIMES modified %
% 1.2 190891 mw %
% - implementation of property PHYSOPNAME for PHYSOPs %
% - procedures SCALOP,VECOP,TENSOP,STATE,INV,ADJ2,INVADJ modified %
% - procedure ORDOP recoded, NCMPCHK and PHYSOPMEMBER modified %
% - procedure PHYSOPSM!* enhanced %
% 1.3 test implementation of a new ordering scheme 260891 mw %
% - Procedure OPNUM!* and RESET!_OPNUMS added %
% - procedure ORDOP recoded %
% - procedure SCALOP,VECOP,TENSOP,STATE,OPORDER modified %
% - procedure !*XADD added %
% - procedure PHYSOPSIMP corrected %
% 1.4 181291 mw %
% - bug in procedures SCALOPP, PHYSOPSIMP and TENSOP corrected %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
module physop;
%-------------------------------------------------------%
% This part has to be modified by the user if required %
%-------------------------------------------------------%
% input the NONCOM2 package here for a compiled version
% input noncom2;
load_package noncom2;
% Modify the infix character for the OPAPPLY function if needed
newtok'((|) opapply);
flag('(opapply), 'spaced);
%-------------------------------------------------------%
% E N D of user modifiable part %
%-------------------------------------------------------%
%**************** the following is needed for REDUCE 3.4 *************
fluid '(!*nosq); % controls proper use of !*q2a
!*nosq := t;
% ************** end of 3.4 modifications **************************
fluid '(frlis!*);
newtok '((d o t) dot);
flag ('(dot), 'spaced);
% define some global variables from REDUCE needed in the package
fluid '(alglist!*);
Global '(tstack!* mul!*);
% ---define global variables needed for the package---
FLuid '(oporder!* defoporder!* physopindices!* physopvarind!*);
Fluid '(physoplist!*);
Global '(specoplist!*);
% define global flags
fluid '(!*anticom !*anticommchk !*contract !*contract2 !*hardstop);
fluid '(!*indflg indcnt!*);
indcnt!* := 0;
% additional flag needed for contraction
!*contract2 := nil;
% flag indicating that one elementary comm or opapply has not
% been found --> print warning message
!*hardstop := nil;
% this are algebraic mode switches
switch contract;
switch anticom;
% reserved operators and variables;
% idx is the basic identifier for system created indices
Global '(idx);
% ----- link new data type PHYSOP in REDUCE ------
% physop is the new datatype containing all subtypes
put('physop,'name,'physop); %datatype name
put('physop,'evfn,'!*physopsm!*); % basic simplification routine
put('physop,'typeletfn,'physoptypelet); % routine for type assignements
% ----RLISP procedures which have been modified -----
% procedure for extended error handling
symbolic procedure rederr2(u,v);
begin
msgpri("Error in procedure ",u, ": ", nil,nil);
rederr v
end ;
% procedures multf and multfnc have to be redefined to avoid
% contraction of terms after exptexpand
symbolic procedure multf(u,v); % changed
%U and V are standard forms.
%Value is standard form for U*V;
begin scalar ncmp,x,y;
a: if null u or null v then return nil
else if u=1 then return v % ONEP
else if v=1 then return u % ONEP
else if domainp u then return multd(u,v)
else if domainp v then return multd(v,u)
else if not(!*exp or ncmp!* or wtl!* or x)
then <<u := mkprod u; v := mkprod v; x := t; go to a>>;
x := mvar u;
y := mvar v;
% if (ncmp := noncomp y) and noncomp x then return multfnc(u,v)
if noncommuting(x,y) then return multfnc(u,v)
% we have to put this clause here to prevent evaluation in case
% of equal main vars
else if noncommutingf(y, lc u) or (ordop(x,y) and (x neq y))
then << x := multf(lc u,v);
y := multf(red u,v);
return if null x then y else lpow u .* x .+ y>>
else if x = y and (not physopp x or !*contract2)
% two forms have the same mvars
% switch contract added here to inhibit power contraction
% if not wanted (used by PHYSOP)
then << x := mkspm(x,ldeg u+ldeg v);
y := addf(multf(red u,v),multf(!*t2f lt u,red v));
return if null x or null(u := multf(lc u,lc v))
then <<!*asymp!* := t; y>>
else if x=1 then addf(u,y)
else if null !*mcd then addf(!*t2f(x .* u),y)
else x .* u .+ y>>;
x := multf(u,lc v);
y := multf(u,red v);
return if null x then y else lpow v .* x .+ y
end;
symbolic procedure multfnc(u,v);
%returns canonical product of U and V, with both main vars non-
%commutative;
begin scalar x,y;
x := multf(lc u,!*t2f lt v);
if null x
then return addf(multf(red u,v),multf(!*t2f lt u,red v));
% switch contract added here to avoid contraction of equal powers
% used by PHYSOP
return addf((if not domainp x and (mvar x = mvar u) and
((not physopp mvar x) or !*contract2)
then addf(if null (y := mkspm(mvar u,ldeg u+ldeg v))
then nil
else if y = 1 then lc x
else !*t2f(y .* lc x),
multf(!*p2f lpow u,red x))
else !*t2f(lpow u .* x)),
addf(multf(red u,v),multf(!*t2f lt u,red v)))
end;
symbolic procedure opmtch!* u;
% same as opmtch but turns subfg!* on
begin scalar x,flg;
flg:= subfg!*; subfg!* := t;x:= opmtch u; subfg!* := flg;
return x
end;
symbolic procedure reval3 u;
% this is just a redefinition of reval2(u,nil)
% which call simp instead of simp!*
% it saves at lot of writing in some procedures
mk!*sq x where x := simp u;
% ---- procedure related to ordering of physops in epxr -------
symbolic procedure oporder u;
% define a new oporder list
begin
if not listp u then rederr2('oporder, "invalid argument to oporder");
if (u = '(nil)) then oporder!* := defoporder!* % default list
else for each x in reverse u do <<
if not physopp x then rederr2('oporder,
list(x," is not a PHYSOP"));
oporder!* := nconc(list(x),physopdelete(x,oporder!*)) >>; %1.01
% write "oporder!* set to: ",oporder!*;terpri();
reset!_opnums(); %1.03
rmsubs()
end;
rlistat '(oporder);
symbolic procedure physopdelete(u,v);
% u is a physop, v is a list of physops
% deletes u from v
if atom u then delete(u,v)
else
delete(u,delete(car u,delete(removeindices(u,collectindices u),v)));
symbolic procedure opnum!* u; % new 1.03
begin scalar op,arglist;
if not idp u then u := removeindices(u,collectindices u);
if idp u then op := u
else << op := car u; arglist := cdr u;>>;
return
if null (u:= assoc(arglist,get(op,'opnum))) then
cdr assoc(nil,get(op,'opnum))
else cdr u
end;
symbolic procedure reset!_opnums();
begin scalar x,lst,n,op,arglist;
lst := oporder!*;
n := 1;
a: if null lst then return;
x := car lst; lst := cdr lst;
if idp x then <<op := x; arglist := nil>>
else <<op := car x; arglist := cdr x>>;
put(op,'opnum,!*xadd((arglist . n),get(op,'opnum)));
n:= n+1;
go to a
end;
symbolic procedure !*xadd(u,v); % new 1.03
% u is assignement , v is a table
% returns updated table
begin scalar x;
x := v;
while x and not (car u = caar x) do x := cdr x;
if x then v := delete(car x,v);
v := u . v;
return v
end;
symbolic procedure ordop(u,v);
% recoded ordering procedure of operators
% checks new list oporder!* for physops or calls ordop2
% default is to put anything ahead of physops
% we use !*physopp instead of physopp in order to use
% ordop even if we hide the physop rtype
begin scalar x,y,z,nx,ny;
% this are the trivial cases
if not (!*physopp u and !*physopp v) then return
if !*physopp u then nil
else if !*physopp v then t
else ordop2(u,v);
% now come the cases with 2 physops
% following section modified 1.02
if idp u then x:= get(u,'physopname)
else
<<
x:=get(car u,'physopname);
x:= x . cdr u;
u := car u;
>>;
if member(u,specoplist!*) then return t;
if idp v then y:= get(v,'physopname)
else
<<
y:= get(car v, 'physopname);
y := y . cdr v;
v := car v;
>>;
if member(v,specoplist!*) then return t;
% end of modifications 1.02
% from here it is 1.03
nx := opnum!* x;
ny := opnum!* y;
return
if nx < ny then t
else if nx > ny then nil
else if idp x then t
else if idp y then nil
else ordop(cdr x, cdr y);
end;
symbolic procedure ordop2(u,v);
% this is nothing but the standard ordop procedure
begin scalar x;
x := kord!*;
a: if null x then return ordp(u,v)
else if u eq car x then return t
else if v eq car x then return;
x := cdr x;
go to a
end;
% obsolete in 1.03
%symbolic procedure physopmember(u,v); % 1.02 order modified
% u is a physop, v is a list
% return part of v starting with u
%member(u,v) or ((not atom u) and (member(car u,v)
% or member(removeindices(u,collectindices u),v)));
symbolic procedure physopordchk(u,v); % new version 080591
% u and v are physopexpr
% builds up a list of physops of u and v
% checks if there is a pair of wrong ordered noncommuting operators
% in these lists
begin scalar x,y,z,oplist,lst;
x := deletemult!* !*collectphysops u; %1.01
y := deletemult!* !*collectphysops v; % 1.01
return
if null x then t
else if null y then nil
else if member('unit,x) or member('unit,y) then nil %further eval needed
else physopordchk!*(x,y);
end;
symbolic procedure ncmpchk(x,y); % order changed 1.02
% x and y are physops
% checks for correct ordering in noncommuting case
(not noncommuting(x,y)) or ordop(x,y);
symbolic procedure physopordchk!*(u,v);
% u and v are lists of physops
% checks if there is a pair of wrong ordered noncommuting operators
% in this list
begin scalar x,y,lst;
x:= car u; u := cdr U;
if null u then
if null cdr v then
return (ncmpchk(x,car v) and not (invp x = car v))
else
<<
lst := for each y in v collect ncmpchk(x,y);
if member(nil,lst) then return nil
else return t
>>
else return (physopordchk!*(list(x),v) and physopordchk!*(u,v));
end;
% ---general testing functions for PHYSOP expressions----
symbolic procedure physopp u;
if atom u then (idp u and (get(u,'rtype) eq 'physop))
else (idp car u and (get(car u,'rtype) eq 'physop));
% slightly more general
symbolic procedure !*physopp u;
% used to determine physops when physop rtype is hidden
if atom u then (idp u and get(u,'phystype))
else (idp car u and get(car u,'phystype));
symbolic procedure physopp!* u;
physopp u or (not atom u and (flagp(car u,'physopfn) or
(flagp(car u,'physoparith) and
hasonephysop cdr u) or (flagp(car u,'physopmapping) and
hasonephysop cdr u)));
symbolic procedure !*physopp!* u;
physopp!* u or getphystype u;
symbolic procedure hasonephysop u;
if null u then nil
else (physopp!* car u) or hasonephysop cdr u;
symbolic procedure areallphysops u;
if null u then nil
else if null cdr u then !*physopp!* car u
else (!*physopp!* car u) and areallphysops cdr u;
% *****defining functions for different data subtypes******
% scalar operator
symbolic procedure scalop u;
begin scalar y;
for each x in u do
if not idp x then
msgpri("cannot declare",x,"a scalar operator",nil,nil)
else if physopp x then
msgpri(x,"already declared as",get(x,'phystype),nil,nil)
else <<y :=gettype x;
if y memq '(matrix operator array procedure) then
msgpri(x,"already defined as",y,nil,nil)
else <<
put(x,'rtype,'physop);
put(x,'phystype,'scalar);
put(x,'psimpfn,'physopsimp);
put(x,'physopname,x); % 1.02
defoporder!* := nconc(defoporder!*,list(x));
oporder!* := nconc(oporder!*,list(x));
physoplist!* := nconc(physoplist!*,list(x));
inv x; adj2 x; invadj x; %1.01
reset!_opnums(); %1.03
>>;
>>;
return nil
end;
symbolic procedure scalopp u;
(idp u and get(u,'phystype) = 'scalar) or (not atom u and (
(get(car u,'phystype) = 'scalar) or
((get(car u,'phystype) = 'vector) and isanindex cadr u)
or ((get(car u,'phystype) = 'tensor) and
(length(cdr u) >= get(car u,'tensdimen)) and
areallindices(for k:=1 :get(car u,'tensdimen) collect nth(cdr u,k)))));
symbolic procedure vecop u;
begin scalar y;
for each x in u do
if not idp x then
msgpri("cannot declare",x,"a vector operator",nil,nil)
else if physopp x then
msgpri(x,"already declared as",get(x,'phystype),nil,nil)
else <<y :=gettype x;
if y memq '(matrix operator array procedure) then
msgpri(x,"already defined as",y,nil,nil)
else <<put(x,'rtype,'physop);
put(x,'phystype,'vector);
put(x,'psimpfn,'physopsimp);
put(x,'physopname,x); % 1.02
defoporder!* := nconc(defoporder!*,list(x));
oporder!* := nconc(oporder!*,list(x));
physoplist!* := nconc(physoplist!*,list(x));
inv x; adj2 x; invadj x; %1.01
reset!_opnums();
>>;
>>;
return nil
end;
symbolic procedure vecopp u;
(idp u and (get(u,'phystype) = 'vector)) or (not atom u and
((get(car u,'phystype) ='vector) and not isanindex cadr u));
symbolic procedure tensop u;
begin scalar y,n;
% write "car u=",car u;terpri();
for each x in u do
<<
if idp x or not numberp cadr x then
msgpri("Tensor operator",x,"declared without dimension",nil,nil)
else
<<
n:= cadr x; x:= car x;
if not idp x then
msgpri("cannot declare",x,"a tensor operator",nil,nil)
else if physopp x then
msgpri(x,"already declared as",get(x,'phystype),nil,nil)
else
<<
y :=gettype x;
if y memq '(matrix operator array procedure) then
msgpri(x,"already defined as",y,nil,nil)
else
<<
put(x,'rtype,'physop);
put(x,'phystype,'tensor);
put(x,'psimpfn,'physopsimp);
put(x,'physopname,x); % 1.02
put(x,'tensdimen,n);
defoporder!* := nconc(defoporder!*,list(x));
oporder!* := nconc(oporder!*,list(x));
physoplist!* := nconc(physoplist!*,list(x));
inv x; adj2 x; invadj x; %1.01
reset!_opnums();
>>
>>
>>
>>;
return nil
end;
symbolic procedure tensopp u;
(idp u and (get(u,'phystype) = 'tensor)) or (not atom u and
((get(car u,'phystype) ='tensor) and not isanindex cadr u));
symbolic procedure state u;
begin scalar y;
for each x in u do
if not idp x then msgpri("cannot declare",x,"a state",nil,nil)
else if physopp x then
msgpri(x,"already declared as",get(x,'phystype),nil,nil)
else <<y :=gettype x;
if y memq '(matrix operator array procedure) then
msgpri(x,"already defined as",y,nil,nil)
else <<put(x,'rtype,'physop);
put(x,'phystype,'state);
put(x,'psimpfn,'physopsimp);
put(x,'physopname,x); % 1.02
defoporder!* := nconc(defoporder!*,list(x));
oporder!* := nconc(oporder!*,list(x));
physoplist!* := nconc(physoplist!*,list(x));
adj2 x;
reset_opnums(); % 1.03
>>
>>;
return nil
end;
symbolic procedure statep u;
(idp u and get(u,'phystype) = 'state) or (not atom u and
(idp car u and get(car u,'phystype) = 'state));
symbolic procedure statep!* u;
% slightly more general since state may be hidden in another operator
(getphystype u = 'state);
% some procedures for vecop and tensop indices
symbolic procedure physindex u;
begin scalar y;
for each x in u do <<
if not idp x then msgpri("cannot declare",x,"an index",nil,nil)
else if physopp x then
msgpri(x,"already declared as",get(x,'phystype),nil,nil)
else <<y :=gettype x;
if y memq '(matrix operator array procedure) then
msgpri(y,"already defined as",y,nil,nil)
else putanewindex x >>
>>;
return nil
end;
symbolic procedure physindexp u;
% boolean function to test if an id is a physindex
% in algebraic mode
if idp u and isanindex u then t
else if idlistp u and areallindices u then t
else nil;
flag ('(physindexp),'opfn);
flag ('(physindexp),'boolean);
deflist('((scalop rlis) (vecop rlis) (tensop rlis)
(state rlis) (physindex rlis)),'stat);
symbolic procedure isanindex u; %recoded 1.01
idp u and (memq(u,physopindices!*) or member(u,physopvarind!*)
or (memq(u,frlis!*) and member(revassoc(u,frasc!*),
physopindices!*)));
symbolic procedure isavarindex u; % recoded 1.01
member(u,physopvarind!*);
symbolic procedure areallindices u;
isanindex car u and (null cdr u or areallindices cdr u);
symbolic procedure putanewindex u;
% makes a new index available to the system
begin scalar indx;
indx := u;
if isanindex indx then nil
else if (not atom indx) or getrtype indx then
rederr2('putanewindex,list(indx,"is not an index"))
else physopindices!* := nconc(physopindices!*,list(indx));
return nil
end;
symbolic procedure putanewindex!* u;
% used by ISANINDEX to recognize unresolved IDXn indices
begin scalar x;
if not idp u then return;
x:= explode u;
if length(x) < 4 then return;
x := for j:= 1 : 3 collect nth(x,j);
if not(x='(I D X ) or x='(!i !d !x)) then return; % check both cases.
physopindices!* := nconc(physopindices!*,list(u));
return t
end;
symbolic procedure makeanewindex();
% generates a new index
begin scalar x,n;
n:=0;
a: n:=n+1;
x:= mkid('idx,n);
if isanindex x then go to a
else putanewindex x;
return x
end;
symbolic procedure makeanewvarindex();
% generates a new variable index
% for patterm matching
% physopvarind!* keeps var indices to avoid inflation
begin scalar x,y,n;
n:=0;
y:= makeanewindex();
x := intern compress append(explode '!=,explode y);
nconc(frlis!*,list(x));
physopvarind!*:= nconc(physopvarind!*,list(x));
frasc!* := nconc(frasc!*,list((y . x)));
return x
end;
symbolic procedure getvarindex n;
begin scalar ilen;
if not numberp n then rederr2 ('getvarindex,
"invalid argument to getvarindex");
ilen := length(physopvarind!*);
return
if n > ilen then makeanewvarindex()
else nth(physopvarind!*,n);
end;
symbolic procedure transformvarindex u;
% u is a free index
% looks for the corresponding index on the frasc
% or creates a new one
begin scalar x;
x := explode u;
if length(x) < 3 or not nth(x,2) eq '= then return u;
x := intern compress pnth(x,3);
putanewindex x;
if not atsoc(x,frasc!*) then
frasc!* := nconc(frasc!*,list((x . u)));
return x
end;
symbolic procedure insertindices(u,x);
% u is a vecopp or tensopp
% x is an index or a list of indices
if (idp x and not isanindex x) or (idlistp x and not areallindices x)
then rederr2('insertindices, "invalid indices to insertindex")
else if vecopp u then if idp u then list(u,x)
else car u . ( x . cdr u)
else if tensopp u then if idp u then u . x
else car u . nconc(x,cdr u)
% do not insert any index in states or scalops
else u;
symbolic procedure insertfreeindices(u,flg);
% procedure to transform vecop and tensop into scalops
% by inserting free indices taken from the varindlist
% flg is set to t if variable indices are requested
begin scalar n,x;
if vecopp u then <<x:= if flg then
transformvarindex getvarindex(indcnt!* + 1)
else getvarindex(indcnt!* + 1);
return insertindices(u,x)>>
else if tensopp u then <<n:= get(u,'tensdimen);
x:= for k:=1 :n collect if flg then
transformvarindex getvarindex(indcnt!* +k)
else getvarindex(indcnt!* +k);
return insertindices(u,x) >>
else rederr2('insertfreeindices,
"invalid argument to insertfreeindices");
end;
symbolic procedure collectindices u;
% makes a list of all indices in a.e. u
begin scalar v,x;
if atom u then
if isanindex u then return list(u)
else return nil;
a: v := car u;
u := cdr u;
x :=nconc(x,collectindices v);
if null u then return x;
go to a
end;
symbolic procedure removeindices(u,x);
% u is physop (scalop) containing physindices
% x is an index or a list of indices
begin scalar op;
trwrite('removeindices,"u= ",u," x= ",x);
if null x or idp u or not !*physopp u then return u;
if (idp x and not isanindex x) or (idlistp x and not areallindices x)
then rederr2('removeindices, "invalid arguments to removeindices");
op:=car u;u := cdr u;
if null u then return op;
if idp x then u := delete(x,u)
else for each y in x do u:= delete(y,u);
return if null u then op else op . u
end;
symbolic procedure deadindices u;
% checks an a.e. u to see if there are dead indices
% i.e. indices appearing twice or more
%returns the list of dead indices in u
begin scalar x,res;
if null u or atom u then return nil;
x := collectindices u;
for each y in x do
if memq(y,memq(y,x)) then res :=nconc(res,list(y));
return res
end;
symbolic procedure collectphysops u;
% makes a list of all physops in a.e. u
begin scalar v,x;
if atom u then
if physopp u then return list(u)
else return nil
else if physopp u then return list(removeindices(u,collectindices u));
a: v := car u;
u := cdr u;
x :=nconc(x,collectphysops v);
if null u then return x;
go to a
end;
symbolic procedure !*collectphysops u;
% makes a list of all physops in a.e. u
% with ALL indices
begin scalar v,x;
if physopp u then return list(u);
if atom u then return nil;
a: v := car u;
u := cdr u;
x :=nconc(x,!*collectphysops v);
if null u then return x;
go to a
end;
symbolic procedure collectphysops!* u;
begin scalar x;
x:= for each y in collectphysops u collect if idp y then y
else car y;
return x
end;
symbolic procedure collectphystype u; % new 1.01
% makes a list of all physops in u
% with ALL indices
if physopp u then list(getphystype u)
else if atom u then nil
else deletemult!* (for each v in u collect getphystype v);
% ---- PHYSOP procedures for type check and assignement ----
% modify the REDUCE GETRTYPE routine to get control over PHYSOP
% expressions
symbolic procedure getrtype u; %modified
% Returns overall algebraic type of u (or NIL if expression is a
% scalar). Analysis is incomplete for efficiency reasons.
% Type conflicts will later be resolved when expression is evaluated.
begin scalar x,y;
return
if atom u
then if not idp u then nil
else if flagp(u,'share) then getrtype eval u
else if x := get(u,'rtype)
then if y := get(x,'rtypefn) then apply1(y,nil)
else x
else nil
else if not idp car u then nil
else if physopp!* u then 'physop % added
else if (x := get(car u,'rtype)) and (x := get(x,'rtypefn))
then apply1(x, cdr u)
else if x := get(car u,'rtypefn) then apply1(x, cdr u)
else nil
end;
symbolic procedure getrtypecadr u;
not atom u and getrtype cadr u;
symbolic procedure getnewtype u;
not atom u and get(car u,'newtype);
symbolic procedure getphystype u;
% to get the type of a PHYSOP object
begin scalar x;
return
if physopp u then
if scalopp u then 'scalar
else if vecopp u then 'vector
else if tensopp u then 'tensor
else if statep u then 'state
else nil
else if atom u then nil
% following line suppressed 1.01
% else if car u = '!*sq then return getphystype physopaeval u
else if (x:=get(car u,'phystype)) then x
else if (x:=get(car u,'phystypefn)) then
apply1(x,cdr u)
% from here it is 1.01
else if null (x := collectphystype u) % 1.01
then nil
else if null cdr x then car x
else if member('state,x) then 'state
else rederr2('getphystype,list(
"PHYSOP type conflict in",u));
end;
symbolic procedure getphystypecar u;
not atom u and getphystype car u;
symbolic procedure getphystypeor u;
not atom u and (getphystype car u or getphystypeor cdr u);
symbolic procedure getphystypeall args; % new 1.01
begin scalar x;
if null (x := collectphystype deleteall(0,args)) then
return nil
else if cdr x then
rederr2('getphystypeall,
list("PHYSOP type mismatch in",args))
else return car x
end;
% ***** dirty trick *****
% we introduce a rtypefn for !*sq expressions to get
% proper type checking in assignements
symbolic procedure physop!*sq U;
% u is a !*sq expressions
% checks if u contains physops
begin scalar x;
x:= !*collectphysops !*q2a car u;
return
if null x then nil
else 'physop
end;
deflist('((!*sq physop!*sq)), 'rtypefn);
% 1.01 we add also a phystypefn for !*sq
symbolic procedure getphystype!*sq u; % new 1.01
getphystypesf caar u;
deflist('((!*sq getphystype!*sq)), 'phystypefn);
symbolic procedure getphystypesf u; % new 1.01
% u is a s.f.
% returns the phystype of u
if null u or domain!*p u then nil
else getphystype mvar u or getphystypesf lc u;
%-----end of 1.01 modifications -----------------
% we have also to modify the simp!*sq routine since
% there is no type checking included
symbolic procedure physopsimp!*sq u;
if cadr u then car u
else if physop!*sq u then physop2sq physopsm!* !*q2a car u
else resimp car u;
put('!*sq,'simpfn,'physopsimp!*sq);
% ***** end of dirty trick ******
% ----PHYSOP evaluation and simplification procedures----
symbolic procedure !*physopsm!* (u,v);
% u is the PHYSOP expression to simplify
begin scalar x,contractflg;
% if contract is set to on we keep its value at the top level
% (first call to physopsm) and set it to nil;
contractflg:=!*contract;!*contract := nil;
!*hardstop := nil;
if physopp u then
if (x:= get(u,'rvalue)) then u := physopaeval x
else if idp u then return u
else if x:=get(car u,'psimpfn) then u:= apply1(x,u)
else return physopsimp u;
u:= physopsm!* u;
if !*hardstop then <<
write " *************** WARNING: ***************";terpri();
write "Evaluation incomplete due to missing elementary relations";
terpri();
return u>>;
% the next step is to do substitutions if there are someones on
% the matching lists
if !*match or powlis1!* then <<
u := physopsubs u;
% now eval u with the substitutions
u := physopsim!* u; >>;
if not contractflg then return u
else <<
!*contract:=contractflg;
return physopcontract u >>
end;
symbolic procedure physopsim!* u;
if !*physopp!* u then physopsm!* u else u;
symbolic procedure physop2sq u; %modified 1.01
% u is a physop expr
% returns standard quotient of evaluated u
begin scalar x;
return
if physopp u then if (x:= get(u,'rvalue)) then physop2sq x
else if idp u then !*k2q u
else if (x:= get(car u,'psimpfn)) then
if physopp (x:=apply1(x,u)) then
!*k2q x
else cadr physopsm!* x
else if get(car u,'opmtch) and
(x:= opmtch!* u) then physop2sq x
else !*k2q u
else if atom u then simp u % added 1.01
else if car u eq '!*sq then cadr u
else if null getphystype u then simp u % moved from top 1.01
else physop2sq physopsm!* u
end;
symbolic procedure physopsm!* u;
% basic simplification routine
begin scalar oper,args,y,v,physopflg;
% the following is 1.02
if (null u or numberp u) then v := u
else if physopp u then v:= if (y:= get(u,'rvalue)) then physopaeval y
else if idp u then u
else if (y:=get(car u,'psimpfn)) then
apply1(y,u)
else if get(car u,'opmtch) and
(y:=opmtch!* u) then y
else u
else if atom u then v := aeval u
else <<
oper := car u;
args := cdr u;
if y:= get(oper,'physopfunction) then
% this is a function which may also have normal scalar arguments
% eg TIMES so we must check if args contain PHYSOP objects
% or if it is an already evaluated expression of physops
if flagp(oper,'physoparith) then
if hasonephysop args then v:= apply(y,list args)
else v := reval3 (oper . args)
else if flagp(oper,'physopfn) then
if areallphysops args then v:= apply(y,list args)
else
rederr2('physopsm!*,
list("invalid call of ",oper," with args: ",args))
else rederr2('physopsm!*,list(oper,
" has been flagged Physopfunction"," but is not defined"))
% this is for fns having a physop argument and no evaluation procedure
else if flagp(oper,'physopmapping) and !*physopp!* args then
v := mk!*sq !*k2q (oper . args)
% special hack for handling of PROG constructs
else if oper = 'PROG then v := physopprog args
else v := aeval u
>>;
return v
end;
symbolic procedure physopsubs u;
% general substitution routine for physop expressions
% corresponds to subs2
% u is a !*sq
% result is u in a.e. form with all substitutions of
% !*MATCH and POWLIS1!* applied
% we use a quite dirty trick here which allows to use
% the pattern matcher of standard REDUCE by hiding the
% PHYSOP rtype temporarily
begin scalar ulist,kord,alglist!*;
% step 1: convert u back to an a.e.
% u := physopaeval u; % 1.01 this line replaced
u := physop2sq u;
% step 2: transform all physops on physoplist in normal ops
for each x in physoplist!* do << remprop(x,'rtype);
put(x,'simpfn,'simpiden)>>;
% since we need it here as a prefix op
remflag('(dot),'physopfn);
put('dot,'simpfn,'simpiden);
% step 3: call simp!* on u
% u := simp!* u; % 1.01 this line replaced
u := subs2 u;
% step 4: transform u back in an a.e.
u := !*q2a u;
% step 5: transform ops in physoplist back to physops
for each x in physoplist!* do <<remprop(x,'simpfn);
put(x,'rtype,'physop)>>;
remprop('dot,'simpfn);
flag('(dot),'physopfn);
% final step return u
return u
end;
symbolic procedure physopaeval u;
% transformation of physop expression in a.e.
begin scalar x;
return
if physopp u then
if (x:=get(u,'rvalue)) then
if car x eq '!*sq then !*q2a cadr x
else x
else if atom u then u
else if (x:= get(car u,'psimpfn)) then apply1(x,u)
else if get(car u,'opmtch) and (x:= opmtch!* u) then x
else u
else if (not atom u) and car u eq '!*sq then !*q2a cadr u
else u
end;
symbolic procedure physopcontract u;
% procedure to contract over dead indices
begin scalar x,x1,w,y,z,ulist,veclist,tenslist,oldmatch,oldpowers,
alglist!*,ncmplist;
u := physopaeval u;
if physopp u then return mk!*sq physop2sq u
else if not getphystype u then return aeval u;
% now came the tricky cases
!*contract2 := t;
% step1 : collect all physops in u
ulist := collectphysops u;
veclist := for each x in ulist collect if vecopp x then x else nil;
tenslist := for each x in ulist collect if tensopp x then x else nil;
veclist:= deletemult!* deleteall(nil,veclist);
tenslist:=deletemult!* deleteall(nil,tenslist);
% step2: we now modify powlis1!* and !*match
oldmatch := !*match; !*match := nil;
oldpowers := powlis1!*; powlis1!* := nil;
% step3: transform all physops on physoplist in normal ops
for each x in physoplist!* do
<<
remprop(x,'rtype);
put(x,'simpfn,'simpiden);
if noncomp!* x then ncmplist := x . ncmplist;
>>;
% we have to declare the ops in the specoplist as noncom to avoid
% spurious simplifications during contraction
remflag('(dot opapply),'physopfn); % needed here as a normal op
flag(specoplist!*,'noncom);
for each x in specoplist!* do
<<
put(x,'simpfn,'simpiden);
put(x,'noncommutes,ncmplist)
>>;
% step4: put new matching for each vecop on the list
y := getvarindex(1);
frlis!* := nconc(frlis!*,list('!=nv));
frasc!* := nconc(frasc!*,list('nv . '!=nv));
for each x in veclist do <<
let2(list('expt,insertindices(x,transformvarindex y),'nv),
list('expt,x,'nv),nil,t);
x1:=delete(x,veclist);
for each w in x1 do
<< z := list(list((insertindices(x,y) . 1),
(insertindices(w,y) . 1)),(nil . t),
list('dot,x,w),nil);
!*match :=append(list(z),!*match) >>
>>;
% step4: put new matching for each tensop on the list
frlis!* := nconc(frlis!*,list('!=nt));
frasc!* := nconc(frasc!*,list('nt . '!=nt));
for each x in tenslist do
let2(list('expt,insertfreeindices(x,t),'nt),
list('expt,x,'nt),nil,t);
% step 6: call simp on u
u := simp!* u;
% step 7: restore previous settings
powlis1!* := oldpowers;!*match := oldmatch;
for each x in physoplist!* do
<<
remprop(x,'simpfn);
put(x,'rtype,'physop)
>>;
flag('(dot opapply),'physopfn);
remflag(specoplist!*,'noncom);
for each x in specoplist!* do
<<
remprop(x,'noncommutes);
remprop(x,'simpfn)
>>;
!*contract2 := nil;
return mk!*sq u
end;
symbolic procedure physopsimp u; % 1 line deleted 1.03
% procedure to simplify the arguments of a physop
% inspired from SIMPIDEN
begin scalar opname,w,x,y,flg;
if idp u then return u;
opname := car u;
x := for each j in cdr u collect
if idp j and (isanindex j or isavarindex j) then j %added 1.01
else physopsm!* j;
u := opname . for each j in x collect
if eqcar(j,'!*sq) then prepsqxx cadr J
else j;
if x := opmtch!* u then return x;
% special hack introduced here to check for
% symmetric and antisymmetric tensor operators
if scalopp u and tensopp opname then
<< y := get(opname,'tensdimen);
% x is the list of physopsindices
x:= for k:=1 :y collect nth(cdr u,k);
% y contains the remaining indices
if length(cdr u) > y then y := pnth(cdr u,y+1)
else y := nil;
if flagp(opname,'symmetric) then u:= opname . ordn x
else if flagp(opname,'antisymmetric) then
<< if repeats x then return 0
else if not permp(w := ordn x, x) then flg := t;
x := w;
u := opname . x >>
else u := opname . x;
if y then u:= append(u,y);
return if flg then list('minus,u) else u
>>
% special hack to introduce unrecognized IDXn indices
else if vecopp u then << if listp u then putanewindex!* cadr u;
return u >>
else if tensopp u then << if listp u then
for j:= 1 : length(cdr u) do
putanewindex!* nth(cdr u,j);
return u >>
else return u
end;
% ---- different procedures for arithmetic in phsyop expressions ----
flag('(quotient times expt difference minus plus opapply),'physoparith);
flag('(adj recip dot commute anticommute),'physopfn);
flag ('(sub),'physoparith);
flag('(sin cos tan asin acos atan sqrt int df log exp sinh cosh tanh),
'physopmapping);
% the following is needed for correct type checking 101290 mw
symbolic procedure checkphysopmap u;
% checks an expression u for unresolved physopmapping operators
begin scalar x;
a: if null u or domain!*p u or atom u or null cdr u then return nil;
x:= car u; u:= cdr u;
if listp x and flagp(car x,'physopmapping) and hasonephysop cdr x
then return t;
go to a;
end;
symbolic procedure physopfn(oper,proc);
begin
put(oper,'physopfunction,proc);
end;
physopfn('difference,'physopdiff);
symbolic procedure physopdiff args;
begin scalar lht,rht,lhtype,rhtype;
lht := physopsim!* car args;
for each v in cdr args do <<
rht := physopsim!* v;
lhtype := getphystype lht;
rhtype := getphystype rht;
if (rhtype and lhtype) and not lhtype eq rhtype then
rederr2('physopdiff,"type mismatch in diff");
lht :=
mk!*sq addsq(physop2sq lht,negsq(physop2sq rht))
>>;
return lht
end;
put('difference,'phystypefn,'getphystypeall); % changed 1.01
physopfn('minus,'physopminus);
symbolic procedure physopminus arg;
begin scalar rht,rhtype;
rht := physopsim!* car arg;
rht :=
mk!*sq negsq(physop2sq rht);
return rht
end;
put('minus,'phystypefn,'getphystypecar);
physopfn('plus,'physopplus);
symbolic procedure physopplus args;
begin scalar lht,rht,lhtype,rhtype;
lht := physopsim!* car args;
for each v in cdr args do <<
rht := physopsim!* v;
lhtype := getphystype lht;
rhtype := getphystype rht;
if (rhtype and lhtype) and not (lhtype eq rhtype) then
rederr2 ('physopplus,"type mismatch in plus ");
lht :=
mk!*sq addsq(physop2sq lht,physop2sq rht)
>>;
return lht
end;
put('plus,'phystypefn,'getphystypeall); % changed 1.01
physopfn('times,'physoptimes);
symbolic procedure physoptimes args;
begin scalar lht, rht,lhtype,rhtype,x,mul;
if (tstack!* = 0) and mul!* then << mul:= mul!*; mul!* := nil; >>;
tstack!* := tstack!* + 1;
lht := physopsim!* car args;
for each v in cdr args do <<
rht :=physopsim!* v;
lhtype := getphystype lht;
rhtype := getphystype rht;
if not lhtype then
if not rhtype then lht := mk!*sq multsq(physop2sq lht,physop2sq rht)
else if zerop lht then lht := mk!*sq (nil . 1)
else if onep lht then lht:= mk!*sq physop2sq rht
else lht:= mk!*sq multsq(physop2sq lht,physop2sq rht)
else if not rhtype then lht:=
if zerop rht then mk!*sq (nil . 1)
else if onep rht then mk!*sq physop2sq lht
else mk!*sq multsq(physop2sq rht,physop2sq lht)
else if physopordchk(physopaeval lht,physopaeval rht)
and (lhtype = rhtype) and (lhtype = 'scalar)
then lht := mk!*sq multsq(physop2sq lht,physop2sq rht)
else lht:= multopop!*(lht,rht)
>>;
b: if null mul!* or tstack!* > 1 then go to c;
lht := apply1(car mul!*,lht);
mul!* := cdr mul!*;
go to b;
c: tstack!* := tstack!* - 1;
if tstack!* = 0 then mul!* := mul;
return lht
end;
put('times,'phystypefn,'getphystypetimes);
symbolic procedure getphystypetimes args; % modified 1.01
begin scalar x;
if null (x := deleteall(nil,collectphystype args)) then
return nil
else if null cdr x then return car x
else rederr2('getphystypetimes,
list("PHYSOP type mismatch in",args))
end;
symbolic procedure multopop!*(u,v);
% u and v are physop exprs in a.e. form
% value is the product of u and v + commutators if needed
begin scalar x,y,u1,v1,stac!*,res;
% if there is no need for additional computations of commutators
% return the product as a standard quotient
u1:= physopaeval u;
v1:= physopaeval v;
if physopp u1 and physopp v1 then res := multopop(u1,v1)
else if physopp v1 then
if car u1 memq '(plus difference minus) then <<
x:= for each y in cdr u1 collect physoptimes list(y,v);
res:= reval3 (car u1 . x) >>
else if car u1 eq 'times then <<
stac!*:= reverse cdr u1; % begin with the last el
y:= v;
while stac!* do <<
x := car stac!*;
y := physoptimes list(x,y);
stac!* := cdr stac!*;
>>;
res:= y >>
else if car u1 eq 'quotient then res:= mk!*sq
quotsq(physop2sq physoptimes list(cadr u1,v),
physop2sq caddr u1)
else res:= physoptimes list(u1,v1)
else if car v1 memq '(plus difference minus) then <<
x:= for each y in cdr v1 collect physoptimes list(u,y);
res:= reval3 (car v1 . x) >>
else if car v1 eq 'times then <<
stac!*:= cdr v1;
y:= u;
while stac!* do <<
x := car stac!*;
y := physoptimes list(y,x);
stac!* := cdr stac!*;
% write "y= ",y," stac= ",stac!*;terpri();
>>;
res:= y >>
else if car v1 eq 'quotient then res:= mk!*sq
quotsq(physop2sq physoptimes list(u,cadr v1),
physop2sq caddr v1)
else res:= physoptimes list(u1,v1);
return res
end;
symbolic procedure multopop(u,v);
% u and v are physops (kernels)
% value is the product of physops + commutators if necessary
begin scalar res,x,ltype,rtype;
ltype := getphystype u;
rtype := getphystype v;
if ltype neq rtype then
rederr2('multopop,"type conflict in TIMES")
else if (invp u = v) then res := mk!*sq !*k2q 'unit
else if u = 'unit then res := mk!*sq !*k2q v
else if v = 'unit then res := mk!*sq !*k2q u
else if ordop(u,v) then
res := mk!*sq !*f2q multfnc(!*k2f u,!*k2f v)
else if noncommuting(u,v) then <<x:= comm2(u,v);
res:= if !*anticommchk then physopplus
list(list('minus,list('times,v,u)),x)
else physopplus
list(list('times,v,u),x) >>
else res := mk!*sq !*f2q multfnc(!*k2f v,!*k2f u);
return res
end;
physopfn('expt,'physopexpt);
symbolic procedure physopexpt args;
begin scalar n1,n2,lht,rht,lhtype,rhtype,x,y,z;
% we have to add a special bootstrap to avoid too much simplification
% in case of dot products raise to a power
lht := physopsm!* car args;
rht := physopsm!* cadr args;
lhtype := physopp lht ;
rhtype := physopp rht;
if rhtype then
rederr2('physopexpt,"operators in the exponent cannot be handled");
if not getphystype lht then lht := reval3 list('expt,lht,rht);
if not lhtype then
if numberp rht then <<
n1 := car divide(rht,2);
n2 := cdr divide(rht,2);
lhtype := getphystype lht;
if (lhtype and zerop rht) then lht := mk!*sq !*k2q 'unit %1.01
else if lhtype = 'vector then <<
x:= for k:= 1 : n1 collect physopdot list(lht,lht);
if onep n1 then x := 1 . x;
lht:= if zerop n2 then physoptimes x
else physoptimes append(x,list(lht));>>
else if lhtype = 'tensor then <<
x:= for k:= 1 : n1 collect physoptens list(lht,lht);
if onep n1 then x := 1 . x;
lht:= if zerop n2 then physoptimes x
else physoptimes append(x,list(lht));>>
else if lhtype = 'state then
rederr2('physopexpt,
"expressions involving states cannot be exponentiated")
else << lht := physopaeval lht;
x := deletemult!* collectindices lht;
z := lht;
for k :=2 :rht do <<
for each x1 in x do
if isavarindex x1 then
lht:= subst(makeanewvarindex(),x1,lht)
else lht:=subst(makeanewindex(),x1,lht);
y := append(y,list(lht));
lht := z; >>;
lht := physoptimes (z . y); >>;
>>
else lht := mk!*sq simpx1(physopaeval lht,physopaeval rht,1)
else if lht = 'unit then lht := mk!*sq !*k2q 'unit
else if numberp rht then lht := exptexpand(lht,rht)
else lht := mk!*sq !*P2q (lht . physopaeval rht); %0.99c
return lht
end;
put('expt,'phystypefn,'getphystypeexpt);
symbolic procedure getphystypeexpt args; % recoeded 0.99c
begin scalar x;
x := getphystypecar args;
return
if null x then nil
else if numberp cadr args and evenp cadr args then 'scalar
else x;
end;
symbolic procedure exptexpand(u,n);
begin scalar bool,x,y,v,n1,n2,res,flg;
if not numberp n then
rederr2('exptexpand,list("invalid argument ",n," to EXPT"));
if zerop n then return mk!*sq !*k2q 'unit; %1.01
bool := if n < 0 then t else nil;
n := if bool then abs(n) else n;
n1 := car divide(n,2);
n2 := cdr divide(n,2);
if zerop n1 then return mk!*sq !*k2q
if bool then invp u else u;
res := (1 . 1);
for k := 1 : n1 do <<
if scalopp u then
if bool then x := multf(!*k2f invp u, !*k2f invp u) . 1
else x := multf(!*k2f u, !*k2f u) . 1
% if bool then x:= list(list((invp u . 1),((invp u . 1) . 1))) . 1
% else x:= list(list((u . 1),((u . 1) . 1))) . 1
else if vecopp u then
if bool then x:= quotsq((1 . 1),physop2sq physopdot list(u,u))
else x:= physop2sq physopdot list(u,u)
else if tensopp u then <<
if bool then x:= quotsq((1 . 1),
physop2sq physoptens list(u,u))
else x:= physop2sq physoptens list(u,u) >>
else rederr2('exptexpand, "cannot raise a state to a power");
res := multsq(res,x)
>>;
b:
if zerop n2 then return mk!*sq res;
u:= if bool then invp u else u;
return mk!*sq multsq(res,!*k2q u)
end;
physopfn('quotient,'physopquotient);
symbolic procedure physopquotient args;
begin scalar lht, rht,y,lhtype,rhtype;
lht := physopsim!* car args;
rht := physopsim!* cadr args;
lhtype := getphystype car args;
rhtype := getphystype cadr args;
if rhtype memq '(vector state tensor) then
rederr2('physopquotient, "invalid quotient")
else if not rhtype then return
mk!*sq quotsq(physop2sq lht,physop2sq rht);
lhtype := physopp lht;
rht := physopaeval rht;
rhtype := physopp rht;
if rhtype then
if not lhtype then lht:= mk!*sq multsq(physop2sq lht,!*k2q invp rht)
else lht:= physoptimes list(lht,invp rht)
else if car rht eq 'times and null deadindices rht then
<< rht := reverse cdr rht;
rht := for each x in rht collect
physopquotient list(1,x);
lht := physoptimes append(list(lht),rht) >>
else lht:= mk!*sq quotsq(physop2sq lht,physop2sq rht);
return lht
end;
put('quotient,'phystypefn,'getphystypeor);
physopfn('recip,'physoprecip);
symbolic procedure physoprecip args;
physopquotient list(1,args);
put('recip,'phystypefn,'getphystypecar);
symbolic procedure inv u;
% inverse of physops
begin scalar x,y;
if not physopp u then rederr2('inv, "invalid argument to INVERSE");
if u = 'unit then return u;
y:= if idp u then u else car u;
x := reversip explode y;
x := intern compress nconc(reversip x,list('!!,'!-,'!1));
put(y,'inverse,x); % 1.01
put(x,'inverse,y); % 1.01
put(x,'physopname,y); % 1.02
if not physopp x then << put(x,'rtype,'physop);
put(x,'phystype,get(y,'phystype));
put(x,'psimpfn,'physopsimp);
put(x,'tensdimen,get(y,'tensdimen));
physoplist!* := nconc(physoplist!*,list(x));
>>;
if idp u then return x
else return nconc(list(x),cdr u)
end;
symbolic procedure invp u; % recoded 1.01
% special cases
if u = 'unit then u
else if atom u then get(u,'inverse)
else if member(car u,'(comm anticomm)) then
list('quotient,1,u)
else get(car u,'inverse) . cdr u;
physopfn('sub,'physopsub); %subcommand;
% ********* redefinition of SUB handling is necessary in 3.4 **********
remprop('sub,'physopfunction);
put('sub,'physopfunction,'subeval);
put('physop,'subfn,'physopsub);
symbolic procedure physopsub(u,v); %redefined
% u is a list of substitutions as an a--list
% v is a simplified physop in prefix form
begin scalar res;
if null u or null v then return v;
v := physopaeval v;
for each x in u do v := subst(cdr x,car x,v);
return physopsm!* V
end;
% *********** end of 3.4 modifications ******************
symbolic procedure physopprog u;
% procedure to handle prog expressions (i.e. loops) containing physops
begin scalar x;
% we use basically the same trick as in physopsubs
% step 1: transform all physops on physoplist in normal ops
for each x in physoplist!* do <<remprop(x,'rtype);
put(x,'simpfn,'simpiden)>>;
% step 2: call normal prog on u
u := aeval ('prog . u);
% step 3: transform u back in an a.e.
u := physopaeval u;
% step 4: transform ops in physoplist back to physops
for each x in physoplist!* do <<remprop(x,'simpfn);
put(x,'rtype,'physop)>>;
% final step return u
return physopsm!* u
end;
% ****** procedures for physopfns ***********
physopfn('dot,'physopdot);
infix dot;
precedence dot,*;
symbolic procedure physopdot args;
begin scalar lht,rht,lhtype,rhtype,x,n,res;
lht := physopaeval physopsim!* car args;
rht := physopaeval physopsim!* cadr args;
lhtype := getphystype lht;
rhtype := getphystype rht;
if not( (lhtype and rhtype) and (lhtype eq 'vector) ) then
rederr2 ('physopdot,"invalid arguments to dotproduct");
lhtype := physopp lht;
rhtype := physopp rht;
if rhtype then
if lhtype then << if !*indflg then<<
lht := insertfreeindices(lht,nil);
rht := insertfreeindices(rht,nil);
indcnt!* := indcnt!* + 1; >>
else <<x := makeanewindex();
lht := insertindices(lht,x);
rht := insertindices(rht,x);>>;
res := physoptimes list(lht,rht)>>
else <<
if car lht eq 'minus then
res := mk!*sq negsq(physop2sq physopdot list(cadr lht,rht))
else if car lht eq 'difference then res := mk!*sq addsq(
physop2sq physopdot list(cadr lht,rht),negsq(physop2sq
physopdot list(caddr lht,rht)))
else if car lht eq 'plus then <<
x := for each y in cdr lht collect physopdot list(y,rht);
res := reval3 append(list('plus),x) >>
else if car lht eq 'quotient then <<
if not vecopp cadr lht then
rederr2('physopdot,"argument to DOT")
else res := mk!*sq quotsq(physop2sq
physopdot list(cadr lht,rht),physop2sq caddr lht) >>
else if car lht eq 'times then <<
for each y in cdr lht do
if getphystype y eq 'vector then x:=y;
lht :=delete(x,cdr lht);
res := physoptimes
nconc(lht,list(physopdot list(x,rht))) >>
else rederr2('physopdot, "invalid arguments to DOT") >>;
if not rhtype then <<
if car rht eq 'minus then
res := mk!*sq negsq(physop2sq physopdot list(lht,cadr rht))
else if car rht eq 'difference then res := mk!*sq addsq(
physop2sq physopdot list(lht,cadr rht),negsq(physop2sq
physopdot list(lht, caddr rht)))
else if car rht eq 'plus then <<
x := for each y in cdr rht collect physopdot list(lht,y);
res := reval3 append(list('plus),x) >>
else if car rht eq 'quotient then <<
if not vecopp cadr rht then
rederr2 ('physopdot,"invalid argument to DOT")
else res := mk!*sq quotsq(physop2sq physopdot
list(lht,cadr rht),physop2sq caddr rht) >>
else if car rht eq 'times then <<
for each y in cdr rht do if getphystype y eq 'vector then x:=y;
rht :=delete(x,cdr rht);
res := physoptimes
nconc(rht,list(physopdot list(lht,x))) >>
else rederr2 ('physopdot,"invalid arguments to DOT") >>;
return res
end;
put('dot,'phystype,'scalar);
symbolic procedure physoptens args;
% procedure for products of tensor expressions
begin scalar lht,rht,lhtype,rhtype,x,n,res;
lht := physopaeval physopsim!* car args;
rht := physopaeval physopsim!* cadr args;
lhtype := getphystype lht;
rhtype := getphystype rht;
if not( (lhtype and rhtype) and (lhtype eq 'tensor) ) then
rederr2 ('physoptens,"invalid arguments to tensproduct");
lhtype := physopp lht;
rhtype := physopp rht;
if rhtype then
if lhtype then << n:= get(lht,'tensdimen);
if (n neq get(rht,'tensdimen)) then
rederr2('physoptens,
"tensors must have the same dimension to be multiplied");
if !*indflg then<<
lht := insertfreeindices(lht,nil);
rht := insertfreeindices(rht,nil);
indcnt!* := indcnt!* + n; >>
else <<x := for k:= 1 : n collect makeanewindex();
lht := insertindices(lht,x);
rht := insertindices(rht,x);>>;
res := physoptimes list(lht,rht)>>
else <<
if car lht eq 'minus then
res := mk!*sq negsq(physop2sq physoptens list(cadr lht,rht))
else if car lht eq 'difference then res := mk!*sq addsq(
physop2sq physoptens list(cadr lht,rht),negsq(physop2sq
physoptens list(caddr lht,rht)))
else if car lht eq 'plus then <<
x := for each y in cdr lht collect physoptens list(y,rht);
res := reval3 append(list('plus),x) >>
else if car lht eq 'quotient then <<
if not tensopp cadr lht then
rederr2 ('physoptens,"invalid argument to TENS")
else res := mk!*sq quotsq(physop2sq
physoptens list(cadr lht,rht),physop2sq caddr lht) >>
else if car lht eq 'times then <<
for each y in cdr lht do
if getphystype y eq 'tensor then x:=y;
lht :=delete(x,cdr lht);
res := physoptimes
nconc(lht,list(physoptens list(x,rht))) >>
else rederr2('physoptens, "invalid arguments to TENS") >>;
if not rhtype then <<
if car rht eq 'minus then
res := mk!*sq negsq(physop2sq physoptens list(lht,cadr rht))
else if car rht eq 'difference then res := mk!*sq addsq(
physop2sq physoptens list(lht,cadr rht),negsq(physop2sq
physoptens list(lht, caddr rht)))
else if car rht eq 'plus then <<
x := for each y in cdr rht collect physoptens list(lht,y);
res := reval3 append(list('plus),x) >>
else if car rht eq 'quotient then <<
if not tensopp cadr rht then
rederr2 ('physoptens,"invalid argument to TENS")
else res := mk!*sq quotsq(physop2sq physoptens
list(lht,cadr rht),physop2sq caddr rht) >>
else if car rht eq 'times then <<
for each y in cdr rht do if getphystype y eq 'tensor then x:=y;
rht :=delete(x,cdr rht);
res := physoptimes
nconc(rht,list(physoptens list(lht,x))) >>
else rederr2('physoptens, "invalid arguments to TENS") >>;
return res
end;
put('tens,'phystype,'scalar);
% -------- procedures for commutator handling -------------
symbolic procedure comm2(u,v);
% general procedure for getting commutators
begin scalar x,utype,vtype,y,z,z1,res;
if not (physopp u and physopp v) then rederr2('comm2,
"invalid arguments to COMM");
utype := getphystype u;
vtype := getphystype v;
if not (utype eq 'scalar) and (vtype eq 'scalar) then
rederr2('comm2, "comm2 can only handle scalar operators");
!*anticommchk:= nil;
if not noncommuting(u,v) then return
if !*anticom then mk!*sq !*f2q multf(!*n2f 2,multfnc(!*k2f v,!*k2f u))
else mk!*sq (nil . 1);
x := list(u,v);
z := opmtch!* ('comm . x);
if null z then z:= if (y:= opmtch!* ('comm . reverse x)) then
physopsim!* list('minus,y)
else nil;
if z and null !*anticom then res:= physopsim!* z
else << z1 := opmtch!* ('anticomm . x);
if null z1 then
z1 := if (y:=opmtch!* ('anticomm . reverse x)) then y
else nil;
if z1 then << !*anticommchk := T;
res:= physopsim!* z1>>
>>;
if null res then
<< !*hardstop:= T;
if null !*anticom then res := mk!*sq !*k2q ('comm . x)
else << !*anticommchk := T;
res := mk!*sq !*k2q ('anticomm . x) >>
>>;
return res
end;
physopfn('commute,'physopcommute);
symbolic procedure physopcommute args;
begin scalar lht,rht,lhtype,rhtype,x,n,res,flg;
lht := physopaeval physopsim!* car args;
rht := physopaeval physopsim!* cadr args;
lhtype := getphystype lht;
rhtype := getphystype rht;
if not (lhtype and rhtype) then return mk!*sq !*d2q 0
else if not(rhtype = lhtype) then
rederr2('physopcommute,
"physops of different types cannot be commuted")
else if not lhtype eq 'scalar then
rederr2 ('physopcommute,
"commutators only implemented for scalar physop expressions");
% flg := !*anticom; !*anticom := nil;
lhtype := physopp lht;
rhtype := physopp rht;
% write "lht= ",lht," rht= ",rht;terpri();
if rhtype then
if lhtype then << res := comm2(lht,rht);
if !*anticommchk then
res := physopdiff list(res,
physoptimes list(2,rht,lht)); >>
else res := mk!*sq negsq(physop2sq physopcommute list(rht,lht))
else <<
if car rht eq 'minus then res:= mk!*sq negsq(physop2sq
physopcommute list(lht, cadr rht));
if car rht eq 'difference then res := mk!*sq addsq(
physop2sq physopcommute list(lht,cadr rht),negsq(physop2sq
physopcommute list(lht,caddr rht)));
if car rht eq 'plus then <<
x:= for each y in cdr rht collect
physopcommute list(lht,y);
res:= reval3 append(list('plus),x) >>;
if car rht memq '(expt dot commute) then
res := physopcommute list(lht,physopsim!* rht);
if car rht eq 'quotient then
if physopp caddr rht then
res:= physopcommute list(lht,physopsim!* rht)
else
res := mk!*sq quotsq(physop2sq physopcommute list(lht,cadr rht),
physop2sq caddr rht);
if car rht eq 'times then <<
n := length cdr rht;
if (n = 2) then res := reval3 list('plus, physopsim!*
list('times,cadr rht,physopcommute list(lht, caddr rht)),
physopsim!* list('times,physopcommute list(lht, cadr rht),caddr rht))
else res := reval3 list('plus, physopsim!*
list('times,cadr rht,physopcommute list(lht,
append('(times),cddr rht))), physopsim!* append(
list('times,physopcommute list(lht, cadr rht)), cddr rht)) >>
>>;
% !*anticom := flg;
return res
end;
put('commute,'phystype,'scalar);
physopfn('anticommute,'physopanticommute);
symbolic procedure physopanticommute args;
begin scalar lht,rht,lhtype,rhtype,x,n,res,flg;
lht := physopaeval physopsim!* car args;
rht := physopaeval physopsim!* cadr args;
lhtype := getphystype lht;
rhtype := getphystype rht;
if not (lhtype and rhtype) then
return mk!*sq aeval list('plus,list('times,lht,rht),
list('times,rht,lht))
else if not(rhtype = lhtype) then
rederr2('physopanticommute,
"physops of different types cannot be commuted")
else if not lhtype eq 'scalar then
rederr2 ('physopanticommute,
"commutators only implemented for scalar physop expressions");
% flg := !*anticom;!*anticom :=t;
lhtype := physopp lht;
rhtype := physopp rht;
% write "lht= ",lht," rht= ",rht;terpri();
if rhtype then
if lhtype then
<<
x := comm2(lht,rht);
if null !*anticommchk then
If !*hardstop then res := mk!*sq !*k2q list('anticomm,lht,rht)
else res := reval3 list('plus,x,physoptimes list(2,rht,lht))
else res := x;
>>
else res := physopsim!* physopanticommute list(rht,lht)
else <<
if car rht eq 'minus then res:= mk!*sq negsq(physop2sq
physopanticommute list(lht, cadr rht));
if car rht eq 'difference then mk!*sq addsq(physop2sq
physopanticommute list(lht,cadr rht),negsq(physop2sq
physopanticommute list(lht,caddr rht)));
if car rht eq 'plus then <<
x:= for each y in cdr rht collect
physopanticommute list(lht,y);
res:= reval3 append(list('plus),x) >>
else res := physopplus list(physoptimes list(lht,rht),
physoptimes list(rht,lht));
>>;
% !*anticom := flg;
return res
end;
put('anticommute,'phystype,'scalar);
symbolic procedure commsimp u;
% procedure to simplify the arguments of COMM or ANTICOMM
% if they are not simple physops
begin scalar opname,x,y,flg,res;
opname := car u;
x := physopsim!* cadr u;
y := physopsim!* caddr u;
% write "op= ",opname," x= ",x," y= ",y;terpri();
flg := !*anticom;
if opname = 'anticomm then !*anticom := t;
res := if physopp x and physopp y then physopaeval comm2(x,y)
else if opname eq 'comm then list('commute,physopaeval x,
physopaeval y)
else list('anticommute,physopaeval x,physopaeval y);
!*anticom := flg;
return res
end;
% -------------- application of ops on states ----------------
physopfn('opapply,'physopapply);
infix opapply;
precedence opapply,-;
symbolic procedure physopapply args; % changed 0.99b
begin scalar lhtype,rhtype,wave,op,wavefct,res,x,y,flg;
lhtype := statep!* car args;
rhtype := statep!* cadr args;
if rhtype and lhtype then
return statemult(car args,cadr args)
else if rhtype then
<<wave := physopaeval physopsim!* cadr args;
op := physopaeval physopsim!* car args >>
else if lhtype then
<<wave := physopaeval physopadj list(car args);
op := physopaeval physopadj cdr args >>
% a previous application of physopapply may have annihilated the
% state
else if zerop car args or zerop cadr args then return mk!*sq (nil . 1)
else rederr2('opapply, "invalid arguments to opapply");
if null getphystype op then
res:= mk!*sq multsq(physop2sq op,physop2sq wave)
else if not physopp op then
if car op eq 'minus then
res := mk!*sq negsq(physop2sq physopapply list(cadr op,wave))
else if car op memq '(plus difference) then <<
for each y in cdr op do <<
res:= nconc(res,list(physopapply list(y,wave)));
if !*hardstop then flg:= t;
!*hardstop := nil;>>;
if flg then !*hardstop := t;
res := reval3 ((car op) . res) >>
else if car op memq '(dot commute anticommute expt) then
res := physopapply list(physopsim!* op,wave)
else if car op eq 'quotient then
if physopp caddr op then
res := physopapply list(physopsim!* op,wave)
else res := mk!*sq quotsq(physop2sq
physopapply list(cadr op,wave),physop2sq caddr op)
else if car op eq 'times then
<<op := reverse cdr op;
while op and not !*hardstop do
<< x := car op; op := cdr op;
wave := physopapply list(x,wave) >>;
if !*hardstop then
if null op then res := wave
else << x:= physopaeval wave;
op := 'times . reverse op;
while x do
<< y := car x; x := cdr x;
if listp y and
(y := assoc('opapply,y)) then
<< wavefct := list('opapply,
nconc(op,list(cadr y)),
caddr y);
wave := subst(wavefct,y,wave);
>>;
>>;
res := wave;
>>
else res := wave;
>>
else rederr2('opapply, "invalid operator to opapply")
% special hack here for unit operator 0.99c
else if op = 'unit then res := mk!*sq physop2sq wave
else if physopp wave or
(flagp(car wave,'physopmapping) and statep!* cdr wave) then
<<x := opmtch!* list('opapply,op,wave);
if null x then x := physopadj list(
opmtch!* list('opapply,adjp wave,adjp op));
if null x then <<!*hardstop := t;
res := mk!*sq !*k2q
list('opapply,op,wave); >>
else res := mk!*sq physop2sq x;
>>
else << x := wave; wave := nil;
while x do <<
wavefct := car x; x := cdr x;
if statep!* wavefct then wave := nconc(wave,
list(physopaeval physopapply list(op,wavefct)))
else wave := nconc(wave,list(wavefct));
if !*hardstop then flg := t;
!*hardstop := nil >>;
if flg then !*hardstop := t;
res := mk!*sq physop2sq wave;
>>;
return res
end;
put('opapply,'phystypefn,'getphystypestate);
symbolic procedure getphystypestate args;
if statep!* car args and statep!* cadr args then nil
else 'state;
symbolic procedure statemult(u,v); % recoded 0.99c
% u and v are states
% returns product of these
begin scalar x,y,res,flg;
if not (statep!* u or statep!* v) then
rederr2 ('statemult,"invalid args to statemult");
if (not atom u and car v eq 'opapply) then
return expectval(u,cadr v,caddr v);
if (not atom u and car u eq 'opapply) then
return expectval(cadr u,caddr u,v);
u := physopaeval physopsim!* u;
v := physopaeval physopsim!* v;
if physopp u then
if physopp v then
<<
x := opmtch!* list('opapply,u,v);
if x then res := physop2sq aeval x
else
<<
x:= opmtch!* list('opapply,v,u);
if null x then
<<
!*hardstop := t;
res:= !*k2q list('opapply,u,v)
>>
else res := physop2sq aeval compconj x
>>;
>>
else
<<
x := deletemult!* !*collectphysops v;
for each y in x do
<<
v := subst(physopaeval statemult(u,y),y,v);
if !*hardstop then flg := t;
!*hardstop := nil;
>>;
if flg then !*hardstop := t;
res := physop2sq v;
>>
else
<<
x := deletemult!* !*collectphysops u;
for each y in x do
<<
u := subst(physopaeval statemult(y,v),y,u);
if !*hardstop then flg := t;
!*hardstop := nil;
>>;
if flg then !*hardstop := t;
res := physop2sq u;
>>;
return mk!*sq res
end;
symbolic procedure expectval(u,op,v);
% u and v are states
% calculates the expectation value < u ! op ! v >
% tries to apply op first on v, then on u
% PHYSOPAPPLY is used rather than STATEMULT to multiply
% resulting states together because of more general definition
begin scalar x,y,z,flg,res;
op := physopaeval physopsim!* op;
if null getphystype op then
return mk!*sq multsq(physop2sq op,physop2sq physopapply list(u,v));
if physopp op then
<<x := physopapply list(op,v);
if !*hardstop then
<< !*hardstop := nil;
x:= physopapply list(u,op);
res := if !*hardstop then mk!*sq
!*k2q list('opapply,list('opapply,u,op),v)
else physopapply list(x,v) >>
else res:= physopapply list(u,x)
>>
else if car op eq 'minus then
res := mk!*sq negsq(physop2sq expectval(u,cadr op,v))
else if car op eq 'quotient then
if physopp caddr op then res := expectval(u,physopsm!* op,v)
else res := mk!*sq quotsq(physop2sq expectval(u,cadr op,v),
physop2sq caddr op)
else if car op memq '(dot commute anticommute expt) then
res := expectval (u,physopsm!* op,v)
else if car op memq '(plus difference) then
<< for each y in cdr op do
<< x:=nconc(x,list(expectval(u,y,v)));
if !*hardstop then flg:= !*hardstop ;
!*hardstop := nil >>;
if flg then !*hardstop := t;
res := reval3 ((car op) . x);
>>
else if car op eq 'times then
<< x := physopapply list(op,v);
if not !*hardstop then return physopapply list(u,x);
x := cdr op;
while (x and !*hardstop and not flg) do
<< y:=car x; x := cdr x;
if not getphystype y then << v:= physopapply list(y,v);
y := v;>>
else << !*hardstop := nil;
z:= physopapply list(u,y);
if !*hardstop then
<< flg := T; x := y . x;
y := if null cdr x then list('opapply,car x,
physopaeval v)
else list('opapply,('times . x),
physopaeval v); >>
else << u:= z;
y:= if null x then v
else if null cdr x then
physopapply list(car x, physopaeval v)
else physopapply
list(('times . x),physopaeval v) >>
>>
>>;
res := if !*hardstop then mk!*sq !*k2q list('opapply,
physopaeval u,physopaeval y)
else physopapply list(u,y);
>>
else rederr2('expectval, "invalid args to expectval");
return res
end;
symbolic procedure compconj u;
% dirty and trivial implementation of
% complex conjugation of everything (hopefully);
% not yet tested for arrays
begin scalar x;
if null u or numberp u then return u
else if idp u and (x:=get(u,'rvalue)) then <<
x:=subst(list('minus,'I),'I,x);
put(u,'rvalue,x); return u >>
else return subst(list('minus,'I),'I,u)
end;
% -------------- adjoint of operators ---------------------
physopfn('adj, 'physopadj);
symbolic procedure physopadj arg;
begin scalar rht,rhtype,x,n,res;
rht := physopaeval physopsim!* car arg;
rhtype := physopp rht;
if rhtype then return mk!*sq !*k2q physopsm!* adjp rht
else <<
if not getphystype rht then res := aeval compconj rht
else if car rht eq 'minus then
res := mk!*sq negsq(physop2sq physopadj list(cadr rht))
else if car rht eq 'difference then res := mk!*sq addsq(
physop2sq physopadj list(cadr rht),negsq(physop2sq
physopadj list(caddr rht)))
else if car rht eq 'plus then <<
x := for each y in cdr rht collect physopadj list(y);
res := reval3 ('plus . x) >>
else if car rht eq 'quotient then <<
if not getphystype cadr rht then
rederr2('physopadj, "invalid argument to ADJ")
else res := mk!*sq quotsq(physop2sq
physopadj list(cadr rht),physop2sq caddr rht) >>
else if car rht eq 'times then <<
x:= for each y in cdr rht collect physopadj list(y);
res := physoptimes reverse x >>
else if flagp(car rht,'physopmapping) then
res := mk!*sq !*k2q list(car rht, physopaeval physopadj cdr rht)
else res :=physopadj list(physopsim!* rht) >>;
return res
end;
Put('adj,'phystypefn,'getphystypecar);
symbolic procedure adj2 u;
begin scalar x,y;
if not physopp u then rederr2('adj2, "invalid argument to adj2");
if u = 'unit then return u;
y:= if idp u then u else car u;
x := reverse explode y;
x := intern compress nconc(reverse x,list('!!,'!+));
put(y,'adjoint,x); %1.01
put(x,'adjoint,y); %1.01
put(x,'physopname,x); % 1.02
if not physopp x then << put(x,'rtype,'physop);
put(x,'phystype,get(y,'phystype));
put(x,'psimpfn,'physopsimp);
put(x,'tensdimen,get(y,'tensdimen));
defoporder!* := nconc(defoporder!*,list(x));
oporder!* := nconc(oporder!*,list(x));
physoplist!* := nconc(physoplist!*,list(x));
>>;
if idp u then return x
else return x . cdr u
end;
symbolic procedure invadj u; %new 1.01
% create the inverse adjoint op
begin scalar x,y;
if not physopp u then rederr2('invadj, "invalid argument to invadj");
if u = 'unit then return u;
y:= if idp u then u else car u;
x := reverse explode y;
x := intern compress nconc(reverse x,list('!!,'!+,'!!,'!-,'!1));
put(x,'adjoint,get(y,'inverse));
put(x,'inverse,get(y,'adjoint));
put(get(y,'inverse),'adjoint,x);
put(get(y,'adjoint),'inverse,x);
put(x,'physopname,get(y,'adjoint)); % 1.02
if not physopp x then << put(x,'rtype,'physop);
put(x,'phystype,get(y,'phystype));
put(x,'psimpfn,'physopsimp);
put(x,'tensdimen,get(y,'tensdimen));
physoplist!* := nconc(physoplist!*,list(x));
>>;
if idp u then return x
else return x . cdr u
end;
symbolic procedure adjp u; %recoded 1.01
% special cases
if u = 'unit then u
else if atom u then get(u,'adjoint)
else if (car u = 'comm) then
list('comm,adjp caddr u,adjp cadr u)
else if (car u = 'anticomm) then
list('anticomm,adjp cadr u,adjp caddr u)
else get(car u,'adjoint) . cdr u;
% --- end of arithmetic routines ---------------------
% ---- procedure for handling let assignements ------
symbolic procedure physoptypelet(u,v,ltype,b,rtype);
% modified version of original typelet
% General function for setting up rules for PHYSOP expressions.
% LTYPE is the type of the left hand side U, RTYPE, that of RHS V.
% B is a flag that is true if this is an update, nil for a removal.
% updated 101290 mw
%do not check physop type in prog exprs on the rhs
begin scalar x,y,n,u1,v1,z,contract;
if not physopp u and getphystype u then goto c; % physop expr
u1 := if atom u then u else car u;
if ltype then
if rtype = ltype then go to a
ELSE IF NULL B OR ZEROP V OR (LISTP V AND ((CAR V = 'PROG)
OR (CAR V = 'COND))) %1.01
or ((not atom u) and (car u = 'opapply)) then
return physopset(u,v,b)
else rederr2('physoptypelet,
list("physop type mismatch in assignement ",
u," := ",v))
else if null (x:= getphystype v) then return physopset(u,v,b)
else << if x = 'scalar then scalop u1;
if x = 'vector then vecop u1;
if x = 'state then state u1;
if x = 'tensor then tensop list(u1,get(v,'tensdimen));
ltype := rtype >>;
A: if b and (not atom u or flagp(u,'used!*)) then rmsubs();
% perform the assignement
physopset(u,v,b);
% phystype checking added 1.01
if b and (getphystype u neq getphystype v) then
rederr2('physoptypelet,
list("physop type mismatch in assignement ",
u," <=> ",v));
% special hack for commutators here
if (not atom u) and (car u = 'comm) then
physopset(list('comm,adjp caddr u,adjp cadr u),list('adj,v),b);
if (not atom u) and (car u = 'anticomm) then
physopset(list(car u,adjp cadr u,adjp caddr u),list('adj,v),b);
if null (x := getphystype u) or (x = 'state) or (x = 'scalar)
then return;
% we have here to add additional scalar let rules for vector
% and tensor operators with arbitrary indices
u1:=u;v1:=v;
if (x eq 'vector) or (x eq 'tensor) then
<< x := collectphysops u;
for each z in x do
u1:= subst(insertfreeindices(z,nil),z,u1);
x := collectphysops v;
for each z in x do
v1:= subst(insertfreeindices(z,nil),z,v1) >>;
physoptypelet(u1,v1,ltype,b,rtype);
return;
C:
% this is for more complicated let rules involving more than
% one term on the lhs
% special hack here to handle let rules involving elementary
% OPAPPLY relations
if car u = 'opapply then return physopset(u,v,b);
% step 1: do all physop simplifications on lhs
% we set indflg!* for dot product simplifications on the lhs
!*indflg:= T; indcnt!* := 0;
contract := !*contract2; !*contract2 := T;
u := physopsm!* u;
!*indflg := nil; indcnt!* := 0; !*contract2 := contract;
% check correct phystype
x := getphystype u;
y := getphystype v;
if b and ((not (y or zerop v)) or (y and (x neq y))) then
rederr2 ('physoptypelet,"phystype mismatch in LET");
% step 2 : transform back in ae
u := physopaeval u;
% write "u= ",u; terpri();
% ab hier neu
% step3 : do some modifications in case of a sum or difference on the lh
if car u = 'PLUS then
<<
u1 := cddr u;
u := cadr u;
v := list('plus,v);
for each x in u1 do
<<
x := list('minus,x);
v := append(v,list(x));
>>;
>>;
if car u = 'DIFFERENCE then
<<
u1:= cddr u;
u:= cadr u;
v := append(list('plus,v),list(u1));
>>;
if car U = 'MINUS then
<<
u := cadr u;
v := list('minus,v);
>>;
% step 4: add the rule to the corresponding list
% expression may still contain quotients and expt
if car u ='EXPT then
<<
u := cadr u . caddr u;
powlis1!* := xadd!*(u .
list(nil . (if mcond!* then mcond!* else t),
v,nil), powlis1!*,b)
>>
else if car u = 'quotient then
<<
v:= list('times,v,caddr u);
physoptypelet(cadr u,v,ltype,b,rtype);
>>
else % car u = times
<<
u1 := nil;
for each x in cdr u do
<<
if car x= 'expt then
u1 := append(u1,list(cadr x . caddr x))
else if car x = 'quotient then
<<
v:= list('times,v,caddr x);
u1 := append(u1,
list(if cadr x = 'expt then
(caddr x . cadddr x)
else (cadr x . 1)));
>>
else u1 := append(u1,list(x . 1));
>>;
!*match := xadd!*(u1 . list(nil .
(if mcond!* then mcond!* else t),
v,nil), !*match,b);
>>;
return;
end;
symbolic procedure physopset(u,v,b);
% assignement procedure for physops
% special hack for assignement of unresolved physop expressions
begin
if not atom u then put(car u,'opmtch,xadd!*(cdr u .
list(nil . (if mcond!* then mcond!* else t),
v,nil), get(car u,'opmtch),b))
else if b then
if physopp u then put(u,'rvalue,physopsim!* v)
else put(u,'avalue,list('scalar,
list('!*sq,cadr physopsim!* v,not !*hardstop)))
else if not member(u,specoplist!*) then
<< remprop(u,'rvalue);
remprop(u,'opmtch);
>>;
!*hardstop := nil;
end;
symbolic procedure clearphysop u;
% to remove physop type from an id
begin scalar y;
for each x in u do <<
if not (physopp x and idp x) then rederr2('clearphysop,
list("invalid argument ",x," to CLEARPHYSOP"));
y := invp x;
remprop(y,'rtype);
remprop(y,'tensdimen);
remprop(y,'phystype);
remprop(y,'psimpfn);
remprop(y,'inverse); %1.01
remprop(y,'adjoint); %1.01
remprop(y,'rvalue); % 1.01
oporder!* := delete(y,oporder!*);
defoporder!* := delete(y,defoporder!*);
physoplist!* := delete(y,physoplist!*);
y:= adjp x;
remprop(y,'rtype);
remprop(y,'tensdimen);
remprop(y,'phystype);
remprop(y,'psimpfn);
remprop(y,'inverse); %1.01
remprop(y,'adjoint); %1.01
remprop(y,'rvalue); % 1.01
oporder!* := delete(y,oporder!*);
defoporder!* := delete(y,defoporder!*);
physoplist!* := delete(y,physoplist!*);
remprop(x,'rtype);
remprop(x,'tensdimen);
remprop(x,'phystype);
remprop(x,'psimpfn);
remprop(x,'inverse); %1.01
remprop(x,'adjoint); %1.01
remprop(x,'rvalue); % 1.01
oporder!* := delete(x,oporder!*);
defoporder!* := delete(x,defoporder!*);
physoplist!* := delete(x,physoplist!*);
>>;
return nil
end;
Rlistat '(clearphysop);
%------ procedures for printing out physops correctly ---------
% we modify the standard MAPRINT routine to get control
% over the printing of PHYSOPs
%**** This section had to be modified for 3.4 **********************
symbolic procedure physoppri u; % modified 3.4
begin scalar x,y,z,x1;
x := if idp u then u else car u;
y := if idp u then nil else cdr u;
trwrite(physoppri,"x= ",x," y= ",y,"nat= ",!*nat," contract= ",
!*contract);
if !*nat and not !*contract then go to a;
% transform the physop name in a string in order not to loose the
% special characters
x:= compress append('!" . explode x,list('!"));
prin2!* x;
if y then << prin2!* "(";
obrkp!* := nil;
inprint('!*comma!*,0,y);
obrkp!* := t;
prin2!* ")" >>;
return u;
a:
x := reverse explode x;
if length(x) > 2 then
if cadr x = '!- then <<z := compress list('!-,'!1);
x := compress reverse pnth(x,4); >>
else if car x = '!+ then << z:='!+;
x:= compress reverse pnth(x,3); >>
else x := compress reverse x
else x := compress reverse x;
x:= compress append('!" . explode x,list('!"));
x1 := if y then x . y else x;
trwrite(physoppri,"x= ",x," z= ",z," x1= ",x1);
% if z then exptpri(get('expt,'infix),list(x1,z))
% the following is 3.4
if z then exptpri(list('expt,x1,z),get('expt,'infix))
else << prin2!* x;
if y then << prin2!* "(";
obrkp!* := nil;
inprint('!*comma!*,0,y);
obrkp!* := t;
prin2!* ")" >>
>>;
return u
end;
symbolic procedure maprint(l,p!*!*); %3.4 version
begin scalar p,x,y;
p := p!*!*; % p!*!* needed for (expt a (quotient ...)) case.
if null l then return nil
else if physopp l then return apply1('physoppri,l)
else if atom l
then <<if not numberp l or (not l<0 or p<=get('minus,'infix))
then prin2!* l
else <<prin2!* "("; prin2!* l; prin2!* ")">>;
return l >>
else if stringp l then return prin2!* l
else if not atom car l then maprint(car l,p)
else if ((x := get(car l,'pprifn)) and
not(apply2(x,l,p) eq 'failed)) or
((x := get(car l,'prifn)) and
not(apply1(x,l) eq 'failed))
then return l
else if x := get(car l,'infix) then <<
p := not x>p;
if p then <<
y := orig!*;
prin2!* "(";
orig!* := if posn!*<18 then posn!* else orig!*+3 >>;
% (expt a b) was dealt with using a pprifn sometime earlier than this
inprint(car l,x,cdr l);
if p then <<
prin2!* ")";
orig!* := y >>;
return l >>
else prin2!* car l;
prin2!* "(";
obrkp!* := nil;
y := orig!*;
orig!* := if posn!*<18 then posn!* else orig!*+3;
if cdr l then inprint('!*comma!*,0,cdr l);
obrkp!* := t;
orig!* := y;
prin2!* ")";
return l
end;
% ******* end of 3.4 modifications ********************
% ------- end of module printout -------------------------
% ------------- some default declarations -------------------
% this list contains operators which when appearing in expressions
% have unknown properties (unresolved expressions)
specoplist!* := list('dot,'comm,'anticomm,'opapply);
% unit,comm and anticomm operators
put('comm,'rtype,'physop);
put('comm,'phystype,'scalar);
put('comm,'psimpfn,'commsimp);
put('anticomm,'rtype,'physop);
put('anticomm,'phystype,'scalar);
put('anticomm,'psimpfn,'commsimp);
physoplist!* := list('comm,'anticomm);
scalop 'unit;
flag ('(unit comm anticomm opapply),'reserved);
endmodule;
end;