%I #20 Jan 25 2023 06:09:41
%S 1,49,1897,69553,2515513,90663937,3264855049,117543378001,
%T 4231639039705,152339702576545,5484235568128681,197432536935184369,
%U 7107571838026381177,255872590744254526273,9211413307971174616393
%N Numerators of the higher order exponential integral constants alpha(k,4).
%C The higher order exponential integrals, see A163931, are defined by E(x,m,n) = x^(n-1)*Integral_{t>=x} E(t,m-1,n)/t^n for m >= 1 and n >= 1, with E(x,m=0,n) = exp(-x).
%C The series expansions of the higher order exponential integrals are dominated by the alpha(k,n) and the gamma(k,n) constants, see A090998.
%C The first Maple program uses the alpha(k,n) formula and the second the GF(z,n) to generate the alpha(k,n) coefficients in each column.
%C Appears to equal the numerator of the multiple harmonic (star) sum Sum_{1 <= k_1 <= ... <= k_n <= 3} 1/(k_1^2*...*k_n^2)). If true, then a(n) = numerator( 3/2 - 3/(5*4^n) + 1/(10*9^n) ). - _Peter Bala_, Jan 31 2019
%H J. W. Meijer and N. H. G. Baken, <a href="https://doi.org/10.1016/0167-7152(87)90041-1">The Exponential Integral Distribution</a>, Statistics and Probability Letters, Volume 5, No. 3, April 1987, pp. 209-211.
%F alpha(k,n) = (1/k) * Sum_{i=0..k-1} (Sum_{p=0..n-1}(p^(2*i-2*k))*alpha(i, n)) with alpha(0,n) = 1, k >= 0 and n >= 1.
%F alpha(k,n) = alpha(k,n+1) -alpha(k-1,n+1)/n^2.
%F GF(z,n) = product((1-(z/k)^2)^(-1), k = 1..n-1) = (Pi*z/sin(Pi*z))/(Beta(n+z,n-z)/Beta(n,n)).
%e a(k=0,n=4) = 1, a(k=1,4) = 49/36, a(k=2,4) = 1897/1296, a(k=3,4) = 69553/46656.
%p coln := 4; nmax := 15; kmax := nmax: k:=0: for n from 1 to nmax do alpha(k, n) := 1 od: for k from 1 to kmax do for n from 1 to nmax do alpha(k, n) := (1/k)*sum(sum(p^(-2*(k-i)), p=0..n-1)*alpha(i, n), i=0..k-1) od; od: seq(alpha(k, coln), k=0..nmax-1);
%p # End program 1
%p coln:=4; nmax1 := 16; for n from 0 to nmax1 do A008955(n, 0):=1 end do: for n from 0 to nmax1 do A008955(n, n) := (n!)^2 end do: for n from 1 to nmax1 do for m from 1 to n-1 do A008955(n, m) := A008955(n-1, m-1)*n^2 + A008955(n-1, m) end do: end do: m:=coln-1: f(m):=0: for n from 0 to m do f(m) := f(m) + (-1)^(n + m)*A008955(m, n)*z^(2*m-2*n) od: GF(z,coln) := m!^2/f(m): GF(z,coln):=series(GF(z,coln), z, nmax1);
%p # End program 2
%Y Cf. A163931 (E(x,m,n)), A090998 (gamma(k,n)).
%Y Cf. A163928 y A163929.
%Y a(k,1) = A000007(k)
%Y a(k,2) = A000012(k) = 1^k.
%Y a(k,3) = A002450(k+1)/A000302(k) with A000302(k) = 4^k.
%Y a(k,4) = A163927(k)/A009980(k) with A009980(k) = 36^k.
%Y The GF(z,n) lead to A008955.
%Y The denominators of a(1,n), n >= 2, lead to A007407.
%K easy,frac,nonn
%O 0,2
%A _Johannes W. Meijer_ and _Nico Baken_, Aug 13 2009, Aug 17 2009