Các giải thuật tính nghịch đảo modulo
White

Nguyễn Gia viết ngày 20/02/2017

Intro

Trong các coding contest, khi mà kết quả là một số rất lớn (n! hoặc tổ hợp), không dùng kiểu dữ liệu primitive để lưu được thì đề bài thường yêu cầu sẽ trả về số dư của kết quả khi chia cho một số nguyên tố lớn (thường là 10^9 + 7)

Vì kết quả là một số rất lớn nên ta không chơi được kiểu tính ra kết quả rồi lấy số dư.
(Một số ngôn ngữ như Java hay Python có tính toán cho số lớn, nhưng độ phức tạp tính toán sẽ tăng cao dẫn đến bị time limit error).
Do vậy ta sẽ phải break kết quả ra từng bước, tính số dư cho từng bước đấy và tổng hợp lại thành kết quả.
Khi tổng hợp kết quả thì ta sẽ sử dụng các công thức tính toán số dư. Với phép cộng, trừ, nhân thì đơn giản như sau:

1. ( a + b ) % c = ( ( a % c ) + ( b % c ) ) % c
2. ( a – b ) % c = ( ( a % c ) – ( b % c ) ) % c
3. ( a * b ) % c = ( ( a % c ) * ( b % c ) ) % c

Vậy với phép chia thì sao?
Phép chia không có tính chất phân phối với phép đồng dư như các phép tính trên. Để tính số dư của một thương cho một số, ta phải dùng nghịch đảo modulo.

Nghịch đảo modulo

Nghịch đảo của một số nguyên là số mà khi nhân với nó sẽ có tích là 1. Để tính thương của 2 số, ta nhân một số với nghịch đảo của số kia.
$ \frac{a}{b} = a * b ^ {-1} $
Tương tự như vậy, nghịch đảo modulo của một số theo modulo M là số mà khi nhân với nó thì được tích chia M dư 1.

$ \text{b là nghịch đảo modulo M của a} \iff (a * b) \mod M = 1 $

hay viết cách khác:
$ \text{b là nghịch đảo modulo M của a} \iff a*b \equiv 1 \mod M $

Ví dụ: 7 là nghịch đảo modulo 11 của 8 vì 7*8 = 56 và 56 % 11 = 1

Và phép chia modulo sẽ tương ứng với phép nhân nghịch đảo modulo.

Các giải thuật tính nghịch đảo modulo

Brute force

Ta chỉ cần quan tâm đến các số nhỏ hơn M.
Các số chia hết M thì không tồn tại nghịch đảo modulo, vì nó nhân với số nào cũng chia hết M.
Các số lớn hơn M thì nghịch đảo modulo của nó (nếu có) sẽ bằng nghịch đảo modulo của số dư của nó khi chia cho M.

int modulo_inverse(int n, int m) {
  for (int i=1; i < m; i++) {
    if (((long long) n * i) % m == 1) {
      return i;
    }
  }
  return -1; // not exist
}

Định lí nhỏ fermat

Định lí nhỏ fermat nói rằng:
Nếu M là số nguyên tố và a không chia hết cho M thì
$ a^{M-1} \equiv 1 \mod M$
hay
$ a * a^{M-2} \equiv 1 \mod M$

Điều đó có nghĩa là nếu M là số nguyên tố và a không chia hết cho M thì $a^{M-2} \text{ là nghịch đảo modulo M của a}$
cũng tương đương với $(a^{M-2} \mod M) \text{ là nghịch đảo modulo M của a}$

Code:

int pow_modulo(long long base, long long exponent, long long m) {
  if (exponent == 0) return 1;
  if (exponent == 1) return base;
  long long half = pow_modulo(base, exponent/2, m);
  if (exponent % 2 == 0) {
    return (half * half) % m;
  } else {
    return (((half * half) % m) * base) % m;
  }
}

int modulo_inverse_fermat(int n, int m) {
  return pow_modulo(n, m-2, m);
}

Euclidean (Ơ-cờ-lít) mở rộng

Giải thuật euclidean

Giải thuật euclidean dùng để tính ước chung lớn nhất của 2 số (Greatest common divisor, viết tắt là gcd).
Xuất phát từ nhận xét:
$gcd(a, b) = gcd(b, r) \text{ với } r = a \mod b $
Chứng minh có thể xem ở đây: http://joequery.me/notes/proof-a-bqr-gcd/
Dựa theo công thức trên ta cứ tính cho đến khi b chia hết cho r thì kết quả là r.

int gcd(int a, int b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

Giải thuật euclidean mở rộng

Giải thuật này dùng để tính x, y thoả mãn phương trình nghiệm nguyên:
$ a*x + b * y = gcd(a, b) $

Giả sử x1, y1 là nghiệm của phương trình
$ b*x1 + (a\mod b) * y1 = gcd(a, b) $

Đây chính là phương trình ở bước tiếp theo của giải thuật euclidean.
(b đóng vai trò a, và a%b đóng vai trò b trong hàm gcd(a, b))
Ta sẽ thay a % b = a - (a / b) * b (phép chia ở đây là phép chia trong programming, tức là phép chia nguyên) để có được phương trình:
$ b*x1 + (a - (a / b) * b)* y1 = gcd(a, b) $

Khai triển ta được:
$ a * y1 + b * ( x1 - (a / b) *y1) = gcd(a, b) $

Như vậy ta được công thức đệ quy:

x = y1;
y = x1 - (a/b) *  y1;

với bước cơ bản (ứng với phương trình cuối cùng là a_last * x = gcd(a, b) = a_last

x = 1;
y = 0;

Code:

int gcd_extend(int a, int b, int *x, int *y) {
  if (b == 0) {
    *x = 1;
    *y = 0;
    return a;
  }
  int x1, y1;
  int gcd = gcd_extend(b, a%b,  &x1, &y1);
  *x = y1;
  *y = x1 - (a / b) * y1;
  return gcd;
}

Dùng euclidean mở rộng để tính nghịch đảo modulo

Nếu gcd (n, M) = 1 thì sử dụng euclidean mở rộng, ta sẽ tính được x, y sao cho:
$ n*x + M * y = 1$

Mà M* y luôn chia hết cho M, nên ta có
$ n*x \equiv 1 (\mod M)$
hay x chính là nghịch đảo modulo của n.

int modulo_inverse_euclidean(int n, int m) {
  int x, y;
  if (gcd_extend(n, m, &x, &y) != 1) return -1; // not exist
  return (x % m  + m) % m; // vì x có thể âm 
}

Bonus code tính số dư của tổ hợp chập k của n khi chia cho số nguyên tố M (Code này là code trâu bò nhưng vẫn đủ nhanh và còn optimize được.)

int modulo_combination(int n, int k, int m) {
  int a = 1, b = 1, c=1;
  for (int i=2; i <= n; i++) {
    a = ((long long)a * i) % m;
  }
  for (int i=2; i <= n-k; i++) {
    b = ((long long)b * i) % m;
  }
  for (int i=2; i <= k; i++) {
    c = ((long long)c * i) % m;
  }
  b = ((long long) b * c) % m;
  return ((long long) a * modulo_inverse_euclidean(b, m)) % m;
}

Nguyen Gia

Bình luận


White
{{ comment.user.name }}
Bỏ hay Hay
{{comment.like_count}}
Male avatar
{{ comment_error }}
Hủy
   

Hiển thị thử

Chỉnh sửa

White

Nguyễn Gia

6 bài viết.
44 người follow
Kipalog
{{userFollowed ? 'Following' : 'Follow'}}
Cùng một tác giả
White
17 12
Cảm xúc khi học javascript y hệt như bài viết này nên mình dịch lại cho vui. Link gốc: https://hackernoon.com/howitfeelstolearnjavascriptin2016d3a7...
Nguyễn Gia viết hơn 1 năm trước
17 12
White
12 1
Thi thoảng la liếm github hay gặp một số từ viết tắt mà chẳng hiểu nghĩa là gì nên tổng hợp vào đây. | Từ viết tắt | Nghĩa | | |::| |AFAICT|As far...
Nguyễn Gia viết hơn 1 năm trước
12 1
White
5 0
20160506 Goroutine Nếu goroutine đang chiếm quyền xử lí ko buông, thì các goroutine khác không được chạy. Cách đơn giản nhất để buông quyền là d...
Nguyễn Gia viết 2 năm trước
5 0
{{like_count}}

kipalog

{{ comment_count }}

bình luận

{{liked ? "Đã kipalog" : "Kipalog"}}


White
{{userFollowed ? 'Following' : 'Follow'}}
6 bài viết.
44 người follow

 Đầu mục bài viết

Vẫn còn nữa! x

Kipalog vẫn còn rất nhiều bài viết hay và chủ đề thú vị chờ bạn khám phá!