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.

Today, we will be solving the Polynomial Sample Shift problem with in $$$O(n \log n)$$$ time using one convolution.

Problem: We are told that $$$f(x)$$$ is a polynomial with degree $$$d$$$ and we are given the samples $$$f(0), f(1), \dots, f(d)$$$, as well as a number $$$m$$$, and are asked to compute $$$f(m), f(m + 1), \dots, f(m + d)$$$.

Since the polynomial is degree $$$d$$$ and we are given $$$d + 1$$$ samples, $$$f(x)$$$ is uniquely defined. Using polynomial interpolation and then multi-point evaluation we can compute this in $$$O(d \log^2 d)$$$ using FFTs, but for this specific problem we can do better.

We start off with the Lagrange Interpolation Formula:

The goal is to plug in $$$x = m + k$$$ and get a formula in terms of $$$f(0), f(1), \dots, f(d)$$$.

We can precompute the product parts of this formula. First, the product in the denominator, letting

We can precompute this for $$$i = 0, 1, \dots, d$$$ by directly computing $$$p[0] = \prod_{j=0,j\neq i}^d (0-j)$$$, then using the formula $$$p[i] = \frac i{i-d-1}\ {p[i-1]}$$$ for the rest.

Next, we can remove the $$$j\neq i$$$ from the product on the right allowing that term to be in the product, but dividing by it, i.e.

We also want to precompute this product. Let

We can similarly precompute this for $$$k = 0, 1, \dots, d$$$ by directly computing $$$q[0] = \prod_{j=0}^d (m - j)$$$ then using the formula $$$q[k] = \frac{m+k}{m+k-d-1}\ q[k -1]$$$ for the rest.

So we have

Now each term sum is the product of two parts: one in terms of $$$i$$$, and one in term of $$$-i$$$, so we can express this as a convolution.

Let

be polynomials with degree $$$d$$$ and $$$2d$$$ respectively.

Let $$$[x^i]p[x]$$$ denote the $$$i$$$-th coefficient of polynomial $$$p[x]$$$. The $$$c$$$-th term of the product is

In particular, the $$$(k + d)$$$-th coefficient of the product is

which matches exactly. Therefore, with $$$O(d)$$$ precomputation of $$$p[i]$$$ and $$$q[k]$$$ plus one FFT for the polynomial product, we have computed

for all $$$k = 0, 1, \dots, d$$$, which is what we wanted.

Summary

The algorithm is as follows:

Input: An array $$$[f(0), f(1), \dots, f(d)]$$$, and a number $$$m$$$.

Output: An array $$$[f(m), f(m + 1), \dots, f(m + d)]$$$.

1. Compute $$$p[0] = \prod_{j = 0, j \neq i}^d (-j)$$$ and $$$p[i]$$$ for all $$$i = 1, \dots, d$$$ with $$$p[i] = \frac i{i - d - 1}p[i - 1]$$$
2. Compute $$$q[0] = \prod_{j=0}^d (m - j)$$$ and $$$q[k]$$$ for all $$$k = 1, \dots, d$$$ with $$$q[k] = \frac{m+k}{m+k-d-1}q[k-1]$$$
3. Compute $$$a[i] = \frac{f(i)}{p[i]}$$$ for $$$i = 0, 1, \dots, d$$$ and $$$b[i] = \frac1{m+i-d}$$$ for $$$i = 0, 1, \dots, 2d$$$.
4. Compute $$$c[i]$$$, the convolution of $$$a[i]$$$ and $$$b[i]$$$.
5. Pointwise multiply the coefficients of $$$q[k]$$$ and the $$$k$$$-th to $$$k + d$$$-th term of $$$c[i]$$$, and output.

Applications

Evenly Spaced Samples

It's worth noting that the samples we're given don't have to be sequential, just evenly spaced apart. Suppose we are told that $$$g(x)$$$ is a polynomial with degree $$$d$$$ and we are given the samples $$$g(a), g(a + b), g(a + 2b), \dots, g(a + db)$$$, and a number $$$m$$$ and are asked to compute $$$g(a + m), g(a + b + m), g(a + 2b + m), \dots, g(a + db + m)$$$.

We define $$$f(x) = g(a + bx)$$$ so that the samples given to us can be interpreted as $$$f(0), f(1), f(2), \dots, f(d)$$$. Performing a shift by $$$\frac mb$$$, we get $$$f(\frac mb), f(1 + \frac mb), f(2 + \frac{m}b), \dots, f(d + \frac mb)$$$ which are $$$g(a + m), g(a + b + m), g(a + 2b + m), \dots, g(a + db + m)$$$, as desired.

Fast Factorial Modulo p

Calculating factorials can be done with multi-point evaluation for a polynomial whose roots and sampling points are evenly spaced. Using this technique, we can calculate $$$n! \bmod p$$$ in $$$O(\sqrt n \log n)$$$ time. You can read about this technique in my blog post, or from here or here.

Problems

• Shift of Sampling Points of Polynomial

Additional Notes

1. Technically the polynomial we are given is promised to have degree at most $$$d$$$, as we do not care whether e.g. the $$$x^d$$$, $$$x^{d-1}, \dots$$$ terms are 0.
2. If we attempt to shift by a value such that the output samples overlaps the given samples, which happens if $$$m$$$ is one of $$$0, 1, \dots, d$$$, then there is a division-by-zero in the above formula. In such a case, we should choose $$$m = d + 1$$$ and discard the extra values.
3. Looking closer at $$$p[i] = \prod_{j = 0, j \neq i}^d (i - j)$$$, this product has $$$d$$$ terms, the $$$i$$$ positive ones span from $$$1$$$ to $$$i$$$ and the $$$d - i$$$ negative ones span from $$$-1$$$ to $$$-(d - i)$$$, so we can alternatively write $$$p[i] = i! (d - i)! (-1)^{d-i}$$$. Either way, it is $$$O(d)$$$ time to compute all $$$d + 1$$$ values.

Sources

• The math and derivations are based on Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator by Alin Bostan, Pierrick Gaudry, and Eric Schost.