Pensieve header: A free-Lie calculator, at long last. (Continues pensieve://2011-02/).

A Lyndon word is a word lexicographically smaller than all of its proper right factors; see http://katlas.math.toronto.edu/drorbn/bbs/show?shot=Chu-071214-182203.jpg

LyndonQ[w_String] := And @@ (
OrderedQ[{w, #}]& /@ Table[StringDrop[w,i], {i,1,StringLength[w]-1}]
);
AllWords[0, _List] = {""};
AllWords[n_ /; n>0, ab_List] := AllWords[n, ab] = Flatten[Outer[
StringJoin[#1, #2]&,
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 /: 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[⟨w⟩];
⟨w_⟩:=LW[w];
Deg[LW[x_]] := StringLength[x];

LyndonQ["abba"]

LyndonQ["ababb"]

AllWords[1, {"1", "2"}]

AllWords[3, {"1", "2"}]

AllLyndonWords[3, {"1", "2"}]

Table[Length[AllLyndonWords[k, {"1", "2"}]], {k,10}]

Table[Length[AllLyndonWords[k, {"1", "2", "3"}]], {k,8}]

TreeForm[LW["12122"] //. w_LW :> LyndonFactorization[w] /. LW[w_] :> w]

b[0,_] =0; b[_,0]=0;
b[c_*x_LW, y_] := Expand[c b[x,y]];
b[x_, c_*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] := Which[
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]],
b[x,b[y,z]]+b[b[x,z],y]
]
]
];

b[LW["112"], LW["122"]]

Outer[b, AllLyndonWords[5, {"1", "2"}], AllLyndonWords[5, {"1", "2"}]]

⟨"1122221222"⟩

⟨"1122221222"⟩//FullForm

Length /@ Flatten[
Outer[b, AllLyndonWords[5, {"1", "2"}], AllLyndonWords[5, {"1", "2"}]]
]

Outer[(b[#1,#2]+b[#2,#1])&,
AllLyndonWords[{6}, {"1", "2"}], AllLyndonWords[{6}, {"1", "2"}]
] // MatrixForm

Outer[(b[#1,b[#2,#3]]+b[#2,b[#3,#1]]+b[#3,b[#1,#2]])&,
AllLyndonWords[{5}, {"1", "2"}], AllLyndonWords[{5}, {"1", "2"}], AllLyndonWords[{5}, {"1", "2"}]
] // Flatten // Union

{x=LW["x"], y=LW["y"]}

ad[x_][y_]:=b[x,y];
MakeLieSeries[expr_, d_] := LieSeries @@ Table[
expr /. w_LW /; Deg[w] ≠ k -> 0,
{k,d}
];
LieSeries /: s1_LieSeries + s2_LieSeries := Module[
{l=Min[Length /@ {s1, s2}]},
LieSeries@@(Take[List@@s1, l]+Take[List@@s2,l])
];
LieSeries /: c_ * s_LieSeries := Expand[c*#]& /@ s;
b[s1_LieSeries, s2_LieSeries] := Module[
{d,k},
LieSeries @@ Table[
Sum[b[s1[[k]], s2[[d-k]]], {k,1,d-1}],
{d, 1+Min[Length /@ {s1, s2}]}
]
];
LieSeries /: EulerE[s_LieSeries] := LieSeries @@ Expand[Range[Length[s]]*(List@@s)];
OperatorSeries[f_, var_ -> op_, d_][ψ_] := Module[
{ser,as,ni,nf,t},
ser=List@@Series[f, {var,0, d}];
{as, ni, nf}=ser[[{3,4,5}]];
t=Nest[op, ψ, ni];
Expand[as.NestList[op,t,nf-ni-1]]
]

OperatorSeries[E^(-ad), ad->ad[y], 3][x]

MakeLieSeries[OperatorSeries[E^(-ad), ad->ad[x], 3][y], 5]

BCH[1]= LieSeries[⟨"x"⟩+⟨"y"⟩]

BCH[n_]:= Module[
{bch,t1,t2},
bch=Append[BCH[n-1],0];
t1=MakeLieSeries[⟨"y"⟩+OperatorSeries[E^(-ad), ad->ad[⟨"y"⟩], n-1][⟨"x"⟩], n];
t2=OperatorSeries[(1-E^(-ad))/ad, ad->ad[bch], n-1][EulerE[bch]];
bch+(t1-t2)/n
]

t1=MakeLieSeries[⟨"y"⟩+OperatorSeries[E^(-ad), ad->ad[⟨"y"⟩], 1][⟨"x"⟩], 2]

bch=Append[BCH[1],0];
t2=OperatorSeries[(1-E^(-ad))/ad, ad->ad[bch], 1][EulerE[bch]]

t1-t2

BCH[2]

BCH[3]

BCH[4]

BCH[5]

BCH[6]

BCH[7]

BCH[8]

BCH[15]

TimeUsed[]