Tim Peters <t...@python.org> added the comment:
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 ;-)
----------
_______________________________________
Python tracker <rep...@bugs.python.org>
<https://bugs.python.org/issue35431>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe:
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
