I still think that the given implementation could be improved. An alternative would be to work in the Symbolic Ring. def Bell_partial_polynomial(n,k,v): @cached_function def T(n, k): if k==0: return v[0]^n return sum(binomial(n-1,j-1)*v[j]*T(n-j,k-1) for j in (0..n-k+1)) return T(n,k) If one wants to see the representation as a multivariate polynomial one can define the variables x_0,.._x_{n+1} R = range(5) v = var(['x_'+str(i) for i in (0..n+1)]) for n in R: [Bell_partial_polynomial(n,k,v).expand() for k in (0..n)] [1] [x_0, x_1] [x_0^2, x_0*x_1 + x_2, x_1^2] [x_0^3, x_0^2*x_1 + 2*x_0*x_2 + x_3, x_0*x_1^2 + 3*x_1*x_2, x_1^3] [x_0^4, x_0^3*x_1 + 3*x_0^2*x_2 + 3*x_0*x_3 + x_4, x_0^2*x_1^2 + 5*x_0*x_1*x_2 + 3*x_2^2 + 4*x_1*x_3, x_0*x_1^3 + 6*x_1^2*x_2, x_1^4] Substituting x_0=0 one gets the more familiar form starting with x_1. for n in R: [Bell_partial_polynomial(n,k,v).subs(x_0=0) for k in (0..n)] [1] [0, x_1] [0, x_2, x_1^2] [0, x_3, 3*x_1*x_2, x_1^3] [0, x_4, 3*x_2^2 + 4*x_1*x_3, 6*x_1^2*x_2, x_1^4] This is the form best suited for the Stirling/Lah family of combinatorial numbers. For example to compute the unsigned Stirling numbers first kind A132393 one can write R = range(6) s = [0]+[factorial(i) for i in R] for n in R: [Bell_partial_polynomial(n,k,s) for k in (0..n)] [1] [0, 1] [0, 1, 1] [0, 2, 3, 1] [0, 6, 11, 6, 1] [0, 24, 50, 35, 10, 1] In any case this implementation seems to be more efficient than the loop in the current implementation for p in Partitions(n, length=k): for part, count in p.to_exp_dict().iteritems(): We give a last example which shows that there are other important examples which do /not/ reduce to integers. R = range(9) t = [2^(1-j) if is_odd(j) else 0 for j in (0..n+2)] for n in R: [Bell_partial_polynomial(n,k,t) for k in (0..n)] [1] [0, 1] [0, 0, 1] [0, 1/4, 0, 1] [0, 0, 1, 0, 1] [0, 1/16, 0, 5/2, 0, 1] [0, 0, 1, 0, 5, 0, 1] [0, 1/64, 0, 91/16, 0, 35/4, 0, 1] [0, 0, 1, 0, 21, 0, 14, 0, 1] These are the coefficients of the central factorials. http://mathworld.wolfram.com/CentralFactorial.html /////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////// Bell polynomials reconsidered. ============================= Consider the transformation T(n, k, x) = Sum_{j=0..n-k+1} K(n,j)*x[j]*T(n-j,k-1,f) if k != 0 else T(n, k, x) = x[0]^n. As kernels K(n,j) we consider here the two cases K(n,j) = binomial(n-1,j-1) and K(n,j) = number_of_partitions(n,j), where number_of_partitions(n,j) is the number of j-partition of n. If the kernel K is fixed than T can be seen as a sequence-to-triangle transformation, mapping the sequence x to the triangle T(x), defined by T(x)_{n,k} = T(n, k, x). Let us look at the first few rows of the generated triangles. Case K(n,j) = binomial(n-1,j-1): x[0]^0 x[0]^1, x[1] x[0]^2, x[1]*x[0]+x[2], x[1]^2 x[0]^3, x[1]*x[0]^2+2*x[2]*x[0]+x[3], x[1]^2*x[0]+3*x[2]*x[1], x[1]^3 x[0]^4, x[1]*x[0]^3+3*x[2]*x[0]^2+3*x[3]*x[0]+x[4], x[1]^2*x[0]^2+5*x[1]*x[2]*x[0] +4*x[3]*x[1]+3*x[2]^2, x[1]^3*x[0]+6*x[2]*x[1]^2, x[1]^4 Case K(n,j) = number_of_partitions(n,j): x[0]^0 x[0]^1, x[1] x[0]^2, x[1]*x[0]+2*x[2], x[1]^2 x[0]^3, x[1]*x[0]^2+2*x[2]*x[0]+3*x[3], x[1]^2*x[0]+4*x[2]*x[1], x[1]^3 x[0]^4, x[1]*x[0]^3+3*x[2]*x[0]^2+4*x[3]*x[0]+5*x[4], x[1]^2*x[0]^2 +5*x[1]*x[2]*x[0]+7*x[3]*x[1]+6*x[2]^2, x[1]^3*x[0]+7*x[2]*x[1]^2, x[1]^4 Note that these triangles have the offset (0,0), i. e. the enumeration is T(0,0) T(1,0) T(1,1) T(2,0) T(2,1) T(2,2) T(3,0) T(3,1) T(3,2) T(3,3) T(4,0) T(4,1) T(4,2) T(4,3) T(4,4) This is important to note because in the older literature the first column T(n,0) often is ignored and the triangle starts at T(1,1). This happens in particular in the important special cases which arise when x[0] = 0. Then the above two triangles reduce to, in the case K(n,j) = binomial(n-1,j-1): 1 0, x[1] 0, x[2], x[1]^2 0, x[3], 3*x[2]*x[1], x[1]^3 0, x[4], 4*x[3]*x[1] + 3*x[2]^2 , 6*x[2]*x[1]^2 , x[1]^4 and in the case K(n,j) = number_of_partitions(n,j) to: 1 0, x[1] 0, 2*x[2], x[1]^2 0, 3*x[3], 4*x[2]*x[1], x[1]^3 0, 5*x[4], 7*x[3]*x[1] + 6*x[2]^2 , 7*x[2]*x[1]^2 , x[1]^4 Let's get numerical now and look for instance at the indicator sequence of the positive numbers, 1 - 0^n (in other words, let's set x[0] = 0 and x[k] = 1 for all k>0). Then the first triangle becomes 1; 0, 1; 0, 1, 1; 0, 1, 3, 1; 0, 1, 7, 6, 1; 0, 1, 15, 25, 10, 1; 0, 1, 31, 90, 65, 15, 1; which is A048993, the triangle of the Stirling set numbers. The second triangle becomes 1; 0, 1; 0, 2, 1; 0, 3, 4, 1; 0, 5, 13, 7, 1; 0, 7, 30, 30, 10, 1; 0, 11, 76, 119, 65, 14, 1; a triangle which has not yet been explored. But let's now switch our point of view. Until now we understood x as a sequence x[0], x[1], x[2],.... Now we fix the kernel binomial(n-1,j-1), write x_n instead of x[n] and understand x_n as a variable in some ring. Then an expression like x_1*x_0^3+3*x_2*x_0^2+3*x_3*x_0+x_4 (see the fifth line in the first triangle above) becomes an element of a polynomial ring in n-k+1 variables. These polynomials are well known and famous: they are called the Bell polynomials in the special case where x_0 is substituted by 0 and the first column is disregarded. x_1 x_2, x_1^2 x_3, 3*x_2*x_1, x_1^3 x_4, 4*x_3*x_1 + 3*x_2^2, 6*x_2*x_1^2, x_1^4 This is the definition as given for example by L. Comtet in Advanced Combinatorics, page 135. However this definition is not satisfactory. For example substituting 1 for the x_n in this beheaded triangle we get 1; 1, 1; 1, 3, 1; 1, 7, 6, 1; 1, 15, 25, 10, 1; 1, 31, 90, 65, 15, 1; which is A008277, the triangle of the trimmed Stirling numbers of the second kind. (In passing note that the first comment in A008277 is not correct. It says: "Also known as Stirling set numbers and written {n,k}." However the notation {n,k} and the term 'Stirling set numbers' were introduced by D. E. Knuth explicitly to denote the full triangle A048993, not the trimmed Stirling numbers.) Of course there are many other instances where the reduced Bell polynomials as defined by Comtet do not comply anymore to the way important numbers nowadays are derived from the Bell polynomials. The reader might check this for example for the case x[n] = 1!,2!,3!,.. (Lah numbers) or x[n] = 0!,1!,2!,3!,.. (signless Stirling numbers of the first kind). To sum up: We strongly suggest to use the polynomials with the full set of variables x_0, x_1, ... and the enumeration of the generated number triangles starting at (0,0). This gives, among other important numbers, the Stirling set numbers, the Stirling cycle numbers and the signless Lah numbers as used today by the majority of writers. --- With Maple one can nicely generate the above triangles with four lines of code: T := proc(n,k,x) option remember; `if`(k=0,x[0]^n, add(binomial(n-1,j-1)*x[j]*T(n-j,k-1,x),j=0..n-k+1)) end: for n from 0 to 6 do seq(simplify(T(n,k,x)),k=0..n) od; for n from 0 to 6 do seq(simplify(subs(x[0]=0,T(n,k,x))), k=0..n) od;