owned this note
owned this note
Published
Linked with GitHub
**Người viết:**
- Nguyễn Minh Hiển - Trường Đại học Công nghệ, ĐHQGHN
**Reviewer:**
- Phạm Công Minh - Trường THPT chuyên Khoa học Tự nhiên, ĐHQGHN
- Phạm Hoàng Hiệp - University of Georgia
| [TOC] |
| ----- |
# Giai thừa modulo p
## **Lời nói đầu**
Khi giải các bài toán, đôi khi chúng ta sẽ gặp phải những công thức liên quan đến giai thừa, có thể ở trên tử, cũng có thể ở dưới mẫu như hệ số nhị thức. $n!$ là một hàm số nguyên tăng khá nhanh, nên thường chúng ta được yêu cầu tính toán nó theo một modulo $p$ nguyên tố nào đó.
Để giải quyết vấn đề đó trong bài viết này, trước hết chúng ta đi đến một số ký hiệu sau:
- $n!$ ($n$ giai thừa) là tích các số từ $1$ đến $n$:
$$n! = 1 \cdot 2 \cdot 3 \cdots n$$
- $v_p(n)$ là số nguyên lớn nhất để $p^{v_p(n)}$ là ước của $n$.
<!-- $n!_{\%p}$ bằng $n!$ sau khi chia tất cả thừa số $p$:
$$n!_{\%p} = \dfrac{n!}{v_p(n!)}$$
- $n!_p$ là tích các số từ $1$ đến $n$, không tính các số chia hết cho $p$: $\require{cancel}$
$$n!_p = 1 \cdot 2 \cdots (p-1) \cancel{p} (p+1) \cdots (2p-1) \cancel{2p} (2p+1) \cdots n$$
- $$n!_{\%p}=\left(n\right)!_p \left(\left\lfloor\frac{n}{p}\right\rfloor\right)!_p \left(\left\lfloor\frac{n}{p^2}\right\rfloor\right)!_p\ldots$$
Như vậy, khi tính được $n!_p, \ldots$ ta hoàn toàn tính được $n!_{\%p}$.-->
## **Thuật toán "ngây thơ"**
Cách đơn giản nhất để tính $n! \bmod p$ là ta lần lượt nhân các số từ $1$ đến $n$:
```cpp=
factorial[0] = 1;
for (int i = 1; i <= n; i++) {
factorial[i] = factorial[i - 1] * i % p;
}
```
Chú ý rằng nếu $n \ge p$ thì $n!$ luôn chia hết cho $p$.
Độ phức tạp: $O(n)$.
<!--Tuy nhiên, nếu $n \ge p$ thì $n!$ luôn chia hết cho $p$. Điều này sẽ gây bất lợi khi $n!$ ở dưới mẫu. Vì vậy, chúng ta sẽ tính $n!_p$ thay thế.
Ở đây, xin giới thiệu với bạn đọc **định lý Wilson**:
$$(p - 1)! \equiv -1 \pmod p$$
Áp dụng, ta có:
$$n!_p \equiv (-1)^{\left\lfloor\frac{n}{p}\right\rfloor} \cdot (n \bmod p)! \pmod p$$
Nếu chuẩn bị trước $0!,~ 1!,~ 2!,~ \dots,~ (p-1)!$ ta có thể tính $n!_p$ trong thời gian $O(1)$ và $n!_{\%p}$ trong thời gian $O(\log_pn)$.
*Dưới đây là code C++ minh họa việc tính $n!_{\%p} \bmod p$*
```cpp=
int factmod(int n, int p) {
vector<int> f(p);
f[0] = 1;
for (int i = 1; i < p; i++)
f[i] = f[i-1] * i % p;
int res = 1;
while (n > 1) {
if ((n/p) % 2)
res = p - res;
res = res * f[n%p] % p;
n /= p;
}
return res;
}
```
-->
## **Công thức Legendre**
Phát biểu [công thức Legendre](https://vi.wikipedia.org/wiki/C%C3%B4ng_th%E1%BB%A9c_Legendre#):
$$v_p(n!) = \left\lfloor\frac{n}{p}\right\rfloor + \left\lfloor\frac{n}{p^2}\right\rfloor + \left\lfloor\frac{n}{p^3}\right\rfloor + \ldots$$
>Chú ý rằng khi $p^k > n$ thì $\left\lfloor\frac{n}{p^k}\right\rfloor = 0$, ta không cần tính tiếp
:::spoiler Chứng minh công thức Legendre
$v_p(n!)$ là số mũ của $p$ khi phân tích $n!$ ra thừa số nguyên tố, vì thế ta cần tính xem từ $1$ đến $n$ cho ta bao nhiêu thừa số $p$.
- Trong các số từ $1$ đến $n$, có $\left\lfloor\frac{n}{p}\right\rfloor$ số chia hết cho $p$. Mỗi số cho ta $1$ thừa số $p$.
- Trong các số chia hết cho $p$, lại có $\left\lfloor\frac{n}{p^2}\right\rfloor$ số chia hết cho $p^2$. Mỗi số cho ta thêm $1$ thừa số $p$ nữa.
- Trong các số chia hết cho $p^2$, lại có $\left\lfloor\frac{n}{p^3}\right\rfloor$ số chia hết cho $p^3$. Mỗi số cho ta thêm $1$ thừa số $p$ nữa.
- ... Cứ như vậy cho đến khi $\left\lfloor\frac{n}{p^k}\right\rfloor = 0$ thì ta kết thúc.
:::
*Code minh họa công thức Legendre:*
```cpp=
int vp(int n, int p) {
int cnt = 0;
while (n) {
n /= p;
cnt += n;
}
return cnt;
}
```
Độ phức tạp: $O(\log_p n)$
Bạn đọc có thể sử dụng công thức Legendre cùng sàng nguyên tố (đã chuẩn bị trước) để phân tích $n! = p_1^{k_1} \cdot p_2^{k_2} \cdots$ rồi tính $n!$ theo modulo bất kỳ với thuật toán lũy thừa nhanh.
Dưới đây là code được tham khảo từ [Geeksforgeeks](https://www.geeksforgeeks.org/compute-n-under-modulo-p/):
```cpp=
#include<bits/stdc++.h>
using namespace std;
const int N = 1000000; // 10^6
vector<int> primes;
// Công thức Legendre
int largestPower(int n, int p){
int x = 0;
while (n) {
n /= p;
x += n;
}
return x;
}
int binpow(int x, int y, int mod){
int res = 1;
x = x % mod;
while (y) {
if (y & 1)
res = (1LL * res * x) % mod;
y >>= 1;
x = (1LL* x * x) % mod;
}
return res;
}
// Sàng nguyên tố
void prepare(){
bool isPrime[N + 1];
memset(isPrime, 1, sizeof(isPrime));
isPrime[0] = isPrime[1] = 0;
for (int i = 2; i * i <= N; i++) {
if (isPrime[i]) {
for (int j = i * i; j <= N; j += i) {
isPrime[j] = 0;
}
}
}
for (int i = 2; i <= N; i++) {
if (isPrime[i]) {
primes.push_back(i);
}
}
}
int modFact(int n, int mod){
int res = 1;
for (int p : primes) {
if (p <= n) {
int k = largestPower(n, p);
res = (1LL * res * binpow(p, k, mod)) % mod;
}
else {
break;
}
}
return res;
}
int main(){
prepare();
int n, mod;
cin >> n >> mod;
cout << modFact(n, mod);
return 0;
}
```
Độ phức tạp thời gian:
- Chuẩn bị: $O(N \log\log N)$
- Truy vấn: $O(N)$
## **Chia căn**
- *Điều kiện sử dụng:* $n < p \le 2 \cdot 10^9$ ($p$ không cần phải nguyên tố).
Có thể chạy đến $5 \cdot 10^{10}$ tùy thuộc vào lượng hằng số có thể sử dụng.
- *Ý tưởng:*
- Chia đoạn $[1;p]$ thành các đoạn riêng biệt có kích thước $S$, ta sẽ có tổng cộng $\text{CNT} = \left\lfloor \frac{p}{S} \right\rfloor + 1$ đoạn.
- Ta sẽ tính sẵn một mảng gồm $\text{CNT}$ giá trị sau: $0!,\, S!,\, (2S)!,\, \ldots,\, \left(\left\lfloor \frac{p}{S} \right\rfloor \cdot S\right)!$ với độ phức tạp $O(p)$. Điều này không có ý nghĩa lắm. Vì thế, ta chỉ cần chạy code ra chương trình khác rồi lấy kết quả vào bài, bước này sẽ chỉ còn độ phức tạp $O\left(\frac{p}{S}\right)$ nhưng lại tăng độ dài của code!
- Khi này, chọn $k = \left\lfloor \frac{n}{S} \right\rfloor$, ta có: $n! = (k\cdot S)! \times \underbrace{(k\cdot S + 1) \ldots n}_{\style{font-family: Cambria Math}{\text{Duyệt mất tối đa }} 10^6 \style{font-family: Cambria Math}{\text{ phép tính}}}$
- *Code C++ minh họa:*
- Phần chuẩn bị:
```cpp=
const int MOD = 1e9 + 7;
const int S = 1e7;
// Phần này thực hiện ở code khác
// và lấy kết quả vào mảng fact[] trong code chính.
void prepare(){
int tmp = 1;
for (int i = 1; i < MOD; i++){
tmp = 1LL * tmp * i % MOD;
if (i % S == 0) {
cout << tmp << ", ";
}
}
}
```
Kết quả in ra màn hình là một dãy: $1, 641102369, 578095319, 5832229, \ldots$
Độ phức tạp phần chuẩn bị này là :$O(p)$.
- Phần code chính:
```cpp=
const int MOD = 1e9 + 7;
const int S = 1e7;
// Dãy được copy từ kết quả in ra của code bên trên
int fact[] = {1, 641102369, 578095319, 5832229, ...};
int get_fact(int n){
int t = n / S;
int res = fact[t];
for (int i = t * S + 1; i <= n; i++)
res = 1LL * res * i % MOD;
return res;
}
```
- *Đánh giá*:
- Độ phức tạp không gian: $O\left(\dfrac{p}{S}\right)$
- Độ phức tạp thời gian:
- Tiền xử lý (để lưu mảng `fact[]`): $O\left(\dfrac{p}{S}\right)$
- Truy vấn $O\left(S\right)$
Tùy thuộc vào giới hạn bài toán, giới hạn độ dài của code, ta sẽ chọn $S$ sao cho phù hợp.
## **Biến đổi Fourier nhanh (FFT)**
>Trong bài này, ta sẽ sử dụng dạng Biến đổi số học - Number-theoretic transform (NTT) của FFT. Bạn đọc tham khảo tại [đây](https://en.wikipedia.org/wiki/Discrete_Fourier_transform_over_a_ring#Number-theoretic_transform).
Đây là thuật toán nhân hai đa thức bậc không vượt quá $n$ với độ phức tạp $O(n \log n)$.
- *Điều kiện sử dụng:* $n < p \le 10^{12}$
- *Ý tưởng:* Ta sẽ tính $n! \bmod p$ như sau:
Xét $m = \left\lfloor \sqrt n \right\rfloor$ và đa thức $P(x) = (x+1)(x+2)\ldots (x + m) = \prod\limits_{i=1}^m (x+i)$
Khi này:
$$n! = \underbrace{P(0) \cdot P(m) \cdot P(2m) \ldots P\Big( (m-1)m \Big)}_{\style{font-family: Cambria Math}{\large\text{chính là }} 1\cdot 2\cdot 3 \ldots m ^ 2} \cdot \underbrace{(m^2 + 1)\ldots n}_{\style{font-family: Cambria Math}{\large\text{tính được trong }} O\left(\sqrt n\right)} \pmod p$$
### Cách 1: Sử dụng FFT đa điểm (FFT multipoint evaluation)
> Cách này sử dụng rất nhiều kiến thức khó và mang giá trị về nghiên cứu là chính. Vì thế, để tham khảo chi tiết, bạn đọc tham khảo tại [đây](https://www.cecm.sfu.ca/CAG/theses/justine.pdf). Tuy nhiên bạn đọc có thể tham khảo các bài toán nhỏ của nó.
- *Kiến thức sử dụng:* Biến đổi Fourier nhanh (FFT), chia để trị, chia đa thức.
- Đầu tiên, ta sử dụng **FFT/NTT** để khai triển đa thức $P(x)$ bên trên. Bạn đọc tham khảo bài tương tự của bước này: [960G - Bandit Blues - Codeforces](https://codeforces.com/contest/960/problem/G)
- Sau đó, ta tiếp tục sử dụng **FFT/NTT** để tính giá trị $P(x) \bmod M$ tại $m$ điểm: $k\cdot m$ với $k=0..(m-1)$. Đó cũng chính là bài toán [POLYEVAL - Codechef](https://discuss.codechef.com/t/polyeval-editorial/12771) hoặc [Library Checker - Multipoint Evaluation](https://judge.yosupo.jp/problem/multipoint_evaluation).
- Chú ý rằng [POLYEVAL - Codechef](https://discuss.codechef.com/t/polyeval-editorial/12771) là bài toán dễ hơn, và bạn có thể sử dụng chia để trị với không gian lưu trữ $O(p)$. Để làm được bài [Library Checker - Multipoint Evaluation](https://judge.yosupo.jp/problem/multipoint_evaluation), ta sử dụng thuật toán chia để trị và chia đa thức để đạt được không gian lưu trữ $O\left(\sqrt p\right)$:
- Ta thấy, $P(x) = Q(x) \cdot (x-k) + P(k) \implies P(k) = P(x) \mod (x - k)$.
- Giả sử cần tính $P(x)$ tại $x_1, x_2, \ldots, x_n$. Mỗi lần, ta chia $P(x)$ thành $2$ đa thức (*Chú ý rằng, phép chia đa thức sẽ được chuyển thành nhân với đa thức nghịch đảo để sử dụng FFT*):
$$\begin{align*}
P_0(x) &= P(x) \mod \Big((x-x_1)(x-x_2)\cdots(x-x_{\lfloor n/2 \rfloor})\Big)\\
P_1(x) &= P(x) \mod \Big((x-x_{n/2+1})\cdots(x-x_n)\Big)
\end{align*}$$$$\implies\begin{align*}
P(x) &= \Big((x-x_1)(x-x_2)\cdots(x-x_{\lfloor n/2 \rfloor})\Big) \cdot Q_0(x) + P_0(x)\\
P(x) &= \Big((x-x_{\lfloor n/2 \rfloor+1})\cdots(x-x_n)\Big) \cdot Q_1(x) + P_1(x)
\end{align*}$$
$\implies P(x_i) = \begin{cases}
P_0(x_i) &\style{font-family: Cambria Math}{\large\text{nếu }}i\le \lfloor n/2 \rfloor\\
P_1(x_i) &\style{font-family: Cambria Math}{\large\text{còn lại}}\end{cases}$
- Bài toán trở thành: Tính $P_0(x)$ tại $x_1, \cdots, x_{\lfloor n/2 \rfloor}$ và tính $P_1(x)$ tại $x_{\lfloor n/2 \rfloor+1}, \cdots, x_n$.
- *Đánh giá:*
- Độ phức tạp không gian: $O\left(\sqrt p\right)$
- Độ phức tạp thời gian: $O\left(\sqrt p \log^2 p\right)$
### Cách 2: Kết hợp sử dụng nội suy Lagrange
- *Kiến thức sử dụng:* Biến đổi Fourier nhanh (FFT), nội suy Lagrange (Lagrange Interpolation).
- *Ý tưởng*:
- Đặt $P_d(x) = \prod\limits_{i=1}^d (x+i)$ là đa thức bậc $d$.
- Đặt tập $G_d(x) = \bigg\{ P_d(x), P_d(m +x), P_d(2m+x), \ldots, P_d\Big( (m-1)m+x \Big) \bigg\}$
Mục tiêu của ta là tính tập $G_m(0)$.
- Nhận thấy: $P_1(0) = 1$ và
$$P_{d}(x) = \begin{align*}
\begin{cases}
P_{d/2}(x) \cdot P_{d/2}(x+d/2) &\text{nếu } d \text{ chẵn}\\
P_{\lfloor d/2 \rfloor}(x) \cdot P_{\lfloor d/2 \rfloor}\big(x+\lfloor d/2 \rfloor\big)\cdot (x + d + 1) &\text{nếu } d \text{ lẻ}
\end{cases}
\end{align*}$$
- Như vậy để tính tập $G_{d}(0)$, ta cần tính tập $G_{d/2}(0)$ và $G_{d/2}(d/2)$. Để tính tập $G_{d/2}(d/2)$ từ tập $G_{d/2}(0)$, ta cần bài toán ngay dưới đây:
- Cho đa thức $h(x)$ có bậc là $d$ và biết $h(0), h(1), \ldots, h(d)$. Ta cần tính $h(x), h(x+1), \dots, h(x+d')$ với nội suy Lagrange. $\forall k = 1,2,\ldots d'$, ta có:
$$ \begin{align*} h(x + k) &= \sum_{i = 0}^d h(i) \prod_{j=0, i \neq j}^d \frac{x+k-j}{i-j} \\\ &= \left(\prod_{j=0}^d (x + k - j) \right) \underbrace{\left( \sum_{i=0}^d \underbrace{\frac{h(i)}{i! (d-i)! (-1)^{d-i}}}_{\text{Hệ số bậc } i} \cdot \underbrace{\frac{1}{x + k - i}}_{\text{Hệ số bậc } x+k-i} \right)}_{\text{Sử dụng FFT/NTT ở đây, hệ số bậc } x+k} \end{align*} $$
- Áp dụng với $h(i) = P_d(i \cdot m)$ và $x = d\cdot m^{-1}$.
Trước đó, ta đã tính tập
$$G_d(0) = \Big\{h(0), h(1), \ldots, h(d)\Big\}$$
Ta cũng sẽ tính tiếp được
$$h(x + i) = P_d((x + i) \cdot m) = P_d(i\cdot m + d)$$
Và đây cũng chính là tập $G_d(d)$ cần tính!
- *Code C++ minh họa (lược bỏ phần FFT)*
```cpp=
int MOD;
// hàm trả về nghịch đảo x modulo MOD
long long inv(long long x);
// hàm trả về đa thức a * b
vector<long long> NTT(vector<long long> a, vector<long long> b);
void mul(long long &x, long long y){
x = __int128(x) * y % MOD;
}
// Biết h(0), h(1), ..., h(d)
// Hàm trả về h(m), h(m + 1), ..., h(m + cnt - 1)
// m > d
vector<long long> Lagrange(vector<long long> h, long long m, int cnt){
int d = h.size() - 1;
// tính h[i] = (-1)^(d-i) h(i)/(i! (d-i)!)
// hệ số x^i
for (int i = 0; i <= d; i++){
mul(h[i], (ifact[i] * ifact[d - i]) % MOD);
if ((d - i) & 1)
h[i] = (MOD - h[i]) % MOD;
}
//tính f[i] = 1/(m+i-d)
// hệ số x^(m + cnt - 1 - i)
vector<long long> f(d + cnt);
long long now = m - d;
for (int i = 0; i < d + cnt; i++)
f[i] = inv(now+i);
// Nhân 2 đa thức và hệ số được lấy mod p
h = NTT(f, h);
h.resize(d + cnt);
// chỉ lấy hệ số x^(m) đến x^(m+cnt-1)
h = vector<long long>(h.begin() + d, h.end());
now = 1;
for (int i = 0; i <= d; i++)
mul(now, m - i);
mul(h[0], now);
for (int i = 1; i < cnt; i++){
mul(now, m + i);
mul(now, inv(m + i - d - 1));
mul(h[i], now);
}
return h;
}
long long factorial(long long n, long long p){
if (n >= p) return 0;
if (n < 2) return 1;
int s = __builtin_sqrtl(n);
MOD = p;
vector<long long> h{1, s + 1};
for (int bit = __lg(s) - 1, d = 1; bit >= 0; bit--){
// Hiện tại h(i) = (i * s + 1) * (i * s + 1) ... (i * s + d)
// Tính h(d+1), ..., h(2d)
auto nh1 = Lagrange(h, d + 1, d);
// Tính h(d.inv(s)), ..., h(d.inv(s) + 2d)
// Như vậy, nh2(i) = (i * s + d + 1) * (i * s + d + 2) ... (i * s + d * 2)
auto nh2 = Lagrange(h, 1LL * inv(s) * d % mod, 2 * d + 1);
// h giờ đây là h(0), h(1), ..., h(2d)
h.insert(h.end(), nh1.begin(), nh1.end());
// d --> d * 2
d <<= 1;
// Hiện tại h(i) = (i * s + 1) * (i * s + 2) ... (i * s + d/2)
// Còn nh2(i) = (i * s + d/2 + 1) * (i * s + d/2 + 2) ... (i * s + d)
for (int i = 0 ; i <= d; i++)
h[i] *= nh2[i];
// Tại đây h(i) = (i * s + 1) * (i * s + 1) ... (i * s + d)
// Nếu bit hiện tại của s là 1
if (s >> bit & 1){
d |= 1;
long long tmp = d;
// h[i] *= (d + i * s)
for (int i = 0; i < d; i++, tmp += s)
mul(h[i], tmp);
long long last = 1, tj = 1LL * s * d;
// last = (d*s+1)(d*s+2)...(d*s+d)
for (int i = 1; i <= d; i++)
tj++, last *= tj;
// Thêm biến last vào h
h.emplace_back(last);
}
}
long long ans = 1;
for (int i = 0; i < s; i++)
mul(ans, h[i]);
for (long long i = 1LL * s * s + 1; i <= n; i++)
mul(ans, i);
return ans;
}
```
- *Đánh giá:*
- Độ phức tạp không gian: $O(\sqrt n)$
- Độ phức tạp thời gian:
- Hàm Lagrange ở phần sử dụng sử dụng FFT có độ phức tạp $O(N \log N)$ với $N = \deg(h) + \text{cnt}$.
Các vòng lặp khác chỉ có độ phức tạp tối đa đến $O(N)$.
- Hàm factorial có độ phức tạp lớn nhất ở phần sử dụng hàm Lagrange với tổng cộng: $O\left(s\log s + \frac{s}{2}\log\frac{s}{2} + \frac{s}{4}\log\frac{s}{4} + \ldots \right)$ với $s = \lfloor \sqrt n \rfloor$. Chú ý rằng:
\begin{align}
s\log s + \frac{s}{2}\log\frac{s}{2} + \frac{s}{4}\log\frac{s}{4} + \ldots &< s\log s + \frac{s}{2}\log s + \frac{s}{4}\log s + \ldots\\
&= \log s \cdot \left(s+ \frac{s}{2} + \frac{s}{4} + \ldots\right)\\
&< 2s\log s
\end{align}
- Điều này có nghĩa là bài toán có tổng độ phức tạp: $O\left(\sqrt n \log n\right)$
## Bài tập luyện tập
[Hackerearth - Army Parade](https://www.hackerearth.com/fr/problem/algorithm/army-parade-7bcfea8e/)
[SPOJ - FACTMODP](https://www.spoj.com/problems/FACTMODP/)
## Tài liệu tham khảo
[GeeksforGeeks - Compute n! under modulo p](https://www.geeksforgeeks.org/compute-n-under-modulo-p/)
[FFT - multipoint evaluation](https://cs.stackexchange.com/questions/60239/multi-point-evaluations-of-a-polynomial-mod-p)
[Factorial mod prime - Prabowo Djonatan](https://blog.prabowodjonatan.id/posts/factorial-mod-p/)
[Codeforces - Fastest way to get factorial modulo a prime](https://codeforces.com/blog/entry/63491?#comment-474051)
[Codeforces - How to calculate n!%p for very large 'n' ?](https://codeforces.com/blog/entry/69002)