################################################################
#                                                              #
#  Discrete Light-Cone Quantized QED Algebraic Reduction Code  #
#                                                              #
#            Frank E. Harris        February 1999              #
#                                                              #
#         Department of Physics,  University of Utah           #
#                                                              #
#                             and                              #
#                                                              #
#        Quantum Theory Project, University of Florida         #
#             P. O. Box 118435, Gainesville, FL                #
#                                                              #
#  This version for Maple V, Release 5.  It will not run       #
#  properly on other versions.                                 #
#                                                              #
#  Copyright (c) by Frank E. Harris, 1999.                     #
#                                                              #
#  May be used, distributed, or modified without limitation,   #
#  but whether modified or not, may not be sold.               #
#                                                              #
#  This file consists of four modules.  They will work         #
#  properly with the enclosed tutorial program if they are     #
#  broken into separate files with the respective names        #
#  "qed", "module1", "module2", and "module3".                 #
#                                                              #
################################################################

### Module "qed" ###

read module1;
read module2;
read module3;
`QED program loading complete.  Input can now be accepted.`;

################################################################

### module1 ###

#  MODULE 1.  Type definitions, procedures "build" and "show" for
#  construction and display of special-type quantities, and functions
#  for defining dynamical operators.

#  Type definitions for vector-matrix operations.  Note "#" is part of
#  the internal symbol for spinor types, "@" for polarization-vector types.
#  Whenever these symbols are used in input, they must be placed in
#  back-quotes.

#  These are for convenience in using the procedure "build" and for
#  type testing

`type/row_spinor` := proc(x) options remember; type(x,`r#`) end:
`type/column_spinor` := proc(x) options remember; type(x,`c#`) end:
`type/spinor_matrix` := proc(x) options remember; type(x,`m#`) end:
`type/row_pol` := proc(x) options remember; type(x,`r@`) end:
`type/column_pol` := proc(x) options remember; type(x,`c@`) end:
`type/pol_matrix` := proc(x) options remember; type(x,`m@`) end:

#  These are the internal representations 

`type/r#`:=proc(x) options remember; type(x,[name,`ts$qty`,`ts$qty`,
           `ts$qty`,`ts$qty`]) and x[1]=`r#` end:
`type/r#z`:=proc(x) options remember; type(x,{`r#`,0}) end:
`type/c#`:=proc(x) options remember; type(x,[name,`ts$qty`,`ts$qty`,
           `ts$qty`,`ts$qty`]) and x[1]=`c#` end:
`type/m#`:=proc(x) options remember;
           type(x,[name,`r#z`,`r#z`,`r#z`,`r#z`]) and x[1]=`m#` end:
`type/r@`:=proc(x) options remember;
           type(x,[name,`ts$qty`,`ts$qty`]) and x[1]=`r@` end:
`type/r@z`:=proc(x) options remember; type(x,{`r@`,0}) end:
`type/c@`:=proc(x) options remember;
           type(x,[name,`ts$qty`,`ts$qty`]) and x[1]=`c@` end:
`type/c@z`:=proc(x) options remember; type(x,{`c@`,0}) end:
`type/m@`:=proc(x) options remember;
           type(x,[name,`r@z`,`r@z`]) and x[1]=`m@` end:
`type/m@z`:=proc(x) options remember; type(x,{`m@`,0}) end:
`type/r#r@`:=proc(x) options remember;
           type(x,[name,`r@z`,`r@z`,`r@z`,`r@z`]) and x[1]=`r#r@` end:
`type/r#r@z`:=proc(x) options remember; type(x,{`r#r@`,0}) end:
`type/r#c@`:=proc(x) options remember;
           type(x,[name,`c@z`,`c@z`,`c@z`,`c@z`]) and x[1]=`r#c@` end:
`type/r#c@z`:=proc(x) options remember; type(x,{`r#c@`,0}) end:
`type/c#r@`:=proc(x) options remember;
           type(x,[name,`r@z`,`r@z`,`r@z`,`r@z`]) and x[1]=`c#r@` end:
`type/c#c@`:=proc(x) options remember;
           type(x,[name,`c@z`,`c@z`,`c@z`,`c@z`]) and x[1]=`c#c@` end:
`type/m#r@`:=proc(x) options remember;
           type(x,[name,`r#r@z`,`r#r@z`,`r#r@z`,`r#r@z`]) and x[1]=`m#r@` end:
`type/m#c@`:=proc(x) options remember;
           type(x,[name,`r#c@z`,`r#c@z`,`r#c@z`,`r#c@z`]) and x[1]=`m#c@` end:
`type/c#m@`:=proc(x) options remember;
           type(x,[name,`m@z`,`m@z`,`m@z`,`m@z`]) and x[1]=`c#m@` end:
`type/r#m@`:=proc(x) options remember;
           type(x,[name,`m@z`,`m@z`,`m@z`,`m@z`]) and x[1]=`r#m@` end:
`type/r#m@z`:=proc(x) options remember; type(x,{`r#m@`,0}) end:
`type/m#m@`:=proc(x) options remember;
           type(x,[name,`r#m@z`,`r#m@z`,`r#m@z`,`r#m@z`]) and x[1]=`m#m@`end:


#  Type definitions for dynamical (creation/annihilation) operators
#  Operators of these types may be constructed using the function
#  definitions b_, b*, etc (the functions containing an asterisk must
#  be back-quoted in input).

#  These are for convenience in type testing

`type/elec_annih` := proc(x) options remember; type(x,`b_$`) end:
`type/elec_creat` := proc(x) options remember; type(x,`b*$`) end:
`type/pos_annih` := proc(x) options remember; type(x,`d_$`) end:
`type/pos_creat` := proc(x) options remember; type(x,`d*$`) end:
`type/phot_annih` := proc(x) options remember; type(x,`a_$`) end:
`type/phot_creat` := proc(x) options remember; type(x,`a*$`) end:

#  These are the internal representations

`type/b_$`:= proc(x) options remember;
            type(x,function) and op(0,x)=b_ end:
`type/d_$`:= proc(x) options remember;
            type(x,function) and op(0,x)=d_ end:
`type/b*$`:= proc(x) options remember;
            type(x,function) and op(0,x)=`b*` end:
`type/d*$`:= proc(x) options remember;
            type(x,function) and op(0,x)=`d*` end:
`type/a_$`:= proc(x) options remember;
            type(x,function) and op(0,x)=a_ end:
`type/a*$`:= proc(x) options remember;
            type(x,function) and op(0,x)=`a*` end:

#  Various more general types.  "ts" stands for "time/space";  ts
#  quantities have no spinor, polarization-vector, or dynamical-
#  operator components.

`type/&*` := proc(x) options remember; type(x,function) and op(0,x)=`&*` end:
`type/&m` := proc(x) options remember; type(x,function) and op(0,x)=`&m` end:
`type/&w` := proc(x) options remember; type(x,function) and op(0,x)=`&w` end:

`type/fermi$`:= proc(x) options remember;
	   type(x,{`b_$`,`b*$`,`d_$`,`d*$`}) end:
`type/bose$`:= proc(x) options remember;
	   type(x,{`a_$`,`a*$`}) end:
`type/matvec$`:= proc(x) options remember;
	   type(x,{`r#`,`c#`,`m#`,`r@`,`c@`,`m@`,`r#r@`,`r#c@`,`c#r@`,`c#c@`,
                `m#r@`,`m#c@`,`c#m@`,`r#m@`,`m#m@`}) end:
`type/dynop$`:= proc(x) options remember;
	   type(x,{`fermi$`,`bose$`}) end:
`type/gen$op` := proc(x) options remember;
         type(x,function) and substring(op(0,x),1..1)=`&` end:
`type/ts$qty`:= proc(x) options remember;
	   type(x,{`ts$op`,`ts$nop`}) end:
`type/adj$`:= proc(x) options remember;
   type(x,function) and (op(0,x)=adj or op(0,x)=conj) end:
`type/ts$nop`:= proc(x) local i,a; options remember;
   if type(x,{`+`,`*`,`^`,`gen$op`,`adj$`}) then a:= true;
      for i from 1 to nops(x) do a := a and type(op(i,x),`ts$nop`) od;
      a
   else type(x,{name,constant}) or type(x,function)
                                   and not type(x,{`gen$op`,`dynop$`})
   fi end:
`type/ts$op^`:= proc(x) options remember;
        type(x,[name]) and member(x[1],{`D1`,`D2`,`D+`}) end:
`type/ts$op`:= proc(x) local i,a,b; options remember;
   if type(x,{`+`,`*`,`^`,`gen$op`,`adj$`}) then a := true; b := false;
      for i from 1 to nops(x) do a := a and type(op(i,x),`ts$qty`);
                                 b := b or type(op(i,x),`ts$op`) od;
      a and b
   else type(x,`ts$op^`) or type(x,[name]) and (x[1]=`I+` or x[1]=`I++`)
   fi end:

`type/scalar` := proc(x) local i,a; options remember;
   if type(x,{`+`,`*`}) then a := true; for i from 1 to nops(x) do
      a := a and type(op(i,x),scalar) od; a
   else
      type(x,{`ts$qty`,`dynop$`})
      or type(x,`gen$op`) and type(op(1,x),scalar) and type(op(2,x),scalar)
   fi end:

#  Function definitions for dynamical operators.  Those
#  containing asterisks must be back-quoted in input.

#  In the following, s is restricted to 1/2 or -1/2 for b_, b*, d_, d*,
#  1 or -1 for a_, a*.

b_ := proc(s,n) options remember; 'procname(args)' end:
`b*` := proc(s,n) options remember; 'procname(args)' end:
d_ := proc(s,n) options remember; 'procname(args)' end:
`d*` := proc(s,n) options remember; 'procname(args)' end:
a_ := proc(s,n) options remember; 'procname(args)' end:
`a*` := proc(s,n) options remember; 'procname(args)' end:

#  General vector-matrix build facility.  First argument indicates the
#  arrangement of the remaining arguments;  it may be in an easily readable
#  format (row_spinor, column_pol, spinor_matrix, ...) or in the internal
#  format in back-quotes (`r#`, `c@`, `m#`, ...).  A matrix in which both
#  rows and columns refer to spinors, or to polarization-vectors, may be
#  built (with first argument `m#` or `m@`) as a list of 16 (for m#) or
#  4 (for m@) elements in row-major order (e.g. 11, 12, 13, 14, 21, 22,...),
#  as a column of rows--as in build(`c@`,R1,R2), where R1 and R2 are each
#  of the form [`r@`,X,X], or as a row of columns--as in build(`r@`,C1,C2),
#  where C1 and C2 are column vectors.  Compound quantities involving both
#  spinor and polarization-vector quantities may be constructed either as a
#  spinor (or spinor-matrix) with polarization-vector (or matrix) elements
#  or the reverse.  The programs presently only permit the construction of
#  a mixed (spinor+polarization) quantity as indicated above--it is not
#  possible, for example, to use a column of row-spinors with row polarization-
#  vector components to make a spinor matrix with row pol-vector components;
#  one must instead either construct directly the spinor matrix, each of whose
#  components is a row pol-vector, or one may construct a row pol-vector
#  with spinor-matrix components.

#  Irrespective of the mode of construction, all compound types are brought
#  to a standard organization;  for the example of the preceding paragraph,
#  the organization is as a spinor-matrix with row-pol components.

build := proc() local x,G11,G12,G21,G22,i; options remember;
x := args[1];
if x=row_spinor then x := `r#`
elif x=column_spinor then x := `c#`
elif x=spinor_matrix then x := `m#`
elif x=row_pol then x := `r@`
elif x=column_pol then x := `c@`
elif x=pol_matrix then x := `m@` fi;
if (x=`m@` or x=`r#` or x=`c#`) and nargs=5 then
      G11:=args[2]; G12:=args[3]; G21:=args[4]; G22:=args[5];
elif (x=`r@` or x=`c@`) and nargs=3 then G11:=args[2]; G12:=args[3];
   G21:=G11; G22:=G12
elif x=`m#` and nargs=17 then G11:=args[2]
else RETURN('procname(args)') fi;
if nargs=17 and type(G11,`ts$qty`) then for i from 3 to 17 do
   if not type(args[i],`ts$qty`) then RETURN('procname(args)') fi od
elif nargs=17 and type(G11,{`r@`,`c@`,`m@`}) then for i from 3 to 17 do
   if not args[i][1] = G11[1] then RETURN('procname(args)') fi od
elif nargs=17 then RETURN('procname(args)') fi;
if nargs < 17 and 
   not(type(G11,`r#`) and type(G12,`r#`) and type(G21,`r#`) and type(G22,`r#`)
    or type(G11,`c#`) and type(G12,`c#`) and type(G21,`c#`) and type(G22,`c#`)
    or type(G11,`m#`) and type(G12,`m#`) and type(G21,`m#`) and type(G22,`m#`)
    or type(G11,`r@`) and type(G12,`r@`) and type(G21,`r@`) and type(G22,`r@`)
    or type(G11,`c@`) and type(G12,`c@`) and type(G21,`c@`) and type(G22,`c@`)
    or type(G11,`m@`) and type(G12,`m@`) and type(G21,`m@`) and type(G22,`m@`)
       or type(G11,`ts$qty`) and type(G12,`ts$qty`) and
          type(G21,`ts$qty`) and type(G22,`ts$qty`))
then RETURN('procname(args)') fi;
if type(G11,`ts$qty`) and x=`m@` then
                      RETURN([`m@`,[`r@`,G11,G12],[`r@`,G21,G22]])
elif type(G11,`ts$qty`) and (x=`r@` or x=`c@`) then RETURN([x,G11,G12])
elif type(G11,`ts$qty`) and (x=`r#` or x=`c#`) then
                                 RETURN([x,G11,G12,G21,G22])
elif type(G11,`ts$qty`) and x=`m#` then RETURN([x,[`r#`,args[2..5]],
              [`r#`,args[6..9]],[`r#`,args[10..13]],[`r#`,args[14..17]]])
elif type(G11,`ts$qty`) then RETURN('procname(args)') fi;
if type(G11,{`r@`,`c@`,`m@`}) and x=`m#` then
    i:=cat(`r#`,G11[1]); RETURN([cat(x,G11[1]),
    [i,args[2..5]],[i,args[6..9]],[i,args[10..13]],[i,args[14..17]]]) fi;
if G11[1]=`m#` and x=`m@` then
                RETURN([`m#`.x,build(x,G11[2],G12[2],G21[2],G22[2]),
                               build(x,G11[3],G12[3],G21[3],G22[3]),
                               build(x,G11[4],G12[4],G21[4],G22[4]),
                               build(x,G11[5],G12[5],G21[5],G22[5])])
elif G11[1]=`m#` and (x=`r@` or x=`c@`) then
                RETURN([`m#`.x,build(x,G11[2],G12[2]),
                               build(x,G11[3],G12[3]),
                               build(x,G11[4],G12[4]),
                               build(x,G11[5],G12[5])])
elif (G11[1]=`r#` or G11[1]=`c#`) and (x=`r@` or x=`c@`) then
     [cat(G11[1],x),[x,G11[2],G12[2]],[x,G11[3],G12[3]],
                    [x,G11[4],G12[4]],[x,G11[5],G12[5]]]
elif (G11[1]=`r#` or G11[1]=`c#`) and x=`m@` then
     [cat(G11[1],x),[x,[`r@`,G11[2],G12[2]],[`r@`,G21[2],G22[2]]],
                    [x,[`r@`,G11[3],G12[3]],[`r@`,G21[3],G22[3]]],
                    [x,[`r@`,G11[4],G12[4]],[`r@`,G21[4],G22[4]]],
                    [x,[`r@`,G11[5],G12[5]],[`r@`,G21[5],G22[5]]]]
elif type(G11,{`r@`,`c@`,`m@`}) and (x=`r#` or x=`c#`) then
     [cat(x,G11[1]),G11,G12,G21,G22]
elif x=`c#` and type(G11,`r#`) then [`m#`,G11,G12,G21,G22]
elif x=`r#` and type(G11,`c#`) then [`m#`,
                    [`r#`,G11[2],G12[2],G21[2],G22[2]],
                    [`r#`,G11[3],G12[3],G21[3],G22[3]],
                    [`r#`,G11[4],G12[4],G21[4],G22[4]],
                    [`r#`,G11[5],G12[5],G21[5],G22[5]]]
elif x=`c@` and type(G11,`r@`) then [`m@`,G11,G12]
elif x=`r@` and type(G11,`c@`) then [`m@`,[`r@`,G11[2],G12[2]],
                    [`r@`,G11[3],G12[3]]]
else 'procname(args)' fi
end:

#  The procedure "show" displays matrix or vector quantities of the types
#  that can be generated by the "build" procedure in conventional matrix
#  tableaux.  The output of show is for display only, and cannot be used
#  as input to the other procedures of this program set.

show := proc(x) local i,j,T; options operator,remember;
if type(x,{`&*`,`&m`}) then
       T:= show(op(1,x))&*show(op(2,x))
elif type(x,function) and op(0,x)=`&++` then T:= show(op(1,x)) + show(op(2,x))
elif type(x,function) and op(0,x)=`&--` then
   if nops(x)=2 then T:= show(op(1,x)) - show(op(2,x))
   elif nops(x)=1 then T := - show(op(1,x)) fi
elif type(x,`+`) then T := x; for i from 1 to nops(x) do
   T := subsop(i=show(op(i,x)),T) od
elif type(x,{`dynop$`,`ts$qty`}) then T:=x;
elif type(x,{`r#`,`r#r@`,`r#c@`,`r#m@`}) then T:=array(1..1,1..4);
                   for i from 1 to 4 do T[1,i]:=show(op(i+1,x)) od
elif type(x,{`c#`,`c#r@`,`c#c@`,`c#m@`}) then T:=array(1..4,1..1);
	           for i from 1 to 4 do T[i,1]:=show(op(i+1,x)) od
elif type(x,{`m#`,`m#r@`,`m#c@`,`m#m@`}) then T:=array(1..4,1..4);
  for i from 1 to 4 do for j from 1 to 4 do
    if op(i+1,x)=0 then T[i,j]:=0 else T[i,j]:=show(op(j+1,op(i+1,x))) fi od od
elif type(x,`r@`) then T:=array(1..1,1..2);
	           for i from 1 to 2 do T[1,i]:=op(i+1,x) od
elif type(x,`c@`) then T:=array(1..2,1..1);
	           for i from 1 to 2 do T[i,1]:=op(i+1,x) od
elif type(x,`m@`) then T:=array(1..2,1..2);
  for i from 1 to 2 do for j from 1 to 2 do 
    if op(i+1,x)=0 then T[i,j]:=0 else T[i,j]:=op(j+1,op(i+1,x)) fi od od
else RETURN('procname(args)') fi;
eval(T)
end:

`QED Module 1 loaded`;

################################################################

### module2 ###

#  MODULE 2.  The general processing procedure, "evalqed", and
#  the special operations needed to implement it (general multiplication
#  &m, general addition/subtraction &++ and &--).  These special
#  operations are not input by the user (read the next paragraph), but
#  are introduced by the expression processor and may occur in output if
#  the program set cannot simplify them away.

#  Expressions should be written using ordinary + and - for
#  addition and subtraction and &* for multiplication (this ensures that
#  factors stay in correct order when multiplication is noncommutative).
#  After an expression EXPRESSION has been built in this way, carrying
#  out the operation evalqed(EXPRESSION) will evaluate it or reduce its
#  complexity.

#  If EXPRESSION is not reduced to a scalar, addition and
#  subtraction will be denoted by &++ and &--, and multiplication
#  of non-scalar objects will be denoted by &m.

#  If EXPRESSION is reduced to a scalar, derivative application is then
#  expanded to individual factors--right now the derivatives are
#  only indicated symbolically as D1, D2, or D+.

#  Note: Do NOT activate Maple's definitions of matrix operations
#  when using this set of procedures.

#  When using +, -, and &* the normal rules for precedence of
#  operations will apply, so that parentheses are needed only if
#  required in conventional algebraic notation.  If "evalqed" cannot
#  reduce an expression to scalar form, the precedence of operations
#  will be explicitly indicated and controlled by parentheses in the output.

evalqed := proc() local x,i,xx,t; options remember;
   x := ckz(args[1]);
   if nargs=2 then i := args[2] else i := 0 fi;
    if type(x,`+`) and (has(x,D_sum) or has(x,d_sum) or has(x,S_sum)
       or has(x,s_sum)) then
  RETURN(map(evalqed,x,i)) fi;
   if type(x,function) and op(0,x)=s_sum then
      xx := evalqed(op(1,x),i); t := op(2,x);
      RETURN(combine(Sum(combine(sub_spin(xx,t)),t))) fi;
   if type(x,function) and op(0,x)=D_sum then
     if nops(x)=2 then RETURN(d_sum(evalqed(op(1,x),i),op(2,x)))
      else xx := d_sum(evalqed(D_sum(op(1,x),op(3..nops(x),x)),i),op(2,x));
         RETURN(combine(simplify(xx,Sum)))
      fi fi;
	xx := _evalqed(x);       # basic reduction of operator sequences
xx := d_reduce(xx);      # applies derivatives and inverse derivatives
xx := eval(subs(exp=exp_R,xx));
xx := expand(xx,`&*`);
xx := simplify(xx,exp);
xx := eval(subs(exp=Delta,xx)); # Integrates to Kronecker Delta functions
                                # Drops those which cannot be satisfied
if i=1 then xx := subs(Delta=delta,xx) fi;
xx := expand(xx,`&*`)  ;        # Make a sum of individual terms
 end:

_evalqed := proc(x) local n,xx; options remember;
if type(x,`ts$qty`) then xx := x
elif type(x,`*`) and nops(x)=2 and op(1,x)=-1 then
    xx := &-- _evalqed(op(2,x))
elif type(x,`+`) then
  if type(op(1,x),`*`) and nops(op(1,x))=2 and op(1,op(1,x))=-1 then
    xx := &-- _evalqed(op(2,op(1,x)))
  else
    xx := _evalqed(op(1,x)) fi;
  for n from 2 to nops(x) do
   if type(op(n,x),`*`) and nops(op(n,x))=2 and op(1,op(n,x))=-1 then
       xx := xx &-- _evalqed(op(2,op(n,x)))
   else
       xx := xx &++ _evalqed(op(n,x)) fi od
  else xx := map(_evalqed,x) fi;
  xx := subs(`&*`=`&m`,xx);
  xx := a_left(xx); a_left(xx) end:

#  This is the basic (necessarily complicated) procedure for
#  carrying out general multiplications involving the vector and
#  matrix quantities defined for this set of procedures.

`&m`:=proc(x,y) local i,t,z; options remember;
    if ckz(x)=0 or ckz(y)=0 then RETURN(0) fi;
    if type(x,function) and op(0,x)=`&++` then
       RETURN(ckz(a_left(op(1,x) &m y) &++ a_left(op(2,x) &m y))) fi;
    if type(y,function) and op(0,y)=`&++` then
       RETURN(ckz(a_left(x &m op(1,y)) &++ a_left(x &m op(2,y)))) fi;
    if type(x,`dynop$`) and type(y,{`ts$qty`,`matvec$`}) then
       RETURN(y &m x) fi;
    if type(x,`&m`) and
       (type(op(2,x),`dynop$`) and type(y,{`ts$qty`,`matvec$`}) or
       not type(op(2,x),constant) and type(y,constant))
       then RETURN((op(1,x) &m y) &m op(2,x)) fi;
    if type(x,`&m`) then 
       RETURN(op(1,x) &m (op(2,x) &m y)) fi;
    if type(x,constant) and type(y,`ts$op`) or
       type(x,`ts$op`) and type(y,constant) or
       type(x,`ts$nop`) and type(y,`ts$nop`) then RETURN(x*y) fi;
    if type(x,`ts$qty`) and type(y,{`r#`,`c#`,`m#`,`r#r@`,`c#r@`,`m#r@`,
         `r#c@`,`c#c@`,`m#c@`,`r#m@`,`c#m@`,`m#m@`}) then z:=y[2..5];
       RETURN(ckz([y[1],op(map2(`&m`,x,z))])) fi;
    if type(x,`ts$qty`) and type(y,{`r@`,`c@`,`m@`}) then z:=y[2..3];
       RETURN(ckz([y[1],op(map2(`&m`,x,z))])) fi;
    if type(x,{`c#`,`r#`,`m#`,`c#r@`,`r#r@`,`m#r@`,`c#c@`,
       `r#c@`,`m#c@`,`c#m@`,`r#m@`,`m#m@`}) and type(y,`ts$qty`) then
       z:=x[2..5]; RETURN(ckz([x[1],op(map(`&m`,z,y))])) fi;
    if type(x,{`r@`,`c@`,`m@`}) and type(y,`ts$qty`) then z:=x[2..3];
       RETURN(ckz([x[1],op(map(`&m`,z,y))])) fi;
    if type(x,{`c#`,`c#r@`,`c#c@`,`c#m@`}) and
         type(y,{`c#`,`m#`,`c#r@`,`m#r@`,`c#c@`,`m#c@`,`c#m@`,`m#m@`})
       or type(x,{`r#`,`m#`,`r#r@`,`m#r@`,`r#c@`,`m#c@`,`r#m@`,`m#m@`}) and
         type(y,{`r#`,`r#r@`,`r#c@`,`r#m@`})
       or type(x,{`r@`,`r#r@`,`c#r@`,`m#r@`,`m@`,`r#m@`,`c#m@`,`m#m@`}) and
         type(y,{`r@`,`r#r@`,`c#r@`,`m#r@`})
       or type(x,{`c@`,`r#c@`,`c#c@`,`m#c@`}) and
         type(y,{`c@`,`m@`,`r#c@`,`r#m@`,`c#c@`,`c#m@`,`m#c@`,`m#m@`})
   or type(x,`ts$qty`) and type(y,`ts$qty`)
       or type(y,{`dynop$`,`ts$op`})       
       then RETURN('procname(args)') fi;
    if type(y,function) and op(0,y)=`&m` then 
       RETURN('procname(args)') fi;
    if type(x,{`m#`,`m#r@`,`m#c@`,`m#m@`}) and type(y,{`c#`,`c#r@`,`c#c@`,
         `c#m@`,`m#`,`m#r@`,`m#c@`,`m#m@`}) then z:=x[2..5];
         RETURN(ckz([suffix2(y,x,y),op(map(`&m`,z,y))])) fi;
    if type(x,{`r#`,`r#r@`,`r#c@`,`r#m@`}) and
      type(y,{`m#`,`m#r@`,`m#c@`,`m#m@`}) then z:=xms(y)[1..4];
      RETURN(ckz([suffix2(x,x,y),op(map2(`&m`,x,z))])) fi;
    if type(x,{`r@`,`c@`,`m@`}) and type(y,{`r#`,`c#`,`m#`}) then
      RETURN(y &m x) fi;
    if type(x,{`r#`,`c#`,`m#`}) and type(y,{`r@`,`c@`,`m@`}) then
      z:=x[2..5]; RETURN(ckz([cat(x[1],y[1]),op(map(`&m`,z,y))])) fi;
    if type(x,{`r#`,`r#r@`,`r#c@`,`r#m@`}) and
       type(y,{`c#`,`c#r@`,`c#c@`,`c#m@`}) then
       RETURN(ckz((x[2] &m y[2])&++(x[3] &m y[3])&++(x[4] &m y[4])
                  &++(x[5] &m y[5]))) fi;
    if type(x,`r@`) and type(y,`c@`) then
       RETURN(ckz((x[2] &m y[2])&++(x[3] &m y[3]))) fi;
    if type(x,`r@`) and type(y,`m@`) then
       z:=xms(y)[1..2]; RETURN(ckz([`r@`,op(map2(`&m`,x,z))])) fi;
    if type(x,`m@`) and type(y,{`c@`,`m@`}) then
       z:=x[2..3]; RETURN(ckz([y[1],op(map(`&m`,z,y))])) fi;
    if type(x,`c@`) and type(y,`r@`) then
       RETURN(ckz([`m@`,x[2] &m y,x[3] &m y])) fi;
    if type(x,{`c#`,`c#r@`,`c#c@`,`c#m@`}) and
       type(y,{`r#`,`r#r@`,`r#c@`,`r#m@`}) then
       RETURN(ckz([suffix2([`m#`],x,y),x[2] &m y,x[3] &m y,
                                    x[4] &m y,x[5] &m y])) fi;
    if type(x,`r@`) and type(y,{`r#c@`,`c#c@`,`m#c@`}) then z:=y[2..5];
       RETURN(ckz([substring(y[1],1..2),op(map2(`&m`,x,z))])) fi;
    if type(x,{`r#r@`,`c#r@`,`m#r@`,`r#m@`,`c#m@`,`m#m@`}) and type(y,`c@`)
       then z:=x[2..5];
       RETURN(ckz([suffix2(x,x,[`r#c@`]),op(map(`&m`,z,y))])) fi;
    if type(x,`r@`) and type(y,{`r#m@`,`c#m@`,`m#m@`}) then z:=y[2..5];
      RETURN(ckz([cat(substring(y[1],1..2),`r@`),op(map2(`&m`,x,z))])) fi;
    if type(x,`m@`) and type(y,{`r#c@`,`c#c@`,`m#c@`,`r#m@`,`c#m@`,`m#m@`})
       then z:=y[2..5]; RETURN(ckz([y[1],op(map2(`&m`,x,z))])) fi;
    if type(x,{`r#r@`,`c#r@`,`m#r@`,`r#m@`,`c#m@`,`m#m@`}) and type(y,`m@`)
       then z:=x[2..5]; RETURN(ckz([x[1],op(map(`&m`,z,y))])) fi;
    if type(x,`c@`) and type(y,{`r#r@`,`c#r@`,`m#r@`}) then z:=y[2..5];
      RETURN(ckz([cat(substring(y[1],1..2),`m@`),op(map2(`&m`,x,z))])) fi;
    if type(x,{`r#c@`,`c#c@`,`m#c@`}) and type(y,`r@`) then z:=x[2..5];
       RETURN(ckz([cat(substring(x[1],1..2),`m@`),op(map(`&m`,z,y))])) fi;
  print(`Being processed`,'procname(args)');
    ERROR(`Invalid type in &*`)
end:

#  Generalization of "+".  Does not combine addends until they have
#  been individually reduced to a standard scalar, vector, or matrix
#  form, and returns ERROR if the reduced forms exhibit a type
#  mismatch.    

`&++`:=proc(x,y) options remember;
   if x=0 then RETURN(y) fi;
   if y=0 then RETURN(x) fi;
   if type(x,`ts$qty`) and type(y,`ts$qty`) then RETURN(x+y) fi;
   if type(x,`matvec$`) and type(y,`matvec$`) then
     if evalb(x[1]<>y[1]) then ERROR(`Type mismatch in + or -`) fi;
     if type(x,{`r#`,`c#`,`m#`,`r#r@`,`c#r@`,`m#r@`,`r#c@`,`c#c@`,`m#c@`,
             `r#m@`,`c#m@`,`m#m@`})
       then RETURN([x[1],x[2]&++y[2],x[3]&++y[3],x[4]&++y[4],x[5]&++y[5]]) fi;
     if type(x,{`r@`,`c@`,`m@`}) then RETURN([x[1],x[2]&++y[2],x[3]&++y[3]]) fi;
   fi;
   'procname(args)' end:

#  Generalization of "-".  Recognizes both unary and binary "-"
#  situations.

`&--` := proc() options remember; if nargs=1 then (-1) &* args[1]
        else args[1] &++ ((-1) &* args[2]) fi end:

#####  Utility procedures used in the above:

#  Interchanges order of factors for convenient use of "map" function

#####`rev&m`:=proc(x,y); y &m x end:

#  Builds compound type indicators

`suffix2` := proc(ww,xx,yy) local x,y,w;
  x := substring(xx[1],3..4);
  y := substring(yy[1],3..4);
  w := substring(ww[1],1..2);
  if length(x)=0 then RETURN(cat(w,y)) fi;
  if length(y)=0 or y=`m@` then RETURN(cat(w,x)) fi;
  if x=`m@` then  RETURN(cat(w,`c@`)) fi;
  if x=`r@` then RETURN(w) fi;
  cat(w,`m@`) end:

#  Replaces matrix or vector quantities all of whose elements are
#  zero by the single number zero.

ckz := proc(x) options remember;
   if type(x,`matvec$`) and op(2,x)=0 and op(3,x)=0 and
      (nops(x)=3 or op(4,x)=0 and op(5,x)=0) then 0 else x fi end:

#  Converts matrix to an array of columns to facilitate multiplication
#  from the left.  Output of "xms" does not contain an overall
#  type designator and is not a standard-form matrix quantity.

xms:=proc(x) local i,j,T,y,z,result; options operator,remember;
    if x=0 then RETURN(0) fi;
    if substring(x[1],1..2)=`m#` then
      y:=`c#`; z:=cat(y,substring(x[1],3..4)); 
      T:=array(1..4,1..4);  for i from 1 to 4 do
        if x[i+1]=0 then for j from 1 to 4 do T[i,j]:=0 od
        else for j from 1 to 4 do T[i,j]:= x[i+1][j+1] od fi od;
      result:=[ckz([z,T[1,1],T[2,1],T[3,1],T[4,1]]),
                    ckz([z,T[1,2],T[2,2],T[3,2],T[4,2]]),
                    ckz([z,T[1,3],T[2,3],T[3,3],T[4,3]]),
                    ckz([z,T[1,4],T[2,4],T[3,4],T[4,4]])];
    elif type(x,`m@`) then
      z:=`c@`; T:=array(1..2,1..2);
      for i from 1 to 2 do 
        if x[i+1]=0 then for j from 1 to 2 do T[i,j]:=0 od
        else for j from 1 to 2 do T[i,j]=x[i+1][j+1] od fi od;
      result:=[ckz([z,T[1,1],T[2,1]]),ckz([z,T[1,2],T[2,2]])];
    fi end:

#  Forces left-association of operator &t

`&t` := proc(x,y) options remember;
       if type(y,function) and op(0,y)=`&t` then
              RETURN(x &t op(1,y) &t op(2,y))
       else RETURN('procname(args)')
       fi end:

#  Substitutes &t for &m temporarily to force left-association of &m

a_left := proc(x) local xx; options remember;
   xx := eval(subs(`&m`=`&t`,x));
   xx := eval(subs(`&t`=`&m`,xx));
end:

`&v` := proc(x,y); x * y end:

`&p` := proc(x,y); x + y end:

#  Applies derivatives

#  f, D1, D2, and D+ must be undefined functions;  they correspond
#  to the defined functions F, DD1, DD2, and DD+.

d_reduce := proc(x) local xx; options remember;
   if not type(x,scalar) then RETURN(x) fi;
   xx := eval(subs(`&++`=`&p`,x));
   xx := eval(subs(`&m`=`&w`,xx));
   xx := aa_right(xx);
   xx := eval(subs(`&w`=`&u`,xx));
   xx := eval(subs(f=F,D1=DD1,D2=DD2,`D+`=`DD+`,xx));
   xx := subs(`&u`=`&w`,xx);
   xx := aa_right(xx);
   xx := do_I(xx);
   xx := eval(fix_mult(xx));
   subs(`&w`=`&*`,xx)
end:

fix_mult := proc(x) local T; global `S^`;
if type(x,{`+`,`*`}) then map(fix_mult,x)
elif type(x,`&w`) then
   if not type(op(1,x),`dynop$`) then op(1,x)*fix_mult(op(2,x))
   else `S^`:=1; T := normal_order(x); `S^`*T fi
else x fi
end:

#  Always called with x of pattern (... &w ...)
normal_order := proc(x) local T; global `S^`;
if type(op(2,x),`&w`) then
   T := normal_order(op(2,x));
   if switch_order(op(1,x),op(1,T)) then 
      op(1,T) &w normal_order(op(1,x) &w op(2,T)) else op(1,x) &w T fi
elif switch_order(op(1,x),op(2,x)) then op(2,x) &w op(1,x)
else x fi end:

switch_order := proc(x,y) global `S^`;
if not type(x,`a*$`) and type(y,`a*$`) or
   type(x,`fermi$`) and type(y,`bose$`) then true
elif type(x,{`b_$`,`d_$`}) and type(y,{`b*$`,`d*$`}) or
   type(x,`d*$`) and type(y,`b*$`) or type(x,`d_$`) and type(y,`b_$`) then
   `S^`:= -`S^`; true
else false fi end:


F := proc(x,n); exp(I*(-x.1*n.1-x.2*n.2+ x.`-`*n.`+`/2)) end:
DD1 := proc(x); diff(I*x,x1) end:
DD2 := proc(x); diff(I*x,x2) end:
`DD+` := proc(x); diff(2*I*x,`x-`) end:   

do_I := proc(x) local i,m,T;
if type(x,`+`) then map(do_I,x)
elif type(x,`&w`) then
    if op(1,x)=[`I+`] or op(1,x)=[`I++`] then
       T := eval(subs(`&w`=`&v`,op(2,x)));
       m := normal(expand(`DD+`(T)/T));
      if m<>0 then
         if op(1,x)=[`I+`] then m := 1/m else m := 1/m^2 fi;
      m &w op(2,x) fi
  else op(1,x) &w do_I(op(2,x)) fi
else x fi
end:

`diff/&u` := proc(x,y,z);
diff(x,z) &u y + x &u diff(y,z) end:

`diff/conj` := proc(x,z);
conj(diff(x,z)) end:

`&w` := proc(x,y) local i; options remember;
if x=0 or y=0 then RETURN(0) fi;
if type(x,`&w`) then op(1,x) &w (op(2,x) &w y)
elif type(x,`+`) then map(`&w`,x,y)
else 'procname(args)' fi end:

aa_right := proc(x) local i,n,t;
if type(x,`*`) then n:=nops(x); t:=aa_right(op(n,x));
   for i from n-1 by -1 to 1 do t:=aa_right(op(i,x)) &w t od
elif type(x,`+`) then map(aa_right,x)
elif type(x,`^`) then aa_right(op(1,x))^op(2,x)
elif type(x,{`&m`,`&w`}) then
   aa_right(op(1,x)) &w aa_right(op(2,x))
else t := x fi end:

`&u` := proc(x,y) local n,T; options remember;
if type(x,`ts$op`) and type(x,[name]) and member(x[1],{D1,D2,`D+`}) then
   x[1](y)
else 'procname(args)' fi end:

exp_R := proc(z) local n,nn,T;
T := sort(collect(z,[`x-`,x1,x2]),[`x-`,x1,x2]);
n := normal(I*op(2,T)/x1);
if type(n,name) then nn := substring(n,1..length(n)-1)
else nn := -n; nn := -substring(nn,1..length(nn)-1) fi;
exp(nn) end:

Delta := proc(x) local T; options remember;
T := fix_names(x);
if is(T <> 0) then 0 else 'procname(args)' fi end:

fix_names := proc(x);
if type(x,string) then assume(D_.x > 0); D_.x
elif type(x,{`+`,`*`}) then map(fix_names,x)
else x fi end:

S_sum := proc() local i; options remember;
if nargs > 2 then
     i := args[nargs];
     S_sum(sum(args[1],i=(-1/2)..(1/2)),args[2..nargs-1])
else s_sum(eval(subs(args[2]=1/2,args[1])),args[2]) fi end:

d_sum := proc(x,m) local n,mm,mm1,mm2,`mm+`; options remember;
if type(x,`+`) then map(d_sum,x,m)
elif not has(x,Delta) then if x<>0 then Sum(x,m) else 0 fi
else n := arg_Delta(x); mm := solve(n=0,m);
if mm=NULL then
   print(`MAPLE can't solve for`,m,` this simple equation`,n=0);
   ERROR('procname(args)') fi;
#  Now change m1, m2, m+ to mm1, mm2, mm+ whatever those are
mm1 := add_suffix(mm,1);
mm2 := add_suffix(mm,2);
`mm+` := add_suffix(mm,`+`);
eval(subs(m=mm,m.1=mm1,m.2=mm2,m.`+`=`mm+`,Delta=Unity,x)) fi end:

add_suffix := proc(mm,n);
if type(mm,name) then mm.n
elif type(mm,{`+`,`*`}) then map(add_suffix,mm,n)
else mm fi end:

arg_Delta := proc(x) local i;
if type(x,function) and op(0,x)=Delta then op(1,x)
else for i from 1 to nops(x) do
   if type(op(i,x),function) and op(0,op(i,x))=Delta then
      RETURN(op(1,op(i,x))) fi od fi end:

Unity := proc(x); 1 end:

#   Changes all spin 1/2 to t, all spin -1/2 to -t
#      and pol 1 to 2*t, pol -1 to -2*t
sub_spin := proc(x,t);
if type(x,{`+`,`*`}) then map(sub_spin,x,t)
#elif type(x,function) and op(0,x)=`&*` then map(sub_spin,x,t)
elif type(x,`&*`) then map(sub_spin,x,t)
elif type(x,function) and op(0,x)=sum or op(0,x)=Sum then
         subsop(1=sub_spin(op(1,x),t),x)
elif type(x,`fermi$`) then 
         if op(1,x)=1/2 then subsop(1=t,x) else subsop(1=-t,x) fi
elif type(x,`bose$`) then
         if op(1,x)=1 then subsop(1=2*t,x) else subsop(1=-2*t,x) fi
else x fi end:

`simplify/Sum` := proc(x) local i,S,T; options remember;
if type(x,{`+`,`*`}) then map(simplify,x,Sum)
elif type(x,function) and op(0,x)=Sum then
   if type(op(1,x),`*`) then S := 1; T := 1;
      for i from 1 to nops(op(1,x)) do
         if type(op(i,op(1,x)),constant) then S := S*op(i,op(1,x));
         else T := T*op(i,op(1,x)) fi od;
         if S <> 1 then S*simplify(Sum(T,op(2,x)),Sum) else x fi
   elif type(op(1,x),function) and op(0,op(1,x))=Sum then
      i := cat(op(2,x),`,`,op(2,op(1,x)));
      simplify(Sum(op(1,op(1,x)),i),Sum)
   else x fi
else x fi end:

`QED Module 2 loaded`;

################################################################

### module3 ###

#  MODULE 3.  Adjoint, transposition, and conjugation operations

#  Caution is needed if a differential operator is within the scope of
#  an adjoint command, as "adjoint" reverses the order of all operands
#  of the &* operator.

#  Takes complex conjugate (assumes simple names refer to real variables)
#  Leaves dynamical operators unaffected

#  Leaves D1^... unaffected in spite of the factor I, so that adjoint
#  will be properly computed

conj := proc(x) local t; options operator,remember;
  if type(x,I) then RETURN(-I) fi;
  if type(x,{numeric,`dynop$`}) then RETURN(x) fi;
  if type(x,function) then
     if op(0,x) = conj and nops(x)=1 then RETURN(op(1,x))
     elif op(0,x) = exp then RETURN(exp(conj(op(1,x)))) fi fi;
   if type(x,{`+`,`*`,`^`,`gen$op`}) then RETURN(map(conj,x)) fi;
  if type(x,{`m#`,`r#`,`c#`,`m#m@`,`m#c@`,`m#r@`,`r#m@`,`r#c@`,`r#r@`,`c#m@`,
             `c#c@`,`c#r@`})
     then RETURN([x[1],conj(x[2]),conj(x[3]),conj(x[4]),conj(x[5])]) fi;
  if type(x,{`m@`,`c@`,`r@`}) then RETURN([x[1],conj(x[2]),conj(x[3])]) fi;
  if type(x,{name,[name]}) and not type(x,indexed) then RETURN(x) fi;
  'procname(args)' end:

#  Transpose spinor or spinor-matrix indices (no action on pol-vectors
#  or pol-matrices; processes only single matrix/vector quantities)

trans4 := proc(x) options operator,remember;
  if x=0 then RETURN(0) fi;
  if type(x,{`m#`,`m#r@`,`m#c@`,`m#m@`}) then RETURN(_trans4(x)) fi;
  if type(x,{`r#`,`r#r@`,`r#c@`,`r#m@`}) then
     RETURN([cat(`c#`,substring(x[1],3..4)),x[2],x[3],x[4],x[5]]) fi;
  if type(x,{`c#`,`c#r@`,`c#c@`,`c#m@`}) then
     RETURN([cat(`r#`,substring(x[1],3..4)),x[2],x[3],x[4],x[5]]) fi;
     x end:

_trans4 := proc(x) local i,j,T,y,z,result; options operator,remember;
     if x=0 then RETURN(0) fi;
     z:=cat(`r#`,substring(x[1],3..4)); T:=array(1..4,1..4);
     for i from 1 to 4 do
       if x[i+1]=0 then for j from 1 to 4 do T[i,j]:=0 od
       else for j from 1 to 4 do T[i,j]:=x[i+1][j+1] od fi od;
     result:=[x[1],ckz([z,T[1,1],T[2,1],T[3,1],T[4,1]]),
                   ckz([z,T[1,2],T[2,2],T[3,2],T[4,2]]),
                   ckz([z,T[1,3],T[2,3],T[3,3],T[4,3]]),
                   ckz([z,T[1,4],T[2,4],T[3,4],T[4,4]])]
     end:

#  Transpose polarization-vector or matrix indices (no action on
#  spinors or spinor-matrices).  Processes single matrix/vector
#  quantities, sums of same, and products containing at most one
#  matrix/vector quantity (designated by either * or &*).

trans2 := proc(x) local t; options operator,remember;
  if x=0 then RETURN(0) fi;
  if type(x,{`+`,`*`}) or (type(x,function) and op(0,x)=`&*`) then
       RETURN(map(trans2,x)) fi;
  if type(x,`r@`) then RETURN([`c@`,x[2],x[3]]) fi;
  if type(x,`c@`) then RETURN([`r@`,x[2],x[3]]) fi;
  if type(x,{`+`,`*`,`&*`}) then
       RETURN(map(trans2,x)) fi;
  if type(x,`r@`) then RETURN([`c@`,x[2],x[3]]) fi;
  if type(x,`c@`) then RETURN([`r@`,x[2],x[3]]) fi;
  if type(x,`m@`) then RETURN(_trans2(x)) fi;
  if type(x,{`r#r@`,`c#r@`,`m#r@`}) then t:=cat(substring(x[1],1..2),`c@`)
  elif type(x,{`r#c@`,`c#c@`,`m#c@`}) then t:=cat(substring(x[1],1..2),`r@`)
  elif type(x,{`r#m@`,`c#m@`,`m#m@`}) then t:=x[1];
  else RETURN(x) fi;
  [t,trans2(x[2]),trans2(x[3]),trans2(x[4]),trans2(x[5])] end:

_trans2 := proc(x) local T; options operator,remember;
  if x=0 then RETURN(0) fi;
  T:=array(1..2,1..2); if x[2]=0 then T[1,1]:=0; T[1,2]:=0
                       else T[1,1]:=x[2][2]; T[1,2]:=x[2][3] fi;
                       if x[3]=0 then T[2,1]:=0; T[2,2]:=0
                       else T[2,1]:=x[3][2]; T[2,2]:=x[3][3] fi;
  [`m@`,ckz([`r@`,T[1,1],T[2,1]]),ckz([`r@`,T[1,2],T[2,2]])] end:

#  General adjoint procedure

adj := proc(x) options operator,remember;
  if x=I then RETURN(-I) fi;
  if type(x,{`+`,`*`,`^`}) then map(adj,x)
  elif type(x,`gen$op`) then op(0,x)(adj(op(2,x)),adj(op(1,x)))
  elif type(x,`dynop$`) then adjd(x)
  elif type(x,`matvec$`) then (conj@trans2@trans4)(x)
  elif type(x,indexed) or (type(x,function) and op(0,x)=conj) then
     'procname(args)'
  else conj(x) fi end:

#  Takes adjoints of dynamical operators

adjd := proc(x) local u; options remember;
if not type(x,`dynop$`) then RETURN('procname(args)') fi;
if substring(op(0,x),2..2)=`_` then u := `*` else u := `_` fi;
subsop(0=cat(substring(op(0,x),1..1),u),x) end:

`QED Module 3 loaded`;

################################################################
####################### end of program #########################