Message 331859 - Python tracker

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 ;-)