OFFSET
0,6
COMMENTS
a(17) is the first term that differs from A001662.
LINKS
Alois P. Heinz, Table of n, a(n) for n = 0..446 (first 151 terms from G. C. Greubel)
R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, On the Lambert W Function, Advances in Computational Mathematics, (5), 1996, pp. 329-359.
R. M. Corless, D. J. Jeffrey and D. E. Knuth, A sequence of series for the Lambert W Function (section 2.2).
FORMULA
From Vladimir Kruchinin, Nov 11 2012: (Start)
a(n) = numerator(1/n!*(Sum_{u=2..n} stirling2(n,u)*(Sum_{k=1..u-1} ((u+k-1)!*Sum_{j=1..k} 2^(-u-j)/(k-j)!*Sum_{l=1..j} (-1)^(l)/((j-l)!)*Sum_{i=0..l} (l^(u+j-i-1))/((l-i)!*i!*(u+j-i-1)!)))+1/2)).
a(n) = numerator(1/n!*Sum_{k=0..n-1} (n+k-1)!*Sum_{j=0..k} ((-1)^(j)/(k-j)!*Sum_{i=0..j} ((1/i!)*Stirling1(n-i+j-1,j-i))/(n-i+j-1)!))*2^(n-j-1)), n>0, a(0)=1. (End)
a(n) = numerator(q(n)/n!) where q(n) = add_{k=0..n-1}(-1)^k*E2(n-1,k) if n>0 and 1 otherwise, E2 the second-order Eulerian numbers. - Peter Luschny, Nov 13 2012
a(n) := numerator(1/n!*Sum_{i=1..n} Stirling2(n,i)*A013703(i)/2^(2*i+1)). - Paolo Bonzini, Jun 23 2016
EXAMPLE
W(exp(x)) = 1 + (x-1)/2 + (x-1)^2/16 - (x-1)^3/192 - ...
MAPLE
a:= n-> numer(coeftayl(LambertW(exp(x)), x=1, n)):
seq(a(n), n=0..30); # Alois P. Heinz, Nov 08 2012
# For large n much faster is:
q := proc(n) if n=0 then 1 else add((-1)^k*combinat[eulerian2](n-1, k), k=0..n-1) fi end: A001662 := n -> numer(q(n)/n!):
seq(A001662(n), n=0..100): # Peter Luschny, Nov 13 2012
MATHEMATICA
CoefficientList[ Series[ ProductLog[ Exp[1+x] ], {x, 0, 22}], x] // Numerator (* Jean-François Alcover, Oct 15 2012 *)
a[0] = 1; a[n_] := 1/n!*Sum[(n+k-1)!*Sum[(-1)^(j)/(k-j)!*Sum[1/i!* StirlingS1[n-i+j-1, j-i]/(n-i+j-1)!, {i, 0, j}]*2^(n-j-1), {j, 0, k}], {k, 0, n-1}] // Numerator; Array[a, 30, 0] (* Jean-François Alcover, Feb 13 2016, after Vladimir Kruchinin *)
PROG
(Sage)
@CachedFunction
def eulerian2(n, k):
if k==0: return 1
if k==n: return 0
return eulerian2(n-1, k)*(k+1)+eulerian2(n-1, k-1)*(2*n-k-1)
def q(n):
return add((-1)^k*eulerian2(n-1, k) for k in (0..n-1)) if n>0 else 1
A001662 = lambda n: (q(n)/factorial(n)).numerator()
[A001662(n) for n in (0..22)] # Peter Luschny, Nov 13 2012
(Maxima)
a(n):=num(if n=0 then 1 else 1/n!*(sum((n+k-1)!*sum(((-1)^(j)/(k-j)!*sum((1/i!*stirling1(n-i+j-1, j-i))/(n-i+j-1)!, i, 0, j))*2^(n-j-1), j, 0, k), k, 0, n-1))); /* Vladimir Kruchinin, Nov 11 2012 */
CROSSREFS
KEYWORD
sign,easy,frac
AUTHOR
Paolo Bonzini, Jun 23 2016
STATUS
approved