Try   HackMD

Randomizer & các ứng dụng trong lập trình thi đấu

Author: Gã coder si tình

Máy tính vốn được sinh ra để tính toán các phép tính phức tạp với độ chính xác có thể nói là tuyệt đối mà con người có siêu việt đến cỡ nào cũng không làm được. Thế nhưng, ta vẫn có thể sử dụng những cỗ máy nổi tiếng vì sự chính xác ấy để tạo ra số ngẫu nhiên. Chúng có thể sinh hàng tỉ số (thậm chí hơn) chỉ trong vài giây với khả năng đoán trước cực thấp và sự phân bố của mẫu số cũng rất đều.

Vậy làm thế nào mà máy tính có thể sinh ra số ngẫu nhiên? Và việc sinh ngẫu nhiên có ứng dụng gì trong Tin học nói chung và lập trình thi đấu nói riêng?

Chúng ta sẽ cùng tìm hiểu trong bài viết này nhé.

Thế nào là ngẫu nhiên

In common usage, randomness is the apparent or actual lack of pattern or predictability in information.

Wikipedia

Tạm dịch theo định nghĩa của Wikipedia, sự ngẫu nhiên là một sự kiện thiếu đi quy luật hoặc khả năng dự đoán trước được.

Thật ra theo mình thì không có gì là hoàn toàn ngẫu nhiên, không có quy luật. Ví dụ như việc tung xúc sắc, theo lý thuyết, nếu đó là

1 viên xúc sắc bình thường thì cả
6
mặt của nó đều có xác suất là
1/6
. Tuy nhiên, giá trị nhận được từ việc tung xúc sắc thực chất không phải ngẫu nhiên mà nó phụ thuộc vào các yếu tố: trạng thái của nó trước khi được tung, trọng lượng của xúc sắc, vận tốc đầu, bề mặt mà nó rơi xuống, rất nhiều yếu tố khác. Có thể thấy, nó cũng có quy luật, nhưng có rất nhiều yếu tố khác nhau quyết định đến kết quả cuối cùng. Nếu viết ra thành công thức thì cực kỳ phức tạp, do đó, người ta coi như xúc sắc sẽ cho ra giá trị "ngẫu nhiên".

Tương tự, với máy tính, những "yếu tố" đó sẽ là tham số truyền vào hàm sinh ngẫu nhiên hoặc những biến số khác. Máy tính tuy là một thứ được sinh ra để đảm bảo độ chính xác nhưng nó vẫn có thể sinh ra "số ngẫu nhiên", nếu ta có thể tạo ra một thuật toán biến đổi các tham số đưa vào một cách đủ phức tạp để cho ra đáp án có thể coi là ngẫu nhiên.

Phân loại randomizer

Pseudo-random Number Generator (PRNG) hoạt động như ý tưởng được trình bày ở trên, máy tính sẽ biến đổi các tham số đưa vào (gọi là seed) một cách đủ phức tạp, để khả năng dự đoán trước được gần như là con số

0. Đây là loại randomizer rất phổ biến trong Tin học và cũng sẽ là chủ đề chính của bài viết này. Một số ví dụ của PRNG là hàm rand(), mt19937 trong C++, thư viện random trong Python, hay randomizer của trang CalculatorSoup.

True Random Number Generator (TRNG) sinh ra số ngẫu nhiên dựa trên một sự kiện không liên quan gì đến máy tính, ví dụ như tiếng ồn trong khí quyển, các thông số liên quan đến chất phóng xạ, Ví dụ của TRNG là trang RANDOM.ORG.

Sinh số ngẫu nhiên trong một số ngôn ngữ lập trình

C++

Hàm rand()

Cách sinh số ngẫu nhiên đơn giản nhất trong ngôn ngữ lập trình C++ đó là sử dụng hàm rand() [1]. Hàm này sẽ sinh

1 số nguyên ngẫu nhiên từ
0
đến RAND_MAX, với RAND_MAX là một giá trị tùy thuộc vào máy, tối thiểu là
2151=32767
.

Ưu điểm của hàm rand() là cú pháp dễ nhớ và dễ sử dụng. Tuy nhiên, khuyết điểm của nó là rất lớn, mà neal (một Codeforces International Grandmaster) đã viết hẳn một bài blog để nói về điều này [2]. Mình sẽ tóm tắt 2 vấn đề chính mà neal đã đề cập:

  • Thứ nhất là giới hạn RAND_MAX của hàm này, tùy thuộc vào máy, trên Codeforces thì là đúng
    32767
    . Nếu ta cần sinh một dãy số với các giá trị trong đoạn
    [1;109]
    thì các giá trị được hàm sinh ra sẽ bị lệch giá trị kỳ vọng rất nhiều.
  • Thứ hai, hàm rand() được cài đặt bằng một thuật toán khá đơn giản, chỉ dựa trên một hàm tuyến tính [3] nên khả năng dự đoán trước khá cao.

MT19937

Randomizer thứ 2 có sẵn trong STL của C++ đó là mt19937. Đây là cài đặt của thuật toán Mersenne Twister [4] thuật toán sinh ngẫu nhiên mạnh nhất trong Tin học tính đến thời điểm hiện tại, được phát triển dựa trên số nguyên tố Mersenne

2199371. Hai người phát minh ra Mersenne Twister tên là Makoto Matsumoto và Takuji Nishimura.

Để sử dụng thuật toán này, ta khai báo một hàm rng() như sau:

mt19937 rng(SEED);

Các bạn có thể chọn seed tùy ý, trong bài viết này mình sẽ sử dụng thời gian. Ở đây, có thể gọi hàm time(0), nhưng với seed như thế thì khi tham gia Codeforces sẽ dễ bị hack. Do đó, mình sẽ dùng đến thư viện chrono.

Thư viện chrono cung cấp

3 loại đồng hồ chính đó là system_clock, steady_clockhigh_resolution_clock với độ chia nhỏ nhất khác nhau và đều có điểm bắt đầu (epoch) là
0
giờ, ngày
1/1/1970
(UTC). Mình sẽ sử dụng thời gian từ epoch đến khoảnh khắc mà hàm được khai báo, tính bằng high_resolution_clock làm seed như sau:

mt19937 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

Sau khi khai báo dòng này thì hàm rng() sẽ trả về một số nguyên ngẫu nhiên trong khoảng

[231;2311]. Nếu muốn tăng giới hạn lên
64
-bit thì ta sử dụng mt19937_64.

Để sinh 1 số ngẫu nhiên trong khoảng

[L;R], ta dùng uniform_int_distribution. Để tiện thì mình sẽ khai báo thêm một hàm khác như sau:

int randRange (int L, int R) {
    return uniform_int_distribution<int>(L, R)(rng);
}

Để kiểm tra chất lượng của mt19937, mình đã sinh thử một mẫu dữ liệu

105 số trong khoảng
[0;109]
và tính toán vài con số trong thống kê để thấy được distribution của mẫu dữ liệu này. Kết quả của lần chạy đầu tiên như sau:

  • Trung bình cộng (mean):
    500 045,548
  • Khoảng tứ phân vị (interquartile range):
    250 085,5
    //
    500 415,5
    //
    750 410
    .

Có thể thấy thuật toán này sinh ra số với distribution khá đều. Bên cạnh đó, mt19937 có thể sinh số ngẫu nhiên nhanh hơn hàm rand().

Python

Ngôn ngữ lập trình Python cũng có randomizer tương tự với thư viện random và hàm randrange sử dụng thuật toán Mersenne Twister. Cách truyền tham số vào hàm randrange khá giống hàm range, cụ thể như sau:

  • randrange(R) trả về một số ngẫu nhiên trong đoạn
    [0;R)
    .
  • randrange(L, R) trả về một số ngẫu nhiên trong đoạn
    [L;R)
    .
  • randrange(L, R, x) trả về một số ngẫu nhiên trong đoạn
    [L;R)
    với bậc là
    x
    .

Ví dụ đoạn code sau in ra

10 số ngẫu nhiên được chọn từ tập
Ω={2;5;8;11;14;17;20;23;26;29}
.

import random

for i in range(10):
    print(random.randrange(2, 32, 3))

Random pick with weights

Khi đã có thể sinh số ngẫu nhiên với xác suất đều, bây giờ ta sẽ tính đến bài toán sinh số ngẫu nhiên với xác suất không đều. Đây cũng là bài tập xuất hiện trên Leetcode [5] có tên Random pick with weights.

Việc sinh số ngẫu nhiên với xác suất không đều không có quá nhiều ứng dụng trong lập trình thi đấu, tuy nhiên ta có thể cần đến nó trong một số bài toán heuristic, hoặc các lĩnh vực Tin học khác như Machine Learning,

Trong bài viết này, mình sẽ trình bày nhanh thuật toán sử dụng binary search, cách giải này cũng đã được đề cập đến trong bài viết của page Code cùng RR [6].

Tóm tắt đề

Viết hàm sinh ngẫu nhiên với distribution không đều, cụ thể giá trị

i (với
0i<n
) sẽ có xác suất xuất hiện là
pip
.

Constraints:

  • 1n104
    .
  • 1pi105
    .
  • Hàm ngẫu nhiên sẽ được gọi không quá
    104
    lần.

Ví dụ cho mảng

p=[10,25,15,50],
p=100
, ta cần viết một hàm sinh ngẫu nhiên với kết quả là:

  • Giá trị
    0
    với xác suất
    10%
    .
  • Giá trị
    1
    với xác suất
    25%
    .
  • Giá trị
    2
    với xác suất
    15%
    .
  • Giá trị
    3
    với xác suất
    50%
    .

Ý tưởng

Ý tưởng của bài toán này đó là ta sẽ tạo ra

p "slots" cho các giá trị, rồi cho giá trị thứ
i
chiếm
pi
slots. Với ví dụ ở trên thì giá trị
0
sẽ chiếm các slots từ
110
, giá trị
1
chiếm các slots từ
1135
, giá trị
2
chiếm các slots từ
3650
,

Sau đó, ta sẽ sinh ra một số ngẫu nhiên trong khoảng

[1;p] rồi tính xem slot tương ứng là slot của giá trị nào. Để làm điều đó ta binary search trên mảng prefix sum của
p
.

Như vậy, độ phức tạp thời gian và bộ nhớ của phần tiền xử lý là

O(n), một lần gọi hàm tốn thời gian xử lý là
O(logn)
.

Cài đặt

Mình sẽ cài đặt trong class Solution, giống với yêu cầu trên Leetcode.

#include<bits/stdc++.h>

mt19937 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

class Solution {
    vector<int> pre;
public:
    Solution(vector<int>& w) : pre(w.size()) {
        pre[0] = w[0];
        for (int i = 1; i < w.size(); i++) pre[i] = pre[i - 1] + w[i];
        return;
    }

    int randRange (int l, int r) {
        return uniform_int_distribution<int>(l, r)(rng);
    }
    
    int pickIndex() {
        int val = randRange(1, pre.back());
        int k = -1;
        for (int b = pre.size() / 2; b >= 1; b /= 2) {
            while (k + b < pre.size() && pre[k + b] < val) k += b;
        }
        return k + 1;
    }
};

XOR Hash

XOR Hash là một trong những hàm Hash được sử dụng phổ biến trong lập trình thi đấu. Các bạn có thể đọc thêm về hàm này cũng như một số hàm Hash khác trong bài viết của page Page này được lập ra để Dế Mèn được gáy [7].

XOR Hash yêu cầu một hàm sinh ngẫu nhiên đủ chất lượng để có thể giảm thiểu sai số và cũng như khả năng bị hack khi tham gia các contests có open hack như trên Codeforces hoặc các dạng tương tự.

Ý tưởng chính

Ý tưởng chính của XOR Hash đó là để so sánh hai dãy

a
b
gồm các phần tử phân biệt và có cùng độ dài
m
không quan tâm thứ tự, ta sẽ tính tổng XOR các giá trị hash của các phần tử trong dãy. Tức là, ta sẽ kiểm tra:

ha1ha2...ham=hb1hb2...hbm

Với

là bitwise XOR và
h
là một dãy số (
k
-bit) được sinh ngẫu nhiên.

Dễ dàng thấy được, nếu hai dãy

a
b
giống nhau (không cần đúng thứ tự) thì biểu thức trên sẽ luôn đúng. Ngược lại, biểu thức trên có thể vẫn đúng (trường hợp xảy ra sai số) với tỉ lệ
12k
, nếu
k
đủ lớn, ta có thể coi như xác suất xảy ra sai số là
0%
.

Bài tập minh họa

Ở phần này, mình sẽ sử dụng bài tập minh họa là bài Perfect [8] (Bedao Grand Contest 10).

Tóm tắt đề

Cho hoán vị

p độ dài
n
q
truy vấn có dạng
(l,r)
, ở mỗi truy vấn, hãy cho biết khi sắp xếp dãy con
pl,pl+1,pl+2,...,pr
theo thứ tự tăng dần thì dãy con này có tạo thành
1
dãy số nguyên liên tiếp hay không?

Constraints:

  • 1n,q105
    .
  • 1pin
    .

Ý tưởng

Nhiệm vụ của chúng ta là so sánh dãy

pl,pl+1,pl+2,...,pr với dãy
v,v+1,v+2,...,v+sz1
với
v
sz
lần lượt là phần tử bé nhất và độ dài của dãy con.

Do

p là hoán vị nên ta có thể áp dụng XOR Hash, ta sẽ tính:

hplhpl+1hpl+2...hpr

rồi đem so sánh với:

hvhv+1hv+2...hv+sz1

Giá trị

v có thể được tính bằng Sparse table, và các giá trị Hash cũng có thể tính trong
O(1)
với prefix XOR. Ngoài ra, để giảm thiểu khả năng sai số, ta sẽ tạo ra dãy
h
là các số ngẫu nhiên
64
-bit, khi đó, xác suất xảy ra sai số là cực thấp (một số có
19
số
0
ngay sau dấu chấm thập phân).

Cài đặt

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int mn = 2e5 + 1;
ll preXor[mn], hashVal[mn], preArr[mn];
int spt[mn][18];

mt19937_64 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

ll queryXor (int l, int r) {
    return preXor[r] ^ preXor[l - 1];
}

ll queryArr (int l, int r) {
    return preArr[r] ^ preArr[l - 1];
}

int query (int l, int r) {
    int p = 31 - __builtin_clz(r - l + 1);
    return min(spt[l][p], spt[r - (1 << p) + 1][p]);
}

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);

    for (int i = 1; i < mn; i++) {
        hashVal[i] = rng();
        preXor[i] = preXor[i - 1] ^ hashVal[i];
    }

    int n, q; cin >> n >> q;
    for (int i = 1; i <= n; i++) {
        int x; cin >> x;
        spt[i][0] = x;
        preArr[i] = preArr[i - 1] ^ hashVal[x];
    }
    for (int s = 1; (1 << s) <= n; s++) {
        for (int i = 1; i + (1 << s) - 1 <= n; i++) {
            int p = s - 1;
            spt[i][s] = min(spt[i][p], spt[i + (1 << p)][p]);
        }
    }

    while (q--) {
        int l, r; cin >> l >> r;
        int val = query(l, r), sz = r - l + 1;
        cout << (queryArr(l, r) == queryXor(val, val + sz - 1)? "YES" : "NO") << "\n";
    }

    return 0;
}

Monte Carlo Algorithm

Monte Carlo là tên gọi chung của những thuật toán hoạt động dựa trên việc sinh nghẫu nhiên, có khả năng đưa ra kết quả sai với tỉ lệ sai số khá thấp (lưu ý Hash không phải là thuật toán Monte Carlo).

Cái tên Monte Carlo cũng chính là tên của một casino nổi tiếng ở Monaco, xuất phát từ việc rủi ro của những thuật toán này khá giống với rủi ro của đặt cược trong cờ bạc. Tuy nhiên, trong Tin học, ta có cách để giảm khả năng sai số của chúng xuống nhiều lần để đảm bảo độ chính xác của chương trình.

Trong bài viết này, chúng ta sẽ cùng tìm hiểu

2 thuật toán Monte Carlo có tính ứng dụng khá cao trong lập trình thi đấu.

Freivalds' algorithm

Freivalds' algorithm là thuật toán kiểm tra phép nhân ma trận

A×B có bằng một ma trận
C
cho trước hay không một cách hiệu quả mà không cần phải tính trực tiếp
A×B
.

Một chút về phép nhân ma trận

Khi nhân ma trận

A có kích thước
n×m
cho ma trận
B
có kích thước
m×q
, ta sẽ được
1
ma trận
C
có kích thước
n×q
sao cho:

Ci,j=k=1mAi,k×Bk,j

Để tính toàn bộ ma trận

C ta cần
O(n×m×q)
phép toán.

Do công thức trên hơi khó tưởng tượng nên các bạn có thể xem ảnh động dưới đây để dễ hình dung cách hoạt động của phép nhân ma trận:

Một tính chất quan trọng của phép nhân ma trận là tính chất kết hợp, tức là:

(A×B)×C=A×(B×C)

Tùy vào cách kết hợp mà phép nhân ma trận sẽ tốn độ phức tạp thời gian khác nhau.

Ý tưởng

Cho

3 ma trận
A
,
B
, và
C
có cùng kích thước
n×n
. Nhiệm vụ là kiểm tra xem
A×B
có bằng
C
hay không. Dĩ nhiên, ta sẽ không tính trực tiếp
A×B
trong
O(n3)
.

Thuật toán này sẽ sinh ngẫu nhiên một ma trận

r có kích thước
1×n
. Mỗi phần tử của
r
có giá trị là
0
hoặc
1
(để tránh tràn số).

Bây giờ ta kiểm tra:

r×A×B=r×C  ()

Để kiểm tra

(), ta chỉ cần thực hiện
3
phép nhân ma trận, cả
3
đều có độ phức tạp thời gian
O(n2)
. Hai ma trận nhận được sau khi tính toán có kích thước
1×n
và chỉ tốn thời gian
O(n)
để so sánh.

Sẽ có

2 trường hợp xảy ra:

  • Trường hợp
    A×B=C
    chắc chắn
    ()
    đúng.
  • Trường hợp
    A×BC
    ()
    có thể sai hoặc đúng. Người ta chứng minh được khả năng
    ()
    vẫn đúng trong trường hợp này là
    1/2
    [9] một xác suất sai số còn rất cao.

Nếu chỉ dừng lại ở đây với xác suất sai số lên đến

50% thì thuật toán này chẳng còn tính ứng dụng gì cả. Do đó ta chạy thuật toán
t
lần, mỗi lần sinh một ma trận
r
mới.

Ta kết luận

A×B=C nếu cả
t
lần chạy đều cho ra
()
đúng. Ngược lại, nếu có ít nhất
1
lần chạy cho ra
()
sai thì ta kết luận
A×BC
.

Lúc này, vẫn sẽ có trường hợp xảy ra sai số, là khi

A×BC nhưng cả
t
lần chạy đều cho ra
()
đúng. Dễ thấy, xác suất để trường hợp này xảy ra là
12t
. Nếu ta chọn
t
đủ lớn thì con số này sẽ trở nên rất nhỏ. Ví dụ, với
t=40
, xác suất sai số giảm xuống còn:

12400,00000000009095%

Như vậy, ta đã có một thuật toán chạy trong

O(t×n2) thay vì
O(n3)
mà còn đạt độ chính xác rất cao.

Cài đặt

Mình sẽ áp dụng thuật toán Freivalds vào bài Ma trận (vòng sơ khảo bảng C1, Hội thi Tin học trẻ toàn quốc năm 2022) [10]. Sau đây là phần cài đặt:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

mt19937 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

bool rnd() {
    return uniform_int_distribution<int>(0, 1)(rng);
}

ll a[501][501], b[501][501];
int n;

vector<ll> multi (vector<ll> vec, ll matr[501][501]) {
    vector<ll> nvec(n);
    for (int i = 0; i < n; i++) {
        for (int k = 0; k < n; k++) nvec[i] += matr[i][k] * vec[k];
    }
    return nvec;
}

bool freivalds() {
    vector<ll> vec(n), vec_a(n), vec_b(n);
    for (int i = 0; i < n; i++) vec[i] = rnd();

    vec_a = multi(vec, a);
    vec_a = multi(vec_a, a);
    vec_a = multi(vec_a, a);
    vec_b = multi(vec, b);

    for (int i = 0; i < n; i++) {
        if (vec_a[i] != vec_b[i]) return 0;
    }
    return 1;
}

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);

    int t; cin >> t;
    while (t--) {
        cin >> n;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) cin >> a[i][j];
        }
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) cin >> b[i][j];
        }

        bool ans = 1;
        for (int test = 0; test < 50; test++) ans &= freivalds();
        cout << (ans? "YES" : "NO") << "\n";
    }

    return 0;
}

MillerRabin primality test

MillerRabin primality test (primality test tạm dịch là kiểm tra tính nguyên tố) là thuật toán kiểm tra số nguyên tố được phát triển bởi Garry L. Miller và hoàn thiện bởi Michael O. Rabin vào năm 1980. Miller đã phát triển ý tưởng cho thuật toán này dựa trên định lý Fermat nhỏ và những tính chất toán học khác.

Thông thường, để kiểm tra xem một số nguyên dương

n có phải số nguyên tố hay không, ta cần thời gian
O(n)
, hoặc
O(1)
sau khi chạy sàng nguyên tố trong
O(nloglogn)
. Với thuật toán của MillerRabin thì ta chỉ tốn độ phức tạp
O(t×logn)
mà không cần tiền xử lý (mình sẽ giải thích biến
t
sau).

Tuy nhiên, là một thuật toán Monte Carlo, MillerRabin primality test có thể xảy ra sai số. Song, như đã nói, ta có thể giảm xác suất đó xuống mức cực thấp.

Định lý Fermat nhỏ & Fermat primality test

Trước khi nói đến ý tưởng của thuật MillerRabin thì mình sẽ nói một chút về định lý Fermat nhỏ [11] và một thuật toán kiểm tra số nguyên tố được phát triển trực tiếp trên định lý này [12].

Theo định lý Fermat nhỏ, ta có:

ap11modp  ()

Với

p là một số nguyên tố và
a
là một số nguyên không phải bội của
p
.

Dựa vào định lý này, nếu ta chọn một cơ số

a không phải bội của
p
mà khi thế vào
()
lại không thỏa, tức là
ap11modp
, thì ta có thể kết luận ngay
p
là một hợp số (không phải số nguyên tố). Khi đó, cơ số
a
đã được chọn để kiểm tra gọi là Fermat witness (tạm dịch là minh chứng Fermat cho việc
p
là hợp số).

Tuy nhiên, với một số hợp số

p, tồn tại
a
không phải bội của
p
()
vẫn thỏa. Khi đó ta gọi
a
Fermat liar (tạm dịch lời nói dối của Fermat cho việc
p
là hợp số?). Đây là lúc mà Fermat primality test xảy ra sai số.

Để khắc phục điều này, ta sẽ thử với nhiều cơ số

a khác nhau bằng cách chọn ngẫu nhiên. Nếu không tìm được cơ số
a
nào làm Fermat witness thì ta kết luận
p
có khả năng rất cao là số nguyên tố.

Để dễ hình dung hơn, các bạn có thể xem decision tree cho Fermat primality test với số lần thử là

t=4. Mũi tên hồng biểu diễn trường hợp
()
không thỏa (đồng nghĩa với việc tìm được Fermat witness cho
p
), mũi tên xanh biểu diễn trường hợp
()
thỏa.

Hợp số Carmichael

Hợp số Carmichael [13] có thể được coi là "khắc tinh" của Fermat primality test bởi tính chất đặc biệt của nó. Với một hợp số Carmichael

n, mọi cơ số
a
nguyên tố cùng nhau với
n
đều là Fermat liar. Tức là
an11modn
luôn đúng với
n
là hợp số Carmichael và
gcd(a,n)=1
.

Một ví dụ cho hợp số Carmichael là

172 081=7×13×31×61. Các bạn có thể xem thêm dãy số Carmichael trên OEIS có mã A002997 [14].

Để Fermat primality test phát hiện được một số Carmichael là hợp số, nó phải sinh ngẫu nhiên ra một cơ số

a thỏa
gcd(a,n)>1
, mà khả năng điều này xảy ra là rất thấp.

Tuy nhiên, hợp số Carmichael khá hiếm (số càng lớn lại càng hiếm). Người ta đã tính được rằng có

646 số Carmichael không quá
109
, và khoảng
20 138 200
số Carmichael không quá
1021
. Bên cạnh đó ta có thể tính và lưu vài số Carmichael đầu tiên vào mảng hằng. Vì vậy, Fermat primality test vẫn được ứng dụng vào thực tế với độ chính xác khá cao.

Ý tưởng của MillerRabin primality test

MillerRabin primality test cũng được phát triển dựa trên định lý Fermat nhỏ

ap11modp (
p
là số nguyên tố,
a
không phải bội của
p
). Nhưng thay vì sử dụng trực tiếp định lý này, MillerRabin đưa ra nhận xét là ta chỉ quan tâm
p
lẻ (vì chỉ có
2
là số chẵn nguyên tố)
p1
chẵn. Từ đó, ta có thể biểu diễn
p1
dưới dạng
2s×d
(với
s
d
là số nguyên
1
). Nếu biểu diễn
p1
dưới dạng nhị phân thì
s
là số bit
0
liên tiếp bên phải
p1
d
p1
sau khi right shift
s
lần.

Thay

p1=2s×d và chuyển số
1
sang vế trái:

a2sd10modp

Ta thấy

a2sd1 có dạng hằng đẳng thức
x2y2
với
x=a2s1d
y=1
. Sau khi phân tích xong, ta lại áp dụng hằng đẳng thức vào
a2s1d1
, cứ như vậy, ta có:

(a2s1d+1)(a2s1d1)0modp

(a2s1d+1)(a2s2d+1)(a2s2d1)0modp

(a2s1d+1)(a2s2d+1)(a2s3d+1)(a2s3d1)0modp

(a2s1d+1)(a2s2d+1)(a2s3d+1)(a20d+1)(a20d1)0modp  ()

Cộng với một số phân tích toán học, MillerRabin kết luận một số có khả năng cao là số nguyên tố khi

1 trong
2
điều kiện này đúng:

  • ad1modp
    .
  • Tồn tại
    1
    số nguyên
    r
    (
    0r<s
    ) sao cho
    a2rd1modp
    , hay
    a2rdp1modp
    (do
    1p1modp
    ).

Nếu cả

2 điều này đều không đúng thì ta kết luận ngay
p
là hợp số và cơ số
a
vừa được dùng để kiểm tra là minh chứng cho việc
p
là hợp số.

Lưu ý,

1 trong
2
điều này đúng
()
đúng
()
đúng. Điều ngược lại không xảy ra vì vẫn có trường hợp các thừa số đều không chia hết cho
p
nhưng tích lại chia hết cho
p
. Điều này khiến MillerRabin mạnh hơn Fermat primality test.

Với một số hợp số

p, vẫn tồn tại cơ số
a
sao cho
1
trong
2
điều kiện nêu trên đúng. Lúc này, ta gọi
a
strong liar.

Khác với Fermat, thuật toán của MillerRabin không có bất kỳ "khắc tinh" nào. Hơn nữa, người ta chứng minh được, với một cơ số

a ngẫu nhiên (
2ap2
), khả năng
a
là một strong liar
1/4
[15]. Đây là một xác suất không quá thấp nhưng tốt hơn rất nhiều so với Fermat primality test khi gặp số Carmichael.

Tương tự các thuật toán Monte Carlo khác, ta có thể chạy MillerRabin primality test

t lần, mỗi lần với một cơ số
a
ngẫu nhiên khác nhau. Có thể thấy, xác suất sai số lúc này chính là xác suất cả
t
lần cho ra cơ số
a
strong liar giảm xuống còn
14t
. Ta chỉ cần chọn
t20
thì xác suất sai số sẽ còn là một con số không đáng kể.

Cài đặt

Để cài đặt MillerRabin primality test, mình cần cài đặt sẵn hàm sinh ngẫu nhiên trong đoạn

[L;R] với xác suất đều và hàm tính mũ nhanh bằng chia để trị. Lưu ý, nếu muốn nhân số
64
-bit với nhau thì cần lưu vào biến
128
-bit để tránh tràn số.

typedef __uint128_t u128;
typedef long long ll;

mt19937_64 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

ll randRange (ll L, ll R) {
    return uniform_int_distribution<ll>(L, R)(rng);
}

u128 binpow (u128 a, u128 b, u128 MOD) {
    if (b == 0) return 1;
    u128 val = binpow(a, b / 2, MOD);
    if (b % 2) return val * val % MOD * a % MOD;
    return val * val % MOD;
}

Tiếp đến, mình sẽ cài hàm kiểm tra điều kiện MillerRabin, chỉ cần

1 điều kiện thỏa thì hàm sẽ trả về true ngay:

bool MillerRabin (ll p, ll a, ll d, int s) {
    u128 val = binpow(a, d, p);
    if (val == 1) return 1;
    for (int r = 0; r < s; r++) {
        if (val == p - 1) return 1;
        val = val * val % p;
    }
    return 0;
}

Ở đây, mình dùng biến val ban đầu lưu giá trị

admodp để kiểm tra điều kiện thứ nhất. Sau đó, mình sẽ chạy vòng lặp, cứ mỗi lần kiểm tra xong, val sẽ lưu giá trị
a2rdmodp
, để tăng
r
lên, ta chỉ việc bình phương biến val lên là xong, vì
a2rd×a2rd=a2rd+2rd=a2r+1d
.

Cuối cùng, mình sẽ cài hàm kiểm tra số nguyên tố, trước khi dùng đến hàm MillerRabin, mình sẽ loại hết các edge cases (

n<4, có thể loại luôn trường hợp
n
chẵn luôn để tối ưu), tính
s
d
bằng cách tận dùng hàm __builtin_ctzll(x) và right shift. Mình sẽ chọn
t=20
, đủ để xác suất sai số gần như là không có.

bool isPrime (ll n) {
    if (n < 4) return n == 2 || n == 3;

    int s = __builtin_ctzll(n - 1);
    ll d = (n - 1) >> s;

    for (int i = 0; i < 20; i++) {
        ll base = randRange(2, n - 2);
        if (MillerRabin(n, base, d, s) == 0) return 0;
    }
    return 1;
}

Dễ thấy,

slogn nên MillerRabin primality test có độ phức tạp thời gian
O(t×logn)
với
t
là số lần thử và
n
là số cần kiểm tra.

Nếu tự cài đặt, các bạn có thể kiểm tra code trên SPOJ [16].

Las Vegas algorithm

Song song với Monte Carlo, ta có Las Vegas algorithm, đây là tên gọi chung của những thuật toán hoạt động dựa trên việc sinh ngẫu nhiên, luôn đưa ra kết quả đúng, nhưng có lúc sẽ rơi vào độ phức tạp xấu nhất với xác suất khá thấp.

Trong phần này, ta sẽ tìm hiểu

2 thuật toán Las Vegas khá hay mà mình biết.

QuickSelect

QuickSelect là thuật toán tìm kiếm được phát triển bởi C. A. R. Hoare và xuất bản năm 1961 [17], thuật toán "anh em" của QuickSort một thuật toán Las Vegas khác rất nổi tiếng thông dụng (được sử dụng cho hàm sort trong C++) cũng do Hoare tạo ra.

Thuật toán này được sinh ra để giải bài toán

k-th order statistic, tức là tìm phần tử lớn thứ
k
trong một dãy số có
n
phần tử chưa được sắp xếp, được đánh số từ
0n1
. Ta có thể giả sử
0k<n
.

Ý tưởng

Thông thường, ta sẽ chạy một thuật toán sắp xếp nào đó trong

O(nlogn) rồi tìm phần tử thứ
k
trong
O(1)
. Tuy nhiên, liệu ta có cần sắp xếp cả dãy số chỉ để tìm
1
phần tử?

Ý tưởng chính của QuickSelect cũng khá giống QuickSort, là một thuật toán chia để trị như sau:

  • Với dãy có
    1
    phần tử, phần tử duy nhất đó chắc chắn sẽ là phần tử lớn thứ
    k
    (base case).
  • Để tìm phần tử lớn thứ
    k
    trong dãy có nhiều hơn
    1
    phần tử, ta chọn ngẫu nhiên một phần tử nào đó trong dãy làm pivot (tạm dịch là trục), gọi là
    x
    . Sau đó, tách dãy cần sắp xếp ra làm
    2
    dãy con: dãy đầu tiên gồm những phần tử
    x
    , dãy thứ hai gồm những phần tử
    >x
    . Gọi
    a
    b
    lần lượt là số phần tử của mỗi dãy.
    • Nếu
      k<a
      , ta tiếp tục tìm phần tử lớn thứ
      k
      trong dãy con đầu tiên.
    • Nếu
      ka
      , ta tìm phần tử lớn thứ
      ka
      trong dãy con thứ hai.

Đây là hình ảnh minh họa quá trình tìm phần tử lớn thứ

3 (đếm từ
0
) trong mảng ngẫu nhiên có
10
phần tử:

Với việc chọn pivot ngẫu nhiên, khả năng cao

2 dãy con được tạo ra sẽ tương đối bằng nhau. Khi đó, không gian tìm kiếm sẽ được chia đôi sau mỗi bước, dẫn đến độ phức tạp thời gian là:

n+n/2+n/4+n/8++n/2j<2n=O(n)

Độ phức tạp tệ nhất của QuickSelect là

O(n2), tuy nhiên xác suất là rất thấp, trong thực tế, chuyện này gần như không thể xảy ra.

Cài đặt

Mình sẽ cài hàm quickSelect tìm phần tử lớn thứ

k trong không gian tìm kiếm là vector<int> vec theo kiểu đệ quy.

mt19937 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

int randRange (int L, int R) {
    return uniform_int_distribution<int>(L, R)(rng);
}

int quickSelect (vector<int> vec, int k) {
    if (vec.size() == 1) return vec[0];

    // pick random pivot
    int pivot = vec[randRange(0, vec.size() - 1)];

    // partition vec into part1 and part2
    vector<int> part1, part2;
    for (auto u : vec) {
        if (u <= pivot) part1.push_back(u);
        else part2.push_back(u);
    }

    if (k < part1.size()) return quickSelect(part1, k);
    return quickSelect(part2, k - part1.size());
}

Mình đã chạy thử đoạn code này với vector<int>

2×106 phần tử, kết quả là Quickselect chạy nhanh hơn thuật toán sắp xếp khoảng
2
-
3
lần.

Eight Queens Problem

Eight Queens [18] là một bài toán kinh điển và rất phổ biến trong Tin học, thường được dùng làm ví dụ minh họa cho backtracking. Nhắc lại, bài toán yêu cầu tìm cách đặt

8 quân Hậu vào một bàn cờ
8×8
sao cho không có
2
quân Hậu nào tấn công nhau (theo luật cờ vua). Dưới đây là ảnh minh họa cho một cách đặt hợp lệ:

Với cách làm backtrack cổ điển, ta sẽ viết hàm

btr(i) với
0i8
, tức là đã có
i
quân Hậu, bây giờ ta sẽ thử các phương án đặt con Hậu thứ
i+1
vào dòng thứ
i+1
rồi gọi đệ quy tiếp
btr(i+1)
. Nếu tất cả các phương án đều không thành công thì quay lui lại
btr(i1)
và cứ tiếp tục thử. Đến khi gọi
btr(8)
thành công nghĩa là ta đã đặt đủ
8
quân Hậu và tìm được
1
đáp án hợp lệ.

Người ta đã tìm ra cách tối ưu thuật toán này như sau: thay vì lần lượt thử các phương án để đặt quân Hậu vào dòng thứ

i+1, ta chọn ngẫu nhiên một phương án để thử, nếu không thành công, loại phương án này rồi chọn ngẫu nhiên và thử tiếp. Trên thực tế, cách này có thể tìm ra đáp án hợp lệ nhanh hơn thuật toán backtrack thông thường.