12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483 |
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % 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;
|