#### 15.7.2 Factorial

Factorials are calculated by a combination of two algorithms. An idea is
shared among them: to compute the odd part of the factorial; a final step
takes account of the power of 2 term, by shifting.

For small n, the odd factor of n! is computed with the simple
observation that it is equal to the product of all positive odd numbers
smaller than n times the odd factor of [n/2]!,
where [x] is the integer part of x, and so on
recursively. The procedure can be best illustrated with an example,

23! = (23.21.19.17.15.13.11.9.7.5.3)(11.9.7.5.3)(5.3)2^19

Current code collects all the factors in a single list, with a loop and no
recursion, and compute the product, with no special care for repeated chunks.

When n is larger, computation pass trough prime sieving. An helper
function is used, as suggested by Peter Luschny:

n
-----
n! | | L(p,n)
msf(n) = -------------- = | | p
[n/2]!^2.2^k p=3

Where p ranges on odd prime numbers. The exponent k is chosen to
obtain an odd integer number: k is the number of 1 bits in the binary
representation of [n/2]. The function L(p,n)
can be defined as zero when p is composite, and, for any prime
p, it is computed with:

---
\ n
L(p,n) = / [---] mod 2 <= log (n) .
--- p^i p
i>0

With this helper function, we are able to compute the odd part of n!
using the recursion implied by n!=[n/2]!^2*msf(n)*2^k. The recursion stops using the
small-n algorithm on some [n/2^i].

Both the above algorithms use binary splitting to compute the product of many
small factors. At first as many products as possible are accumulated in a
single register, generating a list of factors that fit in a machine word. This
list is then split into halves, and the product is computed recursively.

Such splitting is more efficient than repeated Nx1 multiplies since it
forms big multiplies, allowing Karatsuba and higher algorithms to be used.
And even below the Karatsuba threshold a big block of work can be more
efficient for the basecase algorithm.