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

*Author: [Gã coder si tình](https://www.facebook.com/profile.php?id=100086269423732)*

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](https://www.calculatorsoup.com/calculators/statistics/random-number-generator.php).

**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](https://www.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]^](https://en.cppreference.com/w/cpp/numeric/random/rand). 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à $2^{15} - 1 = 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]^](https://codeforces.com/blog/entry/61587). 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; 10^9]$ 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]^](https://en.wikipedia.org/wiki/Linear_congruential_generator) 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]^](https://en.wikipedia.org/wiki/Mersenne_Twister) -- 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 $2^{19937} - 1$. 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_clock` và `high_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 $[-2^{31}; 2^{31} - 1]$. 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 $10^5$ số trong khoảng $[0; 10^9]$ 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 \space 045,548$
- Khoảng tứ phân vị (interquartile range): $250 \space 085,5$ // $500 \space 415,5$ // $750 \space 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 $\Omega = \{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]^](https://leetcode.com/problems/random-pick-with-weight/?fbclid=IwAR3qXcR34m1wA_gtAXMqLyhcyEi5ax9ZQbyvVbWVYGr3Pp-nk3-VehR3vjI) 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]^](https://www.facebook.com/code.cung.rr/posts/pfbid02XYYWRdteDMZrmUkD39S5Lj5fGML39UALvn8XEd8TzCnD4cDxSQkN8EJh3VHP1Ntvl).

### 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 $0 \leq i < n$) sẽ có xác suất xuất hiện là $\frac{p_i}{\sum{p}}$.

Constraints:
- $1 \leq n \leq 10^4$.
- $1 \leq p_i \leq 10^5$.
- Hàm ngẫu nhiên sẽ được gọi không quá $10^4$ lần.

Ví dụ cho mảng $p = [10, 25, 15, 50]$, $\sum{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 $\sum{p}$ "slots" cho các giá trị, rồi cho giá trị thứ $i$ chiếm $p_i$ slots. Với ví dụ ở trên thì giá trị $0$ sẽ chiếm các slots từ $1 \rightarrow 10$, giá trị $1$ chiếm các slots từ $11 \rightarrow 35$, giá trị $2$ chiếm các slots từ $36 \rightarrow 50$,...

Sau đó, ta sẽ sinh ra một số ngẫu nhiên trong khoảng $[1; \sum{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(\log n)$.

### 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]^](https://www.facebook.com/permalink.php?story_fbid=pfbid0RdJGcLwMTZVvKZBmQutuGdBQGxfWECqCt9vuSH8kyrUWwqs1ETF5jb2Y3VkbzL5ul&id=100071054620272).

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$ và $b$ gồm các **phần tử phân biệt** và có cùng độ dài $m$ 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:

$$ h_{a_1} \oplus h_{a_2} \oplus ... \oplus h_{a_m} = h_{b_1} \oplus h_{b_2} \oplus ... \oplus h_{b_m} $$

Với $\oplus$ 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$ và $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ệ $\frac{1}{2^k}$, 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]^](https://oj.vnoi.info/problem/bedao_g10_perfect) (Bedao Grand Contest 10).

#### Tóm tắt đề

Cho hoán vị $p$ độ dài $n$ và $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 $p_l, p_{l+1}, p_{l + 2},..., p_r$ 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:
- $1 \leq n, q \leq 10^5$.
- $1 \leq p_i \leq n$.

#### Ý tưởng

Nhiệm vụ của chúng ta là so sánh dãy $p_l, p_{l + 1}, p_{l + 2},..., p_r$ với dãy $v, v + 1, v + 2,..., v + sz - 1$ với $v$ 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:

$$ h_{p_l} \oplus h_{p_{l+1}} \oplus h_{p_{l+2}} \oplus ... \oplus h_{p_r} $$

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

$$ h_v \oplus h_{v + 1} \oplus h_{v + 2} \oplus ... \oplus h_{v + sz - 1}$$

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 \times 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 \times 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 \times m$ cho ma trận $B$ có kích thước $m \times q$, ta sẽ được $1$ ma trận $C$ có kích thước $n \times q$ sao cho:

$$ C_{i, j} = \sum_{k = 1}^{m}{A_{i, k} \times B_{k, j}}$$

Để tính toàn bộ ma trận $C$ ta cần $O(n \times m \times 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 \times B) \times C = A \times (B \times 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 \times n$. Nhiệm vụ là kiểm tra xem $A \times B$ có bằng $C$ hay không. Dĩ nhiên, ta sẽ không tính trực tiếp $A \times B$ trong $O(n^3)$.

Thuật toán này sẽ sinh ngẫu nhiên một ma trận $r$ có kích thước $1 \times 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 \times A \times B = r \times C \space \space (*) $$

Để 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(n^2)$. Hai ma trận nhận được sau khi tính toán có kích thước $1 \times 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 \times B = C \rightarrow$ chắc chắn $(*)$ đúng.
- Trường hợp $A \times B \neq C \rightarrow$ $(*)$ 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à $\leq 1/2$ [^[9]^](https://en.wikipedia.org/wiki/Freivalds%27_algorithm#Error_analysis) -- 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 \times 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 \times B \neq C$.

Lúc này, vẫn sẽ có trường hợp xảy ra sai số, là khi $A \times B \neq C$ 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à $\frac{1}{2^t}$. 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:

$$ \frac{1}{2^{40}} \approx 0,00000000009095\% $$

Như vậy, ta đã có một thuật toán chạy trong $O(t \times n^2)$ thay vì $O(n^3)$ 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]^](https://tinhoctre.vn/problem/tht2022_matran). 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;
}

### Miller--Rabin primality test

**Miller--Rabin 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(\sqrt{n})$, hoặc $O(1)$ sau khi chạy sàng nguyên tố trong $O(n \log \log n)$. Với thuật toán của Miller--Rabin thì ta chỉ tốn độ phức tạp $O(t \times \log n)$ 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, Miller--Rabin 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 Miller--Rabin thì mình sẽ nói một chút về định lý Fermat nhỏ [^[11]^](https://en.wikipedia.org/wiki/Fermat%27s_little_theorem) 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]^](https://en.wikipedia.org/wiki/Fermat_primality_test).

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

$$ a^{p-1} \equiv 1 \bmod p \space \space (**)$$

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à $a^{p-1} \not\equiv 1 \bmod p$, 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$ mà $(**)$ vẫn thỏa. Khi đó ta gọi $a$ là *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]^](https://en.wikipedia.org/wiki/Carmichael_number) 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à $a^{n-1} \equiv 1 \bmod n$ 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 \space 081 = 7 \times 13 \times 31 \times 61$. Các bạn có thể xem thêm dãy số Carmichael trên OEIS có mã A002997 [^[14]^](https://oeis.org/A002997).

Để 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á $10^9$, và khoảng $20 \space 138 \space 200$ số Carmichael không quá $10^{21}$. 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 Miller--Rabin primality test

Miller--Rabin primality test cũng được phát triển dựa trên định lý Fermat nhỏ $a^{p-1} \equiv 1 \bmod p$ ($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, Miller--Rabin đư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ố) $\Rightarrow p - 1$ chẵn. Từ đó, ta có thể biểu diễn $p-1$ dưới dạng $2^s \times d$ (với $s$ và $d$ là số nguyên $\geq 1$). Nếu biểu diễn $p - 1$ dưới dạng nhị phân thì $s$ là số bit $0$ liên tiếp bên phải $p - 1$ và $d$ là $p - 1$ sau khi right shift $s$ lần.

Thay $p - 1 = 2^s \times d$ và chuyển số $1$ sang vế trái:

$$ a^{2^sd} - 1 \equiv 0 \bmod p$$

Ta thấy $a^{2^sd} - 1$ có dạng hằng đẳng thức $x^2 - y^2$ với $x = a^{2^{s - 1}d}$ và $y = 1$. Sau khi phân tích xong, ta lại áp dụng hằng đẳng thức vào $a^{2^{s-1}d} - 1$, cứ như vậy, ta có:

$$ (a^{2^{s-1}d} + 1)(a^{2^{s-1}d} - 1) \equiv 0 \bmod p $$

$$ \Longleftrightarrow (a^{2^{s-1}d} + 1)(a^{2^{s-2}d} + 1)(a^{2^{s-2}d} - 1) \equiv 0 \bmod p $$

$$ \Longleftrightarrow (a^{2^{s-1}d} + 1)(a^{2^{s-2}d} + 1)(a^{2^{s-3}d} + 1)(a^{2^{s-3}d} - 1) \equiv 0 \bmod p $$

$$ \vdots $$

$$ \Longleftrightarrow (a^{2^{s-1}d} + 1)(a^{2^{s-2}d} + 1)(a^{2^{s-3}d} + 1) \cdots (a^{2^0d} + 1)(a^{2^0d} -1) \equiv 0 \bmod p \space \space (***) $$

Cộng với một số phân tích toán học, Miller--Rabin 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:

- $a^d \equiv 1 \bmod p$.
- Tồn tại $1$ số nguyên $r$ ($0 \leq r < s$) sao cho $a^{2^rd} \equiv -1 \bmod p$, hay $a^{2^rd} \equiv p - 1 \bmod p$ (do $-1 \equiv p - 1 \bmod p$).

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 $\Rightarrow (***)$ đúng $\Leftrightarrow (**)$ đú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 Miller--Rabin 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$ là *strong liar*.

Khác với Fermat, thuật toán của Miller--Rabin 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 ($2 \leq a \leq p - 2$), khả năng $a$ là một *strong liar* là $\leq 1/4$ [^[15]^](https://www.sciencedirect.com/science/article/pii/0022314X80900840?via%3Dihub). Đâ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 Miller--Rabin 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$ là *strong liar* -- giảm xuống còn $\frac{1}{4^t}$. Ta chỉ cần chọn $t \geq 20$ 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 Miller--Rabin 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 Miller--Rabin, 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ị $a^d \bmod p$ để 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ị $a^{2^rd} \bmod p$, để tăng $r$ lên, ta chỉ việc bình phương biến `val` lên là xong, vì $a^{2^rd} \times a^{2^rd} = a^{2^rd + 2^rd} = a^{2^{r+1}d}$.

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$ và $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, $s \leq \log n$ nên Miller--Rabin primality test có độ phức tạp thời gian $O(t \times \log n)$ 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]^](https://www.spoj.com/problems/PON/).

## 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]^](https://dl.acm.org/doi/pdf/10.1145/366622.366644), 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ừ $0 \rightarrow n - 1$. Ta có thể giả sử $0 \leq k < 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(n \log n)$ 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ử $\leq x$, dãy thứ hai gồm những phần tử $> x$. Gọi $a$ và $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 $k \geq a$, ta tìm phần tử lớn thứ $k - a$ 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 + \cdots + n/2^j < 2n = O(n) $$

Độ phức tạp tệ nhất của QuickSelect là $O(n^2)$, 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>` có $2 \times 10^6$ 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]^](https://en.wikipedia.org/wiki/Eight_queens_puzzle) 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 \times 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 $0 \leq i \leq 8$, 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(i - 1)$ 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.