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 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), nội suy Lagrange (Lagrange Interpolation).
Ý tưởng:
Như vậy để tính tập , ta cần tính tập và . Để tính tập từ tập , ta cần bài toán ngay dưới đây:
Ta cũng 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;
}
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' ?