Artifact aaff83206ae75086be91da7ce23e9b1438ebc4e3d704b6c4bf29eb6988efc963:
- Executable file
r36/src/physop.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 88336) [annotate] [blame] [check-ins using] [more...]
module physop; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % 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.) % % <UNP008@DBNRHRZ1.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 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% create!-package('(physop),'(contrib physics)); %-------------------------------------------------------% % 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!* obrkp!*); newtok '((d o t) dot); flag ('(dot), 'spaced); % define some global variables from REDUCE needed in the package fluid '(alglist!* subfg!* wtl!*); 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!* ncmp!*); 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;