import functools import operator product = functools.partial(functools.reduce, operator.mul) def factorial(n): """Implementation of Binary-Split Factorial algorithm See http://www.luschny.de/math/factorial/binarysplitfact.html >>> f = 1 >>> for n in range(1, 1001): ... f *= n ... assert(factorial(n) == f) """ if n < 3: return [1, 1, 2][n] _, r = loop(n) return r << (n - count_bits(n)) def loop_iter(n): p = r = 1 s = n.bit_length() - 2 i = s - 1 + (n >> s & 1) while i >= 0: m = n >> i p *= partial_product((m >> 2) + (m >> 1 & 1), (m >> 1) + (m & 1)) r *= p i -= 1 return p, r def loop(n): p = r = 1 if n > 2: p, r = loop(n >> 1) p *= partial_product((n >> 2) + (n >> 1 & 1), (n >> 1) + (n & 1)) r *= p assert(n < 3 or (p, r) == loop_iter(n)) return p, r def partial_product(start, stop): length = stop - start if length == 1: return start << 1 | 1 if length == 2: x = (start << 1 | 1) return x * (x + 2) middle = start + (length >> 1) return partial_product(start, middle) * partial_product(middle, stop) def count_bits(n): count = 0 while n: n &= n - 1 count += 1 return count if __name__ == '__main__': import doctest doctest.testmod()