﻿ Binary Split Factorial Formula

# The binary-split formulaof the factorial of n.

```  bigint Factorial(int n)
{
bigint p = 1, r = 1;
loop(n, p, r);
return r << nminussumofbits(n);
}

loop(int n, reference bigint p, reference bigint r)
{
if (n <= 2) return;
loop(n / 2, p, r);
p = p * partProduct(n / 2 + 1 + ((n / 2) & 1), n - 1 + (n & 1));
r = r * p;
}

bigint partProduct(int n, int m)
{
if (m <= (n + 1)) return (bigint) n;
if (m == (n + 2)) return (bigint) n * m;
int k =  (n + m) / 2;
if ((k & 1) != 1) k = k - 1;
return partProduct(n, k) * partProduct(k + 2, m);
}

int nminussumofbits(int v)
{
long w = (long)v;
w -= (0xaaaaaaaa & w) >> 1;
w = (w & 0x33333333) + ((w >> 2) & 0x33333333);
w = w + (w >> 4) & 0x0f0f0f0f;
w += w >> 8;
w += w >> 16;
return v - (int)(w & 0xff);
}
```
```   For example let's trace the computation of 30!.

init: p = 1, r = 1, then loop:

p_1          = 3
p = p * p_1  = 3
r = r * p    = 3

p_2          = 5 * 7
p = p * p_2  = 105
r = r * p    = 315

p_3          = 9 * 11 * 13 * 15
p = p * p_3  = 2027025
r = r * p    = 638512875

p_4          = 17 * 19 * 21 * 23 * 25 * 27 * 29
p = p * p_4  = 6190283353629375
r = r * p    = 3952575621190533915703125

return r * 2^26 = 265252859812191058636308480000000

------

The partial products p_n in turn are computed via the standard
recursive product, i.e. evaluated (almost) like binary trees:

|                    |                3053876175
|        19305       |        156009              19575
35   |    99       195    |     323       483       675     29
5 * 7  |  9 * 11 * 13 * 15  |   17 * 19 * 21 * 23 * 25 * 27 * 29

-----

Depending on your temperament you now can start to fiddle with the
implementation, in dependency on your favorite computer language.
```

```# A pure Python version.

# Returns the number of bits necessary to represent an integer in binary,
# excluding the sign and leading zeros.
# Needed only for Python version < 3.0; otherwise use n.bit_length().
def bit_length(self):
s = bin(self)       # binary representation:  bin(-37) --> '-0b100101'
s = s.lstrip('-0b') # remove leading zeros and minus sign
return len(s)       # len('100101') --> 6

def num_of_set_bits(i) :
# assert 0 <= i < 0x100000000
i = i - ((i >> 1) & 0x55555555)
i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
return (((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) & 0xffffffff) >> 24

def rec_product(start, stop):
"""Product of integers in range(start, stop, 2), computed recursively.
start and stop should both be odd, with start <= stop.
"""
numfactors = (stop - start) >> 1
if numfactors == 2 : return start * (start + 2)
if numfactors > 1 :
mid = (start + numfactors) | 1
return rec_product(start, mid) * rec_product(mid, stop)
if numfactors == 1 : return start
return 1

def binsplit_factorial(n):
"""Factorial of nonnegative integer n, via binary split.
"""
inner = outer = 1
for i in range(n.bit_length(), -1, -1):
inner *= rec_product((n >> i + 1) + 1 | 1, (n >> i) + 1 | 1)
outer *= inner
return outer << (n - num_of_set_bits(n))

# Test (from math import factorial).
[[n, binsplit_factorial(n) - factorial(n)] for n in range(99)]
```

Implementations in Java and in C# can be found in the index of the algorithms under the name 'Split Recursive'. They differ from the above pseudo code by incorporating the recursive function 'loop' as a 'while' loop into the factorial function.

The above algorithm was implemented in Python's mathmodule. John Benediktsson writes that it was one of the improvements made during the Python 3.2 development cycle. This is a little surprising as I would have expected a prime-based algorithm for such a library. Here you can see a Python implementation of such an prime algorithm which is shorter than the mathmodule-implementation and much more efficient for large arguments. Here further remarks on related SymPy and Sage implementations.

Back to the Homepage of Factorial Algorithms.