# Factorial algorithm more efficient than naive multiplication

I know how to code for factorials using both iterative and recursive (e.g. n * factorial(n-1) for e.g.). I read in a textbook (without been given any further explanations) that there is an even more efficient way of coding for factorials by dividing them in half recursively.

I understand why that may be the case. However I wanted to try coding it on my own, and I don't think I know where to start though. A friend suggested I write base cases first. and I was thinking of using arrays so that I can keep track of the numbers... but I really can't see any way out to designing such a code.

What kind of techniques should I be researching?

The best algorithm that is known is to express the factorial as a product of prime powers. One can quickly determine the primes as well as the right power for each prime using a sieve approach. Computing each power can be done efficiently using repeated squaring, and then the factors are multiplied together. This was described by Peter B. Borwein, On the Complexity of Calculating Factorials, Journal of Algorithms 6 376–380, 1985. (PDF) In short, $n!$ can be computed in $O(n(\log n)^3\log \log n)$ time, compared to the $\Omega(n^2 \log n)$ time required when using the definition.

What the textbook perhaps meant was the divide-and-conquer method. One can reduce the $n-1$ multiplications by using the regular pattern of the product.

Let $n?$ denote $1 \cdot 3 \cdot 5 \dotsm (2n-1)$ as a convenient notation. Rearrange the factors of $(2n)! = 1 \cdot 2 \cdot 3 \dotsm (2n)$ as $$(2n)! = n! \cdot 2^n \cdot 3 \cdot 5 \cdot 7 \dotsm (2n-1).$$ Now suppose $n = 2^k$ for some integer $k>0$. (This is a useful assumption to avoid complications in the following discussion, and the idea can be extended to general $n$.) Then $(2^k)! = (2^{k-1})!2^{2^{k-1}}(2^{k-1})?$ and by expanding this recurrence, $$(2^k)! = \left(2^{2^{k-1}+2^{k-2}+\dots+2^0}\right) \prod_{i=0}^{k-1} (2^i)? = \left(2^{2^k - 1}\right) \prod_{i=1}^{k-1} (2^i)?.$$ Computing $(2^{k-1})?$ and multiplying the partial products at each stage takes $(k-2) + 2^{k-1} - 2$ multiplications. This is an improvement of a factor of nearly $2$ from $2^k-2$ multiplications just using the definition. Some additional operations are required to compute the power of $2$, but in binary arithmetic this can be done cheaply (depending on what precisely is required, it may just require adding a suffix of $2^k-1$ zeroes).

The following Ruby code implements a simplified version of this. This does not avoid recomputing $n?$ even where it could do so:

def oddprod(l,h)
p = 1
ml = (l%2>0) ? l : (l+1)
mh = (h%2>0) ? h : (h-1)
while ml <= mh do
p = p * ml
ml = ml + 2
end
p
end

def fact(k)
f = 1
for i in 1..k-1
f *= oddprod(3, 2 ** (i + 1) - 1)
end
2 ** (2 ** k - 1) * f
end

print fact(15)


Even this first-pass code improves on the trivial

f = 1; (1..32768).map{ |i| f *= i }; print f


by about 20% in my testing.

With a bit of work, this can be improved further, also removing the requirement that $n$ be a power of $2$ (see the extensive discussion).

• You left out an important factor. The calculation time according to Borwein's paper is not O (n log n log log n). It is O (M (n log n) log log n), where M (n log n) is the time for multiplying two numbers of size n log n. – gnasher729 Jul 30 '16 at 19:01

Keep in mind that the factorial function grows so fast that you'll need arbitrary-sized integers to get any benefit of more efficient techniques than the naive approach. The factorial of 21 is already too big to fit in a 64-bit unsigned long long int.

As far as I know, there is no algorithm to compute $n!$ (factorial of $n$) which is faster than doing the multiplications.¹

However, the order in which you do the multiplications matter. Multiplication on a machine integer is a basic operation that takes the same time no matter what the value of the integer is. But for arbitrary-sized integers, the time it takes to multiply a and b depends on the size of a and b: a naive algorithm operates in time $\Theta(|a| \cdot |b|)$ (where $|x|$ is the number of digits of $x$ — in whatever base you like, as the result is the same up to a multiplicative constant). There are faster multiplication algorithms, but there is an obvious lower bound of $\Omega(|a| + |b|)$ since multiplication has to at least read all the digits. All known multiplication algorithms grow faster than linearly in $\max(|a|,|b|)$.

Armed with this background, the Wikipedia article should make sense.

Since the complexity of multiplications depend on the size of the integers that are being multiplied, you can save time by arranging multiplications in an order that keeps the numbers being multiplied small. It works out better if you arrange for the numbers to be of roughly the same size. The “division in half” that your textbook refers to consists of the following divide-and-conquer approach to multiply a (multi)set of integers:

1. Arrange the numbers to be multiplied (initially, all the integers from $1$ to $n$) in two sets whose product is roughly the same size. This is a lot less expensive than doing the multiplication: $|a \cdot b| \approx |a| + |b|$ (one machine addition).
2. Apply the algorithm recursively on each of the two subsets.
3. Multiply the two intermediate results.

See the GMP manual for more specifics.

There are even faster methods that not only rearrange the factors $1$ to $n$ but split the numbers by decomposing them into their prime factorization and rearranging the resulting very long product of mostly-small integers. I'll just cite the references from the Wikipedia article: “On the Complexity of Calculating Factorials” by Peter Borwein and implementations by Peter Luschny.

¹ There are faster ways of computing approximations of $n!$, but that's not computing the factorial any more, it's computing an approximation of it.

Since the factorial function grows so fast, your computer can only store $n!$ for relatively small $n$. For example, a double can store values up to $171!$. So if you want a really fast algorithm for computing $n!$, just use a table of size $171$.

The question becomes more interesting if you're interested in $\log(n!)$ or in the $\Gamma$ function (or in $\log \Gamma$). In all these cases (including $n!$), I don't really understand the comment in your textbook.

As an aside, your iterative and recursive algorithms are equivalent (up to floating point errors), since you are using tail recursion.

• "your iterative and recursive algorithms are equivalent" you're referring to the their asymptotic complexity, right? as for the comment in the textbook, well i'm translating it from another language so, maybe my translation sucks. – user65165 Sep 20 '13 at 2:13
• The book talks about iterative and recursive, and then comments on how if you use divide and conquer to divide n! in half you can get a way faster solution... – user65165 Sep 20 '13 at 2:14
• My notion of equivalence isn't completely formal, but you could say that the arithmetic operations performed are the same (if you switch the order of the operands in the recursive algorithm). An "inherently" different algorithm will perform a different calculation, perhaps using some "trick". – Yuval Filmus Sep 20 '13 at 2:19
• If you consider the size of the integer as a parameter in the complexity of multiplication, the overall complexity can change even if the arithmetic operations are "the same". – Tpecatte Sep 20 '13 at 7:31
• @CharlesOkwuagwu Right, you could use a table. – Yuval Filmus Feb 8 '19 at 3:04

Is the text book written Fundamentals of Algorithms, or Art of Programming by Knuth?

The following method is very effective:

1. Compute binomial coefficient (n, r), using the formula, by repeated multiplication,

(1+b)^n = 1 + (n,1)b + (n, 2)b^2 + ..... + (n, r)b^r + .... If you choose b large enough, for example, b = 2^n, there is no carry generated between the coefficients, so, you can select (n, r)

1. Suppose you want to compute m!.

Select N, much larger than m, for example, say N = m^m, and compute

K = integer(N^m/(N, m)),. Then K = m!

Here both exponent is computed by repeated multiplication.

Example : 8! = 40320 = integer( N^8/(N, 8)) for any N ≥ 1128977