Today, we will be computing $$$n! \bmod p$$$ in $$$O(\sqrt n \log n)$$$ time.

Assume all math in this article is done $$$\bmod p$$$.

Using a square root decomposition idea, we let $$$v = \text{floor}(\sqrt n)$$$ and $$$f(x) = (x + 1)(x+2)\cdots(x+v)$$$. Then $$$f(0)$$$ gives the the product of first $$$v$$$ terms in the factorial product, $$$f(v)$$$ gives the next $$$v$$$, and so on, so $$$f(0)f(v)\cdots f((v - 1)v) = v!$$$ and we can multiply the remaining terms $$$v + 1, v + 2, \dots, n$$$ (of which there are at most $$$\sqrt n$$$) on to get $$$n!$$$.

Computing the coefficients of $$$f(x)$$$ and doing multi-point evaluation on $$$x = 0, v, 2v, \dots, (v - 1)v$$$ each takes $$$O(v \log^2 v) = O(\sqrt n \log^2 n)$$$, but we can do better and save a lot of of FFTs by skipping computing $$$f(x)$$$ directly.

The Algorithm

Let $$$f_d(x) = (x + 1)(x + 2) \cdots (x + d)$$$. We will show an algorithm to efficient compute $$$f_{2d}(x)$$$ from $$$f_d(x)$$$, then use this algorithm to compute $$$f_1(x)$$$, then $$$f_2(x)$$$, then $$$f_4(x)$$$ and so on, until we compute $$$f(x) = f_v(x)$$$. One small modification we have to make is rounding $$$v$$$ up to the nearest power of 2.

It was mentioned that we won't be computing the coefficients of $$$f(x)$$$. For each $$$f_d(x)$$$, we will instead be computing $$$S_d = [f_d(0), f_d(v), f_d(2v), \dots, f_d(dv)]$$$. This is an array of samples of $$$f_d(x)$$$ of size $$$d + 1$$$, which is enough to uniquely determine $$$f_d(x)$$$, so it works as an alternative representation.

Next, we tackle the problem of computing $$$f_{2d}(x)$$$ from $$$f_d(x)$$$, or more specifically, computing $$$S_{2d}$$$ from $$$S_{d}$$$. The following identity gives us a hint:

Our tool of choice will be the Polynomial Sample Shift. Specifically, let $$$S_{d, a} = [f(a), f(a + v), f(a + 2v), \dots, f(a + dv)]$$$, then the Polynomial Sample Shift is an algorithm that takes in $$$S_d = S_{d, 0}$$$ and $$$a$$$ as input and gives us $$$S_{d, a}$$$ in $$$O(d \log d)$$$ time. To get the first $$$d + 1$$$ elements of $$$S_{2d}$$$, we do one shift by $$$d$$$ to get $$$S_{d, d} = [f_d(d), f_d(d + v), \dots, f_d(d + dv)]$$$ and pointwise multiply with $$$S_d$$$ to get $$$[f_{2d}(0), f_{2d}(v), \dots, f_{2d}(dv)]$$$ since $$$f_{2d}(x) = f_d(x)f_d(x + d)$$$.

To get the rest of the terms, we shift two more times by $$$dv+d$$$ and $$$dv+2d$$$, getting $$$S_{d, dv+d}$$$ and $$$S_{d, dv+2d}$$$ and pointwise multiply the arrays to get $$$[f_{2d}(dv + d), f_{2d}(dv + v), \dots, f_{2d}(2dv), f_{2d}(2dv + v)]$$$. Finally, we concatenate these two arrays together and discard the last value to get $$$[f_{2d}(0), f_{2d}(v), \dots, f_{2d}(2dv)] = S_{2d}$$$ as desired.

After computing $$$S_{1} = [f(0), f(v)] = [1, v + 1]$$$, then $$$S_2, S_4, \dots, S_v$$$, we now have $$$S_v = [f(0), f(v), \dots, f(v^2)]$$$. Since $$$v$$$ can be up to twice as large as $$$\sqrt n$$$, we don't necessarily stop at $$$f((v - 1)v)$$$; instead we compute $$$f(0)f(v)\cdots f(l v)$$$ where $$$l$$$ is the largest number such that $$$lv + v \leq n$$$, then multiply on the rest of the numbers $$$(lv + v + 1)(lv + v + 2)\cdots n$$$ to get $$$n!$$$.

Our time complexity is $$$T(v) = T(v / 2) + O(v \log v)$$$ which is $$$O(v \log v)$$$ or $$$O(\sqrt n \log n)$$$.

Summary

The algorithm is as follows:

Input: $$$n$$$ and $$$p$$$.

Output: $$$n! \bmod p$$$.

1. Choose $$$v$$$ to be the smallest power of 2 at least $$$\sqrt n$$$.
2. Compute $$$S_1 = [f_1(0), f_1(v)] = [1, v + 1]$$$.
3. Let $$$d = 1$$$. While $$$d \lt v$$$, compute $$$S_{2d}$$$ from $$$S_{d}$$$ by doing the following:
   • Perform the Polynomial Sample Shift on $$$S_{d}$$$ 3 times, computing $$$S_{d, d}$$$, $$$S_{d, dv + d}$$$, and $$$S_{d, dv + 2d}$$$
   • Pointwise multiply $$$S_d$$$ and $$$S_{d, d}$$$ to get the first $$$d+1$$$ terms of $$$S_{2d}$$$
   • Pointwise multiply $$$S_{dv + d}$$$ and $$$S_{d, dv + 2d}$$$ and discard the last term to get the last $$$d$$$ terms of $$$S_{2d}$$$
   • Concatenate the two to get $$$S_{2d}$$$
   • Multiply $$$d$$$ by 2
4. We now have $$$S_{v} = [f(0), f(v), \dots, f(v^2)]$$$. Initialize $$$\text{prod} = 1, m = 0$$$.
5. While $$$m + v \leq n$$$:
   • Multiply $$$\text{prod}$$$ by $$$f(m)$$$
   • Increment $$$m$$$ by $$$v$$$
6. Multiply $$$m + 1, m + 2, \dots, n$$$ onto $$$\text{prod}$$$ to get $$$n!$$$

Practice Problems

• [Library Checker] Factorial
• [DMOJ] Fast Factorial Calculator 3
• [SPOJ] Factorial Modulo Prime

Extensions

It is known that computing the $$$n$$$-th term of a linear recurrence that uses the last $$$k$$$ terms can be computed in $$$O(k^3 \log n)$$$ using matrix exponentiation or $$$O(k \log k \log n)$$$ using the Kitamasa Method.

On the other hand, if the recurrence has coefficients that are not constant, but also depend on $$$i$$$, we can instead apply this technique to achieve $$$O(k^3 \sqrt n + k^2 \sqrt n \log n)$$$.

We use the example of the subfactorial $$$!n$$$ which can be used to count derangements, the number of permutations of $$$n$$$ objects that leaves no object where it started. We have the recurrence $$$!(n + 1) = n \times (!n\ +\ !(n - 1))$$$. In matrix form, this is

Instead of defining an $$$f(x) = (x + 1)\cdots(x+v)$$$, we define

as a $$$k\times k$$$ matrix of polynomials of degree $$$d$$$, we can apply this technique directly, with two main differences:

1. $$$S_{d, a}$$$ is an array of $$$k\times k$$$ matrices, so pointwise multiplication is now doing matrix multiplication for each element, which is why we have a $$$O(k^3 \sqrt n)$$$ term in the time complexity.
2. In polynomial sample shift, the $$$f(a)$$$ samples never multiply by other matrices, so we can perform it independently for each of the $$$k^2$$$ matrix entries of the samples, which is the $$$O(k^2 \sqrt n \log n)$$$ term.

You can find a list of things you can do with this from this article's source.

Sources / Additional Reading

This article is based on this blog post by Min25.

You can also read about this technique on this blog post by ho94949.