Người viết:
Reviewer:
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. 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 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:
Cách đơn giản nhất để tính là ta lần lượt nhân các số từ đến :
factorial[0] = 1;
for (int i = 1; i <= n; i++) {
factorial[i] = factorial[i - 1] * i % p;
}
Chú ý rằng nếu thì luôn chia hết cho .
Độ phức tạp: .
Phát biểu công thức Legendre:
Chú ý rằng khi thì , ta không cần tính tiếp
là số mũ của khi phân tích ra thừa số nguyên tố, vì thế ta cần tính xem từ đến cho ta bao nhiêu thừa số .
Code minh họa công thức Legendre:
int vp(int n, int p) {
int cnt = 0;
while (n) {
n /= p;
cnt += n;
}
return cnt;
}
Độ phức tạp:
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 rồi tính 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:
#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:
Điều kiện sử dụng: ( không cần phải nguyên tố).
Có thể chạy đến tùy thuộc vào lượng hằng số có thể sử dụng.
Ý tưởng:
Code C++ minh họa:
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:
Độ phức tạp phần chuẩn bị này là :.
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á:
fact[]): 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 sao cho phù hợp.
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.
Đây là thuật toán nhân hai đa thức bậc không vượt quá với độ phức tạp .
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. Tuy nhiên, các bài toán nhỏ mà nó sử dụng thật sự đáng để tham khảo và luyện tập.
Kiến thức sử dụng: Biến đổi Fourier nhanh (FFT), nội suy Lagrange (Lagrange Interpolation).
Ý tưởng:
Chi tiết thuật toán:
Đặt đa thức giai thừa là một đa thức bậc có dạng:
Xét tập hợp các đa thức:
Mục tiêu của ta là tính hay toàn bộ giá trị của tập .
Tương tự thuật toán lũy thừa nhanh, có thể xây dựng một cách đệ quy như sau: và
Như vậy để tính tập , ta cần tính tập và . Và cứ đệ quy như vậy đến khi dừng.
Tuy nhiên, điểm đặc biệt của thuật toán này là để tính tập , ta sẽ không đệ quy xuống để tính mà sẽ tận dụng tập đã tính toán:
Bổ đề: Cho đa thức có bậc là và biết . Ta cần tính . Với công thức nội suy Lagrange, , ta có:
Ta sẽ tính tiếp được:
Và đây cũng chính là tập cần tính!
Code C++ minh họa (lược bỏ phần FFT)
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;
}
Bạn đọc có thể tham khảo code đầy đủ của yosupo hoặc của tác giả.
Hackerearth - Army Parade
SPOJ - FACTMODP
GeeksforGeeks - Compute n! under modulo p
FFT - multipoint evaluation
Factorial mod prime - Prabowo Djonatan
Codeforces - Fastest way to get factorial modulo a prime
Codeforces - How to calculate n!%p for very large 'n' ?