Lucas Theorem

Prerequisite knowledge: factorial modulo

Introduction

This article discusses the solution of modulus of large combination numbers. Combinatorial numbers, also known as binomial coefficients, refer to the expression:

When the scale is small, the combination number can be solved by recurrence formula with a time complexity of $O(nk)$; it can also be solved in $O(n)$ time by calculating the factorial of the numerator and denominator under a larger prime modulus $p>n$. But when the problem scale is very large ($n\sim 10^{18}$), these methods are no longer applicable.

Based on Lucas' theorem and its extension, this article discusses a method that can solve the combinatorial number when the modulus is not too large ($m \sim 10^6$). More precisely, this method can be used whenever the unique decomposition of the modulus $m=\prod p_i^{e_i}$ has a sum of prime powers (that is, $\sum p_i^{e_i}$) on the order of $10^6$, since the preprocessing of the algorithm corresponds roughly to this scale.

Lucas Theorem

First, we discuss the case where the modulus is a prime number $p$. In this case, Lucas theorem states:

On the left-hand side of the equation, the coefficient of the term $x^k$ is

$$\binom{n}{k}\bmod p.$$

We now calculate the coefficient of the middle term $x^k$ on the right side of the equation. The degree of each term in the first factor must be a multiple of $p$, and the degree of each term in the second factor must be smaller than $p$. The only way to decompose $k$ into the sum of these two parts is by division with remainder: $k=p\lfloor k/p\rfloor +(k\bmod p)$. Therefore, the first factor can only contribute its $p\lfloor k/p\rfloor$ order terms, and the second factor can only contribute its $k\bmod p$ order terms. Therefore, the $x^k$ coefficient in the equation on the right is the product of the coefficients of the terms contributed by the two factors:

$$\binom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}\binom{n\bmod p}{k\bmod p}\bmod p.$$

Equating the coefficients on both sides are equal, we get Lucas' theorem.

The power $\nu_p(n!)$ and factorial remainder $(n!)_p\bmod p$:

$$\begin{aligned} \nu_p(n!) &= \lfloor n/p\rfloor+\nu_p( \lfloor n/p\rfloor!),\\ (n!)_p &\equiv (-1)^{\lfloor n/p\rfloor}\cdot (n\bmod p)!\cdot (\lfloor n/p\rfloor!)_p\pmod p. \end{aligned}$$

The former is the corollary of Legendre's formula, and the latter is the corollary of Wilson's theorem.

Substituting the recurrence into the expression of the binomial coefficient and sorting it out, we get:

$$\begin{aligned} \binom{n}{k} &\equiv (-p)^{\lfloor n/p\rfloor-\lfloor k/p\rfloor-\lfloor(n-k)/p\rfloor}\cdot\dfrac{(n\bmod p)!}{(k\bmod p)!((n-k)\bmod p)!} \\ &\quad \cdot p^{\nu_p(\lfloor n/p\rfloor!)-\nu_p(\lfloor k/p\rfloor!)-\nu_p(\lfloor(n-k)/p\rfloor!)}\dfrac{(\lfloor n/p\rfloor!)_p}{(\lfloor k/p\rfloor!)_p(\lfloor(n-k)/p\rfloor!)_p} \pmod p. \end{aligned}$$

Now examine the value of $\lfloor n/p\rfloor-\lfloor k/p\rfloor-\lfloor(n-k)/p\rfloor$. because there is

$$\begin{aligned} n &= \lfloor n/p\rfloor p + (n\bmod p),\\ k &= \lfloor k/p\rfloor p + (k\bmod p),\\ n-k &= \lfloor (n-k)/p\rfloor p + ((n-k)\bmod p),\\ \end{aligned}$$

So, by subtracting the last two equations from the first equation, we get

$$(\lfloor n/p\rfloor-\lfloor k/p\rfloor-\lfloor(n-k)/p\rfloor)p = (k\bmod p)+((n-k)\bmod p)-(n\bmod p).$$

On the right-hand side of the equation, the sum of the first two terms is strictly less than $2p$, and the third term $n\bmod p$ is the remainder of the sum of the first two terms, so the right side must be non-negative, but smaller than $2p$, and it needs to be a multiple of $p$, so it can only be $0$ or $p$. This means that $\lfloor n/p\rfloor-\lfloor k/p\rfloor-\lfloor(n-k)/p\rfloor$ can only be $0$ or $1$:

• If it is $0$, then $(n\bmod p) = (k\bmod p)+((n-k)\bmod p)$ is also true. Therefore, the exponent of the first factor in the above formula is $0$, which is equal to one; the second factor is $\dbinom{n\bmod p}{k\bmod p}$; and the third factor, as can be seen from the previous expansion, is equal to $\dbinom{\lfloor n/p\rfloor}{\lfloor k/p\rfloor}$. At this point, Lucas’ formula holds;

• If it is $1$, then the first factor has exponent $1$ and that factor is equal to zero, so the remainder of the binomial coefficient is zero. At the same time, $\dbinom{n\bmod p}{k\bmod p}$ on the right side of the equation to be proved by Lucas' theorem must also be zero, because there must be $(n\bmod p)<(k\bmod p)$ at this time; otherwise, there will be

  $$((n-k)\bmod p) = p + (n\bmod p) - (k\bmod p) \ge p.$$

  This obviously contradicts the definition of remainder.

By combining the two situations, the required Lucas theorem is obtained. This proof shows that when solving combinatorial numbers modulo prime numbers, the results obtained by using the Lucas theorem and the exLucas algorithm are equivalent.

The Lucas theorem points out that when the modulus is a prime number $p$, the calculation of large combinatorial numbers can be converted into the calculation of smaller combinatorial numbers. In the right formula, the first combination number can continue to recurse until $n,k [!example]+ means

long long Lucas(long long n, long long k, long long p) {
  if (k == 0) return 1;
  return (C(n % p, k % p, p) * Lucas(n / p, k / p, p)) % p;
}

Among them, C(n, k, p) is used to calculate the number of small-scale combinations.

The recursion is performed at most $O(\log_p n)$ times, so the complexity of the algorithm is $O(f(p)+g(p)\log_p n)$, where $f(p)$ is the complexity of the preprocessing combination number, and $g(p)$ is the complexity of the single calculation combination number.

Reference Implementation

After preprocessing the factorials and their inverses up to $p$ in $O(p)$ time, the following reference implementation computes each binomial coefficient in $O(1)$ time:

#include <iostream>
#include <vector>

class BinomModPrime {
  int p;
  std::vector<int> fa, ifa;

  // Calculate binom(n, k) mod p for n, k < p.
  int calc(int n, int k) {
    if (n < k) return 0;
    long long res = fa[n];
    res = (res * ifa[k]) % p;
    res = (res * ifa[n - k]) % p;
    return res;
  }

 public:
  BinomModPrime(int p) : p(p), fa(p), ifa(p) {
    // Factorials mod p till p.
    fa[0] = 1;
    for (int i = 1; i < p; ++i) {
      fa[i] = (long long)fa[i - 1] * i % p;
    }
    // Inverse of factorials mod p till p.
    ifa[p - 1] = p - 1;  // Wilson's theorem.
    for (int i = p - 1; i; --i) {
      ifa[i - 1] = (long long)ifa[i] * i % p;
    }
  }

  // Calculate binom(n, k) mod p.
  int binomial(long long n, long long k) {
    long long res = 1;
    while (n || k) {
      res = (res * calc(n % p, k % p)) % p;
      n /= p;
      k /= p;
    }
    return res;
  }
};

int main() {
  int t, p;
  std::cin >> t >> p;
  BinomModPrime bm(p);
  for (; t; --t) {
    long long n, k;
    std::cin >> n >> k;
    std::cout << bm.binomial(n, k) << '\n';
  }
  return 0;
}

The time complexity of this implementation is $O(p+T\log_p n)$, where $T$ is the number of queries.

exLucas algorithm

Lucas' theorem requires that the modulus $p$ must be a prime number. If $p$ is not a prime number, the exLucas algorithm needs to be used. Despite its name, this algorithm does not use Lucas' theorem in actual operation. Its key step is factorials modulo prime powers. The second proof above points out its connection with Lucas' theorem.

Prime-Power Modulus Case

First consider the case where the modulus is a prime power $p^\alpha$. By separating the powers of $p$ from other powers in the factorial $n!$, we can get the decomposition:

Among them, $\nu_p(n!)$ is the power of $p$ in the prime factorization of $n!$, and $(n!)_p$ is obviously coprime with $p$. Therefore, the combinatorial number can be written as:

$\nu_p(n!)$ and so on in the formula can be calculated by Legendre formula, and $(n!)_p$ and so on can be calculated by recurrence relation. Since the latter is relatively prime with $p^\alpha$, the inverse of the product in the denominator can be calculated by extended Euclidean algorithm. The problem will be solved.

Note that if the power is $\nu_p(n!)-\nu_p(k!)-\nu_p((n-k)!)\ge\alpha$, the remainder must be zero and there is no need to do any more calculations.

General Modulus Case

For the case where $m$ is a general composite number, just do prime factorization on it first:

Then, calculate the remainders of the combinations $\dbinom{n}{k}$ modulo $p_i^{\alpha_i}$, and get $s$ congruence equations:

Finally, use Chinese remainder theorem to find the remainder modulo $m$.

Reference Implementation

Finally, the reference implementation of the template title binomial coefficient is given.

#include <iostream>
#include <vector>

// Extended Euclid.
void ex_gcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1;
    y = 0;
  } else {
    ex_gcd(b, a % b, y, x);
    y -= a / b * x;
  }
}

// Inverse of a mod m.
int inverse(int a, int m) {
  int x, y;
  ex_gcd(a, m, x, y);
  return (x % m + m) % m;
}

// Coefficient in CRT.
int crt_coeff(int m_i, int m) {
  long long mm = m / m_i;
  mm *= inverse(mm, m_i);
  return mm % m;
}

// Binominal Coefficient Calculator Modulo Prime Power.
class BinomModPrimePower {
  int p, a, pa;
  std::vector<int> f;

  // Obtain multiplicity of p in n!.
  long long nu(long long n) {
    long long count = 0;
    do {
      n /= p;
      count += n;
    } while (n);
    return count;
  }

  // Calculate (n!)_p mod pa.
  long long fact_mod(long long n) {
    bool neg = p != 2 || pa <= 4;
    long long res = 1;
    while (n > 1) {
      if ((n / pa) & neg) res = pa - res;
      res = res * f[n % pa] % pa;
      n /= p;
    }
    return res;
  }

 public:
  BinomModPrimePower(int p, int a, int pa) : p(p), a(a), pa(pa), f(pa) {
    // Pretreatment.
    f[0] = 1;
    for (int i = 1; i < pa; ++i) {
      f[i] = i % p ? (long long)f[i - 1] * i % pa : f[i - 1];
    }
  }

  // Calculate Binom(n, k) mod pa.
  int binomial(long long n, long long k) {
    long long v = nu(n) - nu(n - k) - nu(k);
    if (v >= a) return 0;
    auto res = fact_mod(n - k) * fact_mod(k) % pa;
    res = fact_mod(n) * inverse(res, pa) % pa;
    for (; v; --v) res *= p;
    return res % pa;
  }
};

// Binominal Coefficient Calculator.
class BinomMod {
  int m;
  std::vector<BinomModPrimePower> bp;
  std::vector<long long> crt_m;

 public:
  BinomMod(int n) : m(n) {
    // Factorize.
    for (int p = 2; p * p <= n; ++p) {
      if (n % p == 0) {
        int a = 0, pa = 1;
        for (; n % p == 0; n /= p, ++a, pa *= p);
        bp.emplace_back(p, a, pa);
        crt_m.emplace_back(crt_coeff(pa, m));
      }
    }
    if (n > 1) {
      bp.emplace_back(n, 1, n);
      crt_m.emplace_back(crt_coeff(n, m));
    }
  }

  // Calculate Binom(n, k) mod m.
  int binomial(long long n, long long k) {
    long long res = 0;
    for (size_t i = 0; i != bp.size(); ++i) {
      res = (bp[i].binomial(n, k) * crt_m[i] + res) % m;
    }
    return res;
  }
};

int main() {
  int t, m;
  std::cin >> t >> m;
  BinomMod bm(m);
  for (; t; --t) {
    long long n, k;
    std::cin >> n >> k;
    std::cout << bm.binomial(n, k) << '\n';
  }
  return 0;
}

This algorithm decomposes the modulus $m$ into prime powers during preprocessing, and then preprocesses all $p^\alpha$ by preprocessing the product of all natural numbers from $1$ to $p^\alpha$ that are not multiples of $p$, as well as its corresponding coefficients when merging answers to the Chinese Remainder Theorem. The time complexity of preprocessing is $O(\sqrt{m}+\sum_ip_i^{\alpha_i})$. Each time you ask, the complexity is $O(\log m+\sum_i\log_{p_i}n)$. The two items in the complexity are the complexity of calculating the inverse element and calculating the power and factorial remainder.

Exercises

• Luogu3807【template】Lucas theorem
• SDOI2010 Ancient Pig Text Lucas theorem
• Luogu4720【template】extended Lucas
• Ceizenpok’s formula