(*  Mathematica program for A000677,
  adapted from N. J. A. Sloane's Maple program, 
  by Jean-François Alcover, Aug 23 2022. *)
  
did[m_, n_]:=If[Mod[m, n]==0, 1, 0];

EULER[seq_List] := Module[{coeff, final = {}},
coeff = Table[Sum[d*did[i, d]*seq[[d]], {d, 1, i}],
   {i, 1, Length[seq]}];
Do[AppendTo[final, (coeff[[i]]+Sum[coeff[[d]]*final[[i-d]],
   {d, 1, i-1}])/i], {i, 1, Length[seq]}];
final];

M = 102; Mh = 102; Mp = 20;

R[_] = 0; R[0] = x; R[1] = x^2/(1-x);
r[1, 1] = 0;
Do[r[1, n] = 1, {n, 2, M}];
Do[s[1, n] = 1, {n, 1, M}];
Do[
   a = Table[s[h-1, n], {n, 1, M}];
   b = Join[{1}, EULER[a]];
   Do[s[h, n] = b[[n]], {n, 1, M}];
   Do[r[h, n] = s[h, n]-s[h-1, n], {n, 1, M}];
   R[h] = Sum[r[h, n]*x^n, {n, 1, M}],
{h, 2, Mh}];

t1 = Sum[(1/2)*(R[h]^2+(R[h]/.x->x^2)), {h, 0, Mh}];
t2 = Series[t1, {x, 0, 2*Mh}];
CoefficientList[t2, x]