Bell numbers calculation in C

Question

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;
}

Answer 1

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.

• Zorgatone: +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?

Answer 2

• 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.

• Zorgatone: I thought about memoizing even with factorials and/or binomials.
• Zorgatone: 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 😂
• mdfst13: 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: @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.
• Zorgatone: @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.