Clear[FestoonBackwardsDivision,FestoonSequence,FestoonTreeSymbolic,FestoonTreeSemiSymbolic,FestoonToTree]

FestoonBackwardsDivision[d_Integer,n_Integer]/;
    d>=2&&n>=1 :=Module[{R,epsilon,q},
    R=Mod[-n-1,d];
    epsilon=(d-1)+(d^2-1)R;
    If[epsilon<n,
      (
        Sow[R,FestoonCoefficient];
        (n-epsilon)/d
        ),
      (
        q=If[n==1,-1,Mod[(n-d)/(d-1),d+1]];
        R=(n-d)/(d^2-1)-q/(d+1);
        Sow[{q,R},FestoonCoefficient];
        0
        )
      ]
    ]

FestoonSequence[d_Integer,n_Integer]/;d>=2&&n>=1:=
  Reap[NestWhile[FestoonBackwardsDivision[d,#]&,
        n,#=!=0&]]/.{0,{coefficients_}}:>coefficients

FestoonTreeSymbolic[d_,n_]:=Festoon@@Module[{seq,l,q,r},
      seq=FestoonSequence[d,(d-1)n+1];
      l=Length[seq]-1;
      {q,r}=seq[[-1]];
      Append[
        MapThread[
          Join[
              Table[CompleteTreeSymbolic[d,#1],{If[#1>0,d-1-#2,0]}],
              Table[CompleteTreeSymbolic[d,#1+2],{#2}]
              ]&,
          {Range[0,l-1],Drop[seq,-1]}],
        If[{q,r}=={-1,0},
          Table[CompleteTreeSymbolic[d,l-1],{d}],
          Join[
            Table[CompleteTreeSymbolic[d,l],{d-q-r}],
            Table[CompleteTreeSymbolic[d,l+1],{q}],
            Table[CompleteTreeSymbolic[d,l+2],{r}]
            ]
          ]
        ]
      ]

FestoonToTree[Festoon[l_List, other__List]]:=
  Tree@@Append[l,FestoonToTree[Festoon[other]]]
FestoonToTree[Festoon[l_List]]:=Tree@@l
Festoon/: u:(sigma0|sigma1|M0|M1)[f_Festoon, other___]:=Hold[u][[1,0]][FestoonToTree[f],other]


FestoonTreeSemiSymbolic[d_,n_]:=FestoonToTree[FestoonTreeSymbolic[d,n]]
FestoonTree[d_,n_]:=
  FestoonTreeSemiSymbolic[d,n]/.CompleteTreeSymbolic->CompleteTree

Clear[CompleteTreeSymbolic,CompleteTree]
CompleteTreeSymbolic/:
  u:(sigma0|sigma1|M0|M1)[CompleteTreeSymbolic[d_,h_/;h>0],other___]:=Module[{head,result},
  head=Hold[u][[1,0]];
  result=StoreCompleteTreeData[head,d,h,other];
  If[Head[result]===StoreCompleteTreeData,
    StoreCompleteTreeData[head,d,h,other]=
    head[Tree@@Table[CompleteTreeSymbolic[d,h-1],{d}],other],
    result]
]
CompleteTreeSymbolic[_,0]:=EmptyGraph
CompleteTree[d_,h_/;h>1]:=Tree@@Table[CompleteTree[d,h-1],{d}]
CompleteTree[_,1]:=Tree[]

Clear[sigma0,sigma1,sigma,rho]
sigma0[Tree[branches___]]:=Times@@Map[sigma,{branches}]
sigma1[Tree[branches___]]:=Times@@Map[sigma0,{branches}]
sigma[u_]:=sigma0[u]+sigma1[u]
rho[u_]:=sigma0[u]/sigma[u]
sigma0[EmptyGraph]=1
sigma1[EmptyGraph]=0

Clear[M0,M1,M,z0,z1,z,tau,Energy]
M0[Tree[branches___],x_]:=Collect[Times@@Map[M[#,x]&,{branches}],x]
M1[Tree[branches___],x_]:=Collect[
  x Together[Plus@@Map[tau[#,x]& ,{branches}] M0[Tree[branches],x]],x]
M[u_,x_]:=M0[u,x]+M1[u,x]
tau[u_,x_]:=M0[u,x]/M[u,x]
M0[EmptyGraph,_]:=0
M1[EmptyGraph,_]:=1 
Energy[u_]:=Plus @@ Abs[z /. Solve[(M[u, x] /. x -> -1/z^2) == 0]]
z0[u_]:=M0[u,1]
z1[u_]:=M1[u,1]
z[u_]:=M[u,1]