For those of you interested, here is the math.comb
CPython implementation in mathmodule.c. The basic recursive formula (no memoization) is
j = k // 2
C(n, k) = C(n, j) * C(n-j, k-j) // C(k, j)
This is a repeated use of the identity C(n, k) = (n/k) * C(n-1, k-1). The rationale given is that Karatsuba multiplication, a divide-and-conquer algorithm, runs faster multiplying numbers about the same size (therefore the choice of j).
/* Calculate P(n, k) or C(n, k) using recursive formulas.
* It is more efficient than sequential multiplication thanks to
* Karatsuba multiplication.
*/
static PyObject *
perm_comb(PyObject *n, unsigned long long k, int iscomb)
{
if (k == 0) {
return PyLong_FromLong(1);
}
if (k == 1) {
return Py_NewRef(n);
}
/* P(n, k) = P(n, j) * P(n-j, k-j) */
/* C(n, k) = C(n, j) * C(n-j, k-j) // C(k, j) */
unsigned long long j = k / 2;
PyObject *a, *b;
a = perm_comb(n, j, iscomb);
if (a == NULL) {
return NULL;
}
PyObject *t = PyLong_FromUnsignedLongLong(j);
if (t == NULL) {
goto error;
}
n = PyNumber_Subtract(n, t);
Py_DECREF(t);
if (n == NULL) {
goto error;
}
b = perm_comb(n, k - j, iscomb);
Py_DECREF(n);
if (b == NULL) {
goto error;
}
Py_SETREF(a, PyNumber_Multiply(a, b));
Py_DECREF(b);
if (iscomb && a != NULL) {
b = perm_comb_small(k, j, 1);
if (b == NULL) {
goto error;
}
Py_SETREF(a, PyNumber_FloorDivide(a, b));
Py_DECREF(b);
}
return a;
error:
Py_DECREF(a);
return NULL;
}
There is also a specialization for small n, using some pre-computed values. Pre-computed factorials in general are useful if computing many binomial coefficients, including mod p.
/* Number of permutations and combinations.
* P(n, k) = n! / (n-k)!
* C(n, k) = P(n, k) / k!
*/
/* Calculate C(n, k) for n in the 63-bit range. */
static PyObject *
perm_comb_small(unsigned long long n, unsigned long long k, int iscomb)
{
assert(k != 0);
/* For small enough n and k the result fits in the 64-bit range and can
* be calculated without allocating intermediate PyLong objects. */
if (iscomb) {
/* Maps k to the maximal n so that 2*k-1 <= n <= 127 and C(n, k)
* fits into a uint64_t. Exclude k = 1, because the second fast
* path is faster for this case.*/
static const unsigned char fast_comb_limits1[] = {
0, 0, 127, 127, 127, 127, 127, 127, // 0-7
127, 127, 127, 127, 127, 127, 127, 127, // 8-15
116, 105, 97, 91, 86, 82, 78, 76, // 16-23
74, 72, 71, 70, 69, 68, 68, 67, // 24-31
67, 67, 67, // 32-34
};
if (k < Py_ARRAY_LENGTH(fast_comb_limits1) && n <= fast_comb_limits1[k]) {
/*
comb(n, k) fits into a uint64_t. We compute it as
comb_odd_part << shift
where 2**shift is the largest power of two dividing comb(n, k)
and comb_odd_part is comb(n, k) >> shift. comb_odd_part can be
calculated efficiently via arithmetic modulo 2**64, using three
lookups and two uint64_t multiplications.
*/
uint64_t comb_odd_part = reduced_factorial_odd_part[n]
* inverted_factorial_odd_part[k]
* inverted_factorial_odd_part[n - k];
int shift = factorial_trailing_zeros[n]
- factorial_trailing_zeros[k]
- factorial_trailing_zeros[n - k];
return PyLong_FromUnsignedLongLong(comb_odd_part << shift);
}
/* Maps k to the maximal n so that 2*k-1 <= n <= 127 and C(n, k)*k
* fits into a long long (which is at least 64 bit). Only contains
* items larger than in fast_comb_limits1. */
static const unsigned long long fast_comb_limits2[] = {
0, ULLONG_MAX, 4294967296ULL, 3329022, 102570, 13467, 3612, 1449, // 0-7
746, 453, 308, 227, 178, 147, // 8-13
};
if (k < Py_ARRAY_LENGTH(fast_comb_limits2) && n <= fast_comb_limits2[k]) {
/* C(n, k) = C(n, k-1) * (n-k+1) / k */
unsigned long long result = n;
for (unsigned long long i = 1; i < k;) {
result *= --n;
result /= ++i;
}
return PyLong_FromUnsignedLongLong(result);
}
}
else {
/* Maps k to the maximal n so that k <= n and P(n, k)
* fits into a long long (which is at least 64 bit). */
static const unsigned long long fast_perm_limits[] = {
0, ULLONG_MAX, 4294967296ULL, 2642246, 65537, 7133, 1627, 568, // 0-7
259, 142, 88, 61, 45, 36, 30, 26, // 8-15
24, 22, 21, 20, 20, // 16-20
};
if (k < Py_ARRAY_LENGTH(fast_perm_limits) && n <= fast_perm_limits[k]) {
if (n <= 127) {
/* P(n, k) fits into a uint64_t. */
uint64_t perm_odd_part = reduced_factorial_odd_part[n]
* inverted_factorial_odd_part[n - k];
int shift = factorial_trailing_zeros[n]
- factorial_trailing_zeros[n - k];
return PyLong_FromUnsignedLongLong(perm_odd_part << shift);
}
/* P(n, k) = P(n, k-1) * (n-k+1) */
unsigned long long result = n;
for (unsigned long long i = 1; i < k;) {
result *= --n;
++i;
}
return PyLong_FromUnsignedLongLong(result);
}
}
/* For larger n use recursive formulas:
*
* P(n, k) = P(n, j) * P(n-j, k-j)
* C(n, k) = C(n, j) * C(n-j, k-j) // C(k, j)
*/
unsigned long long j = k / 2;
PyObject *a, *b;
a = perm_comb_small(n, j, iscomb);
if (a == NULL) {
return NULL;
}
b = perm_comb_small(n - j, k - j, iscomb);
if (b == NULL) {
goto error;
}
Py_SETREF(a, PyNumber_Multiply(a, b));
Py_DECREF(b);
if (iscomb && a != NULL) {
b = perm_comb_small(k, j, 1);
if (b == NULL) {
goto error;
}
Py_SETREF(a, PyNumber_FloorDivide(a, b));
Py_DECREF(b);
}
return a;
error:
Py_DECREF(a);
return NULL;
}
while
? you can merely useif
!!