6
\$\begingroup\$

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;
}
\$\endgroup\$
5
\$\begingroup\$

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.

\$\endgroup\$
1
  • \$\begingroup\$ +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? \$\endgroup\$ – Zorgatone Dec 2 '15 at 7:23
3
\$\begingroup\$
  • 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.

\$\endgroup\$
6
  • \$\begingroup\$ I thought about memoizing even with factorials and/or binomials. \$\endgroup\$ – Zorgatone Dec 1 '15 at 22:05
  • \$\begingroup\$ 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 😂 \$\endgroup\$ – Zorgatone Dec 1 '15 at 22:06
  • \$\begingroup\$ 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. \$\endgroup\$ – mdfst13 Dec 2 '15 at 1:15
  • \$\begingroup\$ @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. \$\endgroup\$ – mdfst13 Dec 2 '15 at 1:21
  • \$\begingroup\$ @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. \$\endgroup\$ – Zorgatone Dec 2 '15 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.