import functools
import operator

product = functools.partial(functools.reduce, operator.mul)

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):
    """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))
    """
    _, r = loop(n, 1, 1)
    return r << nminusnumofbits(n)

def loop(n, p, r):
    if n > 2:
        p, r = loop(n // 2, p, r)
        p *= partial_product(n // 2 + 1 + ((n // 2) & 1),  n - 1 + (n & 1))
        r *= p
    return p, r

def partial_product(n, m):
    if m <= n + 1:
        return n
    if m == n + 2:
        return n * m
    k = (n + m) // 2
    if k & 1 != 1:
        k -= 1
    return partial_product(n, k) * partial_product(k + 2, m)

def nminusnumofbits(n):
    nb = 0
    v = n
    while v:
        nb += v & 1
        v >>= 1
    return n - nb


if __name__ == '__main__':
    import doctest
    doctest.testmod()