BHC formula : convert polynomial representation into a nested commutators

Questions
,
1

Hi,

I’m trying to convert polynomial representation of z=log(exp(x)exp(y)) into nested commutators representation.
Polynomial representation can be obtained by Goldberg’s method.
(The code below is taken from Michael Weyraucha, Daniel Scholz, Computing the Baker–Campbell–Hausdorff series and the Zassenhaus product)

g[1] = 1;

g[s_] := g[s] = Expand[1/s*D[t*(t - 1)*g[s - 1], t]];
c[w_] := c[w] = Module[
    {m, m1, m2, k},
    m = Length[w];
    m1 = Floor[m/2];
    m2 = Floor[(m - 1)/2]; 
    Integrate[t^m1*(t - 1)^m2*Product[g[w[[k]]], {k, m}], {t, 0, 1}]
     ];
BCH[n_Integer, alph_List] := Module[
   {p},
   p = Flatten[Permutations /@ IntegerPartitions[n], 1]; 
   Plus @@ (c[Sort[#]]*(words[#, alph] -
           (-1)^n*words[#, Reverse[alph]]) & /@ p) // Expand];
words[p_List, alph_List] := StringJoin @ (ConstantArray @@@
     Partition[Riffle[p, alph, {1, 2*Length[p], 2}], 2]);

For example:

In[2458]:= BCH[3, {"x", "y"}] 
Out[2458]= ("xxy")/12 - ("xyx")/6 + ("xyy")/12 + ("yxx")/12 -("yxy")/6 + ("yyx")/12

The output is a polynomial of the form: P(x,y) = Sum a(x1,x2,..,xk) x1 x2 … xk with xi = x or y
In the same article (eq.14) the map that converts polynomial representation into commutators is given:
F(P) = Sum a(x,y,x3,x4,…,xk)/(nx(x3,x4,…,xk)+1) [[…[[ [x,y],x3],x4],…,],xk],
here nx(x3,x4,…,xk) is a number of x in a monomial after xy, for e.g. n(xy)=1, n(xyx)=2, …

Here is the algorithm I need to code:

  1. for given n find BCH[n, {“x”, “y”}]
  2. remove all monomials (words) not starting with xy
  3. for xy… monomials find a(x1,x2,..,xk) and nx(x3,x4,…,xk)
  4. build [[…[[ [x,y],x3],x4],…,],xk]

Example for z3:
1) (“xxy”)/12 - (“xyx”)/6 + (“xyy”)/12 + (“yxx”)/12 - (“yxy”)/6 +(“yyx”)/122) - (“xyx”)/6 + (“xyy”)/12
3) - (“xyx”)/6: a(x,y,x)=-1/6 and n(x)=1
        (“xyy”)/12: a(x,y,y)=1/12 and n(y)=0
4) a(x,y,x)/(n(x)+1) [[x,y],x] + a(x,y,y)/(n(y)+1) [[x,y],y] = -1/12[[x,y],x]+1/12[[x,y],y]

Any tips on how to implement this map in Mathematica?
Thanks in advance,
I.M.

2

Here is what I’ve managed to do

(* generate poly *)
poly = BCH[3, {"x", "y"}] ;
(* create list1 and list2 for a(x1,x2,...,xk) and x1x2...xk *)
list1 = list2 = Table[0, {i, 1, Length[poly]}] ;
Do[
  {
   list1[[i]] = poly[[i]] [[1]],
   list2[[i]] = poly[[i]] [[2]]
   },
  {i, 1, Length[poly]}
  ];
(* select words starting with xy *)
tmp = StringCases[list2, StartOfString ~~ "xy" ~~ ___] ;
(* find xy words positions and take corr. a's *)
A = list1[[Take[Position[tmp, x_String] // Flatten, {1, -1, 2}]]];
(* shape arr. for words *)
WORDS = tmp // Flatten;
(* build n(x3x4...)+1 array *)
B = StringCount[WORDS, "x"];
(* build the answer *)
answer = Table[0, {i, 1, Length[WORDS]}];
Do[
 {
  answer[[i]] = A[[i]]/B[[i]] Apply[NonCommutativeMultiply,
     ToExpression[StringCases[WORDS[[i]], Repeated[_, 1]]], {0}]
  }, {i, 1, Length[WORDS]}
 ]
Total[answer]

I’ve used NonCommutativeMultiply to shape the answer into the form -(1/12) x ** y ** x + x ** y ** y/12,
but I still need to build correspondence between ** and commutators, i.e. to change x**y**x into [[[x,y],x],…] or ((x <my operation> y) < my operation>) < my operation> … )

Can I redefine ** to be non associative or define non associative (acting from left to right) <my operation> and replace ** with it, or replace ** with <my operation> but correctly group the answer with parentheses?

I.M.

3

Here is a complete (not optimized) solution to the problem stated.I’ve checked terms up to z5, higher terms also seem to be correct from numerical experiments.

(* M.Weyrauch,D.Scholz/Computer Physics Communications 180 (2009) \
1558\[Dash]1565 *)
(* Goldberg's method *)
(* polynomial generator *)
g[1] = 1;
g[s_] := g[s] = Expand[1/s*D[t*(t - 1)*g[s - 1], t]];
c[w_] := c[w] = Module[
    {m, m1, m2, k},
    m = Length[w];
    m1 = Floor[m/2];
    m2 = Floor[(m - 1)/2]; 
    Integrate[t^m1*(t - 1)^m2*Product[g[w[[k]]], {k, m}], {t, 0, 1}]
     ];
BCH[n_Integer, alph_List] := Module[
   {p},
   p = Flatten[Permutations /@ IntegerPartitions[n], 1]; 
   Plus @@ (c[Sort[#]]*(words[#, alph] -
           (-1)^n*words[#, Reverse[alph]]) & /@ p) // Expand];
words[p_List, alph_List] := StringJoin @ (ConstantArray @@@
     Partition[Riffle[p, alph, {1, 2*Length[p], 2}], 2]);
(* polynomial converter *)
CircleTimes[a_, b_, c__] = CircleTimes[CircleTimes[a, b], c];
PolyToCom[n_Integer] := Module[
   {
    POLY, LISTA, LISTB, TMP, A, B, WORDS, ANSWER, i
    },
   POLY = Apply[List, BCH[n, {"x", "y"}]];
   LISTA = POLY[[All, 1]];
   LISTB = POLY[[All, 2]];
   TMP = StringCases[LISTB, StartOfString ~~ "xy" ~~ ___] ;
   WORDS = TMP // Flatten;
   A = LISTA[[Take[Position[TMP, x_String] // Flatten, {1, -1, 2}]]];
   B = StringCount[WORDS, "x"];
   ANSWER = A;
   Do[
    {
     ANSWER[[i]] = 
      A[[i]]/B[[i]] Apply[CircleTimes, 
        ToExpression[StringCases[WORDS[[i]], Repeated[_, 1]]], {0}]
     }, {i, 1, Length[WORDS]}
    ];
   Total[ANSWER]
   ];
(* EG *)
PolyToCom[2]
PolyToCom[3]
PolyToCom[4]
PolyToCom[5]

I.M.