Particiones y números de Stirling de segunda clase

Un número de Stirling de segunda clase es el número de formas de particionar un conjunto de n elementos en k subconjuntos no vacíos, denotándose como S2(n, k) o {nk}.

Se denota S2(n, k) para diferenciarlo de los números de Stirling de primera clase que, sin signo, lo denotábamos como S(n, k) o [nk], mientras que con signo eran s(n, k). En lo que sigue vemos las relaciones de las números de Stirling de primera y segunda clase con el factorial ascendente y descendente, también llamado símbolo de Pochhammer o rising factorial, denotado como x⁽ⁿ⁾ o xⁿ para el ascendente y falling factorial, denotado como (x)_n o xⁿ, para el descendente:

Los números de Stirling de segunda clase obedecen a la siguiente recurrencia:

que usando la denotación S2(n, k) y agregando los casos bases queda así:

En la aplicación Combine Tools, en la pestaña "Numéricos", implementamos este algoritmo:

// S2(n, k) = S2(n-1, k-1) + k * S2(n-1, k)
function stirlingS2(n, k) {
    if (n===k){
        return 1;
    } else if (n===0 || k===0 || n<k) {
        return 0;
    } else {
        return stirlingS2(n-1, k-1) + k * stirlingS2(n-1, k);
    }
}

obteniéndose esta tabla de primeros valores:

Se definen varias recurrencias para obtener los números de Stirling de segunda clase:

Y también una fórmula cerrada, es decir, que no es recursiva:

Veamos algunas identidades cuya demostración y comentarios se incluyen en cada desplegable.

Los números de Stirling de segunda clase se relacionan con los polinomios de Tourchard T_n(x) para un cierto n fijo de esta forma:

Y con el factorial descendente (o falling factorial) de esta forma:

Dado que un número de Stirling S2(n, k) es una partición de orden k-ésima, podemos contar todas las particiones de un conjunto de n elementos así:

En la aplicación lo implementamos con este simple algoritmo bell2(n) que usa el algoritmo anterior stirlingS2(n, k):

function bell2(n) {
    let sum = 0;
    for (let j=0; j<=n; j++) {
        sum += stirlingS2(n, j);
    }
    return sum;
}

Obtenemos estos valores para el número de particiones con los primeros valores de n:

Función generadora de los números de Stirling de segunda clase

La función generadora de los números de Stirling de segunda clase es con una variable:

y con dos variables:

Algoritmo para generar particiones usando números de Stirling de segunda clase

Usando la expresión [4] que vimos antes B(n) = ∑_j=0..n S2(n, j), vamos a implementar un algoritmo que genere las particiones, observando que el código implementa la recurrencia S2(n, k) = k × S2(n-1, k) + S2(n-1, k-1), donde la variable k está implícita en la longitud del subresultado res:

/* B(n) = sum_(j=0)^(n) S2(n, j)
   S2(n, j) = S2(n-1, j-1) + j × S2(n-1, j) */
function partition4(list=[], n=list.length, res=[], result=[]) {
    if (n===0) {
        //iter += 1;
        //if (withCopyCost) iter += list.length;
        result.push(copy2(res));
    } else {
        //iter += 1;
        let pivot = list[list.length-n];
        for (let j=0; j<res.length; j++) {
            //iter += 1;
            res[j].push(pivot);
            partition4(list, n-1, res, result);
            res[j].pop();
        }
        res.push([pivot]);
        partition4(list, n-1, res, result);
        res.pop();
   }
   return result;
}

Sobre el conjunto con tres elementos {a, b, c} obtenemos estas particiones:

que de forma compacta, tal como explicamos en el tema sobre particiones, también puede representarse así:

En cualquier caso la aplicación nos los devuelve como Arrays, aunque deben entenderse que son conjuntos. La vista de la salida en JSON sería esta:

[[["a","b","c"]],
[["a","b"],["c"]],
[["a","c"],["b"]],
[["a"],["b","c"]],
[["a"],["b"],["c"]]]

Observamos la simplicidad del algoritmo con respecto a los que vimos en los dos temas anteriores. El siguiente esquema expone la evolución del algoritmo para la ejecución de partition4(["a", "b", "c"]),

list=[a,b,c] res=[]
// FIRST CALL
partition(n=3, res=[])
n=3  pivot=a j=0..-1
// EMPTY LOOP
res.push([a]) => [[a]]
partition(n=2, res=[[a]])
    n=2  pivot=b  j=0..0
    j=0
        res[0].push(b) => [[a,b]]
        partition(n=1, res=[[a,b]])
            n=1 pivot=c  j=0..0
            j=0
                res[0].push(c) => [[a,b,c]]
                partition(n=0, res=[[a,b,c]]) ==>> [[a,b,c]]
                res[0].pop() = [[a,b]]
            res.push([c]) => [[a,b],[c]]
            partition(n=0, res=[[a,b],[c]])   ==>> [[a,b],[c]]
            res.pop() => [[a,b]]
        res.pop() => []
    res.push([b]) => [[a],[b]]
    partition(n=1, [[a],[b]])
        n=1 pivot=c j=0..1
        j=0
            res[0].push(c) => [[a,c],[b]]
            partition(n=0, res=[[a,c],[b]])   ==>> [[a,c],[b]]
            res[0].pop() => [[a],[b]]
        j=1
            res[1].push(c) => [[a],[b,c]]
            partition(n=0, res=[[a],[b,c]])   ==>> [[a],[b,c]]
            res[1].pop() => [[a],[b]]]
        res.push([c]) => [[a],[b],[c]]
        partition(n=0, [[a],[b],[c]])         ==>> [[a],[b],[c]]
        res.pop() => [[a],[b]]
    res.pop() => [[a]]
res.pop() => []

Con => indicamos el estado del subresultado res en cada momento, observando que cuando hacemos una llamada con n=0 llegamos al caso base, devolviendo una partición que ponemos en una columna a la derecha precedida por ==>>.

Para calcular el coste observamos el algoritmo, donde el bucle for (let j=0; j<res.length; j++) empieza con el array res=[] vacío con longitud 0 y puede llegar a tener la longitud máxima n, pues para una lista inicial [a, b, c] tendremos 1 partición con longitud 1: [a,b,c], 3 particiones con longitud 2: [[a,b],[c]], [[a,c],[b]], [[a],[b,c]] y 1 partición de longitud 3: [[a],[b],[c]]. Por lo tanto la longitud de res varía entre k=0 inicialmente y un máximo de k=n=3 en este ejemplo. Por lo tanto el primer planteamiento del coste sería el siguiente:

Esto es una recurrencia con dos variables, pues aparte de n tenemos k. Dada la dificultad saber cuánto vale k en cada caso, intentaremos obtener el coste de un forma indirecta. Por un lado si ejecutamos el algoritmo contando iteraciones obtendremos los valores Iter que figuran en la siguiente tabla, usando también los valores B(n) ya vistos en la tabla al final del apartado anterior, obteniendo las diferencias Iter - B(n):

Si prescindimos del primer valor para n=0, la secuencia a(n) = 1, 3, 8, 23, 75, 278, 1155, 5295 se encuentra en OEIS A024716 y obedece a la siguiente fórmula para n≥1

También obedece a la siguiente fórmula alternativa según OEIS, en este caso para n≥0:

Así que tomando [5] tenemos que

de donde obtenemos el coste del algoritmo partition4(list)

Si usamos [6] ajustándolo al mismo offset n≥1, para lo cual usamos C(n, i+1) llegamos a una segunda fórmula alternativa para el coste de partition4(list)

Donde B(n) son los números de Bell que se calcula con esta recurrencia

y S2(n, k) son los números de Stirling de segunda clase que tiene la siguiente recurrencia

El coste calculado puede compararse con el recuento de iteraciones en la aplicación obteniendo una tabla de valores, donde se separan con "#" comprobándose que son iguales:

En caso de coste de copia hay que agregar n B(n) pues hay B(n) subresultados y cada uno tiene una longitud n:

Particiones parciales y primer algoritmo para generarlas

En el apartado anterior vimos que no fue fácil desarrollar la recurrencia del coste debido a que el algoritmo partition4(list) hacía uso de la recurrencia de los números de Stirling de segunda clase S2(n, k) = S2(n-1, k-1) + k × S2(n-1, k). Si esta recurrencia depende de dos variables n, k, cualquier algoritmo que la use también tendrá implícito esas dos variables.

Intentaremos implementar directamente S2(n, k) = S2(n-1, k-1) + k × S2(n-1, k) en algo que podríamos denominar particiones parciales de un conjunto de n elementos tomados en k particiones. Para los primeros valores de S2(n, k) tenemos:

Desarrollando los resultados de esos primeros valores:

Este algoritmo lo denominaremos kPartition(list, k), en este caso la primera versión kPartition1:

// S2(n, k) = S2(n-1, k-1) + k × S2(n-1, k)
function kPartition1(list=[], k=0) {
    let result = [];
    let n = list.length;
    if (n===k) {
        iter += n;
        result.push([...list.map(v => [v])]);
    } else if (n===0 || k===0 || n<k) {
        iter += 1;
        result.push([]);
    } else if (k===1) {
        iter += n;
        result.push([[...list]]);
    } else {
        //iter += 1;
        let pivot = list[n-1];
        //iter += (n-1);
        let sublist = list.slice(0, n-1);
        let res = kPartition1(sublist, k-1);
        /* res.length = S2(n-1, k-1) */
        for (let r of res) {
            //iter += 1;
            r.push([pivot]);
            /* n+k = iter copy2(r) =  iter copyArray(r) */
            //if (withCopyCost) iter += n+k;
            result.push(copy2(r));
            r.pop();
        }
        res = kPartition1(sublist, k);
        /* res.length = S2(n-1, k) */
        for (let r of res) {
            //iter += 1;
            /* res.length = k */
            for (let rr of r) {
                //iter += 1;
                rr.push(pivot);
                /* n+k = iter copy2(r) =  iter copyArray(r) */
                //if (withCopyCost) iter += n+k;
                result.push(copy2(r));
                rr.pop();
            }
        }
    }
    return result;
}

En cuanto a las dos sentencias if (withCopyCost) iter += n+k, es necesario copiar los subresultados para pasarlos al array final result, pues en otro caso estaríamos pasasando la misma referencia. Cada subresultado parcial es como, por ejemplo, el subresultado [[a, b], [c]], por lo que hay que copiar con una profundidad de 2. Podemos usar una función copyArray(arr=[]):

function copyArray(arr=[]) {
    let arrCopy = [];
    for (let i=0; i<arr.length; i++) {
        iter += 1;
        arrCopy[i] = [];
        for (let j=0; j<arr[i].length; j++) {
            iter += 1;
            arrCopy[i].push(arr[i][j]);
        }
    }
    return arrCopy;
}

O bien usar copy2(arr=[[]]) que en la aplicación se utiliza para copiar particiones seleccionando algún método. Por defecto se usa slice. El método loop equivale a la anterior función copyArray():

function copy2(arr=[[]]) {
    if (copyMethod==="slice") {
        return arr.slice(0).map(v => v.slice(0));
    } else if (copyMethod==="map") {
        return arr.map(v => v.map(w => w));
    } else if (copyMethod==="concat") {
        return [].concat(arr).map(v => [].concat(v));
    } else if (copyMethod==="spread"){
        return [...arr].map(v => [...v]);
    } else if (copyMethod==="json") {
        return JSON.parse(JSON.stringify(arr));
    } else if (copyMethod==="loop") {
        let arrCopy = [];
        for (let i=0; i<arr.length; i++){
            arrCopy[i] = [];
            for (let j=0; j<arr[i].length; j++) {
                arrCopy[i].push(arr[i][j]);
            }
        }
        return arrCopy;
    } else {
        return arr;
    }
}

El planteamiento de la recurrencia del coste T'(n, k), considerando el coste de copia, es el siguiente:

Si hacemos n+k = 0 que equivale a ignorar el coste de copiar los subresultados, la recurrencia sin coste de copia queda así:

donde se observa claramente que S2(n-1, k-1) + S2(n-1, k) k = S2(n, k) con lo que nos queda por ahora así el coste de kPartition1(n, k) en la aplicación, pues no he podido encontrar algo mejor:

Podemos usar la fórmula cerrada (no recursiva) [3] para calcular los números de Stirling de segunda clase:

Considerando el coste de copiar los subresultados, que según vemos en el código del algoritmo se produce cuando insertamos los subresultados en el array result que contiene los resultados finales, hemos de suma el coste n+k que es del copiar un subresultado:

...
//if (withCopyCost) iter += n+k;
result.push(copy2(r));
...

La fórmula con coste de copia para kPartition1(list, k) es, entonces, lo mismo que [5] que vimos antes al plantear la recurrencia del coste del algortimo, pues como ya dije, no pude obtener una mejor fórmula:

Segundo algoritmo para generar particiones parciales

Un segundo algoritmo para generar particiones parciales lo basamos en la recurrencia [2] con una ligera modificación:

Hemos agregado el caso base k=1 ⇒ S2(n, k) = 1 siempre que n>0. Podría prescindirse, pero mejora el coste. El código es similar al que vimos en el segundo algoritmo mejorado parar generar particiones, que estaba basado en las combinaciones complementarias.

// S2(n, k) = sum_(j=k-1)^(n-1) S2(j, k-1) C(n-1, j)
function kPartition2(list=[], k=0) {
    let result = [];
    let n = list.length;
    if (n===k) {
        iter += n;
        result.push([...list.map(v => [v])]);
    } else if (n===0 || k===0 || n<k) {
        iter += 1;
        result = [];
    } else if (k===1) {
        iter += n;
        result.push([[...list]]);
    } else {
        iter += 1;
        let pivot = list[0];
        iter += (n-1);
        let subList = list.slice(1);
        for (let j=k-1; j<=n-1; j++) {
            iter += 1;
            let c = combineComplement2(subList, j);
            for (let i=0; i<c.length; i++) {
                iter += 1;
                let [combination, complement] = c[i];
                /* complement.length = n-1-j */
                iter += complement.length;
                let pc = [pivot, ...complement];
                let p = kPartition2(combination, k-1);
                /* p.length = S2(j, k-1) */
                for (let m=0; m<p.length; m++) {
                    iter += 1;
                    p[m].push(pc); //unshift
                    result.push(p[m]);
                }
            }
        }
    }
    return result;
}

En la herramienta Combine Tools podemos comparar los resultados generados por este algoritmo kPartition2(abcd, 3) con el anterior kPartition1(abcd, 3), obteniéndose esta tabla comparativa:

Recuerde que se presenta en la vista compacta, de tal forma que b.c.ad se obtiene como un Array [["b"], ["c"], ["a", "d"]]. Al ser subconjuntos el orden no importa, por lo que se compara correctamente b.c.ad = ad.b.c. Podríamos hacer que kPartition2() tuviese el mismo orden en los subconjuntos que kPartition1() con solo cambiar la sentencia p[k].push(pc) por p[k].unshift(pc). Pero no lo hacemos pues el método push() que inserta al final tiene coste unitario, mientras que unshit() que inserta al principio tiene coste igual a la longitud de p[k].

Analizamos el coste de forma similar a como lo hicimos en segundo algoritmo para generar particiones:

Simplificando la parte recurrente:

Aplicando sumatorios a todos los términos menos los que contienen T_cc y T(j, k-1) tenemos.

• ∑_j=k-1..n-1 1 = n-k+1
• ∑_j=k-1..n-1 C(n-1, j) (n-j)
• ∑_j=k-1..n-1 C(n-1, j) S2(j, k-1) = S2(n, k)

La tercera la vimos en [2], llegando a la solución final para el algoritmo kPartition2(list, k) que genera particiones parciales:

Donde T_cc(n, k) son las combinaciones complementarias del algoritmo combineComplement2(n, k) para lo que no he encontrado una fórmula cerrada:

En caso de coste de copia sólo hay que agregar lo resaltado en amarillo.

Un algoritmo para generar particiones usando particiones parciales

Usando las particiones parciales vistas antes y sabiendo que el total de particiones es el número de Bell, entonces B(n) = ∑_j=0..n S2(n, j), implementamos un quinto algoritmo no recursivo para generar particiones:

// B(n) = sum_(j=0)^(n) S2(n, j)
function partition5(list=[]) {
    let result = [];
    let n = list.length;
    //iter += 1;
    for (let k=1; k<=n; k++){
        //iter += 1;
        let res = kPartition1(list, k);
        /* res.length = S2(n, k) */
        for (let j=0; j<res.length; j++) {
            //iter += 1;
            result.push(res[j]);
        }
    }
    return result;
}

El coste de este algoritmo iterativo se deduce fácilmente:

Como n = ∑_k=1..n 1 y además B(n) = ∑_k=1..n S2(n, k)

siendo T_k el coste de las particiones parciales kPartition1(list, k) que vimos en el apartado anterior. Con copia es lo mismo, usando el coste de copia de las particiones parciales:

Otro algoritmo para generar particiones usando particiones parciales

Usando las particiones parciales kPartition2(list, k) vistas antes y sabiendo que el total de particiones es el número de Bell, entonces B(n) = ∑_j=0..n S2(n, j), implementamos un sexto algoritmo no recursivo para generar particiones:

// B(n) = sum_(j=0)^(n) S2(n, j)
function partition6(list=[]) {
    let result = [];
    let n = list.length;
    //iter += 1;
    for (let k=1; k<=n; k++){
        //iter += 1;
        let res = kPartition2(list, k);
        /* res.length = S2(n, k) */
        for (let j=0; j<res.length; j++) {
            //iter += 1;
            result.push(res[j]);
        }
    }
    return result;
}

El coste es igual que el anterior:

siendo T_k el coste de las particiones parciales kPartition2(list, k) que vimos en el apartado anterior. Con copia es lo mismo, usando el coste de copia de las particiones parciales:

Comparando algoritmos para generar particiones

Comparamos coste de los algoritmos vistos para generar particiones: partition1(list), partition2(list), partition3(list) en los temas anteriores y partition4(list), partition5(list) y partition6(list) en este tema:

Con los dos últimos se consiguen los peores resultados, pues hace uso de las particiones parciales. El mejor resultado es para partition4(list).