\\ PARI progam to compute terms for A338487.
\\ This is mostly a mash-up of other simpler programs with some extras.

\\ Multivariate Euler transform.
\\ For the purposes of this sequence could be replaced by single variable version.
EulerMT(u)={my(n=#u, p=x*Ser(u), vars=variables(p)); Vec(exp( sum(i=1, n, substvec(p + O(x*x^(n\i)), vars, [v^i|v<-vars])/i ))-1)}

\\ Replaces all variables v in p by v^e.
\\ For the purposes of this sequence could be replaced by subst(p, x, x^e)
\\ In this sequence, it is only ever used with e=2.
SubPwr(p,e)={my(vars=variables(p)); substvec(p, vars, [v^e|v<-vars])}


\\ 'G' function from A339065.
\\ Computes the number of multigraphs on n unlabeled vertices by number of edges with a given root.
\\ Returns g.f. in the variable x.
permcount(v) = {my(m=1, s=0, k=0, t); for(i=1, #v, t=v[i]; k=if(i>1&&t==v[i-1], k+1, 1); m*=t*k; s+=t); s!/m}
edges(v, t) = {prod(i=2, #v, prod(j=1, i-1, my(g=gcd(v[i], v[j])); t(v[i]*v[j]/g)^g )) * prod(i=1, #v, my(c=v[i]); t(c)^((c-1)\2)*if(c%2, 1, t(c/2)))}
G(n, x, r)={my(s=0); forpart(p=n, s+=permcount(p)*edges(concat(r, Vec(p)), i->1/(1-x^i))); s/n!}

\\ Combines code from A339230 and A339231.
\\ Computes number of oriented series-parallel networks with given edge and vertex weights.
\\ Edge weights are in vector 'u' and vertex weights given as a g.f. W.
EdgeWeighted(u, W=1)={
  my(Z = x*Ser(u), 
     p = Z+O(x^2)); 
     for(n=2, #u, 
         p=x*Ser(EulerMT(Vec(W*p^2/(1 + W*p) + Z)))); 
  Vec(p)
}

\\ Inverse of EdgeWeighted
\\ The algorithm is suboptimal - but performance is not an issue for this sequence.
\\ Algorithm just iteratively uses EdgeWeighted to grow inverse one term at a time.
\\ Edge weights are in vector 'v' and vertex weights in vector R.
InvertOriented(v, R)={
  my(u=vector(#v)); 
  u[1]=v[1]; 
  for(n=2, #v, u[n]=v[n]-EdgeWeighted(u[1..n], Ser(R[1..n]))[n]);
  u
}

\\ Combines code from A339284 and A339285 (adapted for achiral)
\\ Computes number of achiral series-parallel networks with given edge and vertex weights.
EdgeWeightedAch(ua, uo, vo, W=1)={
  my(n = #ua,
     ZA = x*Ser(ua), 
     ZO = x*Ser(uo[1..(n+1)\2]), 
     q = SubPwr(x*Ser(vo[1..(n+1)\2]), 2),
     W2 = SubPwr(W,2),
     s = SubPwr(ZO,2) + W2*q^2/(1+W2*q), 
     p = ZA + O(x^2),
     t = p); 
  for(n=1, n\2, 
     t = ZA + q*(W + W2*p);
     p = ZA + x*Ser(EulerMT(Vec(t+(s-SubPwr(t, 2))/2))) - t); 
  Vec(p + t - ZA);
}

\\ Inverse of EdgeWeightedAch
InvertAch(va, uo, vo, R)={
  my(ua=vector(#va)); 
  ua[1]=va[1];
  for(n=2, #va, ua[n]=va[n]-EdgeWeightedAch(ua[1..n], uo, vo, Ser(R[1..n]))[n]);
  ua
}


\\ Computes sequence of first n terms of A338487.
A338487seq(n)={
  my( A = O(x*x^n),
      g = G(2*n, x+A, []), \\ multigraphs 
      gr = G(2*n, x+A, [1])/g, \\ rooted multigraph
      rooted = Vec(gr),
      vo = Vec(-1+G(2*n, x+A, [1, 1])/(g*gr^2)),
      va = Vec(-1+G(2*n, x+A, [2])/(g*subst(gr, x, x^2))),
      uo = InvertOriented(vo, rooted),
      ua = InvertAch(va, uo, vo, rooted) );
  ((uo + ua)/2)[5..n];  
}

\\ End