Particiones de un conjunto y número de Bell
Las particiones de un conjunto
Una partición de un conjunto son los agrupamientos de sus elementos en subconjuntos no vacíos, de tal forma que cada elemento está incluido en un único subconjunto.
En la Figura vemos las cinco particiones que se generan con el conjunto A = {a, b, c}. Por ejemplo, observando la segunda partición {{a, b}, {c}} se evidencia que la unión de estos subconjuntos {a, b} ⋃ {c} = A es el conjunto de partida. Y que la intersección es vacía {a, b} ⋂ {c} = ∅, que es lo mismo que decir que los subconjuntos son disjuntos.
En el algoritmo que implementaremos después generamos las particiones del conjunto {a, b, c} con el siguiente resultado:
{{{a, b, c}},
{{b}, {a, c}}, {{c}, {a, b}}, {{b, c}, {a}},
{{c}, {b}, {a}}}
De forma compacta suele escribirse como abc, b.ac, c.ab, bc.a, c.b.a
cuando los elementos del conjunto de partida no tienen más de un caracter, separándose los subconjuntos con un punto y las particiones con una coma.
Observe que el de la Figura de forma compacta es abc, ab.c, ac.b, a.bc, a.b.c
, que es el mismo que el anterior, puesto que en los conjuntos el orden no importa. De hecho esta disposición es la que genera un segundo algoritmo más eficiente que veremos en el tema siguiente.
Relacionado con esto, hace algunos años estuve viendo la estructura de partición o de conjuntos disjuntos, estructura que permite resolver mediante algoritmos voraces los problemas de la planificación de tareas y del árbol de recubrimiento mínimo en un grafo conexo no dirigido, con el algoritmo de Kruskal. La estructura tiene un coste de búsqueda esencialmente lineal, lo que mejora esas aplicaciones.
Números de Bell
El número de Bell Bn representa el total de particiones de un conjunto con n elementos. En lugar de denotarlo como Bn aquí lo escribiremos como B(n). Se calcula con la recurrencia siguiente:
B(n) = ∑j=0..n-1 C(n-1, j) B(j)
Otra forma alternativa es
B(n) = ∑j=1..n C(n-1, j-1) B(n-j)
En la aplicación Combine Tools lo implementamos en la pestaña "Numéricos" con el siguiente algoritmo para la primera expresión:
// B(n) = sum_(j=0)^(n-1) C(n-1, j) B(j) function bell(n=0) { if (n===0) { return 1; } else { let sum = 0; for (let j=0; j<=n-1; j++) { sum += binomial(n-1, j) * bell(j); } return sum; } }
Y con este para la alternativa:
// B(n) = sum_(j=1)^(n) C(n-1, j-1) B(n-j) function bell1(n) { if (n===0) { return 1; } else { let sum = 0; for (let j=1; j<=n; j++) { sum += binomial(n-1, j-1) * bell1(n-j); } return sum; } }
Los primeros valores que se generan con esos algoritmos son:
n | B(n) |
---|---|
0 | 1 |
1 | 1 |
2 | 2 |
3 | 5 |
4 | 15 |
5 | 52 |
6 | 203 |
7 | 877 |
8 | 4140 |
Triángulo de Bell
Otra forma de obtener un número de Bell es con el Triángulo de Bell, también denominado Array de Aitken o Triángulo de Pierce, como puede verse en la secuencia OEIS A011971. El algoritmo kbell(n, k)
no es un recursivo, siendo un par de bucles anidados:
// Bell triangle function kbell(n=0, k=0) { let val = [1]; for (let i=0; i<=n-1; i++) { let prev = [...val]; val = [val.pop()]; for (let j=0; j<=i; j++) { val.push(val[j]+prev[j]); } } return k<val.length ? val[k] : 0; }
En la aplicación se vierte en una tabla como esta para valores n=0 hasta n=8:
n↓k→ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
---|---|---|---|---|---|---|---|---|---|
0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2 | 2 | 3 | 5 | 0 | 0 | 0 | 0 | 0 | 0 |
3 | 5 | 7 | 10 | 15 | 0 | 0 | 0 | 0 | 0 |
4 | 15 | 20 | 27 | 37 | 52 | 0 | 0 | 0 | 0 |
5 | 52 | 67 | 87 | 114 | 151 | 203 | 0 | 0 | 0 |
6 | 203 | 255 | 322 | 409 | 523 | 674 | 877 | 0 | 0 |
7 | 877 | 1080 | 1335 | 1657 | 2066 | 2589 | 3263 | 4140 | 0 |
8 | 4140 | 5017 | 6097 | 7432 | 9089 | 11155 | 13744 | 17007 | 21147 |
Según esa página de OEIS, el Triángulo de Bell también puede definirse recursivamente así:
B(n, 0) = B(n-1, n-1),
B(n, k) = B(n, k-1) + B(n-1, k-1)
Se observa en los valores obtenidos que el último elemento de una fila es igual que el primero de la siguiente, como vemos en la recurrencia anterior B(n, 0) = B(n-1, n-1). Y que un elemento es la suma del anterior en la misma fila y del que están en la fila anterior en la misma columna, con lo que se cumple la fórmula recursiva que vimos antes B(n, k) = B(n, k-1) + B(n-1, k-1). Los valores en la primera posición de cada fila para k=0 con B(n, 0) se corresponden con los números de Bell 1, 1, 2, 5, 15, 52, 203, 877, ... También se repiten en la diagonal B(n-1, n-1). Esto se se observa en las celdas resaltadas de la tabla anterior.
Con esa recurrencia implementamos el algoritmo kbell1(n, k)
para obtener el número del Triangulo de Bell en la fila n y en la columna k, entendiendo que el propio numero de Bell es el que está en la columna cero, es decir, el que se obtiene con kbell1(n, 0)
. O también con la diagonal kbell1(n-1, n-1)
para n≥1:
// Bell triangle recursive function kbell1(n=0, k=0) { if (n<k || k<0) { return 0; } else if (n===0 && k===0) { return 1; } else if (k===0) { return kbell1(n-1, n-1); } else { return kbell1(n, k-1) + kbell1(n-1, k-1); } }
El número de Bell y los números de Stirling de segunda clase
Hay otra identidad que veremos en el tema siguiente y que tiene que ver con los números de Stirling de segunda clase que denotamos como S2(n, k), que adelantamos que son las particiones de un conjunto con n elementos en exactamente k subconjuntos no vacíos:
Lo implementamos con este código bell2(n)
:
// B(n) = sum_(j=0)^(n) S2(n, j) // S2(n, k) = S2(n-1, k-1) + k × S2(n-1, k) function bell2(n) { let sum = 0; for (let j=0; j<=n; j++) { sum += stirlingS2(n, j); } return sum; }
Combinando esta con la que vimos en [1] incrementando un índice
tenemos esta fórmula debida a Spivey:
Lo implementamos con este código, observando que no es un recursivo, pero el primer algoritmo para calcular bell(i)
si lo es.
//B(n+k) = sum_(i=0)^(n) sum_(j=0)^(k) j^(n-) S2(k, j) C(n, i) B(i) function bell3(n=0, k=0) { let sum = 0; for (let i=0; i<=n; i++) { for (let j=0; j<=k; j++) { sum += j**(n-i) * stirlingS2(k, j) * binomial(n, i) * bell(i); } } return sum; }
Los valores obtenidos con ese codigo, que se muestran en la siguiente tabla, son todos números de Bell, tanto horizontal como verticalmente y las diagonales derecha, mientras que las diagonales izquierda incluyen el mismo valor.
n↓k→ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 2 | 5 | 15 | 52 | 203 | 877 | 4140 |
1 | 1 | 2 | 5 | 15 | 52 | 203 | 877 | 4140 | 21147 |
2 | 2 | 5 | 15 | 52 | 203 | 877 | 4140 | 21147 | 115975 |
3 | 5 | 15 | 52 | 203 | 877 | 4140 | 21147 | 115975 | 678570 |
4 | 15 | 52 | 203 | 877 | 4140 | 21147 | 115975 | 678570 | 4213597 |
5 | 52 | 203 | 877 | 4140 | 21147 | 115975 | 678570 | 4213597 | 27644437 |
6 | 203 | 877 | 4140 | 21147 | 115975 | 678570 | 4213597 | 27644437 | 190899322 |
7 | 877 | 4140 | 21147 | 115975 | 678570 | 4213597 | 27644437 | 190899322 | 1382958545 |
8 | 4140 | 21147 | 115975 | 678570 | 4213597 | 27644437 | 190899322 | 1382958545 | 10480142147 |
En el tema Particiones y números de Stirling de segunda clase veremos con más detalle esta relación, implementándose un algoritmo más eficiente para generar las particiones.
La fórmula de Dobinski para calcular el número de Bell
Se puede calcular el número de Bell también con esta fórmula de Dobinski:
B(n) = | 1 | ∑j=0..∞ | jn |
e | j! |
No es un fórmula cerrada pues es una suma infinita. Sin embargo nos puede servir para calcular el número de Bell hasta cierto número, por ejemplo, hasta n=20 solo necesitamos un rango j=0..29 en el sumatorio, obteniendo estos valores en la primera columna con el algoritmo bell4(n)
que veremos a continuación, incluyendo todos los decimales que nos da JavaScript. En la segunda columna obtenemos los valores con el algoritmo bell2(n)
que vimos antes y que se basa en los números de Stirling de segunda clase, siendo correctos hasta n=22, pues en color azul en esta columna vemos los valores que se obtienen en la página OEIS 110:
n | B(n) = bell4(n) | B(n) = bell2(n) & OEIS 110 |
---|---|---|
0 | 1.0000000000000002 | 1 |
1 | 1.0000000000000002 | 1 |
2 | 2 | 2 |
3 | 5.000000000000001 | 5 |
4 | 15.000000000000002 | 15 |
5 | 51.99999999999998 | 52 |
6 | 203 | 203 |
7 | 876.9999999999999 | 877 |
8 | 4140 | 4140 |
9 | 21147.000000000004 | 21147 |
10 | 115975.00000000003 | 115975 |
11 | 678569.9999999999 | 678570 |
12 | 4213596.999999999 | 4213597 |
13 | 27644436.99999999 | 27644437 |
14 | 190899322.00000003 | 190899322 |
15 | 1382958545.0000005 | 1382958545 |
16 | 10480142147 | 10480142147 |
17 | 82864869804 | 82864869804 |
18 | 682076806159.0002 | 682076806159 |
19 | 5832742205057 | 5832742205057 |
20 | 51724158235372 | 51724158235372 |
21 | 474869816156750.8 | 474869816156751 |
22 | 4506715738447322 | 4506715738447323 |
23 | 44152005855084340 | 44152005855084344 44152005855084346 |
24 | 445958869294804860 | 445958869294805300 445958869294805289 |
25 | 4638590332229987000 | 4638590332229999000 4638590332229999353 |
Se observa que hasta n=20 podemos redondear al antero más cercano obteniéndose el valor correcto de B(n). E incluso podríamos redondear hasta n=21, pero a partir de ahí ya no se obtienen valores correctos debido a la precisión de los números decimales en JavaScript, pues se usa Math.E = 2.718281828459045
con 15 decimales para el número e y operaciones de división que involucran decimales.
No debemos olvidar que en los enteros en JavaScript el máximo entero puede obtenerse con la constante Number.MAX_SAFE_INTEGER = 9007199254740991
, por lo que B(23) = 44152005855084346 > 9007199254740991 y JavaScript pondrá ceros al final de cada resultado a partir de n≥23.
Implementamos bell4(n)
usando ese rango j=0..29, usando Math.round()
si quieremos redondearlo a entero:
//B(n) = 1/e sum_(j=0)^(∞) j^n/j! //n≤20 ⇒ B(n) = Round(1/e sum_(j=0)^(29) j^n/j!) function bell4(n=0) { let sum = 0; for (let j=0; j<=29; j++) { sum += j**n / factorial(j); } //return sum/Math.E return Math.round(sum/Math.E); }
Función generadora del número de Bell
El número de Bell B(n) tiene por función generadora la siguiente:
∑n≥0 B(n) | xn | = eex-1 |
n! |
Ya vimos en temas anteriores que podemos usar Wolfram Alpha para ver si encuentra una serie para una función generadora. Si lo hacemos con la función anterior obtenemos esta serie:
eex-1 = ∑n≥0 | (-1 + ex)n |
n! |
Como sabemos que ex = ∑n≥0 xn/n! entonces usando otra variable k
eex-1 = ∑n≥0 | (-1 + ∑k≥0 xk/k!)n |
n! |
Lo anterior parece más díficil de tratar que otras que hemos visto en temas anteriores. Por eso en estos temas sobre particiones no usaré la técnica de la función generadora.
Primer algoritmo para generar las particiones de un conjunto
Veamos ahora el código para generar las particiones:
// n≤2 ⇒ B(n)=1, B(n) = sum_(j=0)^(n-1) C(n-1, j) B(j) function partition1(list=[]){ let n = list.length; let result = [] if (n===0) { //iter += 1; result.push([]); } else if (n===1) { //iter += 1; result.push([[list[0]]]); } else if (n===2) { //iter += 1; result.push([[list[0], list[1]]], [[list[1]], [list[0]]]); } else { //iter += 1; //iter += (n-1); let subList = list.slice(1); let pivot = list[0]; for (let j=0; j<=n-1; j++) { //iter += 1; let c = combine2(subList, j); let cc = combine2(subList, n-1-j); for (let i=0; i<c.length; i++) { //iter += 1; let complement = cc[cc.length-1-i]; //iter += complement.length; let pc = [pivot, ...complement]; let p = partition1(c[i]); for (let k=0; k<p.length; k++) { //iter += 1; p[k].push(pc); result.push(p[k]); } } } } return result; }
El algoritmo está basado en la recurrencia del número de Bell que ya vimos:
B(n) = ∑j=0..n-1 C(n-1, j) B(j)
Para reducir el coste forzamos casos bases para n≤2, componiendo un subresultado cuyo coste es unitario, pues es un coste constante independientemente del tamaño del problema.
Para n≥2 recortamos en una unidad la lista inicial, guardando el pivote (pivo
), e iteramos por ella para aplicar el sumatorio de la recurrencia de Bell anterior. Generamos las combinaciones y sus complementarias, tal como vimos en el tema anterior. Usamos el algoritmo combine2(list, k), de tal forma que para la lista subList
que tiene longitud n-1, las combinaciones serán combine2(subList, j)
y las complementarias combine2(subList, n-1,j)
.
Iteramos por las combinaciones componiendo el array pc = [pivot, ...complement]
y realizamos una nueva llamada, tras lo cual obtenemos las particiones, por las que iteramos para agregar ese array pc
y obtener, finalmente, una partición.
Con el conjunto de elementos list = {a, b, c, d} obtenemos en la aplicación ejecutando el algoritmo partition1(list)
las siguientes 15 particiones, pues B(4) = 15:
{{a, b, c, d}}, {{b}, {a, c, d}}, {{c}, {a, b, d}}, {{d}, {a, b, c}}, {{b, c}, {a, d}}, {{c}, {b}, {a, d}}, {{b, d}, {a, c}}, {{d}, {b}, {a, c}}, {{c, d}, {a, b}}, {{d}, {c}, {a, b}}, {{b, c, d}, {a}}, {{c}, {b, d}, {a}}, {{d}, {b, c}, {a}}, {{c, d}, {b}, {a}}, {{d}, {c}, {b}, {a}}
Realmente en la aplicación se manejan estructuras de array en lugar de conjuntos, devolviendo arrays también.
[[a, b, c, d]], [[b], [a, c, d]], [[c], [a, b, d]], [[d], [a, b, c]], [[b, c], [a, d]], [[c], [b], [a, d]], [[b, d], [a, c]], [[d], [b], [a, c]], [[c, d], [a, b]], [[d], [c], [a, b]], [[b, c, d], [a]], [[c], [b, d], [a]], [[d], [b, c], [a]], [[c, d], [b], [a]], [[d], [c], [b], [a]]
Así una partición como la penúltima [[c, d], [b], [a]]
equivale a {{c, d}, {b}, {a}}
, donde el orden no importa en los conjuntos, pues es lo mismo que, por ejemplo, {{c, d}, {a}, {b}}
o bien {{d, c}, {b}, {a}}
entre otras posibilidades de cambiar de sitio los subconjuntos o los elementos de un subconjunto.
La aplicación permite también una vista compacta como comentamos en un apartado anterior:
abcd, b.acd, c.abd, d.abc, bc.ad, c.b.ad, bd.ac, d.b.ac, cd.ab, d.c.ab, bcd.a, c.bd.a, d.bc.a, cd.b.a, d.c.b.a
Podemos representar esquemáticamente como evoluciona el algoritmo para este ejemplo:
partition([a,b,c,d]) pivot=a subList=[b,c,d] j=0 c=[[]] cc=[[b,c,d]] i=0 pc=[a,b,c,d] partition([]) => [[]] k=0 ==>> [a,b,c,d] abcd j=1 c=[[b],[c],[d]] cc=[[b,c],[b,d],[c,d]] i=0 pc=[a,c,d] partition([b]) => [[b]] k=0 ==>> [[b],[a,c,d]] b.acd i=1 pc=[a,b,d] partition([c]) => [[c]] k=0 ==>> [[c],[a,b,d]] c.abd i=2 pc=[a,b,c] partition([d]) => [[d]] k=0 ==>> [[d],[a,b,c]] d.abc j=2 c=[[b,c],[b,d],[c,d]] cc=[[b],[c],[d]] i=0 pc=[a,d] partition([b,c]) => [[b,c]],[[c],[b]] k=0 ==>> [[b,c],[a,d]] bc.ad k=1 ==>> [[c],[b],[a,d]] c.b.ad i=1 pc=[a,c] partition([b,d]) => [[b,d]],[[d],[b]] k=0 ==>> [[b,d],[a,c]] bd.ac k=1 ==>> [[d],[b],[a,c]] d.b.ac i=2 pc=[a,b] partition([c,d]) => [[c,d]],[[d],[c]] k=0 ==>> [[c,d],[a,b]] cd.ab k=1 ==>> [[d],[c],[a,b]] d.c.ab j=3 c=[[b,c,d]] cc=[] i=0 pc=[a] partition([b,c,d]) pivot=b subList=[c,d] j=0 c=[[]] cc=[[c,d]] i=0 pc=[b,c,d] partition([]) => [[]] k=0 ==>> [[b,c,d]] ==>> [[b,c,d],[a]] bcd.a j=1 c=[[c],[d]] cc=[[d],[c]] i=0 pc=[b,d] partition([c]) => [[c]] k=0 ==>> [[c],[b,d]] ==>> [[c],[b,d],[a]] c.bd.a i=1 pc=[b,c] partition([d]) => [[d]] k=0 ==>> [[d],[b,c]] ==>> [[d],[b,c],[a]] d.bc.a j=2 c=[[c,d]] cc=[[]] i=0 pc=[b] partition([c,d]) => [[c,d]],[[d],[c]] k=0 ==>> [[c,d],[b]] ==>> [[c,d],[b],[a]] cd.b.a k=1 ==>> [[d],[c],[b]] ==>> [[d],[c],[b],[a]] d.c.b.a
Como se observa en el algoritmo, en c
tenemos las combinaciones y en cc
sus complementarias. Un flecha simple =>
indica la salida por un caso base. Una doble ==>>
indica la recomposición para obtener una partición, es decir, el resultado de la sentencia result.push(p[k])
en el código. Disponemos a la derecha una columna con las particiones en forma compacta, para verificar que coinciden con la de la lista anterior que se obtiene en la aplicación.
Coste del primer algoritmo para generar particiones
En primer lugar vemos en el código del algoritmo que expusimos en el apartado anterior que tenemos casos base para n≤2 donde el coste es unitario, pues en cualquiera de esos casos el coste de composición del subresultado es siempre constante respecto al tamaño del problema.
Al entrar en la parte de la alternativa else
sumamos uno y también n-1 que es el coste de list.slice(1)
, acción que recorta la lista en una unidad, pasando de longitud n a n-1 que guardamos en subList
.
El bucle itera en en el rango j=0..n-1 sumando los siguientes costes:
- 1 por cada iteración de este bucle
- 2 Tc(n-1, j) pues ejecutamos el algoritmo combine2(list, k) para obtener las combinaciones C(n-1, j) y sus complementarias C(n-1, n-1-j) usando
let c = combine2(subList, j)
para las combinaciones ylet cc = combine2(subList, n-1-j)
para las complementarias. El coste es el mismo puesto que en general C(n, k) = C(n, n-k). - Viendo que
c.length = cc.length
= C(n-1, j) = C(n-1, n-1-j), iteramos en el buclefor (let i=0; i<c.length; i++)
por cada combinación:- 1 por cada iteración de este bucle
- En este bucle
c[i]
es una combinación ycc[cc.length-1-i]
es su complementaria. Componemos el pivote con las combinaciones complementarias conlet pc = [pivot, ...complement]
que tiene un coste de n-j-1 pues esa es la longitud decomplement
. - Realizamos una nueva llamada
let p = partition1(c[i])
usando como lista las combinaciones, con coste T(j). - El número de particiones
p
que se obtienen es B(j), iterando por este bucle en cuyo interior sumamos uno, pues los dospush()
tienen coste unitario que se acumula en serie con el uno.
Sumando obtenemos esta definición para la recurrencia del coste:
T(n) = 1 + (n-1) + ∑j=0..n-1 ( 1 + 2 Tc(n-1, j) + C(n-1, j)×( 1 + (n-j-1) + T(j) + B(j)×1 ) )
simplificando
T(n) = n + ∑j=0..n-1 ( 1 + 2 Tc(n-1, j) + C(n-1, j) ( n - j + T(j) + B(j) ) )
Aplicamos el sumatorio a cada término
Recordando que el coste de combine2(list, k) era T(n, k) = 2 C(n+1, k) - 1 + k, resolviéndolo en WolframAlpha
El siguiente término del sumatorio también lo resolvemos con WolframAlpha:
Por último tenemos la propia definición del número de Bell como vimos en un [1]:
Tras simplificar 2n+2 + (n+1) 2n-2 = 24 2n-2 + (n+1) 2n-2 = (n+17) 2n-2, llegamos a la solución final recurrente para el coste del algoritmo partition1(list)
:
T(n) = (n+17) 2n-2 + n2 - n - 4 + B(n) + ∑j=0..n-1 C(n-1, j) T(j)
Donde B(n) son los números de Bell, para lo cual ya vimos que no hay una fórmula cerrada, calculándose con la recurrencia ya vista:
B(n) = ∑j=0..n-1 C(n-1, j) B(j)
La solución general no es una fórmula cerrada, es decir, es una fórmula recursiva. Veáse que es como la de B(n), un sumatorio en el rango j=0..n-1 que contiene C(n-1, j) multiplicado por una llamada j-ésima. Si para la recurrencia de B(n) no se conoce una fórmula cerrada, es posible que tampoco se encuentre para la nuestra T(n).
Incorporando el coste de copiar los resultados al primer algoritmo
En este caso hemos de considerar el coste de la copia en el algoritmo combine2(list, k) para las combinaciones Tc y sus complementarias Tcc:
T(n) = 1 + (n-1) + ∑j=0..n-1 1 + Tc(n-1, j) + Tcc(n-1, j) +
C(n-1, j)×( 1 + (n-j-1) + T(j) + B(j)×1 )
Para las combinaciones es el siguiente, como puede calcularse en WolframAlpha
Y para las complementarias es el siguiente, como puede calcularse en WolframAlpha
Vemos que ambas son iguales, así que sumando estas dos más el resto de términos que son iguales que en el apartado anterior, llegamos a la siguiente solución final para el algoritmo partition1(list)
para el caso de coste con copia de subresultados:
T(n) = (3n+15) 2n-2 + n2 - n - 4 + B(n) + ∑j=0..n-1 C(n-1, j) T(j)