-1

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))
lordfoog
  • 9
  • 3
  • Are you using Python 2? – mkrieger1 Jan 19 '23 at 17:58
  • 2
    `print "Bell Number: " + str(bell)` indicates that this is Python 2 code. Python 2 has been completely dead, out of date and unsupported for over 3 years (it was sunsetted January 1, 2020). Division in Python 2 works differently than Python 3, and so you're likely running into errors with using float vs integer division. Python 2 does integer division by default with the `/` operator if the operands are both integers, but Python 3 requires you to use `//`. – Random Davis Jan 19 '23 at 18:02
  • @lordfoog I think you either need to indicate why you must be using a deprecated and unsupported programming language, or convert this to Python 3. There's no reason to be on Python 2, it's obviously just going to cause confusion since everyone's on Python 3 now. – Random Davis Jan 19 '23 at 18:04
  • 1
    @RandomDavis I fixed the code so it's for Python 3 now, it didn't take much. Anyway, the formula I'm using doesn't let there be decimals anyway. I ran a couple of the plugged in formulas into a calculator and it never comes out as a decimal. The point of partitions is that it's an exact number, you can't have .5 of a partition. So I don't believe it's an issue with float vs. integer division, as the result is the same in either. – lordfoog Jan 19 '23 at 19:06
  • @lordfoog Even though the issue was not that this was Python 2, the fact that it was Python 2 made this take a lot longer to solve, so I still recommend making sure everything is on Python 3 when you submit it here in the future, for expedience purposes. – Random Davis Jan 19 '23 at 19:50

1 Answers1

2

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
Tim Roberts
  • 48,973
  • 4
  • 21
  • 30