(* 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]