 Author mark.dickinson PedanticHacker, Stefan Pochmann, mark.dickinson, mcognetta, pablogsal, rhettinger, serhiy.storchaka, tim.peters 2021-12-21.08:42:17 -1.0 Yes <1640076137.31.0.0213484115566.issue37295@roundup.psfhosted.org>
Content
One approach that avoids the use of floating-point arithmetic is to precompute the odd part of the factorial of n modulo 2**64, for all small n. If we also precompute the inverses, then three lookups and two 64x64->64 unsigned integer multiplications gets us the odd part of the combinations modulo 2**64, hence for small enough n and k gets us the actual odd part of the combinations.

Then a shift by a suitable amount gives comb(n, k).

Here's what that looks like in Python. The "% 2**64" operation obviously wouldn't be needed in C: we'd just do the computation with uint64_t and rely on the normal wrapping semantics. We could also precompute the bit_count values if that's faster.

import math

# Max n to compute comb(n, k) for.
Nmax = 67

# Precomputation

def factorial_odd_part(n):
f = math.factorial(n)
return f // (f & -f)

F = [factorial_odd_part(n) % 2**64 for n in range(Nmax+1)]
Finv = [pow(f, -1, 2**64) for f in F]
PC = [n.bit_count() for n in range(Nmax+1)]

# Fast comb for small values.

def comb_small(n, k):
if not 0 <= k <= n <= Nmax:
raise ValueError("k or n out of range")
return (F[n] * Finv[k] * Finv[n-k] % 2**64) << k.bit_count() + (n-k).bit_count() - n.bit_count()

# Exhaustive test

for n in range(Nmax+1):
for k in range(0, n+1):
assert comb_small(n, k) == math.comb(n, k)
History
Date User Action Args
2021-12-21 08:42:17mark.dickinsonsetrecipients: + mark.dickinson, tim.peters, rhettinger, serhiy.storchaka, PedanticHacker, mcognetta, Stefan Pochmann, pablogsal
2021-12-21 08:42:17mark.dickinsonsetmessageid: <1640076137.31.0.0213484115566.issue37295@roundup.psfhosted.org>