(* Mathematica computation of A005142 from A318869 and A318870, 
by J.-F. Alcover, after Andrew Howroyd and Alois P. Heinz. *)

nmax = 30;
mob[m_, n_] := If[Mod[m, n] == 0, MoebiusMu[m/n], 0];
EULERi[b_] := Module[{a, c, i, d}, c = {}; For[i = 1, i <= Length[b],
 i++, c = Append[c, i*b[[i]] -
 Sum[c[[d]]*b[[i - d]], {d, 1, i - 1}]]]; a = {};
 For[i = 1, i <= Length[b], i++, a = Append[a, (1/i)*
 Sum[mob[i, d]*c[[d]], {d, 1, i}]]]; Return[a]];
permcount[v_] := Module[{m = 1, s = 0, k = 0, t},
 For[i = 1, i <= Length[v], i++, t = v[[i]]; k =
 If[i > 1 && t == v[[i - 1]], k + 1, 1]; m *= t*k; s += t]; 
 s!/m];
edges[v_] := Sum[GCD[v[[i]], v[[j]]], {i, 2, Length[v]},
 {j, 1, i - 1}] + Total@Quotient[v + 1, 2]
b[n_] := (s = 0; Do[s += permcount[p]*2^edges[p], {p,
 IntegerPartitions[n]}]; s/n!);
A318869 = Join[{1}, EULERi[Array[b, nmax]]];

b[n_, i_] := b[n, i] =
 If[n == 0, {0}, If[i < 1, {}, Flatten@
 Table[Map[Function[{p}, p + j*x^i],
 b[n - i*j, i - 1]], {j, 0, n/i}]]];
g[n_, k_] := g[n, k] =
 Sum[Sum[2^Sum[Sum[GCD[i, j]*Coefficient[s, x, i]*
 Coefficient[t, x, j], {j, 1, Exponent[t, x]}],
 {i, 1, Exponent[s, x]}]/Product[i^Coefficient[s, x, i]*
 Coefficient[s, x, i]!, {i, 1, Exponent[s, x]}]/
 Product[i^Coefficient[t, x, i]*Coefficient[t, x, i]!,
 {i, 1, Exponent[t, x]}], {t, b[n + k, n + k]}], {s, b[n, n]}];
A[n_, k_] := g[Min[n, k], Abs[n - k]];
b[d_] := Sum[A[n, d - n], {n, 0, d}];
A318870 = Join[{1}, EULERi[Array[b, nmax]]];

a[0] = 1; a[n_] := a[n] = If[OddQ[n], A318870[[n + 1]]/2, (a[n/2] +
 A318869[[n/2 + 1]] + A318870[[n + 1]] - A318870[[n/2 + 1]])/2];
a /@ Range[0, nmax]