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



(* ::Input::Initialization:: *)
BeginPackage["FreeLie`"];
Print["FreeLie` implements / extends ",
Sort@{"*","+","**","$SeriesShowDegree","\[LeftAngleBracket]\[RightAngleBracket]","\[Integral]","\[Congruent]",ad,Ad,adSeries,AllCyclicWords,AllLyndonWords,AllWords,Arbitrator,ASeries,AW,b,BCH,BooleanSequence,BracketForm,BS,CC,Crop,cw,CW,CWS,CWSeries,D,Deg,DegreeScale,DerivationSeries,div,DK,DKS,DKSeries,EulerE,Exp,Inverse,j,J,JA,LieDerivation,LieMorphism,LieSeries,LS,LW,LyndonFactorization,Morphism,New,RandomCWSeries,Randomizer,RandomLieSeries,RC,SeriesSolve,Support,t,tb,TopBracketForm,tr,UndeterminedCoefficients,\[Alpha]Map, \[CapitalGamma],\[Iota],\[CapitalLambda],\[Sigma],\[HBar],\[OverBracket],\[OverParenthesis]},
"."];
Print["FreeLie` is in the public domain. Dror Bar-Natan is committed to support it within reason until July 15, 2022. This is version 150814."];
$SeriesShowDegree=3;
Begin["`Private`"];


(* ::Input::Initialization:: *)
SetAttributes[New, HoldAll];
(*New[type_[srn_], def_] := Module[   (* "srn" for "self-reference name" *)
{un}, (* "un" for "unique name" *)
un=Unique[ToString[type]<>"$"];
ReleaseHold[Hold[def] /. srn \[Rule] un];
type[un]
];*)
New[type_[srn_], def_] := (   (* "srn" for "self-reference name" *)
(*Print["new object of type ", type//ToString, " and srn ", srn];*)
ReleaseHold[Hold[def; type[srn]] /. srn -> Unique[ToString[type]<>"$"]]
);


(* ::Input::Initialization:: *)
AllWords[n_, ab_List] := AW @@@ Tuples[ab, n];
LW /: LW[] < LW[__] = True;
LW /: _LW < LW[] = False;
LW /: LW[x_, xs___] < LW[x_, ys___] := LW[xs] < LW[ys];
LW /: LW[x_, ___] < LW[y_, ___] /; (x=!=y) := OrderedQ[{x,y}];
LW /: x_LW >y_LW := y<x;
LW /: x_LW <= y_LW := !(y<x);
LW /: x_LW >= y_LW := !(x<y);
LyndonQ[w_AW] := And @@ (
(LW@@w<LW@@#)& /@ Table[Drop[w,i], {i,1,Length[w]-1}]);
AllLyndonWords[n_Integer, ab_List] := AllLyndonWords[n, ab] =
LW @@@ Select[AllWords[n, ab /. LW[w_]:>w], LyndonQ];
AllLyndonWords[{n_}, ab_List] := Join@@Table[AllLyndonWords[k,ab], {k,n}];
Deg[w_LW] := Length[w];
LyndonFactorization[w_LW /; Deg[w] == 1] := w;
LyndonFactorization[w_LW /; Deg[w] > 1] := Module[{rf},
rf=First[Sort[Table[Drop[w,i], {i,1,Length[w]-1}], Less]];
LW /@ {Drop[w, -Length[rf]], rf}];
LW[w_LW]:=w;
BracketForm[LW[c_]] := ToString[c];
BracketForm[w_LW] := BracketForm[w] = StringJoin[Flatten[{
"[", BracketForm /@ LyndonFactorization[w],"]"
}]];
BracketForm[expr_] := expr /. w_LW :> BracketForm[w];
topbracketform[LW[c_]] := c;
topbracketform[w_LW] := topbracketform[w] = Overscript[
Row[Riffle[topbracketform /@ LyndonFactorization[w], "\[ThinSpace]"]],\[OverBracket]];
TopBracketForm[LW[c_]] := Overscript[c, \[OverBracket]];
TopBracketForm[w_LW] := topbracketform[w];
TopBracketForm[w_CW] := Overscript[StringJoin@@(ToString /@ w), \[OverParenthesis]];
TopBracketForm[expr_] := expr /. w_LW|w_CW|w_DK :> TopBracketForm[w];
Format[w_LW] := TopBracketForm[w];


(* ::Input::Initialization:: *)
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_LW, z_LW] := LWAdjoint[w][z];
ad[x_][y_]:=b[x,y];


(* ::Input::Initialization:: *)
LWAdjoint[w_] := LWAdjoint[w] = Module[{u},
u=Unique[LWAct];
u[z_] := u[z] = Which[
w===z, 0,
z<w, Expand[-b[z,w]],
Deg[w]==1, Join[w,z],
True, Module[{x,y},
{x,y}=LyndonFactorization[w];
If[y>=z,
Join[w,z],
b[x,LWAdjoint[y][z]]+b[LWAdjoint[x][z],y]
]
]
];
u
];


(* ::Input::Initialization:: *)
BS[v:(True|False)] := BS[v] = New[BooleanSequence[ser], ser[_] = v];
BooleanSequence[bs_Symbol][{dd_Integer}] := TopBracketForm[Append[
BS@@(Table[bs[d], {d, 0, dd}] //. {n_. True, True,rest___} :> {(n+1)True, rest}),
"..."
]];
BooleanSequence[bs_Symbol][e___] := bs[e];
Format[bs_BooleanSequence, StandardForm] := bs[{$SeriesShowDegree}];
BooleanSequence /: And[bss__BooleanSequence] := New[BooleanSequence[bs],
bs[d_Integer] := bs[d] = And @@ ((#[d])& /@ {bss})
];
BooleanSequence /: \[HBar]^k_.*BooleanSequence[bs_Symbol] := New[BooleanSequence[nbs],
nbs[d_Integer] := bs[d-k]
];
BooleanSequence /: TrueQ[bs_BooleanSequence] := New[BooleanSequence[nbs],
nbs[d_Integer] := nbs[d] = TrueQ[bs[d]]
];


(* ::Input::Initialization:: *)
LieSeries[ser_Symbol][{dd_Integer}] := Append[TopBracketForm[LS@@Table[ser[d], {d, dd}]],"..."];
LieSeries[ser_Symbol][e___] := ser[e];
Format[s_LieSeries, StandardForm] := TopBracketForm[s[{$SeriesShowDegree}]];
ShowLieSeries[d_Integer][s_LieSeries] := s[{d}];
LS[s_] := MakeLieSeries[s];
MakeLieSeries[s_LieSeries] := s;
MakeLieSeries[expr_] := MakeLieSeries[expr] = New[LieSeries[ser],
ser[d_Integer] := ser[d] = Expand[expr /. w_LW /; Deg[w] != d -> 0]
];
s1_LieSeries \[Congruent] s2_LieSeries := New[BooleanSequence[bs],
bs[0]=True;
bs[d_Integer] := bs[d-1]&&(s1[d]==s2[d]);
];
Crop[s_LieSeries, d_Integer] := Crop[s,d] = New[LieSeries[ser],
ser[dd_Integer] := If[dd<=d, s[dd], 0]
];
RandomLieSeries[ab_List, opts___Rule] := Module[
{rand = Randomizer /. {opts} /. Randomizer -> (RandomInteger[{-2Deg[#]!, 2Deg[#]!}]/Deg[#]!&)},
New[LieSeries[ser],
ser[d_Integer] := ser[d] = Plus @@ ((rand[#]*#)& /@ AllLyndonWords[d, ab])
]
];


(* ::Input::Initialization:: *)
AddLieSeries[ss___LieSeries] := AddLieSeries[ss] = New[LieSeries[ser],
ser[d_Integer] := ser[d] = Plus @@ ((#[d])& /@ {ss})
];
LieSeries /: Plus[ss__LieSeries] := AddLieSeries[ss];
ScaleLieSeries[c_, s_LieSeries] := ScaleLieSeries[c, s] = New[LieSeries[ser],
ser[d_Integer] := ser[d] = Expand[c*s[d]]
];
LieSeries /: c_*s_LieSeries := ScaleLieSeries[c,s];
LieSeries /: D[ls_LieSeries, vars__] := New[LieSeries[ser],
ser[d_Integer] := ser[d] = D[ls[d], vars]
];
IntegrateLieSeries[ls_LieSeries, {s_, s0_, s1_}] := IntegrateLieSeries[ls, {s, s0, s1}] = New[LieSeries[ser],
ser[d_Integer] := ser[d] = Expand[\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(s0\), \(s1\)]\(ls[d] \[DifferentialD]s\)\)]
];
LieSeries /: Integrate[ls_LieSeries, {s_, s0_, s1_}] := IntegrateLieSeries[ls, {s,s0,s1}];


(* ::Input::Initialization:: *)
b[s1_LieSeries, s2_LieSeries] := b[s1,s2] = New[LieSeries[ser],
ser[d_Integer] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d - 1\)]\(b[s1[k], \ s2[d - k]]\)\)
];
b[s_LieSeries, y_] := b[s, MakeLieSeries[y]];
b[x_, s_LieSeries] := b[MakeLieSeries[x], s];


(* ::Input::Initialization:: *)
LieSeries /: EulerE[s_LieSeries] := New[LieSeries[ser],
ser[d_Integer] := ser[d] = Expand[d*s[d]]
];
DegreeScale[h_][s_LieSeries|s_CWSeries] := New[Head[s][ser],
ser[d_Integer] := ser[d] = Expand[h^d s[d]]
];


(* ::Input::Initialization:: *)
adSeries[1, x_LieSeries][\[Psi]_LieSeries] := adSeries[1,x][\[Psi]] = \[Psi];
adSeries[ad^n_., x_LieSeries][\[Psi]_LieSeries] := adSeries[ad^n,x][\[Psi]] = New[LieSeries[ser],
ser[d_Integer] := ser[d] =  b[x,adSeries[ad^(n-1), x][\[Psi]]][d]
];
adSeries[f_, x_LieSeries][\[Psi]_LieSeries] := adSeries[f, x][\[Psi]] = New[LieSeries[ser],
ser[d_Integer] := ser[d] = Module[{c},
Expand[Sum[
c=SeriesCoefficient[f,{ad, 0, k}];
If[c === 0, 0, c*adSeries[ad^k, x][\[Psi]][d]],
{k,0,d-1}
]]
]
];
adSeries[f_, x_][\[Psi]_] := adSeries[f, MakeLieSeries[x]][MakeLieSeries[\[Psi]]];
Ad[x_] := adSeries[E^ad, x];


(* ::Input::Initialization:: *)
LieDerivation[der_][es___] := der[es];
LieDerivation[rules___Rule] := LieDerivation[{rules}];
LieDerivation[rules_List] := LieDerivation[rules] = New[LieDerivation[der],
der[Support] = First /@ rules;
(der[w_LW] /; Deg[w] == 1) := (der[w] = MakeLieSeries[w /. Append[rules, _LW -> 0]]);
der[w_LW] := der[w] = Module[{x,y},
{x,y}=LyndonFactorization[w];
AddLieSeries[b[der[x],y], b[x,der[y]]]
];
der[s_LieSeries] := der[s] = New[LieSeries[ser],
ser[d_] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d\)]\(\(der[s[k]]\)[d]\)\)];
der[as_ASeries] := der[as] = New[ASeries[ser],
ser[d_] := ser[d] = Sum[
Expand[as[k] /. w_AW :> Sum[
NonCommutativeMultiply[
Take[w, j-1],
\[Iota][der[LW[w[[j]]]][d-k+1]],
Drop[w,j]
],
{j,k}
]],
{k,1, d}
]
];
der[cws_CWSeries] := der[cws] = New[CWSeries[ser],
ser[d_] := ser[d] = Sum[
Expand[cws[k] /. w_CW :> Sum[
tr[NonCommutativeMultiply[
AW@@Take[w, j-1],
\[Iota][der[LW[w[[j]]]][d-k+1]],
AW@@Drop[w,j]
]],
{j,k}
]],
{k,1, d}
]
];
der[expr_][d_] := Expand[expr /. {w_LW :> der[w][d], s_LieSeries :> der[s][d]}];
];


(* ::Input::Initialization:: *)
LieDerivation /: Plus[ders__LieDerivation] := LieDerivation[Table[
u -> Total[#[u]& /@ {ders}],
{u, Union@@(#[Support]& /@ {ders})}
]];
LieDerivation /: c_*der_LieDerivation := LieDerivation[Table[
u -> (c der[u]),
{u, der@Support}
]];


(* ::Input::Initialization:: *)
b[der1_LieDerivation, der2_LieDerivation] := LieDerivation[Table[
u -> der1[der2[u]]-der2[der1[u]],
{u, (der1@Support)\[Union](der2@Support)}
]];


(* ::Input::Initialization:: *)
DerivationPower[0, der_LieDerivation][\[Psi]_LieSeries|\[Psi]_CWSeries] := DerivationPower[0,der][\[Psi]] = New[Head[\[Psi]][ser],
ser[d_Integer] := ser[d] = \[Psi][d]
];
DerivationPower[n_Integer, der_LieDerivation][\[Psi]_LieSeries|\[Psi]_CWSeries] := DerivationPower[n,x][\[Psi]] = New[Head[\[Psi]][ser],
ser[d_Integer] := ser[d] =  der[DerivationPower[n-1,der][\[Psi]]][d]
];
DerivationSeries[__][0] = 0;
DerivationSeries[f_, ld_LieDerivation][\[Psi]_LieSeries|\[Psi]_CWSeries] := DerivationSeries[f, ld][\[Psi]] = New[Head[\[Psi]][ser],
ser[d_Integer] := ser[d] = Module[{c},
Expand[Sum[
c=SeriesCoefficient[f,{der, 0, k}];
If[c==0, 0, c*DerivationPower[k, ld][\[Psi]][d]],
{k,0,d}
]]
]
];
DerivationExp[ld_LieDerivation] := DerivationSeries[E^der, ld];
LieDerivation /: Exp[ld_LieDerivation] := DerivationExp[ld];
LieDerivation /: E^ld_LieDerivation := DerivationExp[ld];


(* ::Input::Initialization:: *)
LieMorphism[mor_][es___] := mor[es];
LieMorphism[rules__Rule] := LieMorphism[{rules}];
LieMorphism[rules_List] := LieMorphism[rules] = New[LieMorphism[mor],
mor[Support] = First /@ rules;
(mor[w_LW] /; Deg[w] == 1) := (mor[w] = MakeLieSeries[w /. rules]);
mor[w_LW] := (mor[w] = b @@ (mor /@ LyndonFactorization[w]));
mor[AW[]] = MakeASeries[AW[]];
mor[AW[w_]] := mor[AW[w]] = \[Iota][MakeLieSeries[LW[w] /. rules]];
mor[w_AW] := mor[w] = Module[{w1,w2},
w1=Take[w, Floor[Length[w]/2]];
w2=Drop[w, Floor[Length[w]/2]];
mor[w1]**mor[w2]
];
mor[w_CW] := tr[mor[AW@@w]];
mor[s_LieSeries] := mor[s] = New[LieSeries[ser],
ser[d_] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d\)]\(\(mor[s[k]]\)[d]\)\)];
mor[cws_CWSeries] := mor[cws] = New[CWSeries[ser],
ser[d_] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d\)]\(\(mor[cws[k]]\)[d]\)\)];
mor[expr_][d_] := Expand[expr /. (w_LW|w_AW|w_CW) :> mor[w][d]];
(* Added 150217, commented out later that day *) (*mor[expr_] := Expand[expr /. (w_LW|w_AW|w_CW) \[RuleDelayed] mor[w]]*)
];


(* ::Input::Initialization:: *)
LieMorphism /: Inverse[mor_LieMorphism] := InvertLieMorphism[mor];
InvertLieMorphism[mor_LieMorphism] := InvertLieMorphism[mor] =
LieMorphism[Table[
LW@u -> New[LieSeries[uimg],
uimg[1] = LW@u;
uimg[d_Integer] /; d>1 := uimg[d] = -\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d - 1\)]\((\(mor[uimg[k]]\)[d])\)\)
],
{u,mor[Support]}
]]


(* ::Input::Initialization:: *)
StableApply[mor_LieMorphism, (type:(LieSeries|ASeries|CWSeries))[s_]] := (
StableApply[mor, type[s]] = New[type[ser],
ser[d_] := ser[d] = Nest[mor, type[s], d][d]
]
);


(* ::Input::Initialization:: *)
BCHBase = New[LieSeries[bch],
bch[1] = LW@"x"+LW@"y";
bch[d_Integer] := bch[d] = Expand[Plus[
adSeries[E^-ad, MakeLieSeries[LW@"y"]][MakeLieSeries[LW@"x"]][d],
-adSeries[(1-E^-ad)/ad-1, LieSeries[bch]][EulerE[LieSeries[bch]]][d]
]/d]
];
BCH[x_, y_] := LieMorphism[{LW["x"] -> x, LW["y"] -> y}][BCHBase];


(* ::Input::Initialization:: *)
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;
Deg[w_AW] := Length[w];
AW[AW[w_]]:=AW[w];
AW[w1___] ** AW[w2___] := AW[w1,w2];
b[w_AW, z_AW] := w**z-z**w;


(* ::Input::Initialization:: *)
ASeries[ser_Symbol][{dd_Integer}] := AS@@Table[ser[d], {d, 0, dd}];
ASeries[as_Symbol][es___] := as[es];
Format[s_ASeries, StandardForm] := s[{$SeriesShowDegree}];
MakeASeries[as_CWSeries] := as;
MakeASeries[expr_] := MakeASeries[expr] = New[ASeries[ser],
ser[d_Integer] := ser[d] = Expand[expr /. w_AW /; Deg[w] != d -> 0]
];
(s1_ASeries ** s2_ASeries) := (s1**s2) = New[ASeries[ser],
ser[d_Integer] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 0\), \(d\)]\(s1[k] ** s2[d - k]\)\)];


(* ::Input::Initialization:: *)
\[Iota][w_LW] /; Deg[w]==1 := AW@@w;
\[Iota][w_LW] := \[Iota][w] = b @@ (\[Iota] /@ LyndonFactorization[w]);
\[Iota][expr_] := Expand[expr /. w_LW :> \[Iota][w]];
\[Iota][ls_LieSeries] := \[Iota][ls] = New[ASeries[as],
as[0]=0;
as[d_] /; d>0 := as[d] = \[Iota][ls[d]]
];


(* ::Input::Initialization:: *)
\[Tau][y_LW, w_LW] /; Deg[y]==1 := \[Tau][y,w] = Which[
y===w, AW[],
Deg[w]===1, 0,
True, Module[{w1,w2},
{w1,w2}=LyndonFactorization[w];
\[Iota][w1]**\[Tau][y,w2]-\[Iota][w2]**\[Tau][y,w1]
]
];
\[Tau][y_, ls_LieSeries] := \[Tau][y, ls] = New[ASeries[as],
as[d_] := as[d] = \[Tau][LW[y], ls[d+1]]
];
\[Tau][y_, expr_] := Expand[expr /. w_LW :> \[Tau][LW[y], w]];


(* ::Input::Initialization:: *)
Subscript[tr, y_][expr_] /; Deg[LW@y]==1 := tr[\[Tau][LW@y, expr]];


(* ::Input::Initialization:: *)
RotateToMinimal[l_]:=Module[
{bestl=l,rotatedl=RotateLeft[l]},
While[rotatedl=!=l,
bestl=First[Sort[{bestl,rotatedl}]];
rotatedl=RotateLeft[rotatedl]
];
bestl
];


(* ::Input::Initialization:: *)
Deg[w_CW] := Length[w];
Format[w_CW] := TopBracketForm[w];
AllCyclicWords[d_Integer, ab_List] := AllCyclicWords[d, ab] = Union[tr[AllWords[d, ab]]];
CWSeries[cws_Symbol][es___] := cws[es];
CWSeries[ser_Symbol][{dd_Integer}] := Append[TopBracketForm[CWS@@Table[ser[d], {d, dd}]], "..."];
Format[s_CWSeries, StandardForm] := TopBracketForm[s[{$SeriesShowDegree}]];
CWS[cws_] := MakeCWSeries[cws];
MakeCWSeries[cws_CWSeries] := cws;
MakeCWSeries[expr_] := MakeCWSeries[expr] = New[CWSeries[ser],
ser[d_Integer] := ser[d] = Expand[
expr /. w_CW :> If[Deg[w] == d , RotateToMinimal@w, 0]
]
];
RandomCWSeries[ab_List, opts___Rule] := New[CWSeries[ser],
Module[
{rand = Randomizer /. {opts} /. Randomizer -> (RandomInteger[{-2Deg[#]!, 2Deg[#]!}]/Deg[#]!&)},
ser[d_Integer] := ser[d] = Plus @@ ((rand[#]*#)& /@ AllCyclicWords[d, ab])
]];
s1_CWSeries \[Congruent] s2_CWSeries := New[BooleanSequence[bs],
bs[0]=True;
bs[d_Integer] := bs[d-1]&&(s1[d]==s2[d]);
];
AddCWSeries[ss___CWSeries] := AddCWSeries[ss] = New[CWSeries[ser],
ser[d_Integer] := ser[d] = Plus @@ ((#[d])& /@ {ss})
];
CWSeries /: Plus[ss__CWSeries] := AddCWSeries[ss];
ScaleCWSeries[c_, s_CWSeries] := ScaleCWSeries[c, s] = New[CWSeries[ser],
ser[d_Integer] := ser[d] = Expand[c*s[d]]
];
CWSeries /: c_*s_CWSeries := ScaleCWSeries[c,s];
IntegrateCWSeries[cws_CWSeries, {s_, s0_, s1_}] := IntegrateCWSeries[cws, {s, s0, s1}] = New[CWSeries[ser],
ser[d_Integer] := ser[d] = Expand[\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(s0\), \(s1\)]\(cws[d] \[DifferentialD]s\)\)]];
CWSeries /: Integrate[cws_CWSeries, {s_, s0_, s1_}] := IntegrateCWSeries[cws, {s,s0,s1}];


(* ::Input::Initialization:: *)
cw[gs___] := CW@@Replace[RotateToMinimal@{gs}, LW[g_]:>g, {1}];


(* ::Input::Initialization:: *)
tr[w_AW] := tr[w] = cw@@w;
tr[expr_] := expr /. aw_AW :> tr[aw];
tr[as_ASeries] := tr[as] = New[CWSeries[cws], cws[d_] := cws[d] = tr[as[d]]];


(* ::Input::Initialization:: *)
div[LW[y_], w_LW]  := div[LW@y,w] = tr[(AW@y)**\[Tau][LW@y,w]];
div[y_, ls_LieSeries] := div[y, ls] = New[CWSeries[cws],
cws[d_] := cws[d] = div[LW[y], ls[d]]
];
div[y_, expr_] := Expand[expr /. w_LW :> div[LW[y], w]];
div[y_, p_Plus] := Total[div[LW[y], List@@p]];
Subscript[div, y_][expr_] := div[y, expr];


(* ::Input::Initialization:: *)
JA[-1, ___] = MakeCWSeries[0];
JA[n_, y_LW, \[Mu]_LieSeries, ss_] := JA[n,y,\[Mu],ss] = Module[
{s, s\[Mu], \[Mu]s},
s\[Mu] = ScaleLieSeries[s, \[Mu]];
\[Mu]s=StableApply[LieMorphism[{y -> Ad[ScaleLieSeries[1,s\[Mu]]][LW[z]]}], \[Mu]];
\[Mu]s = \[Mu]s // LieMorphism[{LW[z] -> y}];
IntegrateCWSeries[
AddCWSeries[
JA[n-1, y, \[Mu], s] // LieDerivation[{y -> b[\[Mu]s, y]}],
div[y, \[Mu]s ]
],
{s,0,ss}
]
];
JA[y_LW, \[Mu]_LieSeries] := JA[y,\[Mu]] = Module[{s}, New[CWSeries[cws],
cws[d_Integer] := cws[d] = JA[d-1, y, \[Mu], s][d] /. s->1
]];


(* ::Input::Initialization:: *)
SetAttributes[AngleBracket, Orderless];
\[LeftAngleBracket]{\[Lambda]___}\[RightAngleBracket] := \[LeftAngleBracket]\[Lambda]\[RightAngleBracket];
(*\[LeftAngleBracket]x_ \[Rule] 0, rest___\[RightAngleBracket] := \[LeftAngleBracket]x\[Rule]LS[0], rest\[RightAngleBracket];*)
Support[\[LeftAngleBracket]\[Lambda]___\[RightAngleBracket]] := First /@ {\[Lambda]};
\[Lambda]_\[Backslash]key_ := DeleteCases[\[Lambda], key -> _];
\[Lambda]_ \[Backslash] keys_List := Fold[#1\[Backslash]#2&, \[Lambda], keys];
Subscript[\[LeftAngleBracket]\[Lambda]___\[RightAngleBracket], s_] := s /. {\[Lambda]};
(\[Lambda]_AngleBracket)[d_] := \[LeftAngleBracket]Table[u -> Subscript[\[Lambda], u][d],{u, Support[\[Lambda]]}]\[RightAngleBracket];
AngleBracket /: EulerE[\[Lambda]_AngleBracket] := MapAt[EulerE, \[Lambda], {All,2}];
Crop[\[Lambda]_AngleBracket, d_] := MapAt[Crop[#,d]&, \[Lambda], {All,2}];
AngleBracket /: \[Lambda]1_AngleBracket \[Congruent] \[Lambda]2_AngleBracket /; Support[\[Lambda]1]==Support[\[Lambda]2] :=
And @@ Table[Subscript[\[Lambda]1, s]\[Congruent]Subscript[\[Lambda]2, s], {s, Support[\[Lambda]1]}];
AngleBracket /: Plus[\[Lambda]s__AngleBracket] := \[LeftAngleBracket]Table[
u -> Total[Subscript[#, u]& /@ {\[Lambda]s}],
{u, Union@@(Support[#]& /@ {\[Lambda]s})}
]\[RightAngleBracket];
AngleBracket /: c_*\[Lambda]_AngleBracket := \[LeftAngleBracket]Table[
u -> Expand[c Subscript[\[Lambda], u]],
{u, Support[\[Lambda]]}
]\[RightAngleBracket];
AngleBracket /: \[Lambda]_AngleBracket // der_LieDerivation := MapAt[der, \[Lambda], {All,2}];
AngleBracket /: \[Lambda]_AngleBracket // DerivationSeries[f_, ld_LieDerivation]:= MapAt[DerivationSeries[f,ld], \[Lambda], {All,2}];
AngleBracket /: \[Lambda]_AngleBracket // mor_LieMorphism := MapAt[mor, \[Lambda], {All,2}];
AngleBracket /: Integrate[\[Lambda]_AngleBracket, {s_, s0_, s1_}] := \[LeftAngleBracket]Table[u -> Integrate[Subscript[\[Lambda], u], {s,s0,s1}],{u, Support[\[Lambda]]}]\[RightAngleBracket];
\[Lambda]_AngleBracket //  DegreeScale[h_] := MapAt[DegreeScale[h], \[Lambda], {All,2}];
b[\[Lambda]1_AngleBracket, \[Lambda]2_AngleBracket] := \[LeftAngleBracket]Table[s -> b[Subscript[\[Lambda]1, s], Subscript[\[Lambda]2, s]],{s, Support[\[Lambda]1]\[Union]Support[\[Lambda]2]}]\[RightAngleBracket];
AngleBracket /: D[\[Lambda]_AngleBracket, vars__] := \[LeftAngleBracket]Table[s->D[Subscript[\[Lambda], s], vars], {s, Support[\[Lambda]]}]\[RightAngleBracket];


(* ::Input::Initialization:: *)
AngleBracket /: D[\[Lambda]_AngleBracket] := LieDerivation[List@@\[Lambda] /. (s_->\[Lambda]s_):>(LW[s]->b[LW[s],\[Lambda]s])];
AngleBracket /: Subscript[D, \[Lambda]_AngleBracket] := D[\[Lambda]];


(* ::Input::Initialization:: *)
tb[\[Lambda]1_AngleBracket, \[Lambda]2_AngleBracket] := \[LeftAngleBracket]Table[
s -> b[Subscript[\[Lambda]1, s], Subscript[\[Lambda]2, s]] + Subscript[D, \[Lambda]1][Subscript[\[Lambda]2, s]] - Subscript[D, \[Lambda]2][Subscript[\[Lambda]1, s]],
{s, Support[\[Lambda]1]\[Union]Support[\[Lambda]2]}
]\[RightAngleBracket];


(* ::Input::Initialization:: *)
AngleBracket /: CC[\[LeftAngleBracket]\[Lambda]___\[RightAngleBracket]] := LieMorphism[{\[Lambda]} /. (s_->\[Lambda]s_):>(LW[s]->Ad[\[Lambda]s][LW[s]])];
CC[u_, \[Gamma]_] := CC[\[LeftAngleBracket]u->\[Gamma]\[RightAngleBracket]];
Subscript[CC, u_][\[Gamma]_] := CC[u,\[Gamma]];
RC[\[Lambda]_AngleBracket] := Inverse[CC[-\[Lambda]]];
RC[u_, \[Gamma]_] := RC[\[LeftAngleBracket]u->\[Gamma]\[RightAngleBracket]];
Subscript[RC, u_][\[Gamma]_] := RC[u,\[Gamma]];
ad[u_, \[Gamma]_LieSeries] := LieDerivation[u->b[\[Gamma],u]];
Subscript[ad, u_][\[Gamma]_] := ad[u,\[Gamma]];


(* ::Input::Initialization:: *)
Subscript[J, u_][\[Gamma]_] := Subscript[J, u][\[Gamma]] = Module[{s},\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(1\)]\(\((\(\(\[Gamma]\  // \ 
\(\*SubscriptBox[\(RC\), \(u\)]\)[s\ \[Gamma]]\)\  // \ 
\*SubscriptBox[\(div\), \(u\)]\)\  // \ 
\(\*SubscriptBox[\(CC\), \(u\)]\)[\(-s\)\ \[Gamma]])\) \[DifferentialD]s\)\)];


(* ::Input::Initialization:: *)
J[u_, \[Gamma]_] := Subscript[J, u][\[Gamma]];


(* ::Input::Initialization:: *)
div[\[Lambda]_AngleBracket] := Sum[Subscript[div, a][Subscript[\[Lambda], a]], {a, Support[\[Lambda]]}];
j[\[Lambda]_AngleBracket] := div[\[Lambda]] // DerivationSeries[(E^der-1)/der, Subscript[D, \[Lambda]]];


(* ::Input::Initialization:: *)
LieMorphism[rules__Rule, keys_AngleBracket, br_] := LieMorphism[{rules}, keys, br];
LieMorphism[rules_List, keys_AngleBracket, br_] := New[LieMorphism[mor],
mor[Support] = First /@ rules;
mor[0] := mor[0] = AngleBracket@@Table[s->LS[0], {s, List@@keys}];
(mor[w_LW] /; Deg[w] == 1) := (mor[w] = w /. rules);
mor[w_LW] := (mor[w] = br @@ (mor /@ LyndonFactorization[w]));
mor[expr_][d_] := Expand[expr /. w_LW :> mor[w][d]];
mor[ls_LieSeries] := mor[ls] = Module[{ser},
\[LeftAngleBracket]Table[
ReleaseHold[Hold[
ser[]=Hold[AngleBracketFromLieMorphism[mor, ls, ss]];
ser[d_Integer] := ser[d] = Subscript[(\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d\)]\(\(mor[ls[k]]\)[d]\)\)), ss];
ss->LieSeries[ser]
] /. {ss->s, ser -> Unique[LieSeriesInAngleBracket]}],
{s, List@@keys}
]\[RightAngleBracket]]
];


(* ::Input::Initialization:: *)
BCH[x_, y_, keys_AngleBracket, br_] := LieMorphism[{LW["x"] -> x, LW["y"] -> y}, keys, br][BCHBase];
Subscript[BCH, br_][\[Lambda]1_AngleBracket, \[Lambda]2_AngleBracket] := BCH[\[Lambda]1, \[Lambda]2, \[LeftAngleBracket]Support[\[Lambda]1]\[Union]Support[\[Lambda]2]\[RightAngleBracket], br];
adSeries[f_,\[Lambda]1_AngleBracket, br_][\[Lambda]2_AngleBracket] := LieMorphism[{LW["x"] -> \[Lambda]1, LW["y"] -> \[Lambda]2}, \[LeftAngleBracket]Support[\[Lambda]1]\[Union]Support[\[Lambda]2]\[RightAngleBracket], br][adSeries[f, LW["x"]][LW["y"]]];
adSeries[f_,\[Lambda]1_AngleBracket][\[Lambda]2_AngleBracket] := adSeries[f,\[Lambda]1,b][\[Lambda]2];


(* ::Input::Initialization:: *)
Subscript[\[CapitalGamma], 0,t_][\[Lambda]_AngleBracket] := t \[Lambda];
Subscript[\[CapitalGamma], n_,t_][\[Lambda]_AngleBracket] := Subscript[\[CapitalGamma], n,t][\[Lambda]] = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(t\)]\(\((\(\[Lambda]\  // \ 
\*SuperscriptBox[\(E\), \(\(-\[Tau]\)\ 
\*SubscriptBox[\(D\), \(\[Lambda]\)]\)]\)\  // \ adSeries[
\*FractionBox[\(ad\), \(
\*SuperscriptBox[\(E\), \(ad\)] - 1\)], \ 
\(\*SubscriptBox[\(\[CapitalGamma]\), \(n - 1, \[Tau]\)]\)[\[Lambda]]])\) \[DifferentialD]\[Tau]\)\);
Subscript[\[CapitalGamma], t_][\[Lambda]_AngleBracket] := Subscript[\[CapitalGamma], t][\[Lambda]] = \[LeftAngleBracket]Table[
s -> New[LieSeries[ser],
With[{s=s}, ser[d_Integer] := ser[d] = (Subscript[Subscript[\[CapitalGamma], d-1,t][\[Lambda]], s])[d]];
],
{s, Support[\[Lambda]]}
]\[RightAngleBracket];
\[CapitalGamma][\[Lambda]_AngleBracket] := Subscript[\[CapitalGamma], 1][\[Lambda]];
(*Subscript[\[CapitalGamma], -1,_][\[Lambda]_AngleBracket] := \[LeftAngleBracket]Table[s\[Rule]MakeLieSeries[0], {s, Support[\[Lambda]]}]\[RightAngleBracket];*)
(*Subscript[\[CapitalGamma], n_,t_][\[Lambda]_AngleBracket] := Subscript[\[CapitalGamma], n,t][\[Lambda]] = Module[{\[Tau]},
\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(t\)]\(\((\(\[Lambda]\  // \ 
\*SuperscriptBox[\(\[ExponentialE]\), \(\(-\[Tau]\)\ 
\*SubscriptBox[\(D\), \(\[Lambda]\)]\)]\)\  // \ adSeries[
\*FractionBox[\(ad\), \(
\*SuperscriptBox[\(\[ExponentialE]\), \(ad\)] - 1\)], \ 
\(\*SubscriptBox[\(\[CapitalGamma]\), \(n - 1, \[Tau]\)]\)[\[Lambda]]])\)\[DifferentialD]\[Tau]\)\)
];*)
(*Subscript[\[CapitalGamma], n_,t_][\[Lambda]_AngleBracket] := Subscript[\[CapitalGamma], n,t][\[Lambda]] = Module[{\[Tau],\[CapitalGamma]0},
\[CapitalGamma]0=Subscript[\[CapitalGamma], n-1,\[Tau]][\[Lambda]];
\[LeftAngleBracket]Table[
s\[Rule]\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(t\)]\(\((\(
\*SubscriptBox[\(\[Lambda]\), \(s\)]\  // \ 
\*SuperscriptBox[\(\[ExponentialE]\), \(\(-\[Tau]\)\ 
\*SubscriptBox[\(D\), \(\[Lambda]\)]\)]\)\  // \ adSeries[
\*FractionBox[\(ad\), \(
\*SuperscriptBox[\(\[ExponentialE]\), \(ad\)] - 1\)], \ 
\*SubscriptBox[\(\[CapitalGamma]0\), \(s\)]])\)\[DifferentialD]\[Tau]\)\),
{s, Support[\[Lambda]]}
]\[RightAngleBracket]
];*)


(* ::Input::Initialization:: *)
Subscript[\[CapitalLambda], 0,t_][\[Lambda]_AngleBracket] := t \[Lambda];
Subscript[\[CapitalLambda], n_,t_][\[Lambda]_AngleBracket] := Subscript[\[CapitalLambda], n,t][\[Lambda]] = Module[{\[Tau],\[CapitalLambda]0},
\[CapitalLambda]0=Subscript[\[CapitalLambda], n-1,\[Tau]][\[Lambda]];
\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(t\)]\(\((\(\[Lambda]\  // \ 
\*SuperscriptBox[\(E\), 
SubscriptBox[\(D\), \(\[CapitalLambda]0\)]]\)\  // \ adSeries[
\*FractionBox[\(ad\), \(
\*SuperscriptBox[\(E\), \(ad\)] - 1\)], \ \[CapitalLambda]0, \ tb])\) \[DifferentialD]\[Tau]\)\)
];
Subscript[\[CapitalLambda], t_][\[Lambda]_AngleBracket] := Subscript[\[CapitalLambda], t][\[Lambda]] = \[LeftAngleBracket]Table[
s -> New[LieSeries[ser],
With[{s=s}, ser[d_Integer] := ser[d] = (Subscript[Subscript[\[CapitalLambda], d-1,t][\[Lambda]], s])[d]];
],
{s, Support[\[Lambda]]}
]\[RightAngleBracket];
\[CapitalLambda][\[Lambda]_AngleBracket] := Subscript[\[CapitalLambda], 1][\[Lambda]];


(* ::Input::Initialization:: *)
LS[ab_List, coefs_] := New[LieSeries[ser],
ser[setter] = Null;
ser[d_Integer, UndeterminedCoefficients] := Cases[coefs @@@ AllLyndonWords[d, ab], _coefs];
ser[d_Integer] := If[ser[setter] =!= Null,
ser[setter][d]; ser[d],
Plus @@ ((#*coefs@@#)& /@ AllLyndonWords[d, ab])
];
];
CWS[ab_List, coefs_] := New[CWSeries[ser],
ser[setter] = Null;
ser[d_Integer, UndeterminedCoefficients] := Cases[coefs @@@ AllCyclicWords[d, ab], _coefs];
ser[d_Integer] := If[ser[setter] =!= Null,
ser[setter][d]; ser[d],
Plus @@ ((#*coefs@@#)& /@ AllCyclicWords[d, ab])
];
];


(* ::Input::Initialization:: *)
Options[SeriesSolve] = {
Arbitrator -> 0 (* should be 0 or a pure function that takes a list of unsettled variables and returns a list of their arbitrated values *)
};
SeriesSolve::ArbitrarilySetting = "In degree `1` arbitrarily setting `2`.";
SeriesSolve::NoSolution = "No solution in degree `1`.";
SeriesSolve[unknown_, eqns_, opts___] /; Head[unknown]=!=List := SeriesSolve[{unknown}, eqns, opts];
SeriesSolve[unknowns_List, eqns_, opts___] := Module[{
arbitrator = Arbitrator /. {opts} /. Options[SeriesSolve],
msgs, MessageStream, solver, lineqs, gens, vars, sol, fvars, avals, d
},
If[arbitrator===0, arbitrator = Replace[#, _->0, 1]&];
msgs={};
MessageStream /: Read[MessageStream] := msgs;
solver[0]=Null;
solver[n_] := (
solver[n-1];
(#[setter] = Null)& /@ (First /@ unknowns);
lineqs=eqns[n];
lineqs=If[Head[lineqs]===And, List @@ lineqs, List@lineqs];
gens=Union[Cases[lineqs, _LW|_CW, \[Infinity]]];
lineqs=Flatten[Replace[lineqs ,
lhs_==rhs_ :> ((Coefficient[lhs, #]==Coefficient[rhs,#])& /@ gens),
{1}]];
vars = Union@@((#[n, UndeterminedCoefficients])& /@ unknowns);
sol= Quiet[Solve[lineqs, vars], {Solve::svars}];
If[sol==={},
AppendTo[msgs, {"NoSolution", n}]; Message[SeriesSolve::NoSolution, n]; Abort[],
(* Else *) sol=sol[[1]];
fvars=Complement[vars, First /@ sol];
avals=arbitrator[fvars];
AppendTo[msgs, {"ArbitrarilySetting", n, Thread[(Hold /@ fvars)->avals]}];
If[fvars =!={}, Message[SeriesSolve::ArbitrarilySetting, n, Thread[fvars->avals]]];
MapThread[(#1=#2)& , {fvars, avals}];
sol /. (Rule -> Set);
(#[n]=#[n])& /@ (First /@ unknowns);
(#[setter] = solver)& /@ (First /@ unknowns);
solver[n]=Null
]
);
(#[setter] = solver)& /@ (First /@ unknowns);
MessageStream
];


(* ::Input::Initialization:: *)
DK[_, 0] = 0;
DK /: DK[k_, x1_] + DK[k_, x2_] := DK[k, x1+x2];
DK /: c_*DK[k_, x_] := DK[k, Expand[c x]];
DK /: b[DK[k1_, c_*x1_LW], DK[k2_, x2_]] := Expand[c b[DK[k1, x1], DK[k2,x2]]];
DK /: b[DK[k1_, x1_LW], DK[k2_, c_*x2_LW]] := Expand[c b[DK[k1, x1], DK[k2,x2]]];
DK /: b[DK[k1_, x1_Plus], DK[k2_, x2_]] := b[DK[k1, #], DK[k2,x2]]& /@ x1;
DK /: b[DK[k1_, x1_], DK[k2_, x2_Plus]] := b[DK[k1, x1], DK[k2, #]]& /@ x2;
DK /: b[DK[k_, x1_], DK[k_, x2_]] := DK[k, b[x1,x2]];
DK /: b[DK[k1_, x1_], DK[k2_, x2_]] /; k1>k2 := b[DK[k2, Expand[-x2]], DK[k1, x1]];
DK /: b[DK[k1_, LW@i1_], DK[k2_, LW@i2_]] /; k1<k2 := b[DK[k1,LW@i1], DK[k2, LW@i2]] = Which[
i1==i2, DK[k2, -b[LW@k1, LW@i2]],
k1==i2, DK[k2, -b[LW@i1, LW@i2]],
True, 0
];
DK /: b[DK[k1_, w1_LW], DK[k2_, w2_LW]] /; Deg[w1]>1 := b[DK[k1, w1], DK[k2, w2]] = Module[{x,y},
{x,y}=LyndonFactorization[w1];
b[b[DK[k1, x], DK[k2,w2]], DK[k1,y]] + b[DK[k1,x], b[DK[k1,y], DK[k2,w2]]]
];
DK /: b[DK[k1_, w1_LW], DK[k2_, w2_LW]] /; Deg[w2]>1 := b[DK[k1, w1], DK[k2, w2]] = Module[{x,y},
{x,y}=LyndonFactorization[w2];
b[b[DK[k1, w1], DK[k2,x]], DK[k2,y]] + b[DK[k2,x], b[DK[k1,w1], DK[k2,y]]]
];
t[i_,j_] := DK[Max[i,j], LW@Min[i,j]];
TopBracketForm[DK[k_, x_]] := x // LieMorphism[Table[LW@i->LW@Subscript[t, 10i+k], {i,k-1}]];
(*Format[dk_DK] := TopBracketForm[dk];*)


(* ::Input::Initialization:: *)
\[Sigma] /: expr_^s_\[Sigma] := expr // s;\[Sigma][lft___, i_Integer, rgt___] := \[Sigma][lft, IntegerDigits[i], rgt];
_\[Sigma][0] = 0;
x_Plus // s_\[Sigma] := s[#]& /@ x;
DK[k_, x_Plus] //s_\[Sigma]:= s[DK[k, #]]& /@ x;
DK[k_, c_*w_LW] // s_\[Sigma] := Expand[c*s[DK[k,w]]];
DK[k_, LW@i_] // s_\[Sigma] := Sum[t[\[Alpha], \[Beta]], {\[Alpha],s[[i]]}, {\[Beta], s[[k]]}];
DK[k_, w_LW] // s_\[Sigma] := b@@(s[DK[k,#]]& /@ LyndonFactorization[w]);


(* ::Input::Initialization:: *)
DKSeries[ser_Symbol][{dd_Integer}] := TopBracketForm[Append[
DKS@@Table[
ser[d] /. DK[k_,x_] :>  (x // LieMorphism[Table[LW@i->LW@Subscript[t, 10i+k], {i,k-1}]])@d,
{d, dd}
],
"..."]];
Format[dks_DKSeries, StandardForm] := dks[{$SeriesShowDegree}];
DKSeries[ser_Symbol][e___] := ser[e];
b[dk1_DKSeries, dk2_DKSeries]:= b[dk1, dk2] = New[DKSeries[ser],
ser[d_Integer] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(d - 1\)]\(b[dk1[j], \ dk2[d - j]]\)\)
];
DKS[dks_DKSeries] := dks;
DKS[expr_] := DKS[expr] = New[DKSeries[ser],
ser[d_Integer] := ser[d] = Expand[expr /. w_LW /; Deg[w] != d -> 0]
];
DKS[k_Integer, coefs_] := New[DKSeries[ser],
ser[setter] = Null;
ser[d_Integer, UndeterminedCoefficients] := Cases[
Join@@Table[coefs[j, Sequence@@#]& /@ AllLyndonWords[d, Range[j-1]], {j,2,k}],
_coefs];
ser[d_Integer] := If[ser[setter] =!= Null,
ser[setter][d]; ser[d],
Sum[
Plus @@ ((DK[j,#]*coefs[j,Sequence@@#])& /@ AllLyndonWords[d, Range[j-1]]),
{j,2,k}
]
];
];
AddDKSeries[ss___DKSeries] := AddDKSeries[ss] = New[DKSeries[ser],
ser[d_Integer] := ser[d] = Plus @@ ((#[d])& /@ {ss})
];
DKSeries /: Plus[ss__DKSeries] := AddDKSeries[ss];
ScaleDKSeries[c_, s_DKSeries] := ScaleDKSeries[c, s] = New[DKSeries[ser],
ser[d_Integer] := ser[d] = Expand[c*s[d]]
];
DKSeries /: c_*s_DKSeries := ScaleDKSeries[c,s];
s1_DKSeries \[Congruent] s2_DKSeries := New[BooleanSequence[bs],
bs[0]=True;
bs[d_Integer] := bs[d-1]&& 
(And @@ Cases[{Expand[s1[d]-s2[d]]}, DK[_,x_] :> (x==0), \[Infinity]]);
];
dks_DKSeries // s_\[Sigma] := dks//s = New[DKSeries[ser],
ser[d_Integer] := ser[d] = dks[d] // s
];
Subscript[(dks_DKSeries), k_Integer] := New[LieSeries[ser],
ser[d_Integer] := ser[d] = dks[d] /. {DK[k, x_] :> x, DK[__] -> 0}
];


(* ::Input::Initialization:: *)
Morphism[mor_][es___] := mor[es];
Morphism[LS, DKS, rules__Rule] := Morphism[LS, DKS, {rules}];
Morphism[LS, DKS, rules_List] := New[Morphism[mor],
mor[Support] = First /@ rules;
(mor[w_LW] /; Deg[w] == 1) := (mor[w] = w /. rules);
mor[w_LW] := (mor[w] = b @@ (mor /@ LyndonFactorization[w]));
mor[expr_][d_] := Expand[expr /. w_LW :> mor[w][d]];
mor[ls_LieSeries] := mor[ls] = New[DKSeries[ser],
ser[d_] := ser[d] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(d\)]\(\(mor[ls[k]]\)[d]\)\)];
];
BCH[x_DKSeries, y_DKSeries] := BCH[LW@"x", LW@"y"] // Morphism[LS, DKS, LW@"x"->x, LW@"y"->y];
DKSeries /: x_DKSeries ** y_DKSeries := BCH[x,y];


(* ::Input::Initialization:: *)
\[Alpha]Map[S_List][dks_DKSeries] := Sum[
LieMorphism[
Table[
LW@S[[i]] -> AngleBracket@@Table[
S[[j]]->LS[Switch[j,
i, LW@S[[k]],
k, LW@S[[i]],
 _, 0
]],
 {j, Length[S]}],
{i,k-1}
],
AngleBracket@@S, tb
][Subscript[dks, S[[k]]]],
{k,2,Length[S]}
];
\[Alpha]Map[S___][dks_DKSeries] := \[Alpha]Map[{S}][dks];


(* ::Input::Initialization:: *)
End[]; EndPackage[];
