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);
}
# 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)]