(* Mathematica program for A303831, coded by J.-F. Alcover 
 after Andrew Howroyd's PARI programs for A000088, 
 A126100 and A126122. *)

nmax = 30;
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, 2]];
a000088[n_] := Module[{s = 0}, Do[s += permcount[p]*2^edges[p],
 {p, IntegerPartitions[n]}]; s/n!];
g[n_, r_] := (s = 0; Do[s += permcount[p]*(2^(r*Length[p] + edges[p])),
 {p, IntegerPartitions[n]}]; s/n!);
seq[m_] := Sum[g[n - 1, 1] x^(n - 1), {n, 0, m}]/
 Sum[g[n - 1, 0] x^(n - 1), {n, 0, m}] + O[x]^m //
 CoefficientList[#, x] & // Prepend[#, 0] &;
A126100 = seq[nmax];
cross[u_, v_] := 
 Sum[GCD[u[[i]], v[[j]]], {i, 1, Length@u}, {j, 1, Length@v}];
a126122[n_] := 
 If[n < 2, 0, s = 0; 
 Do[s += permcount[
 p]*(2^(edges[p])*(2^cross[{1, 1}, p] + 2^cross[{2}, p])),
 {p, IntegerPartitions[n - 2]}]; s/(2 (n - 2)!)];
B[x_] = 2 Sum[a126122[n] x^n, {n, 0, nmax}];
G[x_] = Sum[a000088[n] x^n, {n, 0, nmax}];
c[x_] = Sum[A126100[[n + 1]] x^n, {n, 0, nmax}];
B[x]/G[x] - (c[x^2] + c[x]^2)/2 + O[x]^nmax //
 CoefficientList[#, x] & // Rest