OFFSET
1,6
COMMENTS
Coefficients of (D^k f)(g(t))*(D g(t))^k1*(D^2 g(t))^k2*... in the Faa di Bruno formula for D^m(f(g(t))) where k = k1 + k2 + ..., m = 1*k1 + 2*k2 + ....
Number of set partitions whose block sizes are the prime indices of n (i.e., the integer partition with Heinz number n). - Gus Wiseman, Sep 12 2018
LINKS
Alois P. Heinz, Table of n, a(n) for n = 1..20000
Eric Weisstein's World of Mathematics, Bell Polynomial
Eric Weisstein's World of Mathematics, Faà di Bruno's Formula
FORMULA
For n = p1^k1*p2^k2*... where 2 = p1 < p2 < ... are the sequence of all primes, a(n) = a([k1,k2,...]) = (k1+2*k2+...)!/((k1!*k2!*...)*(1!^k1*2!^k2*...)).
a(2*prime(n)) = n + 1, for n > 1. See A065475. - Bill McEachen, Oct 11 2023
EXAMPLE
The a(6) = 3 set partitions of type (2,1) are {{1},{2,3}}, {{1,3},{2}}, {{1,2},{3}}. - Gus Wiseman, Sep 12 2018
MAPLE
with(numtheory):
a:= n-> (l-> add(i*l[i], i=1..nops(l))!/mul(l[i]!*i!^l[i],
i=1..nops(l)))([seq(padic[ordp](n, ithprime(i)),
i=1..pi(max(1, factorset(n))))]):
seq(a(n), n=1..100); # Alois P. Heinz, Feb 14 2020
MATHEMATICA
numSetPtnsOfType[ptn_]:=Total[ptn]!/Times@@Factorial/@ptn/Times@@Factorial/@Length/@Split[ptn];
Table[numSetPtnsOfType[If[n==1, {}, Flatten[Cases[FactorInteger[n], {p_, k_}:>Table[PrimePi[p], {k}]]]]], {n, 100}] (* Gus Wiseman, Sep 12 2018 *)
PROG
(PARI) a(n) = my(f=factor(n)); sum(k=1, #f~, primepi(f[k, 1])*f[k, 2])!/(prod(k=1, #f~, f[k, 2]!)*prod(k=1, #f~, primepi(f[k, 1])!^f[k, 2])); \\ Michel Marcus, Oct 11 2023
CROSSREFS
KEYWORD
nonn
AUTHOR
Max Alekseyev, Nov 07 2006
STATUS
approved