# Bell numbers calculation in C

This is my implementation of the bell() function to calculate the n-th bell number in C.

Along with it I made the factorial() and binomial() functions.

I have serious doubts about their efficiency, in particular whether my implementation of the factorial() function should be correctly optimized for tail-recursion, and whether it is possible to rewrite the bell() to be more efficient.

Let me know any critics about my code, and if I can rely on this implementation or maybe use different approaches.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h> // For uint64_t
#include <inttypes.h> // For PRIu64

#define factorial(n) fact(n, 1)

static inline uint64_t fact(uint8_t n, uint64_t inc);
static inline uint64_t binomial(uint8_t n, uint8_t k);
static inline uint64_t bell(uint8_t n);

int main(int argc, char **argv) {
printf("Bell: %" PRIu64 ".\n", bell(21));

return EXIT_SUCCESS;
}

// Tail recursion?
static inline uint64_t fact(uint8_t n, uint64_t inc) {
return (n == 0) ? inc : fact(n - 1, inc * (uint64_t)n);
}

static inline uint64_t binomial(uint8_t n, uint8_t k) {
return factorial(n) / (factorial(k) * factorial(n - k));
}

// can calculate up to bell(21)
static inline uint64_t bell(uint8_t n) {
if (n == 0) {
return 1;
}

uint8_t i;
uint64_t sum = 0;

n--;

for (i = 0; i <= n; i++) {
sum += binomial(n, i) * bell(i);
}

return sum;
}


Even without memoization, I cut the calculation time by about 40% when I switched from your version of binomial to the following:

static uint64_t permute(uint8_t n, uint8_t k) {
uint64_t result = 1;

for (; n > k; --n) {
result *= n;
}

return result;
}

static inline uint64_t binomial(uint8_t n, uint8_t k) {
if ( n - k > k ) {
return permute(n, n - k) / permute(k, 1);
} else {
return permute(n, k) / permute(n - k, 1);
}
}


Note that you could use the ternary operator rather than the if in binomial.

static inline uint64_t binomial(uint8_t n, uint8_t k) {
return ( n - k > k )
? permute(n, n - k) / permute(k, 1)
: permute(n, k) / permute(n - k, 1);
}


I find the if easier to follow in a complex statement like that.

The reason this works is that $\binom{n}{k} = \frac{n!}{k!(n-k)!}$ can also be written as

$$\frac{n * (n-1) * ... * (k+1)*k!}{k!(n-k)!} = \frac{n * (n-1) * ... * (k+1)}{(n-k)!}$$

and replacing factorial(k) with permute(k, 1) works because

$$\frac{n!}{1!} = n!$$

It's unclear if memoization of the factorials would help more. It may depend on how many values you need to calculate.

• +1 for the alternative formula. I can see before trying it out it's more efficient. @mdfst13 for the bell numbers I know I get an overflow when calculating bell(22). I even thought I could just return the first 22 bell numbers using a simple switch/case instead. If you think about it, there are not much bell numbers that fit even into a uint64_t. That is if I can just write uint64_t number literals, right? – Zorgatone Dec 2 '15 at 7:23
• Computing binomials

Using a formula $\binom{n}{i} = \frac{n!}{i!(n-i)!}$ results in calculating very large numbers, and computing all the binomials your way you need to recompute factorials over and over again. An identity $\binom{n}{i} =\frac{n-i + 1}{i} \binom{n}{i-1}$ let you compute them on the fly:

binomial = 1;
for (i = 0; i < n; i++) {
sum += binomial * bell(i);
binomial = binomial * (n - i + 1) / i;
}

• Computing Bell numbers

Your code recompute them recursively over and over again. You will gain quite a performance by memoizing them.

• I thought about memoizing even with factorials and/or binomials. – Zorgatone Dec 1 '15 at 22:05
• Anyway there can be 21 numbers that can be saved in a 64-bit integer. It could even be a switch/case that returns one of those 21 😂 – Zorgatone Dec 1 '15 at 22:06
• I tried comparing this method to my version, but I kept getting a runtime error, even after I changed it to start with i = 1 to avoid the divide by zero error. – mdfst13 Dec 2 '15 at 1:15
• @Zorgatone bell(22) is 0x1ADCB3D5, which is only 61 bits. You just can't calculate it in 64-bit integers with your method. You might manage with a different, more complicated method. Or of course you could use a larger integer size to hold the intermediate results. – mdfst13 Dec 2 '15 at 1:21
• @vnp yeah @mdfst13 is right, without even trying it out, I can see there is a division by i in the first loop iteration with i = 0. – Zorgatone Dec 2 '15 at 7:19