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