(*
shows >1 Meg left when Mathematica crashes.
        Anyway,  this program is not well debugged even apart from the
nasty crashes, but maybe someone will find it useful.

        I would also like to get the addresses of user groups on the
Internet,
if you have a list handy.

        Finally, do you wnat to be listed as the Mathematica SIG leader for
the Stanford Mac. Users' Group?  I noticed their latest newsletter did
not
show you in that position.  As a new officer, I can have your name
included
(or explicitly excluded, if you prefer) from future newsletters.

			Regards,

                        Gary Woods
                        gwoods@portia
                        415-857-1882

P.S. I suspect that the tabs are badly out of alignment
        in the program text.  My word processor did not
        have any non-proportional fonts.  Sorry.
*)

(* 
  Commutator & operator algebra for QM 
   Gary Woods, Nov. '89                                 
  This version always crashes (code 28)
   on my SE-30 for reasons unknown...   
*)

(* Conventions:
        op[a]           "operator a", any mathematical 
                                        object that in general does
                                        not commute.
        prod[a,b,...]   non-commutative product of
                                     operators
        com[a,b]        commutator from QM, i.e. [a,b]
                                        defined as  (a b) - (b a)
*)

(*************************)
(* Rules for commutators *)
(*************************)

        (*                                            *)
        (* Put operators in a standard order *)
        (*                                             *)
com[op[a_],op[b_]]              :=      -com[op[b],op[a]] /;
                                                        !OrderedQ[{a,b}]
                                                        
        (*               *)                                             
        (* Linearity *)
        (*               *)
com[a_+b_, c_]                  :=      com[a,c] + com[b,c]
com[a_, b_+c_]                  :=      com[a,b] + com[a,c]

        (*                                   *)
(* Numbers always commute *)
        (*                                   *)
com[a_?NumberQ b_, c_]  :=      prod[a,com[b,c]]
com[a_, b_?NumberQ c_]  :=      prod[b,com[a,c]]
com[a_, b_]                             :=      0       /;( NumberQ[a] ||
                                                                    NumberQ[b]  )
                                                                    
        (*                                *)
        (* Self-commutator is zero *)
        (*                               *)
com[a_, a_]                             :=      0

        (*                                              *)
        (* Expansion rules for products:                *)
        (*              [ab,c] = a[b,c] + [a,c]b        *)
        (*              [a,bc] = [a,b]c + b[a,c]        *)
        (*                                              *)
com[prod[a_,b__],c_]    :=      
                prod[a,  com[prod[b], c] ] +
                prod[com[a,c],  prod[ b] ]
com[a_,prod[b_,c__]]    :=
               prod[ com[a, b], prod[c] ] +
                prod[ b, com[a, prod[c]] ]
                
        (*                                              *)
        (*      Rules for nested commutators    *)
        (*              [a,[b,c] ] = a[b,c] - [b,c]a    *)
        (*              [ [a,b],c] = [a,b]c - c[a,b]    *)
        (*                                              *)
com[a_,com[b_,c_]]              :=      prod[a, com[b,c]] -
                                                        prod[com[b,c], a]
com[com[a_,b_],c_]              :=      prod[com[a,b], c] -
                                                        prod[c, com[a,b]]
                                                        
(* End of properties of commutator              *)

(****************************************)
(* Properties of operators              *)
(****************************************)

        (*                                      *)
        (* Numbers always commute       *)
        (*                                      *)
op[a_]                                        :=      a       /;NumberQ[a]

        (*                                       *)
        (* No nested operators necessary *)
        (*                                       *)
op[op[a_]]                              :=  op[a]

        (*                                      *)
        (* Product of ops. is also on op.       *)
        (*                                      *)
op[prod[a_,b__]]                :=      prod[ op[a], 
                                                        op[ prod[b] ]  ]
                                                        
        (*                                              *)
        (* Coefficients can be pulled out       *)
        (*                                              *)
op[a_?NumberQ b_]               :=      a op[b]

        (*                                                      *)
        (* Commutator is already an op.         *)
        (*                                                      *)
op[com[a_,b_]]                  :=      com[a,b]
(****************************************)
(* Properties of prod, (non-comm. mult.)*)
(****************************************)

        (*                      *)
        (* Associativity        *)
        (*                      *)
Attributes[prod]                =       {Flat}
        (*                                  *)
        (* Can pull numbers outside *)
        (*                                  *)
prod[a__,b_?NumberQ,c___]       :=      b prod[a,c] 
prod[b_?NumberQ,a__]            :=      b prod[a]

        (*                                    *)
        (* Single product is just itself *)
        (*                                    *)
prod[a_]                                        :=      a

        (*                                      *)
        (* Rearrangement identities:    *)
        (*      abc... = bac... + [a,b]c...     *)
	(*      dab        = dba + d[a,b]       *)
        (*      dabc...= dbac...+ d[a,b]c...    *)
        (*                                      *)
prod[op[a_],op[b_],c___]        :=
                prod[op[b],op[a],c] + 
                prod[ com[ op[a], op[b] ], c] /;
                        !OrderedQ[{a,b}]

prod[d_,op[a_],op[b_]]  :=
                prod[d,op[b],op[a]] + 
                prod[d, com[ op[a], op[b] ]] /;
                        !OrderedQ[{a,b}]
                        
prod[d_,op[a_],op[b_],c__]      :=
                prod[d,op[b],op[a],c] + 
                prod[d, com[ op[a], op[b] ], c] /;
                        !OrderedQ[{a,b}]

(******************************************************)
(* Printing out the results in readable form             *)
(*              com[a,b]                ->      "[ a , b ]"      *)
(*              op[a]                   ->      a                *)
(*            prod[a,b,...]   ->      a b ...                  *)
(*                                                               *)
(******************************************************)

Format[com[a_,b_] ]                     :=      StringForm[
                        "[ `` , `` ]",a,b ]
Format[op[a_]]                          :=      a
Format[prod[a_,b__]]            :=      a prod[b]

(* Possible test cases *)
(* com[ op[b], op[ a] ] *)
(* com[ prod[ op[a], op[b]], op[c] ]   *)
(* com[ prod[ op[b],op[c],op[a] ], op[d] ] *)
(* prod[ com[ op[a], op[b] ], com[ op[d], op[c] ]   ]   *)


