(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



AH[h_AH] := h;
AH[0]=0;
AH /: AH[h1_]*AH[h2_] := AH[h1*h2];
AH /: AH[1]*a_Ar := a;
AH /: AH[h1_]*Y[ijk__, AH[h2_]] := Y[ijk, AH[h1*h2]];
AH /: c_*AH[h_] /; FreeQ[c, AH|PH]:= AH[c*h];
AH /: AH[h1_]+AH[h2_] := CanonicalForm[AH[h1+h2]];
AH /: AH[h_]~Mod~x[i_] := AH[Limit[h, x[i] -> 0]];
AH /: Coefficient[AH[h_], eps_, 1] := AH[Coefficient[Expand[h], eps, 1]];
ToAH[expr_] := expr /. PH[h_] :> AH[Normal[h] /. z->1] /. Y[_,_,_, 0] -> 0;


ToPH[n_, expr_] := expr /. {
AH[h_] :> PH[(h /. x[i_] :> z x[i])+O[z]^n],
PH[h_] :> PH[h + O[z]^n]
};
PH[0]=0;
PH /: PH[h1_]*PH[h2_] := CanonicalForm[PH[h1*h2]];
PH /: PH[Literal[SeriesData[z,0,{1,0...},0,__]]]*a_Ar := a;
PH /: PH[h1_]*Y[ijk__, h2_] := Y[ijk, PH[h1]*h2];
PH /: AH[h1_]*PH[h2_] := CanonicalForm[PH[(h1 /. x[i_] :> z*x[i])*h2]];
PH /: c_*PH[h_] /; FreeQ[c, _Ar]:= CanonicalForm[PH[(c /. x[i_] :> z*x[i])*h]];
PH /: PH[h1_]+PH[h2_] := PH[h1+h2];
PH /: AH[h1_]+PH[h2_] := PH[(h1 /. x[i_] :> z*x[i])+h2];
PH /: PH[h_]~Mod~x[i_] := PH[h /. x[i] -> 0];
PH /: Coefficient[h_PH, eps_, 1] := CanonicalForm[
MapAt[Coefficient[#, eps, 1]&, h, {1,3}]
];
DeclareSeries[f_[vars___], deg_] := (
PH[f] = CanonicalForm[PH[
(f[vars] /. Thread[{vars} -> z*{vars}])+O[z]^deg
]];
PH[f] = PH[f] /. Derivative[ds___][f][0...] :> f[ds];
Cases[PH[f], f[___], Infinity]
);


HMatrixExp[mat_] := Module[
{mat1, target},
If[FreeQ[mat, PH], (
mat1 = mat /. AH[h_] :> h;
CanonicalForm[Map[AH, MatrixExp[mat1], {2}]]
), (
target = Min[Cases[
mat,
Literal[SeriesData[z,0,_, _, max_, _]] :> max,
Infinity
]];
mat1 = Map[ToPH[target, #]&, mat, {2}] /. PH[h_] :> h;
CanonicalForm[Map[PH,
Sum[MatrixPower[mat1, k]/k!, {k,0,target}],
{2}] /. {
PH[1] -> AH[1],
PH[O[z]^target] -> 0
}
]
)]
];
HSolve[eqns_List, etc___] := If[
FreeQ[eqns, PH],
Solve[eqns /. AH[h_] :> h, etc],
Module[
{vars, target},
vars = Union[Cases[eqns, x[i_], Infinity]];
If[vars==={}, Solve[eqns, etc],
Solve[
(
target = Min[Cases[
#,
Literal[SeriesData[z,0,_, _, max_, _]] :> max,
Infinity
]];
Series[Normal[ToPH[target, #] /. PH[h_] :> h] /. z-> 1, Sequence @@ ({#, 0, target}& /@ vars)]
)& /@ eqns,
etc
]
]
]
];
HSolve[eqn_Equal, etc___] := HSolve[{eqn}, etc];


CanonicalForm[expr_] := expr /. {
AH[h_] :> AH[Factor[h]],
PH[h_SeriesData] :> PH[MapAt[Expand, h, 3]]
};
CanonicalForm[0] = 0;


ReducePrimitives[prims_] := Module[{l, h0, h1}, prims 
//. {
(* Anti-symmetry, linearity *)
Y[i_,i_,__] -> 0,
Y[i_,j_,k_,h_] /; i>j :> Y[j,i,k, -h],
c_*Y[i_,j_,k_, h_] :> Y[i,j,k, c*h],
Y[i_,j_,k_,h1_]+Y[i_,j_,k_,h2_] :> Y[i,j,k,h1+h2],
(* IHX *)
Y[i_,j_,k_, h_] /; !FreeQ[h, x[l_] /; l<i]:> (
l=Min[Cases[{h}, x[l_]:> l, Infinity]];
h0 = h~Mod~x[l];
h1=(h-h0)/x[l];
Y[i,j,k, h0]
-Y[j,l,k,h1* x[i]]-Y[l,i,k,h1*x[j]]
)
} 
/. Y[__, 0] -> 0
];


(prims_ // S[srules___Rule]) := ReducePrimitives[prims
/. {
Ar[i_,j_] :> Distribute[Ar[Ar[i,0], Ar[0,j]] /. {srules}],
Y[i_,j_,k_,h_] :> Distribute[Y[
Ar[i,0] , Ar[j,0], Ar[0,k], h
] /. {srules}]
}
/. {Ar[i_,0] :> i, Ar[0,j_] :> j}
/. {
Ar[Y[i_,j_,0,h_], k_] :> Y[i,j,k,h],
Y[_Y, _Y, __] :> 0,
Y[i_Integer, Y[j_,k_,0,h_], l_, h1_] :> Y[j,k,l, AH[x[i]]*h*h1],
Y[Y[j_,k_,0,h_], i_Integer, l_, h1_] :> Y[j,k,l, AH[-x[i]]*h*h1]
}
/. {
Ar[i_, Y[0,j_,k_,h_]] :> Y[i,j,k,h],
Ar[i_, Y[j_,0,k_,h_]] :> Y[j,i,k,h],
Y[i_,j_,Y[0,k_,l_,h_], h1_] :> Y[i,j,l, AH[-x[k]]*h*h1],
Y[i_,j_,Y[k_, 0,l_,h_], h1_] :> Y[i,j,l, AH[x[k]]*h*h1]
}
];

(prims_ // ts_TS) := prims // (S @@ ts);

S /: S[sr1___Rule]**S[sr2___Rule] := Module[{elements, rule},
elements = Union[First /@ {sr1, sr2}];
S @@ DeleteCases[ 
(* The introduction of the symbol "rule" is an ugly hack *)
Thread[rule[elements, elements // S[sr1] // S[sr2]]],
rule[e_,e_]
] /. rule -> Rule
]

S[sigma[i_, j_], more___] := S[
Ar[0, j] -> Ar[0,j] + Y[0,i,j, AH[-(Exp[-x[i]]-1)/x[i]]],
Ar[0, i] -> Ar[0,i] + Y[0,i,j, AH[(Exp[-x[i]]-1)/x[i]]],
Ar[j, 0] -> Ar[j,0] + Y[i,j,0, AH[(Exp[x[i]]-1)/x[i]]]
] ** S[more];
S[sigbar[i_, j_], more___] := S[
Ar[0, j] -> Ar[0,j] + Y[0,i,j, AH[-(Exp[x[i]]-1)/x[i]]],
Ar[0, i] -> Ar[0,i] + Y[0,i,j, AH[(Exp[x[i]]-1)/x[i]]],
Ar[j, 0] -> Ar[j,0] + Y[i,j,0, AH[(Exp[-x[i]]-1)/x[i]]]
] ** S[more];


SnG[]:= SnG[S[], 0];
SnG /: SnG[s1_S, g1_] ** SnG[s2_S, g2_] := SnG[s1**s2, ReducePrimitives[(g1 // s2)+g2]];

SnG[sigma[i_, j_], more___] := SnG[S[sigma[i,j]], Ar[i,j]] ** SnG[more];
SnG[sigbar[i_, j_], more___] := SnG[S[sigbar[i,j]], -Ar[i,j]] ** SnG[more];


Der[drules___Rule][expr_] := Module[{s,eps},
s = (S[drules] /. (a_ -> b_) :> ReducePrimitives[a -> Expand[a+eps*b]]);
(expr // s) /. {
_Ar -> 0,
Y[i_,j_,k_,h_] :> Y[i,j,k,Coefficient[h, eps, 1]]
} /. Y[__,0] -> 0
];

Der[a_Plus] := Der /@ a;
Der /: a_*Der[drules___Rule] := Der[drules] /. (b_->c_) :> (b->a*c);
Der /: Der[dr1___Rule] + Der[dr2___Rule] := Module[{elements},
elements = Union[First /@ {dr1, dr2}];
Der @@ DeleteCases[
Thread[Rule[elements, ReducePrimitives[
(elements /. {dr1, _Ar -> 0}) + (elements /. {dr2, _Ar -> 0})
]]],
_ -> 0
]
];

Der[a_.*Ar[i_, i_]] := Der[];
Der[a_.*Ar[i_, j_]] /; i!=j := Der[
Ar[j,0] -> Y[i,j,0,AH[-a]],
Ar[0, i] -> Y[0,i,j,AH[a]],
Ar[0, j] -> Y[i,0,j,AH[a]]
];
Der[a_.*Y[i_,j_,k_, h_]] /; i!=j!=k := Module[{d1, d2, elements},
d1=Der[Ar[i,k]];
d2=Der[Ar[j,k]];
elements = Union[First /@ List @@ d1, First /@ List @@ d2];
DeleteCases[
Der @@ Thread[elements -> ReducePrimitives[d1@d2@elements - d2@d1@elements]] /.
Y[ijk__, h1_] :> Y[ijk, a*h*h1],
_ -> 0
]
];
Der[_.*Y[i_,i_,__]] = Der[];
Der[a_.*Y[i_,j_,j_, h_]] /; i!=j := Der[
Ar[j,0] -> Y[i,j,0,x[j]*a*h],
Ar[0, i] -> Y[i, 0,j,x[j]*a*h],
Ar[0, j] -> Y[i, 0,j,-x[j]*a*h]
];
Der[a_.*Y[j_,i_,j_, h_]] /; i!=j := Der[
Ar[j,0] -> Y[i,j,0,-x[j]*a*h],
Ar[0, i] -> Y[i, 0,j,-x[j]*a*h],
Ar[0, j] -> Y[i, 0,j,x[j]*a*h]
];


Der /: Exp[Der[], _] = {};Der /: Exp[Der[drules___Rule], on_] := Module[
{k0, ins, outs, k, newout, hs, zero, e, mat, expmat},
If[Head[on] === List,
k0=Length[ins=on],
ins = First /@ {drules};
If[on === Tails, ins=Cases[ins, Ar[_,0]]];
k0 = Length[ins]
];
outs={};
For[k=1, k<=Length[ins], ++k,
AppendTo[outs, newout=Der[drules][ins[[k]]]];
ins = ins ~Join~ Complement[
Union[Cases[{newout}, Y[ijk__, _] :> Y[ijk, AH[1]], Infinity]],
ins
]
];
--k;
zero=Table[0,{k}];
e[{{i_}}] := ReplacePart[zero, 1, i];
mat = Replace[
outs /. Y[ijk__, h_] :> -h e[Position[ins, Y[ijk,AH[1]]]],
0 -> zero,
{1}
];
expmat = HMatrixExp[mat];
Sort[Thread[
Take[ins, k0] -> Take[expmat.ins, k0]
]]
];


S[Exp[d_Der]] := S @@ Exp[d, All];
S[Exp[prims_]] := S[Exp[Der[prims]]];


TS[Exp[d_Der]] := TS @@ Exp[d, Tails];
TS[Exp[prims_]] := TS[Exp[Der[prims]]];


Y /: Coefficient[expr_, Y[i_,j_,k_]] := expr /. {
_Ar -> 0,
Y[i,j,k, h_] :> h,
_Y -> 0
}
