BeginPackage["NormalOrder`"]
Unprotect[NormalOrder]
Off[General::spell, General::spell1]

NormalOrderHelper::usage="";

NormalOrder::usage =
"NormalOrder[x] converts x to a normal ordered product. Operators are
defined as objects for which OperatorQ[x] is true.

Operators are ordered based on their priority which must be
specified by Priority[x]. Higher priorities are ordered to the
right (i.e. alphabetical ordering will work as usual). If the
formatting specifiers BraQ[], KetQ[], CreationOperatorQ[] or
AnnihilationOperatorQ[] are defined, these will specify the
formatting of the operators, otherwise the priority will be used as
follows:
1000 <= Priority are displayed as bras: <x|
0 < Priority < 1000 are displayed as creation operators
-1000 < Priority < 0 are displayed as annihilation operators
Priority <= -1000 are displayed as keta |x>

Operators are permuted by looking for matches of the form:
GradedCommutator[a,b] = c
Commutator[a,b] = d
AntiCommutator[a,b] = e

If a GradedCommutator is defined for a and b, and a and b both have
defined grades, then a.b is replaced by
a.b = c + (-1)^(Grade[a]*Grade[b])b.a
otherwise, if a Commutator is defined,
a.b = d + b.a
or if an AntiCommutator is defined
a.b = e - b.a
If both a Commutator and an AntiCommutator is defined
a.b = (d+e)/2

If a.b appears but only a bracket [b,a] is defined, then the
appropriate parity rules are applied to form [a,b]. If both [a,b]
and [b,a] are defined inconsistently the one matching the term is
used and a warning is displayed.

If none of these is defined, terms will just be returned without
normal ordering them. You could, for example, specify that
everything should commute by default by using a global pattern:
Commutator[_,_] = 0
for example.

Note: to affect simplifications, you can use up-values to define for
example
a.b ^= 3
This will of course simplify wherever the expression a.b occurs at
no cost. A good use for this is to define operators that act as
states: an eigenstate x of operator a could be written:
a.x ^= 5*x

NormalOrder first calls Distribute. The order of operators that is
considered normal depends on the priority of the operator.

NormalOrder is not smart: used for sorting it essentially implements
bubble sort. It is quite general though.";

LNormalOrder::usage =
"LNormalOrder[x] performs the same function as NormalOrder but removes all
terms with creation operators to the far left.";

RNormalOrder::usage =
"RNormalOrder[x] performs the same function as NormalOrder but removes all
terms with annihilation operators to the far right.";

LRNormalOrder::usage =
"RNormalOrder[x] performs the same function both LNormalOrder and
LNormalOrder together.";

NormalSimplify::usage =
"NormalSimplify[x] performs some simplifications on x."

Operator::usage =
"Operator[\"name\",priority,grade,format]
is one possible representation of an operator. This is the
representation used internally by the shortcuts CreationOperator[]
AnnihilationOperator etc. Priority[], Name[], Grade[] and Format[] are
defined for these.

Conjugate[Operator[]] will change the sign of priority and the
format from
creation <> annihilation
bra <> ket

Transpose[Operator[]] will not affect the operator.";

CreationOperator::usage =
"CreationOperator[\"name\",<priority>] Creates a creation operator.
An optional priority can be given (default = 10).";

AnnihilationOperator::usage =
"AnnihilationOperator[\"name\",<priority>] Creates an annihilation
operator. An optional priority can be given (default = -10).";

Ket::usage =
"Ket[\"name\",<priority>] Creates a ket state. An optional priority
can be given (default = -2000).";

Bra::usage =
"Bra[\"name\",<priority>] Creates a bra state. An optional priority
can be given (default = 2000).";

Priority::usage =
"Gives the Priority of an operator to determine how to normal
order it.";

Arguments::usage =
"Gives the arguments or an operator for formatting purposes.";

Commutator::usage =
"Specify commutation relations like Commutator[a,b]=1.";

AntiCommutator::usage =
"Specify anti-commutation relations like AntiCommutator[a,b]=1.";

GradedCommutator::usage =
"Specify graded commutation relations. GradedCommutator[a,b]=1. Note,
this will only be used if both Grade[a] and Grade[b] are defined.";

Grade::usage =
"Grade of operator for use with GradedCommutator.";

BraQ::usage =
"If BraQ[x] is true, then x will be displayed as <x|, <Name[x]| or
<Name[x][Arguments[x]]| depending on whether or not Name[x] and
Arguments[x] are defined. Note: Formatting only works for
Operator[] type expressions because Mathematica is inconsistent and
Format requires a name to which to attach the format.";

KetQ::usage =
"If BraQ[x] is true, then x will be displayed as |x>, |Name[x]> or
|Name[x][Arguments[x]]> depending on whether or not Name[x] and
Arguments[x] are defined. Note: Formatting only works for
Operator[] type expressions because Mathematica is inconsistent and
Format requires a name to which to attach the format.";

CreationOperatorQ::usage =
"If CreationOperatorQ[x] is true, then x will be displayed as
SuperDagger[x], SuperDagger[Name[x]] or
SuperDagger[Name[x]][Arguments[x]] depending on whether or not
Name[x] and Arguments[x] are defined. Note: Formatting only works for
Operator[] type expressions because Mathematica is inconsistent and
Format requires a name to which to attach the format.";

AnnihilationOperatorQ::usage =
"If AnnihilationOperatorQ[x] is true, then x will be displayed as
x, Name[x] or Name[x][Arguments[x]] depending on whether or not
Name[x] and Arguments[x] are defined. Note: Formatting only works for
Operator[] type expressions because Mathematica is inconsistent and
Format requires a name to which to attach the format.";

OperatorQ::usage =
"If this is true, then the object will be assumed to not commute and
NormalOrder[] will apply commutation relations.";

NoOperatorsQ::usage =
"If this is true, then the object is assumed to contain no operators
and will be treated as a scalar by NormalOrder[]: i.e. will commute
with everything and will be multiplied by Times rather than NonCommutativeMultiply. It
should be defined something like
NoOperatorsQ[x_]:=FreeQ[x,_?OperatorQ]]
This is the default definition and should not need to be overloaded.";

Comm::usage="";

KronSimplify::usage =
"KronSimplify looks for simple ocurrances of KroneckerDelta and tries
to eliminate occurances of the specified variable by assuming that
the expression occurs under the appropriate sum.";

Begin["`Private`"]
(* NormalSimplify collects terms and distributes. *)

NoOperatorsQ[x_] := FreeQ[x,_?OperatorQ];

NormalSimplify[Plus[x_,y_]] := NormalSimplify[x]+NormalSimplify[y];

NormalSimplify[Times[a_?NoOperatorsQ,x__]] := a*NormalSimplify[x];

NormalSimplify[NonCommutativeMultiply[x__,a_?NoOperatorsQ,y___]] := a*NormalSimplify[x**y];

NormalSimplify[NonCommutativeMultiply[x___,a_?NoOperatorsQ,y__]] := a*NormalSimplify[x**y];

NormalSimplify[NonCommutativeMultiply[x___,Times[a_?NoOperatorsQ,y__],z___]] :=
a*NormalSimplify[x**y**z];

NormalSimplify[x___] := x;

(* Here is the front end for NormalOrder *)

NormalOrder[x_] := NormalOrderHelper[NormalSimplify[MapAll[Distribute,x]]];

(* Now the basic rules *)
NormalOrderHelper[Plus[x_,y_]] :=
NormalOrderHelper[x]+NormalOrderHelper[y];

NormalOrderHelper[Times[a_?NoOperatorsQ,x__]] :=
a*NormalOrderHelper[x]

(* We need two of these because one of x or y must be defined *)
NormalOrderHelper[NonCommutativeMultiply[x__,a_?NoOperatorsQ,y___]] :=
a*NormalOrderHelper[x**y]

NormalOrderHelper[NonCommutativeMultiply[x___,a_?NoOperatorsQ,y__]] :=
a*NormalOrderHelper[x**y]

NormalOrderHelper[NonCommutativeMultiply[x___,Times[a_?NoOperatorsQ,y_],z___]] :=
a*NormalOrderHelper[x**y**z]

NormalOrderHelper[o1_?OperatorQ] := o1;
NormalOrderHelper[NonCommutativeMultiply[x___,Comm[y__,z__],t___]] :=
NormalSimplify[NormalOrderHelper[x**y]**NormalOrderHelper[z**t]];
NormalOrderHelper[Comm[y__,z__]] := NormalSimplify[y**z];

NormalOrderHelper[Commutator[x___]] := Commutator[x];
NormalOrderHelper[AntiCommutator[x___]] := AntiCommutator[x];
NormalOrderHelper[GradedCommutator[x___]] := GradedCommutator[x];

(* This turns products of grassman variables to zero if they are the *)
(* same. *)
NormalOrderHelper[NonCommutativeMultiply[___,
o1_?OperatorQ,
o1_?OperatorQ,
___]] /; Mod[Grade[o1],2]==1 := 0;

(* This is the main function which does normal ordering. It only works on
terms that are of the form NonCommutativeMultiply[...] but this should be ensured by the
previous rules.
*)
NormalOrderHelper[NonCommutativeMultiply[x___,
o1_?OperatorQ,
o2_?OperatorQ,
y___]] /; Priority[o1] > Priority[o2] :=
NormalOrder[x**CheckErrCommute[o1,o2]**y];

(* This is the function that commutes two variables. It first checks
to see if a GradedCommutator is define. If so, it uses that, otherwise
it checks to see if regular commutators etc. are defined.*)

FastCommute[o1_?OperatorQ,
o2_?OperatorQ] :=
Which[
ValueQ[AntiCommutator[o1,
o2]],
AntiCommutator[o1,o2]-o2**o1,
ValueQ[AntiCommutator[o2,
o1]],
AntiCommutator[o2,o1]-o2**o1,
ValueQ[GradedCommutator[o1,
o2]],
GradedCommutator[o1,
o2]+
(-1)^(Grade[o2]*Grade[o1])*o2**o1,
ValueQ[GradedCommutator[o2,
o1]],
(-1)^(Grade[o2]*Grade[o1])*GradedCommutator[o2,
o1]+
o2**o1,
ValueQ[Commutator[o1,
o2]],
Commutator[o1,o2]+o2**o1,
ValueQ[Commutator[o2,
o1]],
-Commutator[o2,o1]+o2**o1,
True,
Comm[o1,o2]];

NormalOrder::CommErr =
"Commutator[`1`,`2`] + Commutator[`2`,`1`] != 0.";

NormalOrder::AntiCommErr =
"AntiCommutator[`1`,`2`] != AntiCommutator[`2`,`1`].";

NormalOrder::GradedCommErr =
"GradedCommutator[`1`,`2`] != -(-1)^(Parity[`1`]*Parity[`2`])*GradedCommutator[`2`,`1`].";

CheckErrCommute[o1_?OperatorQ,
o2_?OperatorQ] :=
Check[
Module[{
vgc1 = ValueQ[GradedCommutator[o1,
o2]],
vgc2 = ValueQ[GradedCommutator[o2,
o1]],
vg1 = ValueQ[Grade[o1]],
vg2 = ValueQ[Grade[o2]],
vc1 = ValueQ[Commutator[o1,
o2]],
va1 = ValueQ[AntiCommutator[o1,
o2]],
vc2 = ValueQ[Commutator[o2,
o1]],
va2 = ValueQ[AntiCommutator[o2,
o1]],
gc1 = GradedCommutator[o1,
o2],
gc2 = GradedCommutator[o2,
o1],
g1 = Grade[o1],
g2 = Grade[o2],
c1 = Commutator[o1,
o2],
a1 = AntiCommutator[o1,
o2],
c2 = Commutator[o2,
o1],
a2 = AntiCommutator[o2,
o1]},
Module[{vgc = And[Or[vgc1,vgc2],vg1,vg2],
vc = Or[vc1,vc2],
va = Or[va1,va2],
gc = If[vgc1,gc1,-(-1)^(g1*g2)*gc2],
c = If[vc1,c1,-c2],
a = If[va1,a1,a2]},
Which[And[vgc1,vgc2,gc1 != -(-1)^(g1*g2)*gc2],
(* Graded Commutation Error *)
Message[NormalOrder::GradedCommErr,
ToString[Format[o1]],
ToString[Format[o2]]],
vgc,
(* Use GradedCommutator *)
gc+(-1)^(g1*g2)*o2**o1,
And[vc1,vc2,c1 != -c2],
(* Commutation Error *)
Message[NormalOrder::CommErr,
ToString[Format[o1]],
ToString[Format[o2]]],
And[va1,va2,a1 != a2],
(* Anticommutation Error *)
Message[NormalOrder::AntiCommErr,
ToString[Format[o1]],
ToString[Format[o2]]],
And[va,vc],
(* Both commutation and anticommutation are defined *)
1/2*(a+c),
vc,
(* Commutation defined *)
c+o2**o1,
va,
(* AntiCommutation defined *)
a-o2**o1,
True,
(* No commutation relation is defined, so return unCommutated. *)
Comm[o1,o2]]]],
err[o1,o2]];

(* Anything not matched by the previous operators is considered normal
ordered. *)

NormalOrderHelper[x_] := x;

(******************************************************)
(* These are here for convienience: they define the Operator[] type *)

CreationOperator[n_] = Operator[n,-10,0,"creation"];
CreationOperator[n_,p_] = Operator[n,p,"creation"];

AnnihilationOperator[n_] = Operator[n,10,0,"annihilation"];
AnnihilationOperator[n_,p_] = Operator[n,p,0,"annihilation"];

Bra[n_] = Operator[n,2000,0,"bra"];
Bra[n_,p_] = Operator[n,p,"bra"];

Ket[n_] = Operator[n,-2000,0,"ket"];
Ket[n_,p_] = Operator[n,p,"ket"];

AnnihilationOperator[n_] = Operator[n,10,0,"annihilation"];
AnnihilationOperator[n_,p_] = Operator[n,p,0,"annihilation"];

CreationOperatorQ[Operator[_,_,_,"creation",___] |
Operator[_,_,_,"creation",___][___]] := True;
CreationOperatorQ[Operator[___] | Operator[___][___]] := False;

AnnihilationOperatorQ[Operator[_,_,_,"annihilation",___] |
Operator[_,_,_,"annihilation",___][___]] := True;
AnnihilationOperatorQ[Operator[___] | Operator[___][___]] := False;

BraQ[Operator[_,_,_,"bra",___] |
Operator[_,_,_,"bra",___][___]] := True;
BraQ[Operator[___] |
Operator[___][___]] := False;

KetQ[Operator[_,_,_,"ket",___] |
Operator[_,_,_,"ket",___][___]] := True;
KetQ[Operator[___] | Operator[___][___]] := False;

OperatorQ[Operator[___] | Operator[___][___]] := True;

PrivateCreationOperatorQ[Operator[_,_,_,"creation",___] |
Operator[_,_,_,"creation",___][___]] := True;

PrivateAnnihilationOperatorQ[Operator[_,_,_,"annihilation",___] |
Operator[_,_,_,"annihilation",___][___]] := True;

PrivateBraQ[Operator[_,_,_,"bra",___] |
Operator[_,_,_,"bra",___][___]] := True;

PrivateKetQ[Operator[_,_,_,"ket",___] |
Operator[_,_,_,"ket",___][___]] := True;

PrivateOperatorQ[Operator[___] | Operator[___][___]] := True;

Priority[Operator[_,p_,___] |
Operator[_,p_,___][___]] := p;

Name[Operator[n_,___] |
Operator[n_,___][___]] := n;

Arguments[Operator[___][a___]] := a;

(* The formatting needs some work for creationoperators (type of *)
(* brakets displayed) and for bra's and kets. *)
CList[x_] := x
CList[x_, y__] := SequenceForm[x, ",", CList[y]]
Format[Operator[o___]] :=
Which[CreationOperatorQ[Operator[o]],
If[ValueQ[Name[Operator[o]]],
If[ValueQ[Arguments[Operator[o]]],
SuperDagger[Name[Operator[o]]][Arguments[Operator[o]]],
SuperDagger[Name[Operator[o]]]],
Operator[o]],
AnnihilationOperatorQ[Operator[o]],
If[ValueQ[Name[Operator[o]]],
If[ValueQ[Arguments[Operator[o]]],
Name[Operator[o]][Arguments[Operator[o]]],
Name[Operator[o]]],
Operator[o]],
BraQ[Operator[o]],
If[ValueQ[Name[Operator[o]]],
If[ValueQ[Arguments[Operator[o]]],
SequenceForm["<", Name[Operator[o]],":",CList[Arguments[Operator[o]]], "|"],
SequenceForm["<", Name[Operator[o]], "|"]],
SequenceForm["<", Operator[o], "|"]],
KetQ[Operator[o]],
If[ValueQ[Name[Operator[o]]],
If[ValueQ[Arguments[Operator[o]]],
SequenceForm["|", Name[Operator[o]],":",CList[Arguments[Operator[o]]], ">"],
SequenceForm["|", Name[Operator[o]], ">"]],
SequenceForm["|", o, ">"]],
True, SequenceForm["Operator"[o]]]

Format[Operator[o___][a___]] :=
Which[CreationOperatorQ[Operator[o][a]],
If[ValueQ[Name[Operator[o][a]]],
If[ValueQ[Arguments[Operator[o][a]]],
SuperDagger[Name[Operator[o][a]]][Arguments[Operator[o][a]]],
SuperDagger[Name[Operator[o][a]]]],
Operator[o][a]],
AnnihilationOperatorQ[Operator[o][a]],
If[ValueQ[Name[Operator[o][a]]],
If[ValueQ[Arguments[Operator[o][a]]],
Name[Operator[o][a]][Arguments[Operator[o][a]]],
Name[Operator[o][a]]],
Operator[o][a]],
BraQ[Operator[o][a]],
If[ValueQ[Name[Operator[o][a]]],
If[ValueQ[Arguments[Operator[o][a]]],
SequenceForm["<", Name[Operator[o][a]],":",CList[Arguments[Operator[o][a]]], "|"],
SequenceForm["<", Name[Operator[o][a]], "|"]],
SequenceForm["<", Operator[o][a], "|"]],
KetQ[Operator[o][a]],
If[ValueQ[Name[Operator[o][a]]],
If[ValueQ[Arguments[Operator[o][a]]],
SequenceForm["|", Name[Operator[o][a]],":",CList[Arguments[Operator[o][a]]], ">"],
SequenceForm["|", Name[Operator[o][a]], ">"]],
SequenceForm["|", o, ">"]],
True, SequenceForm["Operator"[o][a]]]

(* Conjugation up-value *)
Conjugate[Operator[n_,p_,g_,"creation",z___]] ^= Operator[n,-p,g,"annihilation",z];
Conjugate[Operator[n_,p_,g_,"annihilation",z___]] ^= Operator[n,-p,g,"creation",z];
Conjugate[Operator[n_,p_,g_,"bra",z___]] ^= Operator[n,-p,g,"ket",z];
Conjugate[Operator[n_,p_,g_,"ket",z___]] ^= Operator[n,-p,g,"bra",z];
Conjugate[Operator[n_,p_,g_,"creation",z___][s___]] ^= Operator[n,-p,g,"annihilation",z][s];
Conjugate[Operator[n_,p_,g_,"annihilation",z___][s___]] ^= Operator[n,-p,g,"creation",z][s];
Conjugate[Operator[n_,p_,g_,"bra",z___][s___]] ^= Operator[n,-p,g,"ket",z][s];
Conjugate[Operator[n_,p_,g_,"ket",z___][s___]] ^= Operator[n,-p,g,"bra",z][s];
Transpose[Operator[z___]] ^= Operator[z];
Transpose[Operator[z___][s___]] ^= Operator[z][s];

(* To improve performance, undefined these *)

NonCommutativeMultiply[Operator[_,_,_,"annihilation",___],Operator[_,_,_,"ket",___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"bra",___],Operator[_,_,_,"creation",___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"annihilation",___],Operator[_,_,_,"ket",___][___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"bra",___],Operator[_,_,_,"creation",___][___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"annihilation",___][___],Operator[_,_,_,"ket",___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"bra",___][___],Operator[_,_,_,"creation",___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"annihilation",___][___],Operator[_,_,_,"ket",___][___]] ^= 0;
NonCommutativeMultiply[Operator[_,_,_,"bra",___][___],Operator[_,_,_,"creation",___][___]] ^= 0;

(* This might speed things up a bit. *)
NormalOrderHelper[NonCommutativeMultiply[x___,
Operator[n1_,p1_,o1___],
Operator[n2_,p2_,o2___],
y___]] /; p1 > p2 :=
NormalOrder[x**Commute[Operator[n1,p1,o1],Operator[n2,p2,o2]]**y];

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[_?PrivateAnnihilationOperatorQ,_?PrivateKetQ] = 0;
NonCommutativeMultiply[_?PrivateBraQ,_?PrivateCreationOperatorQ] = 0;
Protect[NonCommutativeMultiply];

KronSimplify[ex_,{}]:=ex
KronSimplify[ex_,{q1_,q2___}]:=KronSimplify[KronSimplify[ex,q1],{q2}]
KronSimplify[Times[KroneckerDelta[a_,b_],h___],q_]:=
Block[{sols=Solve[a==b,q][[1]]},
If[Length[sols]==1,
KronSimplify[Times[h]/.sols,q],
Times[KroneckerDelta[a,b],h]]]
KronSimplify[Plus[a_,b___],q_]:=Plus[KronSimplify[a,q],KronSimplify[Plus[b],q]]
KronSimplify[h_,_]:=h
End[]
EndPackage[];