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.