Mark Dickinson <dicki...@gmail.com> added the comment:
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) ---------- _______________________________________ Python tracker <rep...@bugs.python.org> <https://bugs.python.org/issue37295> _______________________________________ _______________________________________________ Python-bugs-list mailing list Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com