Author mark.dickinson
Recipients PedanticHacker, Stefan Pochmann, mark.dickinson, mcognetta, pablogsal, rhettinger, serhiy.storchaka, tim.peters
Date 2021-12-21.08:42:17
SpamBayes Score -1.0
Marked as misclassified Yes
Message-id <>
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)
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: <>
2021-12-21 08:42:17mark.dickinsonlinkissue37295 messages
2021-12-21 08:42:17mark.dickinsoncreate