import functools import operator product = functools.partial(functools.reduce, operator.mul) def partial_product(j, i): if i == j: return j << 1 | 1 if i == j + 1: return (j << 1 | 1) * (i << 1 | 1) l = i + j >> 1 return partial_product(j, l) * partial_product(l + 1, i) def naive_factorial(n): """Naive implementation of factorial: product([1, ..., n]) >>> naive_factorial(4) 24 """ return product(range(1, n+1), 1) def factorial(n, pp=partial_product): """Implementation of Binary-Split Factorial algorithm See http://www.luschny.de/math/factorial/binarysplitfact.html >>> for n in range(20): ... assert(factorial(n) == naive_factorial(n)) >>> import math >>> assert(factorial(100) == math.factorial(100)) """ _, r = loop(n, pp) return r << (n - count_bits(n)) def loop(n, pp): p = r = 1 for i in range(n.bit_length() - 2, -1, -1): m = n >> i if m > 2: p *= pp(((m >> 1) + 1) >> 1, (m - 1) >> 1) r *= p return p, r def partial_product1(j, i): return product((l << 1 | 1 for l in range(j, i + 1)), 1) def partial_product2(j, i): a = [l << 1 | 1 for l in range(j, i + 1)] n = i - j + 1 p = 1 while n > 1: for k in range(n>>1): a[k] *= a[n-k-1] n = (n>>1) + (n&1) return a[0] def partial_product3(j, i): a = [l << 1 | 1 for l in range(j, i + 1)] while 1: n = len(a) if n == 1: return a[0] a = [a[k<<1] * a[k<<1|1] for k in range(n>>1)] + a[(n >> 1)<<1:] def partial_product4(j, i): if i == j: return j << 1 | 1 if i == j + 1: return (j << 1 | 1) * (i << 1 | 1) return (j << 1 | 1) * (i << 1 | 1) * partial_product4(j + 1, i - 1) def count_bits(n): count = 0 while n: n &= n - 1 count += 1 return count def f0(n): return factorial(n) def f1(n): return factorial(n, partial_product1) def f2(n): return factorial(n, partial_product2) def f3(n): return factorial(n, partial_product3) if __name__ == '__main__': for n in list(range(12, 20)) + [100, 101, 200]: assert(f0(n) == f1(n) == f2(n) == f3(n)) import doctest doctest.testmod()