BeginPackage["RKOtrees`"];
Clear["RKOtrees`*"]
RKO::usage = " RKO[a,b,c,e,taxis] finds RKO principal truncation \
error terms of order taxis."
Begin["`Private`"];
Clear["Rooted\[OpenCurlyQuote]Private`*"];
RKO[a0_, b0_, c0_, e0_, taxis_] := 
  If[taxis > 1, 
   1/Sym[taxis]*(Taf[a0, b0, c0, e0, taxis] - Gama[taxis]), {b0 . 
      e0 - 1}];
rle[u_List] := (Through[{First, Length}[#1]] &) /@ Split[u];
Comb[lsta_, n1_] := 
  Module[{j}, 
    Table[Map[Prepend[#, lsta[[j]]] &, 
      Flatten[Comb[lsta, 
         n1 - 1][[Array[Identity, Length[lsta] - j + 1, j]]], 1
       ], {1}
      ], {j, 1, Length[lsta]}
     ]] /; (n1 > 1);
Comb[lsta_, 1] := Comb[lsta, 1] = Map[{{#}} &, lsta];
Comb2[lsta_, n1_] := 
  Apply[Times, Flatten[Comb[lsta, n1], 1], {1}] /; (n1 > 1);
Comb2[lsta_, 1] := lsta;
(*--------------------------------------------------------------------------*)
(*the taxis conditions of taxis i are 1/Sym[i]*(Taf[i]-1/Gama[i])*)
Tau[a_, c_, e_, 1] = {c, a . e}; Gama[1] = {1, 1}; Sym[1] = {1, 1};
Gama[taxis_] := 
  Gama[taxis] = 
   Module[{dum}, 
    dum = Map[Comb2[Gama[#[[1]]], #[[2]]] &, 
      Map[rle[#] &, IntegerPartitions[taxis - 1], {1}], {2}];
    dum = Apply[Times, dum, {3}];
    dum = Map[Prepend[#, Times] &, dum, 1];
    dum = Apply[Outer, dum, {1}];
    dum = Flatten[dum];
    dum = (1/taxis)*dum];
Tau[a_, c_, e_, taxis_] := 
  Tau[a, c, e, taxis] = 
   Module[{dum}, 
    dum = Map[Comb2[Tau[a, c, e, #[[1]]] /. ax -> a, #[[2]]] &, 
      Map[rle[#] &, IntegerPartitions[taxis - 1], {1}], {2}];
    dum = Map[CL[#] &, dum, {3}];
    dum = Apply[Exoter, dum, {1}];
    dum = Flatten[dum, 1];
    dum = dum /. CL[each_] -> each;
    dum = Map[(ax . #) &, dum, {1}]];
Taf[a_, b_, c_, e_, 1] = {b . e};
Taf[a_, b_, c_, e_, taxis_] := Tau[a, c, e, taxis] /. ax -> b;
Sym[taxis_] := 
  Sym[taxis] = 
   Module[{dum}, 
    dum = Map[Comb2[MapIndexed[f0, Sym[#[[1]]]], #[[2]]] &, 
      Map[rle[#] &, IntegerPartitions[taxis - 1], {1}], {2}];
    dum = 
     dum /. {f0[u_, v_]^q_ -> Factorial[q]*u^q, f0[u_, v_] -> u};
    dum = Apply[Exoter, dum, {1}];
    dum = Flatten[dum, 1]];
Exoter[lstas__] := Flatten[Outer[Times, lstas], Length[{lstas}] - 1];
End[];
EndPackage[];