Just for fun, here's a gonzo implementation (without arg-checking) using ideas from the sketch. All factors of 2 are shifted out first, and all divisions are done before any multiplies. For large arguments, this can run much faster than a dumb loop. For example, combp(10**100, 400) takes about a quarter the time of a dumb-loop divide-each-time-thru implementation. # Return number of trailing zeroes in `n`. def tzc(n): result = 0 if n: mask = 1 while n & mask == 0: result += 1 mask <<= 1 return result # Return exponent of prime `p` in prime factorization of # factorial(k). def facexp(k, p): result = 0 k //= p while k: result += k k //= p return result def combp(n, k): from heapq import heappop, heapify, heapreplace if n-k < k: k = n-k if k == 0: return 1 if k == 1: return n firstnum = n - k + 1 nums = list(range(firstnum, n+1)) assert len(nums) == k # Shift all factors of 2 out of numerators. shift2 = 0 for i in range(firstnum & 1, k, 2): val = nums[i] c = tzc(val) assert c nums[i] = val >> c shift2 += c shift2 -= facexp(k, 2) # cancel all 2's in factorial(k) assert shift2 >= 0 # Any prime generator is fine here. `k` can't be # huge, and we only want the primes through `k`. pgen = psieve() p = next(pgen) assert p == 2 for p in pgen: if p > k: break pcount = facexp(k, p) assert pcount # Divide that many p's out of numerators. i = firstnum % p if i: i = p - i for i in range(i, k, p): val, r = divmod(nums[i], p) assert r == 0 pcount -= 1 while pcount: val2, r = divmod(val, p) if r: break else: val = val2 pcount -= 1 nums[i] = val if pcount == 0: break assert pcount == 0 heapify(nums) while len(nums) > 1: a = heappop(nums) heapreplace(nums, a * nums[0]) return nums[0] << shift2 I'm NOT suggesting to adopt this. Just for history in the unlikely case there's worldwide demand for faster `comb` of silly arguments ;-)