I am trying to calculate the Bell number for large sets of elements using integer partitions and Faà di Bruno's formula. The Bell number is the amount of possible partitions that can be made in a set of n elements. I want to be able to do this by calculating the amount of combinations possible in each specific integer partition possible for a set of n elements, then adding them all up. The sum should then be the Bell number.
I used Nicolas Blanc's code on calculating all of the integer partitions in a set. Then I used Faà di Bruno's formula as described in this video, to go through each integer partition and calculate the combinations possible.
The code outputs the formula with the numbers plugged in. The left factorial numbers are the subsets in the integer partitions, and the right factorial numbers are the amounts of each of those subsets. n! is then divided by those numbers. The amount of combinations is then appended to the end of the list.
For smaller sets the code works perfectly fine. This is the output for a set of 4 elements.
4! / (1!1!1!1! * 4!)
[1, 1, 1, 1, [1])
4! / (2!1!1! * 1!2!)
[2, 1, 1, [6]]
4! / (3!1! * 1!1!)
[3, 1, [4]]
4! / (2!2! * 2!)
[2, 2, [3]]
4! / (4! * 1!)
[4, [1]]
Bell Number: 15
It's up until you input numbers 12 and up that the Bell number is incorrect. The entire output is long, but inputting 12 yields 4213663 instead of 4213597, and inputting 13 yields 27645945 instead of 27644437.
After looking into it all, I cannot figure out why this happens only after 11. It could be an issue with the formula I'm using, or I made a mistake in my code or something else.
Other important links: Bell polynomials, Partition of a set, Bell numbers
import math
in_par = []
stack = []
bell = 0
def partitions(remainder, start_number = 1):
if remainder == 0:
in_par.append(list(stack))
#print(stack)
else:
for nb_to_add in range(start_number, remainder+1):
stack.append(nb_to_add)
partitions(remainder - nb_to_add, nb_to_add)
stack.pop()
x = partitions(13) # <------- input element count here
for part in in_par:
combo = 1
n = 13 # <------- input element count here
part.reverse()
refined = []
counts = []
for elem in part:
if str(refined).find(str(elem)) == -1:
refined.append(elem)
#print(str(combo) + " * " + str(math.factorial(elem)) + " a")
combo = combo * math.factorial(elem)
for elem in refined:
#print(str(combo) + " * !" + str(part.count(elem)) + " b")
combo = combo * math.factorial(part.count(elem))
counts.append(str(part.count(elem)))
for i in range(len(part)):
part[i] = str(part[i])
print(str(n) + "! / (" + ("!").join(part) + "! * " + ("!").join(counts) + "!)")
for i in range(len(part)):
part[i] = int(part[i])
combo = math.factorial(n) // combo
bell = bell + combo
part.append([combo])
print(part)
print("")
#print(str(bell))
print("Bell Number: " + str(bell))
I've reorganized things a little bit, but the key point is your string search for integers:
import math
bell = 0
N = 13 # <------- input element count here
stack = []
def partitions(remainder, start_number = 1):
if remainder == 0:
yield list(reversed(stack))
else:
for nb_to_add in range(start_number, remainder+1):
stack.append(nb_to_add)
yield from partitions(remainder - nb_to_add, nb_to_add)
stack.pop()
for part in partitions(N):
combo = 1
n = N
refined = []
counts = []
for elem in part:
if elem not in refined:
refined.append(elem)
combo = combo * math.factorial(elem)
for elem in refined:
combo = combo * math.factorial(part.count(elem))
counts.append(str(part.count(elem)))
for i in range(len(part)):
part[i] = str(part[i])
print(str(n) + "! / (" + ("!").join(part) + "! * " + ("!").join(counts) + "!)")
for i in range(len(part)):
part[i] = int(part[i])
combo = math.factorial(n) // combo
bell = bell + combo
part.append([combo])
print(part)
print("")
print("Bell Number:", bell)
Output:
...
Bell Number: 27644437