(* ::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.                                                         *)
(************************************************************************)



Unprotect[NonCommutativeMultiply];
x_**0=0; 0**y_=0;
(c_*x_AW) ** y_ := Expand[c(x**y)];
x_ ** (c_*y_AW) := Expand[c(x**y)];
x_Plus ** y_ := (#**y)& /@ x;
x_**y_Plus := (x**#)& /@ y;
AW[w1_String] ** AW[w2_String] := AW[w1<>w2];


LyndonQ[AW[w_String]] := And @@ (
OrderedQ[{w, #}]& /@ Table[StringDrop[w,i], {i,1,StringLength[w]-1}]
);
AllWords[0, _List] = {AW[""]};
AllWords[n_ /; n>0, ab_List] := AllWords[n, ab] = AW /@ Flatten[Outer[
StringJoin[#1, #2]&,
First /@ AllWords[n-1, ab],
ab
]];
AllLyndonWords[n_Integer, ab_List] := LW @@@ Select[AllWords[n, ab], LyndonQ];
AllLyndonWords[{n_}, ab_List] := Join@@Table[AllLyndonWords[k,ab], {k,n}];
LyndonFactorization[LW[w_String] /; StringLength[w] == 1] := LW[w];
LyndonFactorization[LW[w_String] /; StringLength[w] > 1] := Module[
{rf},
rf=First[Sort[Table[StringDrop[w,i], {i,1,StringLength[w]-1}]]];
LW /@ {StringDrop[w, -StringLength[rf]], rf}
];
LW[s_Symbol] := LW[ToString[s]];
LW[LW[w_]]:=LW[w];
LW /: LW[x_] <= LW[y_] := OrderedQ[{x,y}];
LW /: x_LW >=y_LW := y<=x;
LW /: x_LW > y_LW := !(x<=y);
LW /: x_LW<y_LW := !(y<=x);
Format[LW[w_], StandardForm] := Defer[\[LeftAngleBracket]w\[RightAngleBracket]];
\[LeftAngleBracket]w__\[RightAngleBracket]:=LW[w];
LW[is__Integer] := LW[StringJoin@@(StringTake["1234567890abcdefghijklmnopqrstuvwxyz", {#}]& /@ {is})];
Deg[LW[x_]] := StringLength[x];


b[0,_] =0; b[_,0]=0;
b[c_*(x_AW|x_LW), y_] := Expand[c b[x,y]];
b[x_, c_*(y_AW|y_LW)] := Expand[c b[x,y]];
b[x_Plus, y_] := b[#,y]& /@ x;
b[x_, y_Plus] := b[x,#]& /@ y;
b[w_AW, z_AW] := w**z-z**w;
b[w_LW, z_LW] := LWBracket[w,z];


LWBracket[w_LW, z_LW] := Which[
(* If[Deg[w]+Deg[z]>4, Dialog[]]; *)
w===z, 0,
z<w, Expand[-b[z,w]],
Deg[w]==1, LW[First[w]<>First[z]],
True, Module[{x,y},
{x,y}=LyndonFactorization[w];
If[y>=z,
LW[First[w]<>First[z]],
LWBracket[w,z]=b[x,LWBracket[y,z]]+b[LWBracket[x,z],y]
]
]
];


ad[x_][y_]:=b[x,y];
MakeLieSeries[d_, l_List] := MakeLieSeries[d, #]& /@ l;
MakeLieSeries[d_, s_LieSeries] /; Length[s]<=d := s;
MakeLieSeries[d_, s_LieSeries] /; Length[s]>d := Take[s,d];
MakeLieSeries[d_, a_ -> b_] := (a -> MakeLieSeries[d, b]);
MakeLieSeries[d_, expr_] := LieSeries @@ Table[
expr /. w_LW /; Deg[w] != k -> 0,
{k,d}
];
MakeLieSeries[w_LW] := Append[LieSeries@@Table[0, {Deg[w]-1}], w];
LieSeries[]=0;
LieSeries /: Expand[s_LieSeries] := Expand /@ s;
LieSeries /: Plus[ss__LieSeries] /; Length[{ss}]>1 := Module[
{l=Min[Length /@ {ss}]},
LieSeries @@ Total[Take[List @@ #, l]& /@ {ss}]
];
LieSeries /: 0*s_LieSeries := 0;
LieSeries /: c_?NumberQ * s_LieSeries := Expand[c*#]& /@ s;
LieSeries /: s_LieSeries + c_.*w_LW /; NumberQ[c] := Module[{d},
d=Deg[w];
If[Length[s]<d, 0,
ReplacePart[s, d -> s[[d]]+Expand[c*w]]
]
];
b[s1_LieSeries, s2_LieSeries] := Module[
{d,k,m1,m2},
m1=1+LengthWhile[s1, #==0&];
m2=1+LengthWhile[s2, #==0&];
LieSeries @@ Table[
Sum[b[s1[[k]], s2[[d-k]]], {k,m1,d-m2}],
{d, Min[m1+Length[s2], m2+Length[s1]]}
]
];
b[w_LW, s_LieSeries] := Join[
LieSeries@@Table[0, {Deg[w]}],
ad[w] /@ s
];
b[s_LieSeries, w_LW] := Expand[-b[w,s]];
LieSeries /: EulerE[s_LieSeries] := LieSeries @@ Expand[Range[Length[s]]*(List@@s)];
adSeries[f_, x_, d_][\[Psi]_] := Module[
{ser,as,ni,nf,t,l,m},
ser=List@@Series[f, {ad,0, d}];
{as, ni, nf}=ser[[{3,4,5}]];
t=Nest[ad[x], \[Psi], ni];
If[Head[t]===LieSeries,
l=Length[t];
Expand[First[as]*t+Sum[
m=1+LengthWhile[t, #==0&];
t=ad[MakeLieSeries[l-m,x]][t];
as[[k+1]]*t,
{k, 1, nf-ni-1}
]],
Expand[as.NestList[ad[x],t,nf-ni-1]]
]
];
Ad[s_LieSeries] := adSeries[E^(-ad), s, Length[s]-1];


BCH[1]= LieSeries[\[LeftAngleBracket]"x"\[RightAngleBracket]+\[LeftAngleBracket]"y"\[RightAngleBracket]];
BCH[n_]:= BCH[n] = Module[
{bch,t1,t2},
bch=Append[BCH[n-1],0];
t1=MakeLieSeries[n, \[LeftAngleBracket]"y"\[RightAngleBracket]+adSeries[E^(-ad), \[LeftAngleBracket]"y"\[RightAngleBracket], n-1][\[LeftAngleBracket]"x"\[RightAngleBracket]]];
t2=adSeries[(1-E^(-ad))/ad, bch,n-1][EulerE[bch]];
bch+(t1-t2)/n
]


LieMorphism[rules_List] := LieMorphism[rules] = LieMorphism[Unique[LieMorphism], rules];
LieMorphism[mor_Symbol, rules_List] := (
mor[]=rules;
(mor[w_LW] /; Deg[w] == 1) := (mor[w] = w /. rules);
mor[w_LW] := (mor[w] = b @@ (mor /@ LyndonFactorization[w]));
mor[s_LieSeries] := Module[
{l=Length[s]},
MakeLieSeries[l, Sum[
LieMorphism[MakeLieSeries[l-k+1, rules]][s[[k]]],
{k,l}
]]
];
mor[expr_] := Expand[expr /. {s_LieSeries :> mor[s], w_LW :> mor[w]}];
mor
);
BCH[n_, x_, y_] := LieMorphism[{LW["x"] -> x, LW["y"] -> y}][BCH[n]];
BCH[s1_LieSeries, s2_LieSeries] := BCH[Min[Length /@ {s1,s2}], s1,s2];


\[Iota][w_LW] /; Deg[w]==1 := AW@@w;
\[Iota][w_LW] := \[Iota][w] = b @@ (\[Iota] /@ LyndonFactorization[w]);
\[Iota][s_LieSeries] := Prepend[\[Iota]/@(ASeries @@ s), 0];
\[Iota][expr_] := Expand[expr /. w_LW :> \[Iota][w]];


ASeries /: Expand[s_ASeries] := Expand /@ s;
ASeries /: Plus[ss__ASeries] := Module[
{l=Min[Length /@ {ss}]},
ASeries @@ Total[Take[List @@ #, l]& /@ {ss}]
];
ASeries /: c_ * s_ASeries := Expand[c*#]& /@ s;
s1_ASeries**s2_ASeries := Module[
{d,k,m1,m2},
m1=LengthWhile[s1, #==0&];
m2=LengthWhile[s2, #==0&];
ASeries @@ Table[
Sum[s1[[k+1]]**s2[[d-k+1]], {k,m1,d-m2}],
{d, 0, Min[m1+Length[s2]-1, m2+Length[s1]-1]}
]
];
ASeries /: EulerE[s_ASeries] := ASeries @@ Expand[Range[{0,1+Length[s]}]*(List@@s)];


\[Sigma][y_LW, w_LW] /; Deg[y]==1 := \[Sigma][y,w] = Which[
y===w, AW[""],
Deg[w]===1, 0,
True, Module[{w1,w2},
{w1,w2}=LyndonFactorization[w];
\[Iota][w1]**\[Sigma][y,w2]-\[Iota][w2]**\[Sigma][y,w1]
]
];
\[Sigma][y_, expr_] := Expand[expr /. w_LW :> \[Sigma][LW[y], w]] /. LieSeries -> ASeries;



