Nguyễn Minh Hiển - Trường Đại học Công nghệ, ĐHQGHN
Reviewer:
Phạm Công Minh - Trường Đại học Công nghệ, ĐHQGHN
Đôi khi, chúng ta sẽ gặp những bài tập như tính xmodp hay thậm chí như tính số Fibonacci Fnmodp. Mà chúng ta biết, công thức tổng quát:
Fn=51[(21+5)n−(21−5)n]
Việc xuất hiện 5 đặt ra nhiều thách thức cho việc tính toán nhanh Fn, nhưng đồng thời cũng mở ra những phương pháp mới để chinh phục được bài toán Fnmodp
Ta sử dụng tiêu chuẩn Euler (Euler's criterion) như sau. Với p nguyên tố lẻ:
(pa)≡a2p−1(modp)
Đến đây, ta sử dụng lũy thừa nhanh để tính.
int pow_mod(long long a, long long n, long long p); // hàm tính lũy thừa nhanh modulo p
int legendre_symbol(int a, int p) {
return pow_mod(a, (p - 1) >> 1, p);
}
Trước hết, ta cần tìm thặng dư "không chính phương" để thực hiện hai thuật toán bên dưới.
Vì một nửa số phần tử trong tập {1,2,⋯,p−1} là thặng dư không chính phương, nên ta sẽ duyệt từng số từ 1 cho đến khi gặp được số thỏa mãn. Để kiểm tra một số thỏa mãn hay không, ta sử dụng cách Tiêu chuẩn Euler bên trên.
Để thuật toán hiệu quả hơn, bạn nên sinh số ngẫu nhiên và kiểm tra đến khi tìm được. Xác suất 1 lần thử tìm được là 21, nên xác suất sau 32 lần thử mà bạn chưa tìm ra là 2321.
Bài viết xin không đề cập phần chứng minh thuật toán. Bạn đọc tham khảo tại Wikipedia.
Bước 1: ta phân tích p=Q⋅2S+1 với Q lẻ
Bước 2: Chọn z là một thặng dư không chính phương bất kỳ.
Bước 3: Gán
xb←a2Q+1←aQ
Bước 4: Lặp
Tìm m nhỏ nhất (0≤m<r) sao cho b2m≡1(modp)
Nếu m=0⟺b≡1(modp) thì x chính là đáp án cần tìm.
Nếu m>0 thì đặt e=2m+1p−1=Q⋅2S−m−1 gán:
xb←x⋅ze←b⋅z2e
Code C++ minh họa:
int pow_mod(long long a, long long k, long long M) {
long long ans = 1;
for (; k > 0; a = a * a % M, k >>= 1) {
if (k & 1) {
ans = ans * a % M;
}
}
return ans;
}
int Tonelli_Shanks(int a, int p) {
if (p == 2) {
return (a & 1);
}
int S = 0, Q = p - 1;
while (Q % 2 == 0) {
S++;
Q /= 2;
}
int z = 2;
while (legendre_symbol(z, p) != p - 1) {
z++;
}
int x = pow_mod(a, (Q + 1) >> 1, p), b = pow_mod(a, Q, p);
int m, v, e, u;
while (b % p != 1) {
m = 0, v = 1; // v = 2^m
while (pow_mod(b, v, p) != 1) {
m++;
v <<= 1;
}
e = Q << (S - m - 1);
u = pow_mod(z, e, p);
x = (1LL * x * u) % p;
b = (((1LL * u * u) % p) * b) % p;
}
return x;
}
Bài viết xin không đề cập phần chứng minh thuật toán. Bạn đọc tham khảo tại Wikipedia.
Bước 1: Tìm b sao cho b2−a là thặng dư không chính phương modulo p
Bước 2: Ta tính x+yb2−a=(b+b2−a)(2p+1).
Khi đó, xmodp tìm được chính là nghiệm của bài toán.
Nói cách khác là ⟨x,y⟩=⟨b,1⟩(2p+1) trên Fp(b2−a)
Code C++ minh họa
Về cài đặt, như đã nói ở trên, ⟨x,y⟩ khá giống số phức nên việc cài đặt cũng tương tự như vậy.
int a, p;
int k; // thặng dư không chính phương mod p
struct Complex {
int re, im;
Complex(int a = 0, int b = 0) {
re = a;
im = b;
}
Complex operator*(const Complex &o) {
Complex res;
res.re = (1LL * re * o.re + (1LL * im * k % p) * o.im) % p;
res.im = (1LL * re * o.im + 1LL * im * o.re) % p;
return res;
}
Complex pow(long long k) {
Complex res = Complex(1, 0), A = *this;
while (k) {
if (k & 1)
res = res * A;
A = A * A;
k >>= 1;
}
return res;
}
};
int Cipolla(long long a, long long p) {
if (p == 2) {
return (a & 1);
}
::p = p; // struct Complex sử dụng biến toàn cục p và k
// Tìm k = b^2 - a, sao cho k không chính phương
int b = 2;
while (true) {
b %= p;
k = (b * b - a) % p;
if (k < 0)
k += p;
if (legendre_symbol(k, p) == p - 1)
break;
b++;
}
// Ta cần tìm <b, 1>^((p+1)/2)
return Complex(b, 1).pow((p + 1) >> 1).re;
}
Ngoài các phương pháp như Nhân ma trận hay Khử nhân ma trận, còn có một phương pháp khác sử dụng
Công thức tổng quát của Fibonacci:
Fn=51[(21+5)n−(21−5)n]
Xét modulo p nguyên tố.
Nếu 5 là thặng dư bình phương modulo p
Ví dụ: Bài Codeforces - DZY Loves Fibonacci Numbers với p=109+9.
Ta tính được: 5=383008016(modp)
Sử dụng nghịch đảo modulo, ta có:
So với việc tính lũy thừa của ma trận, tính lũy thừa của 2 số vẫn nhanh hơn rất nhiều.
Nếu 5 không là thặng dư bình phương modulo p
Ví dụ: Bài VNOI - Fibonacci với p=109+7
Ở bài này, ta sử dụng trường hữu hạn như ở trên.
Ta sẽ viết (21+5)n=⟨u1,v1⟩ và (21−5)n=⟨u2,v2⟩
Trên thực tế, vì Fn nguyên nên u1−u2=0. Từ đó suy ra Fn≡v1−v2(modp).
Do sử dụng công thức tổng quát, cách này có một ưu điểm mà không cách nào có được, thể hiện qua bài toán bên dưới đây.