@peterluschny also, could you explain how you arrived upon your 2011 Maple program in https://oeis.org/A090674? I already knew a bit about the series for W(x) about x=-1/e (or as I prefer, Cayley/Euler T(x) about x=1/e) in https://oeis.org/wiki/User:Natalia_L._Skirrow/Euler's_tree_function#expansion_about_the_vertical_point, but have no idea how it relates to.
Here also is a SymPy translation of said program (I have found that ring_series operations are much more efficient than the regular series ones (taking a matter of centiseconds instead of seconds); its rs_LambertW primitive doesn't allow expanding about other points so I inverted it directly); I will next try likewisely expanding the Laplace transform to get a similar formula in Gregories, but the fact that this process should coincide with that threw me off considerably.

from sympy import fraction,factorial2
from sympy.polys.domains import QQ
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_exp as exp,rs_nth_root,rs_series_reversion as invert

def A090674_A090675(N):
M=2*N+3
R,u,x = ring('u,x', QQ)
f = u*rs_nth_root(2*(1+exp(u,u,M)*(u-1))/u**2, 2,u,M) #sqrt(2*(1+exp(u)*(u-1)))
y = exp(invert(f,u,M,x),x,M).to_dict() #exp(1 + LambertW((x^2/2-1)/e))
return [-y.get((0,2*n+1))*factorial2(2*n+1) for n in range(1,N+1)]