(* Mathematica program for A000676, adapted from N. J. A. Sloane's Maple program for A000677, 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}]; A000677[n_]:=SeriesCoefficient[t2, {x,0,n}]; Clear[a,b]; b[0] = 0; b[1] = 1; b[n_] := b[n] = Sum[d*b[d]*b[n-j], {j, 1, n-1}, {d, Divisors[j]}]/(n-1); A000055[n_] := If[n==0,1,b[n] - (Sum[b[k]*b[n-k], {k, 0, n}] - If[Mod[n, 2] == 0, b[n/2], 0])/2]; a[n_]:=A000055[n]-A000677[n]; Table[a[n],{n,0,200}]