Thuật toán Miller - Kiểm tra số nguyên tố

Miller là một biến thể của thuật toán Rabin-Miller, dùng để kiểm tra xem một số nguyên có phải là số nguyên tố hay không.

Bài này có dùng các kiển thức về đồng dư, nếu bạn chưa biết thì nên đọc bài toán đồng dư.

Thuật toán Miller

Xét số lẻ \(n > 1\), để kiểm tra xem \(n\) có phải là số nguyên tố hay không, ta làm như sau:

• Phân tích \(n - 1\) thành dạng \(2 ^ s \times d\) (\(s\) và \(d\) là số nguyên dương, \(d\) lẻ)
• Xét mỗi số nguyên \(a\) trong khoảng \([2, \min(n-1, \lfloor 2(\ln n) ^ 2 \rfloor)]\)
  • Nếu \(a^d \not\equiv 1 \pmod n\) và \((a ^ d)^{2^r} \not\equiv -1 \pmod n\) với mọi \(r\) từ \(0\) đến \(s-1\) thì \(n\) không phải là số nguyên tố.
• Nếu \(n\) vượt qua được tất cả các lần thử với \(a\) ở trên thì \(n\) là số nguyên tố.

Sử dụng Miller trong thực tế

Một cách tổng quát thì thuật toán này chưa được chứng minh là đúng. Tính đúng đắn của nó phụ thuộc vào giả thuyết Reimann.

Tuy nhiên, với \(n\) nhỏ (\(n < 3,317,044,064,679,887,385,961,981\)) thì thuật toán đã được kiểm tra và chứng minh là đúng. Ngoài ra, ta cũng không cần kiểm tra hết tất các các số nguyên \(a\) như ở trên, mà chỉ cần dùng một vài số là đủ.

Bảng sau được sử dụng để kiểm tra số nguyên tố trong thực tế:

Xem bảng trên, ta thấy chỉ cần dùng 3 số thử: 2, 7, 61 là đủ để kiểm tra tính nguyên tố của mọi số nguyên 32-bit.

Để hiểu rõ hơn cách dùng, bạn đọc xem thêm phần dưới.

Cài đặt thuật toán

Trong phần này mình sẽ dùng Miller để giải 2 bài PNUMBER và PRIME1.

Thuật toán gồm các giai đoạn chính:

• Phân tích \(n-1\) thành dạng \(2^s \times d\).
• Thử lần lượt các số \(a\) trong bảng ở trên.
  • Khi thử, ta cần dùng hàm lũy thừa nhanh để tính \(a^d \bmod n\).

Ta sẽ lần lượt cài đặt các hàm cơ bản trên. Trước hết, ta khai báo các thư viện và kiểu dữ liệu cần sử dụng:

#include <tuple>
#include <iostream>
using namespace std;
typedef unsigned long long ll;

Hàm phân tích thành dạng \(2^s \times d\):

pair<ll, ll> factor(ll n) {
    ll s = 0;
    while ((n & 1) == 0) {
        s++;
        n >>= 1;
    }
    return {s, n};
}

Hàm tính \(a^d \bmod n\):

ll pow(ll a, ll d, ll n) {
    ll result = 1;
    a = a % n;
    while (d > 0) {
        if (d & 1) result = result * a % n;
        d >>= 1;
        a = a * a % n;
    }
    return result;
}

Hàm kiểm tra xem “nếu \(a^d \not\equiv 1 \pmod n\) và \((a ^ d)^{2^r} \not\equiv -1 \pmod n\) với mọi \(r\) từ \(0\) đến \(s-1\) thì \(n\) không phải là số nguyên tố”:

// Trả về false nếu n không phải là số nguyên tố.
bool test_a(ll s, ll d, ll n, ll a) {
    if (n == a) return true;
    ll p = pow(a, d, n);
    if (p == 1) return true;
    for (; s > 0; s--) {
        if (p == n-1) return true;
        p = p * p % n;
    }
    return false;
}

Hàm cuối cùng: Thử kiểm tra \(n\) với các \(a\) khác nhau. Vì ta đang giải bài PNUMBER, giới hạn số trong bài này là \(n < 200000\), tra bảng trên ta thấy chỉ cần kiểm tra với \(a = 2\) và \(a = 3\) là đủ:

// Trả về true nếu n là số nguyên tố,
// false nếu n là hợp số.
bool miller(ll n) {
    if (n < 2) return false;
    if ((n & 1) == 0) return n == 2;
    ll s, d;
    tie(s, d) = factor(n-1);
    return test_a(s, d, n, 2) && test_a(s, d, n, 3);
}

Xem full code giải bài PNUMBER bằng thuật toán Miller: /code/139

Bài tương tự với PNUMBER nhưng có giới hạn cao hơn: PRIME1.

Giới hạn của bài PRIME1 là \(n \leq 10^9\), nên ta cần dùng \(a = 2\), \(a = 7\) và \(a = 61\). Xem code: /code/138.

Xử lý tràn số

Bảng ở trên có thể kiểm tra đến các số lớn hơn 64-bit, tuy nhiên khi cài đặt, nếu không sử dụng kiểu số nguyên lớn thì ta chỉ có thể sử dụng thuật toán với các số 32-bit.

Lý do là thuật toán có dùng phép nhân. Khi nhân 2 số 32-bit kết quả phải được chứa trong số 64-bit, nhân 2 số 64-bit thì kết quả phải được chứa trong số 128-bit. Vì vậy, để không bị tràn số thì thuật toán chỉ có thể dùng với số 32-bit.

Để khắc phục điều này, ta có thể tự viết hàm để vừa thực hiện nhân vừa thực hiện mod.

Đây là code tính \(a \times b \bmod n\) của Dana Jacobsen trong một câu trả lời trên Quora:

#include <stdint.h>

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t n) {
    // if (a >= n) a %= n;   /* Careful attention from the caller */
    // if (b >= n) b %= n;   /* should make these unnecessary.    */
    uint64_t r = 0;
    if ((a|b) < (1ULL << 32)) return (a*b) % n;
    if (a < b) { uint64_t t = a; a = b; b = t; }
    if (n <= (1ULL << 63)) {
        while (b > 0) {
            if (b & 1)  { r += a;  if (r >= n) r -= n; }
            b >>= 1;
            if (b)      { a += a;  if (a >= n) a -= n; }
        }
    } else {
        while (b > 0) {
            if (b & 1)  r = ((n-r) > a) ? r+a : r+a-n;    /* r = (r + a) % n */
            b >>= 1;
            if (b)      a = ((n-a) > a) ? a+a : a+a-n;    /* a = (a + a) % n */
        }
    }
    return r;
}

Đoạn code trên hơi phức tạp, bạn đọc có thể đọc đoạn code sau để hiểu tư tưởng thuật toán:

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t n) {
    uint64_t r = 0;
    while (b > 0) {
        if (b % 2 != 0) r = (r + a) % n;
        a = a * 2 % n;
        b = b / 2;
    }
    return r;
}

Đoạn code ở dưới đơn giản hơn, tuy nhiên nó chạy chậm hơn và chỉ chạy đúng vớn n có tối đa 63 bit.

Thay thế các phép nhân trong code Miller ở phần trên bằng hàm mulmod, ta có thể sử dụng được Miller cho các số 64-bit.

Độ phức tạp

Độ phức tạp của hàm mulmod là \(O(\log n)\).

Mỗi lần kiểm tra tốn \(O(s + \log d) = O(\log n)\) phép nhân, kết hợp với phép ĐPT của phép nhân nữa là \(O(\log ^ 2 n)\).

Vậy độ phức tạp của Miller là \(O(k \log ^ 2 n)\) với \(k\) là số các số \(a\) cần kiểm tra.