(PARI)
dpermcount(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*=2*t*k; s+=2*t); s!/m}
S(n, x)={vector(n, n, if(n>1, sum(k=0, n, binomial(2*n-k, k)*2*n/(2*n-k)*x^k), 0))}
q(n, s)={my(t=0); if(n>1, forpart(p=n, t+=dpermcount(p)*prod(i=1, #p, s[p[i]]), [2, n])); t}
a(n)={my(p=q(n, S(n, x))); sum(i=0, poldegree(p), polcoeff(p, n-i)*(-1)^(n-i)*(2*i)!/(2^i*i!))} \\
Andrew Howroyd, Dec 18 2017
