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

DBTR::usage = " DBTR[a,db,c,e,order]  finds RKN principal \
truncation error of order order+1 for DB. "
BTR::usage = " BTR[a,b,c,e,order]  finds RKN principal truncation \
error of order order+1 for B. "
DBOC::usage = "  DBOC[a,db,c,e,order] finds RKN order conditions of \
orders 1 to order. "
BOC::usage = "  DBOC[a,b,c,e,order] finds RKN order conditions of \
orders 1 to order. "
Begin["`Private`"];
Clear["RKN`Private`*"];

DBTR[aa_, dbb_, cc_, ee_, orderr_] := 
  1/S1[orderr + 1]*(T1[aa, dbb, cc, ee, orderr + 1] - G1[orderr + 1]);
(*--------------------------------------------------------------------------*)


BTR[aa_, bb_, cc_, ee_, orderr_] := 
  1/S0[orderr + 1]*(T0[aa, bb, cc, ee, orderr + 1] - G0[orderr + 1]);
(*--------------------------------------------------------------------------*)


DBOC[aa_, dbb_, cc_, ee_, orderr_] := 
  Table[Map[First, Split[Sort[(T1[aa, dbb, cc, ee, j] - G1[j])]]], {j,
     1, orderr}];
(*--------------------------------------------------------------------------*)


BOC[aa_, bb_, cc_, ee_, orderr_] := 
  Table[Map[First, Split[Sort[(T0[aa, bb, cc, ee, j] - G0[j])]]], {j, 
    1, orderr}];
(*--------------------------------------------------------------------------*)


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;

(*--------------------------------------------------------------------------*)

s = 0;  (* change to s=1 if A.e=c^2/2 holds *)
T[a_, c_, e_, 1] = {c};
T[a_, c_, e_, 2] = {a.e};
G[1] = {1};
G[2] = {1/2};
S[1] = {1}; S[2] = {1};
Switch[s, 1, T[a_, c_, e_, 2] = {c^2}; G[2] = {1}; S[2] = {2};];


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

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

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

S[order_] := 
  S[order] = 
   Module[{temp}, 
    temp = 
     Map[Combinations2[MapIndexed[ff, S[#[[1]]]], #[[2]]] &, 
      Map[RunLengthEncode[#] &, 
       IntegerPartitions[order - 2], {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]];

T1[a_, db_, c_, e_, 1] = {db.e}; G1[1] = {1};
T1[a_, db_, c_, e_, order_] := 
 T1[a, db, c, e, order] = 
  Module[{temp}, 
   temp = 
    Map[Combinations2[T[a, c, e, #[[1]]], #[[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[(db.#) &, temp, {1}]]

G1[order_] := 
  G1[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];
S1[order_] := S[order + 1];

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

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

S0[order_] := S[order];
(*--------------------------------------------------------------------------*)

End[];
EndPackage[];