login
Exponentiation of e.g.f. for primes.
(Formerly M1785)
16

%I M1785 #34 Nov 27 2017 12:41:00

%S 1,2,7,31,162,973,6539,48410,390097,3389877,31534538,312151125,

%T 3271508959,36149187780,419604275375,5100408982825,64743452239424,

%U 856157851884881,11768914560546973,167841252874889898,2479014206472819045,37860543940437797897

%N Exponentiation of e.g.f. for primes.

%C From _Tilman Neumann_, Oct 05 2008: (Start)

%C a(n) is also given by

%C - substituting the primes (A000040) into (the simplest) Faa di Bruno's formula, or

%C - the complete Bell polynomial of the first n prime arguments, or

%C - computing n-th moments from the first n primes as cumulants

%C The examples show that the coefficients of the prime power products are just A036040/A080575 (these are just rearrangements of the same coefficients). Moreover, the prime products of the additional terms span the whole space of natural numbers, thus what we see here is a reordering of the natural numbers! (End)

%D N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).

%H Alois P. Heinz, <a href="/A007446/b007446.txt">Table of n, a(n) for n = 0..500</a>

%F E.g.f.: exp(Sum_{k>=1} prime(k)*x^k/k!). - _Ilya Gutkovskiy_, Nov 26 2017

%e From _Tilman Neumann_, Oct 05 2008: (Start)

%e Let p_i denote the i-th prime A000040(i). Then

%e a(1)=2 = 1*p_1

%e a(2)=7 = 1*p_2 + 1*p_1^2

%e a(3)=31 = 1*p_3 + 3*p_2*p_1 + 1*p_1^3

%e a(4)=162= 1*p_4 + 4*p_3*p_1 + 3*p_2^2 + 6*p_2*p_1^2 + 1*p_1^4

%e a(5)=973= 1*p_5 + 5*p_4*p_1 + 10*p_3*p_2 + 10*p_3*p_1^2 + 15*p_2^2*p_1 + 10*p_2*p_1^3 + 1*p_1^5

%e (End)

%p a:= proc(n) option remember; `if`(n=0, 1,

%p add(binomial(n-1, j-1)*ithprime(j)*a(n-j), j=1..n))

%p end:

%p seq(a(n), n=0..30); # _Alois P. Heinz_, Mar 18 2015

%t a[n_] := a[n] = If[n==0, 1, Sum[Binomial[n-1, j-1]*Prime[j]*a[n-j], {j, 1, n}]]; Table[a[n], {n, 0, 30}] (* _Jean-François Alcover_, Mar 30 2015, after _Alois P. Heinz_ *)

%t Table[Sum[BellY[n, k, Prime[Range[n]]], {k, 0, n}], {n, 0, 20}] (* _Vladimir Reshetnikov_, Nov 09 2016 *)

%o (MuPAD)

%o completeBellMatrix := proc(x,n)

%o // x - vector x[1]...x[m], m>=n

%o local i,j,M;

%o begin

%o M:=matrix(n,n): // zero-initialized

%o for i from 1 to n-1 do

%o M[i,i+1]:=-1:

%o end_for:

%o for i from 1 to n do

%o for j from 1 to i do

%o M[i,j] := binomial(i-1,j-1)*x[i-j+1]:

%o end_for:

%o end_for:

%o return (M):

%o end_proc:

%o completeBellPoly := proc(x, n)

%o begin

%o return (linalg::det(completeBellMatrix(x,n))):

%o end_proc:

%o x:=[2,3,5,7,11,13,17,19,23,29]:

%o for i from 1 to 10 do print(i,completeBellPoly(x,i)): end_for:

%o // _Tilman Neumann_, Oct 05 2008

%Y Cf. A036040, A080575. - _Tilman Neumann_, Oct 05 2008

%K easy,nonn

%O 0,2

%A _N. J. A. Sloane_