Partitioning a set into lists

The Bell numbers count how many ways there are to partition N elements into sets. What if we want to partition the elements into ordered groups instead, as this Quora question asks?

For example, {a,b} {c,d,e} and {b,a} {c,d,e} are the same partitioning into sets, but not into ordered lists.

The Bell numbers obey a recurrence:

B(n+1) = sum k=0..n { binom(n,k) * B(k) }

which is justified in terms of removing the set containing the “first” (lowest-valued, say) element. k
indexes over the number of elements not in this set, there are binom(n,k ways to pick them, and B(k) ways to partition them.

We can easily turn this into a recurrence over ways to order the elements in each partition:

C(n+1) = sum k=0..n { (n+1-k)! * binom(n,k) * B(k) }

The set we removed (that contains the “first” element of the original set) is of size n+1−k, so it can be ordered in (n+k-1)! distinct ways, and each of these orderings can be paired with any of the available choices and orders of the remaining k elements.

I threw together some Python code to implement this recurrence:

# From https://gist.github.com/rougier/ebe734dcc6f4ff450abf
def binomial(n, k):
    """
    A fast way to calculate binomial coefficients by Andrew Dalke.
    See http://stackoverflow.com/questions/3025162/statistics-combinations-in-python
    """
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in xrange(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok // ktok
    else:
        return 0

def factorial( n ):
    if n == 0:
        return 1
    ret = 1
    for i in xrange( 2, n + 1 ):
        ret *= i
    return ret

def partitions_into_total_orders( m ):
    val = [None] * (m+1)
    val[0] = 1
    val[1] = 1
    for n in xrange( 1, m ):
        # Compute n+1
        tmp = 0
        for k in xrange( 0, n+1 ):
            tmp += binomial( n, k ) * factorial(n + 1 - k ) * val[k]
        val[n+1] = tmp
    return val

(The Haskell version would probably be prettier, as a self-recursive array definition.)

>>> pprint.pprint( partitions_into_total_orders(20) )
[1,
 1,
 3,
 13,
 73,
 501,
 4051,
 37633,
 394353,
 4596553,
 58941091,
 824073141,
 12470162233L,
 202976401213L,
 3535017524403L,
 65573803186921L,
 1290434218669921L,
 26846616451246353L,
 588633468315403843L,
 13564373693588558173L,
 327697927886085654441L]

This sequence is A000262 in the Online Encyclopedia of Integer Sequences, and more information can be found about it there. But, unlike the Bell Numbers, it doesn't appear to have a well-known name (though one or more of the references might name it.)

It seems trivial to extend the list beyond the 444 values currently in the OEIS:

>>> a = partitions_into_total_orders(500)
>>> a[445]
3549292370508805159816569987954285344396657510795554322174219200284938021552459305649764453676879331381726227308397349706332625735109322766384742689512837781933552271696452295507430453816017653770782548820800979477740762113630101472674590941070346886789347195166700719170015194143370749361014294622013831855459606619587788493906733145088876951867607064510581907335203726336222112549590988879213985659408433063349187438926065642124233265796247240308820380730636134217026702329240813645548435857866916318250327238530003440951684492429199855392402883015419785783782503556422217719938263845879893216604137311263061583644474420492077874479568291772742379482126859778320972094748833190234383500913296967990387469226679302739115758348855376993890250878022951237465510031943420668964594817338460644504177345168761851434514344256171704617556018164197594391584888906518228567907405660857519991265072822823355837407025410807672466016140913725137763011076901797240961561552844550260012058991623122900932503937295341L

However, one of the other formulas given in the OEIS, such as the recursion a_n=(2 * n - 1) * a_{n-1} - (n-1) * (n-2) * a{n-2} may be more efficient than my code above.