Algorithm: kiểm tra số nguyên tố

Kiểm tra tính nguyên tố của một số luôn là một vấn đề “đau đầu”. Số nguyên tố luôn là một trong số những vấn đề toán học hấp dẫn, cũng vì thế mà các kỹ thuật kiểm tra số nguyên tố luôn luôn được phát minh, cải tiến nhằm đáp ứng nhu cầu thực tế hiện nay là tìm ra các số nguyên tố càng lớn càng tốt.

Trong bài viết này, tôi sẽ trình bày một số kỹ thuật kiểm tra tính nguyên tố của một số cho trước. Lưu ý, một số kỹ thuật trong bài viết này thực ra là kiểm tra tính hợp số, tuy nhiên, cũng dựa vào đó mà chúng ta có thể biết được một số có là số nguyên tố hay không.

Nhắc lại kiến thức phổ thông

Trước hết, xin nhắc lại một số kiến thức về định nghĩa số nguyên tố cũng như hợp số. Đây là những kiến thức chúng ta đều đã được học, nên bạn có thể bỏ qua phần này.

• Số nguyên tố là một số tự nhiên dương có đúng hai ước là 1 và chính nó.
• Hợp số là số chia hết cho các số khác ngoài 1 và chính nó.

Với định nghĩa trên thì số 0, số 1 không phải số nguyên tố, cũng không phải là hợp số. Ngoài ra, các số tự nhiên khác, thì chỉ có thể là số nguyên tố hoặc hợp số. Vì vậy để kiểm tra tính nguyên tố của một số, chúng ta có thể kiểm tra gián tiếp thông qua tính hợp số của nó.

Hơi ngoài lề một chút, nhưng chúng ta có thể hiểu số nguyên tố là những số không thể phân tích thành tích của những số tự nhiên khác, vì chúng chẳng có ước nào như thế. Còn hợp số có thể phân tích thành thích của các số khác, có lẽ cái tên “hợp số” được đặt cũng vì lý do này. Phân tích hợp số thành nhân tử thì cuối cùng, chúng ta sẽ thu được một tích của toàn các số nguyên tố.

Đây chỉ là những kiến thức rất cơ bản, trong phần tiếp theo, tôi sẽ trình bày một số phương pháp kiểm tra tính nguyên tố của một số.

Phương pháp chân phương

Đây là một phương pháp cổ điển nhất và cũng dễ hiểu nhất. Phương pháp này hoàn toàn dựa trên những kiến thức mà chúng ta đã được học ở trường.

Ý tưởng của phương pháp này rất đơn giản: với một số n cho trước, kiểm tra các số trong khoảng [2, n - 1] nếu có số nào là ước của n thì nó là hợp số, còn lại nó là số nguyên tố. Trường hợp số 0, 1 quá rõ ràng nên không xét ở đây.

>>> def is_prime(n):
...     return not (n < 2 or any(n % i == 0 for i in range(2, n - 1)))
...
>>> [i for i in range(100) if is_prime(i)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Phương pháp trên rất dễ hiểu nhưng chưa được tối ưu cho lắm. Để tối ưu hơn, chúng ta có thể giảm số vòng lặp bằng cách phương pháp như sau:

• Không cần duyệt tất cả các số trong khoảng [2, n - 1] mà chỉ cần duyệt các số trong khoảng từ [2, sqrt(n)].
• Không cần duyệt tất cả các số tự nhiên, chỉ cần duyệt các số nguyên tố trong khoảng đó là được.

Những lý thuyết này đã được chứng minh về mặt toán học, chúng ta cũng không cần sa đà vào chuyện đó làm gì. Ở đây chúng ta cần tìm cách vận dụng nó là đủ. Tuy nhiên, vấn đề ở đây là rất khó để có thể duyệt qua chỉ các số nguyên tố. Bởi vì điều đó đòi hỏi chúng ta phải có danh sách các số nguyên tố sẵn rồi.

Chúng ta có thể biến đổi điều này một chút như sau: Các số nguyên tố sẽ có dạng 6k ± 1 (ngoại trừ 2, 3) (không tin ư, bạn có thể tìm một phản ví dụ thử xem). Vì vậy chúng ta có thể kiểm tra việc chia hết cho 2, 3 trước rồi kiểm tra việc chi hết cho các số 6k ± 1 này.

>>> def is_prime(n):
...     if n in (2, 3):
...         return True
...     return not (n < 2 or
...                 n % 2 == 0 or
...                 n % 3 == 0 or
...                 any(n % i == 0 or n % (i + 2) == 0
...                     for i in range(5, int(n ** 0.5) + 1, 6)))
...
>>> [i for i in range(100) if is_prime(i)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Dùng regular expression

Chúng ta có thể dùng regular expression để kiểm tra tính nguyên tố của một số cho trước. Xem thêm ở đây

Ngoại trừ, 0 và 1 là các trường hợp đặc biệt, các trường hợp khác đều là số nguyên tố, hoặc hợp số. Giả sử n là một hợp số, thì sẽ tồn tại hai số tự nhiên sao cho n = a * b (với 1 < a, b < n).

Bây giờ, để áp dụng regular expression kiểm tra các số nguyên tố, chúng ta sẽ biểu diễn một số tự nhiên thành một string toàn chữ số 1 như sau:

0 = ''
1 = '1'
2 = '11'
3 = '111'
...

Các phép tính sẽ như sau:

1 + 2 = '1' + '11' = '111' = 3
2 * 3 = '11' * 3 = '111111' = 6

Với quy tắc như trên, chúng ta có thể kiểm tra một số có phải là hợp số hay không bằng regular expression. Chúng ta cần kiểm tra xem một số có thoả mãn biểu thức a * b hay không.

a sẽ được biểu diễn bằng một string, không rõ độ dài bao nhiêu, nhưng chắc chắn là lớn hơn 2. Vì vậy chúng ta sẽ dùng 11+ trong regex, sau đó chúng ta cần kiểm tra xem số a này có được lặp lại hay không (nghĩa là có tồn tại số b sao cho a * b = n hay không). Vì vậy chúng ta có thể dùng regex như sau ^(11+)\1+$.

Tuy nhiên, regex như vậy chưa tối ưu lắm, chúng ta có thể tìm số a nhỏ nhất có thể bằng cú pháp tìm kiếm non-greedy trong regex như sau ^(11+?)\1+$. Cuối cùng, chúng ta bổ sung thêm regex cho trường hợp đặc biệt là số 0 và 1 là được một regex hoàn hảo để kiểm tra tính nguyên tố của một số.

>>> import re
>>> def is_prime(n):
...     return not re.match(r'^1?$|^(11+?)\1+$', '1' * n)
...
>>> [i for i in range(100) if is_prime(i)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Phương pháp nâng cao

Các phương pháp ở phần trên cho chúng ta độ chính xác rất cao (100%), tuy nhiên, chúng có một nhược điểm là số vòng lặp quá nhiều. Vì vậy phương pháp đó chỉ thích hợp với những số tự nhiên nhỏ, với những số lớn hơn, dài tới hàng trăm chữ số, chúng ta cần một phương pháp hiệu quả hơn.

Trong phần này, tôi sẽ trình bày một số phương pháp nâng cao giúp chúng ta kiểm tra tính nguyên tố của một số. Cần lưu ý rằng, các phương pháp trong phần này không cho kết quả chính xác tuyệt đối mà chỉ cho chúng ta kết luận rằng một số cho trước có thể là số nguyên tố hay không. Tất nhiên, chúng ta sẽ cài đặt thuật toán sao cho tỉ lệ chính xác cao nhất có thể.

Phương pháp Fermat

Thực ra phương pháp này gọi là phương pháp Fermat nhưng không phải do Fermat nghĩ ra, mà là nó dựa vào định lý nhỏ Fermat để thực hiện. Nội dung định lý như sau:

Nếu n là một số nguyên tố thì với mọi số a nguyên tố cùng n, ta luôn có a ^ (n - 1) ≡ 1 (mod n) hoặc a ^ (n - 1) % n == 1

Thực ra Fermat chỉ đưa ra định lý chứ không chứng minh. Tuy nhiên, giờ đây, định lý đã được chứng minh bằng nhiều cách khác nhau. Lưu ý rằng, định lý trên chỉ có một chiều, nghĩa là nếu n là số nguyên tố thì chắc chắn a ^ (n - 1) % n == 1 nhưng nếu a ^ (n - 1) % n == 1 thì n vẫn có thể là hợp số.

Từ định lý này, nếu ta muốn kiểm tra số n có là nguyên tố không, ta lấy ngẫu nhiên các số a < n và kiểm tra xem đẳng thức trên có đúng không. Nếu nó không đúng với một giá trị a nào đó thì n là hợp số. Nếu đẳng thức đúng với nhiều giá trị của a, ta có thể nói rằng n “có khả năng cao” là số nguyên tố. Tất nhiên là chúng ta sẽ tìm cách tăng “khả năng” này lên bằng cách thực hiện nhiều tính toán hơn.

>>> from random import randint
>>> loop_count = 5
>>> def is_prime(n):
...     return not (n < 2 or
...                 any(pow(randint(1, n - 1), n - 1, n) != 1
...                     for _ in range(loop_count)))
...
>>> [i for i in range(100) if is_prime(i)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Lưu ý rằng, kiểm tra bằng phương pháp này vẫn có thể cho kết quả sai (vì định lý nhỏ Fermat là một chiều như đã nói ở trên). Chúng ta có thể sử dụng phương pháp Fermat này khi cần thực hiện các tính toán nhanh với độ chính xác vừa đủ chấp nhận được. Ví dụ chúng ta dùng phương pháp này để lọc bớt một lượng lớn hợp số sau đó thực hiện kiểm tra nâng cao với những số “có thể là nguyên tố” còn lại thì lượng tính toán đã giảm đi đáng kể.

Phương pháp Miller – Rabin

Phương pháp Miller – Rabin là một phương pháp kiểm tra tính nguyên tố của một số, thực chất là kiểm tra tính hợp số của nó. Như những kiến thức chúng ta đã học, một số tự nhiên lớn hơn 1 thì chắc chắn sẽ là số nguyên tố hoặc hợp số. Phương pháp này được cải tiến nhiều lần, và hiện nay nó là một phương pháp mang tính xác suất. Về cài đặt, nó cũng không phức tạp lắm.

Cơ sở lý thuyết của phương pháp này như sau:

• Định lý nhỏ Fermat cho chúng ta biết rằng, nếu n là một số nguyên tố thì với mọi số a sao cho 1 <= a < n thì a ^ (n - 1) % n = 1.
• n chắc chắn là số lẻ (đương nhiên, vì 2 là số nguyên tố chẵn duy nhất mà 2 thì chúng ta không cần dùng đến phương pháp này). Vì n lẻ nên n - 1 là một số chẵn và nó có thể biểu diễn dưới dạng n - 1 = d * 2 ^ s, trong đó d là số lẻ và s > 0.
• Nhờ hai điều trên, với một số bất kỳ trong khoảng [2, n - 2], thì a ^ (d * 2 ^ r) % n phải bằng 1.
• Một tính chất toán học quan trọng khác, đó là nếu x ^ 2 % n = 1 có nghĩa là x ^ 2 - 1 = (x + 1)(x - 1) % n = 0, vậy nếu n là số nguyên tố thì x + 1 hoặc x - 1 phải chia hết cho n.

Dựa trên lý thuyết như vậy, dưới đây là thuật toán của phương pháp này. Lưu ý rằng, phương pháp này cũng tương tự như phương pháp Fermat, chúng ta chỉ có kết quả là một số “có khả năng” là số nguyên tố chứ không thể chắc chắn 100% được.

write n − 1 as 2r·d with d odd by factoring powers of 2 from n − 1
WitnessLoop: repeat k times:
   pick a random integer a in the range [2, n − 2]
   x ← ad mod n
   if x = 1 or x = n − 1 then
      continue WitnessLoop
   repeat r − 1 times:
      x ← x2 mod n
      if x = 1 then
         return composite
      if x = n − 1 then
         continue WitnessLoop
   return composite
return probably prime

Độ chính xác của phương pháp sẽ tăng lên tuỳ thuộc vào số phép tử mà chúng ta thực hiện. Ngoài ra, nó còn phụ thuộc vào việc sinh các số ngẫu nhiên nữa.

>>> import random
>>> _mrpt_num_trials = 5
>>> def is_probable_prime(n):
...     assert n >= 2
...     if n == 2:
...         return True
...     if n % 2 == 0:
...         return False
...     s = 0
...     d = n - 1
...     while True:
...         quotient, remainder = divmod(d, 2)
...         if remainder == 1:
...             break
...         s += 1
...         d = quotient
...     assert 2 ** s * d == n-1
...     def try_composite(a):
...         if pow(a, d, n) == 1:
...             return False
...         for i in range(s):
...             if pow(a, 2**i * d, n) == n-1:
...                 return False
...         return True
...     for i in range(_mrpt_num_trials):
...         a = random.randrange(2, n)
...         if try_composite(a):
...             return False
...     return True
...
>>> [i for i in range(2, 100) if is_probable_prime(i)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
>>> is_probable_prime(743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153)
False
>>> is_probable_prime(643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153)
True

Phương pháp Solovay – Strassen

Phương pháp Solovay – Strassen cũng là một phương pháp kiểm tra tính nguyên tố dựa trên xác suất. Tuy nhiên, phương pháp này tương đối khó hiểu do liên qua đến một số khái niệm sau:

Ký hiệu Legendre

Ký hiệu Legendre được định nghĩa cho một cặp số tự nhiên, a, p trong đó p là một số nguyên tố. Ký hiệu Legendre được viết là a/p là tính như sau:

      = 0    nếu  a % p = 0
(a/p) = 1    nếu tồn tại một số tự nhiên k sao cho k^2 % p = a
      = -1   với các trường hợp khác.

Euler đã chứng minh rằng

(a/p) = a ^ ((p - 1) / 2) % p

Ký hiệu Jacobian Symbol

Ký hiệu Jacobi là tổng quát hóa của Ký hiệu Legendre, trong đó thay p bằng một số tự nhiên n bất kỳ.

Nếu n = p1k1 * .. * pnkn thì ký hiệu Jacobian sẽ như sau:

(a/n) = ((a/p1) ^ k1) * ((a/p2) ^ k2) * ..... * ((a/pn) ^ kn)

Nếu n là một số nguyên tố thì ký hiệu Jacobian và ký hiệu Legendre sẽ bằng nhau. Ký hiệu Jacobian có một số tính chất như sau:

• (a/n) = 0 nếu gcd(a, n) != 1.
• (ab/n) = (a/n) * (b/n).

• Nếu a là một số chẵn thì (a/n) = (2/n) * ((a/2)/n). Chúng ta dễ dàng thấy rằng:

        = 1 if n = 1 ( mod 8 ) or n = 7 ( mod 8 )
  (2/n) = -1 if n = 3 ( mod 8 ) or n = 5 ( mod 8 )
        = 0 otherwise
  

• Nếu cả a, n đều là số lẻ thì (a/n) = (n/a) * (-1)^((a-1) * (n-1) /4)).

Thuật toán

Với một số n cần kiểm tra tính nguyên tố, chúng ta chọn ra một số ngẫu nhiên trong khoảng [2, n-1] và tính Jacobian (a/n). Nếu n là số nguyên tố thì ký hiệu Jacobian phải bằng ký hiệu Legendre và bằng kết quả mà Euler đã chứng minh.

Tương tự như các thuật toán phía trên, thuật toán này dựa nhiều vào xác xuất, độ chính xác của nó phụ thuộc vào số phép thử mà chúng ta thực hiện cũng như việc sinh số ngẫu nhiên.

Cài đặt

>>> from random import randrange
>>> _sspt_num_trials = 10
>>> def is_probable_prime(p):
...     if p < 2:
...         return False
...     if p in (2, 3):
...         return True
...     if p % 2 == 0:
...         return False
...     def jacobian(a, n):
...         if not a:
...             return 0
...         ans = 1
...         if a < 0:
...             a = -a
...             if n % 4 == 3:
...                 ans -= ans
...         if a == 1:
...             return ans
...         while a:
...             if a < 0:
...                 a = -a
...                 if n % 4 == 3:
...                     ans = -ans
...             while a % 2 == 0:
...                 a = a // 2
...                 if n % 8 in (3, 5):
...                     ans = -ans
...             a, n = n, a
...             if a % 4 == 3 and n % 4 == 3:
...                 ans = -ans
...             a = a % n
...             if a > n // 2:
...                 a = a - n
...         if n == 1:
...             return ans
...         return 0
...     for _ in range(_sspt_num_trials):
...         a = randrange(2, p - 1)
...         x = (p + jacobian(a, p)) % p
...         y = pow(a, (p - 1) // 2, p)
...         if not x or y != x:
...             return False
...     return True
...
>>> [i for i in range(2, 100) if is_probable_prime(i)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
>>> is_probable_prime(743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153)
False
>>> is_probable_prime(643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153)
True

Kết luận

Kiểm tra tính nguyên tố của một số tự nhiên thực sự không phải là một công việc dễ dàng, nhưng tôi đảm bảo nó là một công việc thú vị. Hy vọng bài viết cung cấp cho các bạn nhiều phương thức khác nhau để làm điều đó, và một trong số chúng sẽ giúp ích cho các bạn trong công việc.