That is a terrible way to compute a binomial. 18!, as well as 13!, are larger than 231-1, and so overflow the usual int
. In fact, you only need to go to 21! to overflow a uint64_t
. (No, this is not a "compiler error" as suggested in a tag. You are not understanding the limitations of the data types.)
You can compute much, much larger binomials even with just the int
type by never multiplying out n!. There is a lot of cancellation going on. This shows why it is a good idea to do a calculation by hand first, so that you can see what's going on, instead of just jumping in and writing code.
For example 18C2 should not be a problem at all. It is 18!/(16!2!), where 18!/16! is just 18 * 17. You should pick k to be small, e.g. if you had 18C16, you would change that to the equivalent 18C2. Then n!/(n-k)! will throw out at least half of the terms of n!. Then with what's left, and the use of some binomial identities, you can cancel terms in the numerator and denominator (k!) as you go along, avoiding overflows and fractional results.
Even better, at least for smaller n, you could do it with just additions by filling in Pascal's triangle. Along with saving the results in memory, this is super fast.
Here is code that combines Pascal's triangle with direct calculations. This uses fast table lookups up to MAXN
, and calculations after that, with the ability to return any binomial that is less than 264 (except for (264-1)C1, which will appear to be an overflow due to the special return value):
#include <stdint.h>
#include <limits.h>
// Binary GCD algorithm. This assumes that a and b are positive.
uint64_t gcd(uint64_t a, uint64_t b) {
int z = __builtin_ctzll(a | b);
a >>= __builtin_ctzll(a);
do {
b >>= __builtin_ctzll(b);
b -= a;
uint64_t m = (int64_t)b >> 63;
a += b & m;
b = (b + m) ^ m;
} while (b);
return a << z;
}
// Return n C k or (uint64_t)-1 on overflow. This assumes that 1 < k <= n/2.
static uint64_t binom_calc(uint64_t n, uint64_t k) {
uint64_t c = n & 1 ? n * ((n - 1) >> 1) : (n >> 1) * (n - 1);
for (uint64_t i = 3; i <= k; i++) {
uint64_t m = n - i + 1;
uint64_t g = gcd(m, i);
m /= g;
c /= i / g;
if ((uint64_t)-1 / c < m)
return -1;
c *= m;
}
return c;
}
// MAXN is the maximum n for which the table is used. Special MAXN values:
// 967 -- binom_calc() only used for k <= 7
// 1913 -- binom_calc() only used for k <= 6
// 4868 -- binom_calc() only used for k <= 5
// 18580 -- binom_calc() only used for k <= 4
#define MAXN 1913
static uint64_t binom_table[MAXN-3][32] = {0};
// Return n C k or (uint64_t)-1 on overflow. This assumes that 0 < 2*k <= n+1.
static uint64_t binom_recur(uint64_t n, uint64_t k) {
if (n < (k << 1))
k--;
if (k == 1)
return n;
if (binom_table[n - 4][k - 2])
return binom_table[n - 4][k - 2];
uint64_t a = binom_recur(n - 1, k - 1);
uint64_t b = binom_recur(n - 1, k);
a += b;
return binom_table[n - 4][k - 2] = a < b ? (uint64_t)-1 : a;
}
// Return n C k, or 0 if k > n, or (uint64_t)-1 on overflow. This is not
// thread-safe unless and until binom_fill() is called. Otherwise, the table is
// filled in on an as-needed basis.
uint64_t binom(uint64_t n, uint64_t k) {
if (k > n)
return 0;
if (k > (n >> 1))
k = n - k;
if (k == 0)
return 1;
if (k == 1)
return n;
return n > MAXN ? binom_calc(n, k) : binom_recur(n, k);
}
// Fill in the binomial table. If you do this once from the main thread, then
// binom() can safely be used from multiple threads at the same time. This only
// takes about 140 microseconds on an Apple M1 for MAXN = 1913.
void binom_fill(void) {
for (uint64_t n = MAXN, k = 2; n >= 4; n--)
for (; k <= (n << 1); k++)
if (binom_recur(n, k) + 1 == 0)
break;
}
To go further, you can bust through the 64-bit barrier and use the bignum library.