BeginPackage[ "Trees16`"]; 
Clear[ "Trees16`*" ]

RKTrunc::usage = "  RK[a,b,c,e,order] finds RK principal truncation error of order order+1. "
RKCond::usage = "  RKCond[a,b,c,e,order] finds RK order conditions of orders 1 to order. "

Begin["`Private`"];
Clear[ "Trees16`Private`*" ];


RKTrunc[aa_,bb_,cc_,ee_,orderr_]:= 1/S[orderr+1]*(T[aa,bb,cc,ee,orderr+1]-G[orderr+1]);
RKCond[aa_,bb_,cc_,ee_,orderr_]:= Table[T[aa,bb,cc,ee,j]-G[j],{j,orderr,1,-1}];

RunLengthEncode[x_List] := (Through[{First, Length}[#1]] &) /@ Split[x];

Combinations[list_, num_] :=
 Module[{i},
        Table[Map[Prepend[#, list[[i]]]&,
                  Flatten[Combinations[list, num - 1]
                            [[Array[Identity, Length[list] - i + 1, i]]], 1
                         ], {1}
                 ],
              {i, 1, Length[list]}
             ]]/;  (num > 1) ; 
Combinations[list_, 1] := Combinations[list, 1] = Map[{{#}}&, list];

Combinations2[list_, num_] :=
 Apply[Times, Flatten[Combinations[list, num], 1], {1} ]/;  (num > 1); 
Combinations2[list_, 1] := list;

(*--------------------------------------------------------------------------*)
(* the order conditions of order i are    1/S[i]*(T[i]-1/G[i]) *)

TT[a_,c_,e_,1] = {c};  G[1] = {1};  

G[order_] := G[order] =
 Module[{temp},
        temp = Map[Combinations2[G[#[[1]]], #[[2]]]&, Map[RunLengthEncode[#] &, IntegerPartitions[order-1], {1}], {2}];
        temp = Apply[Times, temp, {3}];
        temp = Map[Prepend[#, Times]&, temp, 1];
        temp = Apply[Outer, temp, {1}];
        temp = Flatten[temp];
        temp = (1/order) * temp
       ];

TT[a_,c_,e_,order_] := TT[a,c,e,order] =
Module[{temp},
        temp = Map[Combinations2[TT[a,c,e,#[[1]] ]/.at->a, #[[2]]]&, Map[RunLengthEncode[#] &, IntegerPartitions[order-1], {1}], {2}];
        temp = Map[CoverList[#]&, temp, {3}];
        temp = Apply[MyOuter, temp, {1}];
        temp = Flatten[temp, 1];
        temp = temp /. CoverList[every_] -> every;
        temp = Map[(at . #)&, temp, {1}]
       ];
  T[a_,b_,c_,e_,1]= {b.e};     
  T[a_,b_,c_,e_,order_]:= TT[a,c,e,order] /. at->b;  

  
  S[1] = {1}; 
S[order_] := S[order] =
Module[{temp},
    temp = Map[Combinations2[MapIndexed[ff, S[#[[1]]]], #[[2]]] &, Map[RunLengthEncode[#] &, IntegerPartitions[order-1], {1}], {2}];
    temp = temp /. {ff[a_, b_]^p_ -> Factorial[p]*a^p, ff[a_, b_] -> a};
    temp = Apply[MyOuter, temp, {1}]; 
    temp = Flatten[temp, 1]    ];

MyOuter[lists__] := Flatten[Outer[Times, lists], Length[{lists}] - 1];

End[];
EndPackage[];