HiveBrain v1.2.0
Get Started
← Back to all entries
patterncMinor

Bell numbers calculation in C

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
calculationnumbersbell

Problem

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 
#include 
#include  // For uint64_t
#include  // 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;
}

Solution

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.

Code Snippets

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

Context

StackExchange Code Review Q#112497, answer score: 5

Revisions (0)

No revisions yet.