[issue8692] Use divide-and-conquer for faster factorials

[Alexander Belopolsky]

Alexander Belopolsky <belopolsky at users.sourceforge.net> added the comment:

On Tue, May 11, 2010 at 5:58 PM, Mark Dickinson <report at bugs.python.org> wrote:
> Yes, I'm interested in seeing the pure Python version.

Here is my translation of the reference implementation.

>Â It could go into test_math, and would be a useful form of documentation.

Note that I copied the reference implementation recursive logic rather
than that in the proposed patch.  It may be better for documentation
this way.

If we end up using something like this in documentation, I would
rename nminusnumofbits() to something more readable.  Maybe "ntwos" or
"count_trailing_zeros" with an explanation why number of factors of 2
in factorial(n) is n - popcount(n).

----------
Added file: http://bugs.python.org/file17310/factorial.py

_______________________________________
Python tracker <report at bugs.python.org>
<http://bugs.python.org/issue8692>
_______________________________________
-------------- next part --------------
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()