Extended Lucas theorem

Posted on In Disqus: Word count in article: 2k Reading time ≈ 7 mins.

中文版本请见这里。

\(\newcommand{\stirlingI}[2]{\genfrac{[}{]}{0pt}{}{#1}{#2}}\) Lucas theorem is quite well known for compute \(\binom{n}{m} \pmod{p}\):

Let \(p\) be a prime, \(n = u_1 p + u_2\)，\(m = v_2 p + v_2\) where \(0 \leq u_2, v_2 < p\), then we have \[ \binom{n}{m} \equiv \binom{u_1}{v_1} \binom{u_2}{v_2} \pmod{p}. \]

However, Lucas theorem only works for prime \(p\). For composites CRT can rescue us, but we still have to solve the hard case where \(p\) is a prime power, which this post mainly discusses. The post borrows a lot from a post by prabowo (although I guess he also references min_25's blogpost), while having its own novelties (mostly the upper bound).

Andrew Granville also has a paper on Extended Lucas theroem. However, there are too many cases in the paper and the overall time complexity is also high. I'd never write it again after I did it once.

Analysis on \(n!\)

Define \(E(n) = \max_k \{p^k | n!\}\) as the number of trailing 0's of the base-\(p\) representation of \(n!\). The efficient algorithm for computing \(E(n)\) is trivial. Now we define \[ (n!)_p = \prod_{i=1, p \nmid i}^n i, \] that is, the product of integers in \([n]\) which are not multiples of \(p\). Now, we have \[ n! = p^{E(n)} \prod_{k=0} \left( \left\lfloor n/p^k \right\rfloor !\right)_p, \] and our task is to compute \((n!)_p \bmod p^e\).

Suppose \(n = up + v\) where \(0 \leq v < m\). We may have two terms for \((n!)_p\): \[ (n!)_p = \left( \prod_{i=0}^{u-1} \prod_{j=1}^{p-1} (ip + j) \right) \left( \prod_{i=1}^v (up + j) \right). \] Note that both terms have a structure of \(\prod_{j=1}^r (ip + j)\), which can be seen as a polynomial of \(p\). What's more interesting is that we only care about the polynomial modulo \(p^e\). More precisely, notice that \[ \prod_{i=0}^{k-1} (n + i) = \sum_{i=0}^k \stirlingI{k}{i} n^i, \] where \(\stirlingI{k}{i}\) is the Stirling number of the first kind satisfying \[ \stirlingI{n + 1}{k} = n \stirlingI{n}{k} + \stirlingI{n}{k-1}, \quad \stirlingI{0}{0} = 1, \stirlingI{0}{n} = \stirlingI{n}{0} = 0. \] Therefore, we may represent \((n!)_p\) by \[ \begin{aligned} (n!)_p & = \left( \prod_{i=0}^{u-1} \prod_{j=1}^{p-1} (ip + j) \right) \left( \prod_{i=1}^v (up + j) \right) \\ & = \left( \prod_{i=0}^{u-1} \sum_{k=1}^{p} \stirlingI{p}{k} (ip)^{k-1} \right) \left( \sum_{k=1}^{v} \stirlingI{v+1}{k} (up)^{k-1} \right) \\ & \equiv \left( \prod_{i=0}^{u-1} \sum_{k=1}^{e} \stirlingI{p}{k} (ip)^{k-1} \right) \left( \sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1} \right) \pmod{p^e} \\ & = \stirlingI{p}{1}^u \left( \prod_{i=0}^{u-1} \sum_{k=1}^{e} \frac{\stirlingI{p}{k}}{\stirlingI{p}{1}} (ip)^{k-1} \right) \left( \sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1} \right). \end{aligned} \] We further define \[ f_{p, e}(u) := \prod_{i=0}^{u-1} \sum_{k=1}^{e} \frac{\stirlingI{p}{k}}{\stirlingI{p}{1}} (ip)^{k-1} \bmod p^e = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1} \frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \bmod{p^e}. \] At this point, we have three terms for \((n!)_p\): \(\stirlingI{p}{1}^u\)，\(f_{p, e}(u)\) and \(\sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1}\). The first term and the last term are straightforward to compute, so we'll focus on the middle term (\(f_{p, e}(u)\)).

Here comes a useful observation: \(f_{p, e}(u)\) is polynomial in \(u\).

Prove \(f_{p, e}(u)\) is polynomial in \(u\)

We here prove a more general statement:

Lemma: Suppose for any \(1 \leq k < e\), \(w_k(u)\) is a polynomial in \(u\), and \[ g(x, u) = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1} w_k(u) x^k \right) \mod x^e, \] then \(g\) is a polynomial in \(u\) with a bounded degree \[ \deg_u g(x, u) \leq \lambda (e-1), \quad \lambda := \max_{1 \leq k < e} \frac{\deg w_k + 1}{k}. \]

Proof: For \(0 \leq k < e\), let \(h_k(u)\) be the coefficient of \(x^k\) in \(g(x, u)\). We will prove that \(h_k(u)\) is polynomial in \(u\) with degree \(\leq \lambda k\). With the property of \(h_k(u)\) ahead, we may simple conclude \[ g(x, u) = \sum_{k=0}^{e-1} h_k(u) x^k \quad \Longrightarrow \quad \deg_u g \leq \max_{0 \leq k < e} \deg h_k \leq \lambda (e-1). \]

Here we use induction on \(k\):

  • When \(k = 0\), \(h_k(u) = 1\).
  • Suppose the statement holds for \(0 < k < e\), we have \[ h_k(u) = h_{k-1}(u-1) + \sum_{t=0}^{k-1} h_t(u-1) w_{k-t}(u - 1) = h_{k-1}(0) + \sum_{t=0}^{k-1} \sum_{u' \geq 0}^{u-1} h_t(u') w_{k-t}(u'). \] As \(h_t(u') w_{k-t}(u')\) is polynomial in \(u'\) with degree \(\leq \deg h_t + \deg w_{k-t}\), \(h_k\) is also polynomial in \(u\) with degree \[ \deg h_k \leq \max_{0 \leq t < k} \left( \deg h_t + \deg w_{k-t} + 1 \right), \] where the additional constant of 1 comes from summation over \(u'\). The definition of \(\lambda\) and induction give \[ \deg h_k \leq \max_{0 \leq t < k} \left( \lambda t + \lambda (k - t) \right) = \lambda k. \] Indeed the recursion on \(\deg h_k\) is just a Knapsack problem, so a trivial upper bound is just the item with largest density.

QED.

Bound the degree of \(f_{p, e}(u)\)

We still have one question to answer: given that \(f_{p, e}(u)\) is a polynomial, how do we bound its degree? We'd like to use the lemma above, but the problem, how do we define \(w_k(u)\)?

If we use the original definition, we may have \[ f_{p, e}(u) = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \mod p^e, \] that is, \(1 \leq k < e\), \(w_k(i) = \frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} i^k\). This definition leads to \[ \lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k} = \max_{1 \leq k < e} \frac{k + 1}{k} = 2, \quad \Longrightarrow \quad \deg_u f_{p, e}(u) \leq 2 (e - 1). \]

However, prabowo hints that: when \(p\) is prime and \(1 < k < p\), \(p\) divides \(\stirlingI{p}{k}\).

We first consider the case of \(e < p\), all relevant \(\stirlingI{p}{k+1}\) is a multiple of \(p\), and we can extract the common factor: \[ \begin{aligned} f_{p, e} & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \mod p^e \\ & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1} / p}{\stirlingI{p}{1}} i^k p^{k+1} \right) \mod p^e. \end{aligned} \] Hence, we use \(w_k(i) = \frac{\stirlingI{p}{k} / p}{\stirlingI{p}{1}} i^{k-1}\) for \(k > 1\) and \(w_1(i) = 0\). The \(\lambda\) here can be made smaller: \[ \lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k} = 1, \quad \Longrightarrow \quad \deg_u f_{p, e}(u) \leq e - 1. \] (It's not rigorous here as we can't divide a number by \(p\) under \(\bmod p^e\), although we may do it in \(\mathbb{R}\). The division here is only used for proof and will not be use in code.)

The other case is \(e \geq p\), and we have one more Stirling number to discuss: \(\stirlingI{p}{p} = 1\). We may re-organize the equations: \[ \begin{aligned} f_{p, e} & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \mod p^e \\ & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{p-2}\frac{\stirlingI{p}{k + 1} / p}{\stirlingI{p}{1}} i^k p^{k+1} + \frac{1}{\stirlingI{p}{1}} (ip)^{p-1} \right) \mod p^e. \end{aligned} \] Hence we define for \(1 < k < p\)，\(w_k(i) = \frac{\stirlingI{p}{k} / p}{\stirlingI{p}{1}} i^k\)，and \(w_{p-1} = \frac{\stirlingI{p}{p-1}/p}{\stirlingI{p}{1}} {i^{p-2}}+ \frac{1}{\stirlingI{p}{1}} i^{p-1}\). In this case, \(\lambda\) is \[ \lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k} = \frac{p}{p - 1}, \quad \Longrightarrow \quad \deg_u f_{p, e}(u) \leq \frac{p}{p - 1} (e - 1). \] As \(p \geq 2\), this bound is always better than \(2(e-1)\).

Overall, we have \[ \deg_u f_{p, e}(u) \leq d^* := \left\lfloor \frac{p (e - 1)}{p-1} \right \rfloor. \] How nicely it captures both cases!

The divisibility of the Stirling numbers

The proof above makes use of a special property of Stirling numbers of the first kind, which we'll give a short proof here:

For prime \(p\) and \(1 < k < p\), \(p\) divides \(\stirlingI{p}{k}\).

Proof: Consider a polynomial \(f(x) = \prod\limits_{i=1}^{p-1} (x+i) \in \mathbb{F}_p[x]\). The Stirling numbers of the first kind satisfy \[f(x) = \prod_{i=1}^{p-1} (x + i) = \sum_{i=1}^p \stirlingI{p}{i} x^{i-1}.\] Note that \(f(0) = (p-1)! \equiv -1 \pmod{p}\), as well as \(f(x) = 0\) for any \(x \neq 0\). On the other hand, the polynomial \(g(x) = x^{p-1} - 1\) has identical values as \(f\) in these \(p\) points, and both \(f\) and \(g\) have degree of \(p - 1\), thus \(f = g\). Now we compare the coefficients of \(f\) and \(g\): \[ \stirlingI{p}{i} \equiv \begin{cases} 1, & i = p, \\ 0, & 1 < i < p, \\ -1, & i = 1.\end{cases} \pmod{p} \]

Algorithm for computing \((n!)_p \bmod p^e\)

  • For \(0 \leq v \leq p, 0 \leq k \leq e\), preprocess \(\stirlingI{v}{k}\) with overall complexity \(O(pe)\);
  • Preprocess \(f_{p, e}\). As a polynomial with degree of \(d^*\) requires \(d^* + 1\) points to determine, we may simply brute-force \(f_{p, e}(u)\) for \(0 \leq u \leq d^*\) and use Lagrangian interpolation. The complexity is \(O(e^2)\).
  • For each query, we simply use the following equation: \[ (n!)_p \equiv \stirlingI{p}{1}^u f_{p, e}(u) \left( \sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1} \right) \pmod{p^e}, \] with time complexity of \(O(e)\) per query.

If \(d^* \geq p\), we may be careful in Lagrange interpolation. The following equation can be helpful, \[ f_{p, e}(u) = \sum_{i=0}^{d^*} f_{p, e}(i) \prod_{j \neq i} \frac{u - j}{i - j} = \sum_{i=0}^{d^*} f_{p, e}(i) \binom{u}{j}\binom{u-j-1}{i-j} (-1)^{i-j}, \] which means that it must be an integer. However, \(\prod_{j \neq i} \frac{u - j}{i - j}\) might be a little tricky as \(i - j\) might not be invertible under \(\bmod{p^e}\)... Anyway we have to represent all numbers in the form of \(r \equiv \frac{n}{p^a} \pmod {p^e}\) and do multiplication/division later...

For reference code, prabowo provides a C++ version while I can share my own Python version. Unfortunately, min_25's Python version was lost.

Hexo

When I wrote the blog post for the first time, I use commands like \newcommand{\aaa}[1]{#1} while Hexo recognize # as SwigTags, which leads to a rendering issue. The solution is quite simple: just add 

1
disableNunjucks: true
to the front-matter. Please refer to this issue (there is a typo in disableNunjucks in the comment though) and this issue.