Randomized Algorithms

Lấy mẫu

Nếu đáp án được xác định bởi một (hay hai) phần tử, và các phần tử “tốt” chiếm ít nhất tỉ lệ $q$, thì chọn ngẫu nhiên một phần tử sẽ trúng với xác suất $q$. Thử $T$ lần độc lập, xác suất trượt cả $T$ lần là $(1 - q)^T \le e^{-qT}$. Mỗi ứng viên được kiểm lại chính xác nên không bao giờ in sai “có”; chỉ có thể sót.

Ví dụ: đường thẳng chứa ít nhất 20% số điểm. Một cặp điểm ngẫu nhiên cùng nằm trên đường ấy với xác suất khoảng $0.2^2 = 0.04$; thử 300 cặp thì xác suất trượt là $0.96^{300} \approx 5 \cdot 10^{-6}$.

Với nhiều truy vấn, cộng các xác suất sai (chặn hợp): $q$ truy vấn, mỗi truy vấn sai với xác suất $\varepsilon$, tổng không quá $q\varepsilon$. Chọn $T$ sao cho con số này rất nhỏ, ví dụ $10^{-9}$.

Dùng mt19937 (hoặc mt19937_64) với hạt giống chrono::steady_clock::now(). Đừng dùng rand(): trên một số trình biên dịch nó chỉ cho tối đa 32767, và hạt giống cố định dễ bị dựng test phá.

Hình dưới thử từng cặp điểm: mỗi lần lấy hai điểm ngẫu nhiên rồi đếm chính xác số điểm trên đường qua chúng.

Ví dụĐường Thẳng Đủ Đông

Đề. $n \le 10^5$ điểm phân biệt, tọa độ tới $10^9$, và $p \ge 20$. Có đường thẳng nào đi qua ít nhất $\lceil np/100 \rceil$ điểm không?

Lời giải. Nếu có đường đi qua ít nhất 20% số điểm, một cặp điểm ngẫu nhiên nằm trên nó với xác suất khoảng $0{,}2^2 = 0{,}04$. Thử 300 cặp, mỗi cặp đếm số điểm thẳng hàng; cả 300 lần đều trượt với xác suất $0{,}96^{300} \approx 5 \cdot 10^{-6}$.

Bước 1: trường hợp nhỏ. need = (n * p + 99) / 100 là $\lceil np/100 \rceil$ tính bằng số nguyên. Cần không quá 2 điểm thì luôn có đường (qua hai điểm bất kỳ, hoặc qua điểm duy nhất).

ll need = (n * p + 99) / 100;
if (need <= 2) {
    cout << "possible\n";
    return 0;
}

Bước 2: thử một cặp. Đếm mọi điểm $k$ thẳng hàng với $i$, $j$ bằng cách so hai tích chéo thay cho chia lấy hệ số góc. Hiệu tọa độ tới $2 \cdot 10^9$ nên mỗi tích tới $4 \cdot 10^{18}$, vẫn vừa long long.

ll i = rng() % n, j = rng() % n;
if (i == j) continue;
ll cnt = 0;
for (ll k = 0; k < n; k++)
    if ((x[j] - x[i]) * (y[k] - y[i]) == (y[j] - y[i]) * (x[k] - x[i])) cnt++;

mt19937 gieo bằng đồng hồ, để test cố định không thể nhắm vào một dãy số ngẫu nhiên biết trước.

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

const ll N = 1e5 + 5;
ll n, p, x[N], y[N];
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

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

    cin >> n >> p;
    for (ll i = 0; i < n; i++) cin >> x[i] >> y[i];
    ll need = (n * p + 99) / 100;
    if (need <= 2) {
        cout << "possible\n";
        return 0;
    }
    for (ll it = 0; it < 300; it++) {
        ll i = rng() % n, j = rng() % n;
        if (i == j) continue;
        ll cnt = 0;
        for (ll k = 0; k < n; k++)
            if ((x[j] - x[i]) * (y[k] - y[i]) == (y[j] - y[i]) * (x[k] - x[i])) cnt++;
        if (cnt >= need) {
            cout << "possible\n";
            return 0;
        }
    }
    cout << "impossible\n";
}

ChạyInput6
50
0 1
2 5
4 9
0 0
1 0
5 2

Outputpossible

Ước lượng phép tính: $300n = 3 \cdot 10^7$ phép.

Ví dụGiá Trị Nổi Bật Trên Đoạn

Đề. Mảng $n \le 2 \cdot 10^5$, $q \le 2 \cdot 10^5$ truy vấn $[l, r]$: in tăng dần mọi giá trị xuất hiện nhiều hơn $\lfloor len/3 \rfloor$ lần, hoặc $-1$.

Lời giải. Có nhiều nhất hai giá trị như vậy, mỗi giá trị chiếm hơn $1/3$ đoạn. Lấy 60 vị trí ngẫu nhiên làm ứng viên: mỗi giá trị đáp án bị sót với xác suất dưới $(2/3)^{60} \approx 3 \cdot 10^{-11}$. Đếm chính xác từng ứng viên, nên không bao giờ in sai, chỉ có thể sót với xác suất rất nhỏ đó.

Bước 1: nén giá trị, danh sách vị trí. Mỗi bộ test nén lại và dựng pos[v], các vị trí của giá trị $v$ theo thứ tự tăng.

for (int i = 1; i <= n; i++) {
    id[i] = lower_bound(vals.begin(), vals.end(), a[i]) - vals.begin();
    pos[id[i]].push_back(i);
}

Bước 2: chọn ứng viên. Đoạn không dài hơn 60 thì lấy hết, khỏi phụ thuộc may rủi; dài hơn thì lấy 60 vị trí ngẫu nhiên.

if (len <= S)
    for (int i = l; i <= r; i++) c.push_back(id[i]);
else
    for (int k = 0; k < S; k++) c.push_back(id[l + rng() % len]);

Bước 3: đếm chính xác. Số lần xuất hiện trong $[l, r]$ là upper_bound(r) - lower_bound(l) trên pos[v]. Ứng viên đã sắp theo chỉ số nén, mà chỉ số nén theo giá trị tăng, nên ans ra đúng thứ tự cần in.

for (int v : c) {
    int cnt = upper_bound(pos[v].begin(), pos[v].end(), r) - lower_bound(pos[v].begin(), pos[v].end(), l);
    if (cnt > len / 3) ans.push_back(vals[v]);
}

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

const int N = 2e5 + 5, S = 60;
int t, n, q, a[N], id[N];
vector<int> vals, pos[N];
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

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

    cin >> t;
    while (t--) {
        cin >> n >> q;
        vals.clear();
        for (int i = 1; i <= n; i++) cin >> a[i], vals.push_back(a[i]);
        sort(vals.begin(), vals.end());
        vals.erase(unique(vals.begin(), vals.end()), vals.end());
        for (int v = 0; v < (int)vals.size(); v++) pos[v].clear();
        for (int i = 1; i <= n; i++) {
            id[i] = lower_bound(vals.begin(), vals.end(), a[i]) - vals.begin();
            pos[id[i]].push_back(i);
        }
        while (q--) {
            int l, r;
            cin >> l >> r;
            int len = r - l + 1;
            vector<int> c;
            if (len <= S)
                for (int i = l; i <= r; i++) c.push_back(id[i]);
            else
                for (int k = 0; k < S; k++) c.push_back(id[l + rng() % len]);
            sort(c.begin(), c.end());
            c.erase(unique(c.begin(), c.end()), c.end());
            vector<int> ans;
            for (int v : c) {
                int cnt = upper_bound(pos[v].begin(), pos[v].end(), r) - lower_bound(pos[v].begin(), pos[v].end(), l);
                if (cnt > len / 3) ans.push_back(vals[v]);
            }
            if (ans.empty()) cout << -1;
            for (int x : ans) cout << x << ' ';
            cout << '\n';
        }
    }
}

ChạyInput1
9 5
4 7 4 9 7 4 7 8 10
1 9
1 7
2 6
4 9
8 8

Output-1
4 7
4 7
-1
8

Độ phức tạp: mỗi truy vấn $O(60 \log n)$.

Freivalds

Kiểm $AB = C$ bằng cách nhân thật tốn $O(n^3)$. Freivalds chọn vector ngẫu nhiên $v$ rồi so $A(Bv)$ với $Cv$: ba phép nhân ma trận với vector, $O(n^2)$.

Nếu $AB = C$ thì hai vế luôn bằng nhau. Nếu $D = AB - C \ne 0$, có hàng $d$ của $D$ khác 0, và $d \cdot v = 0$ là một phương trình tuyến tính không tầm thường theo $v$. Với $v$ chọn đều trong $\mathbb{Z}_M^n$, xác suất thỏa chỉ là $1/M$.

Cùng ý tưởng: so hai đa thức bằng giá trị tại một điểm ngẫu nhiên (đa thức bậc $d$ khác 0 có không quá $d$ nghiệm), so hai xâu bằng băm.

Ví dụKiểm Tra Phép Nhân Ma Trận

Đề. Ba ma trận $n \times n$, $n \le 1000$, phần tử dưới $M = 10^9 + 7$. Kiểm $AB \equiv C \pmod M$.

Lời giải. Tính $AB$ tốn $O(n^3) = 10^9$ phép. Freivalds: chọn vector ngẫu nhiên $v$ và so $A(Bv)$ với $Cv$, chỉ cần ba phép nhân ma trận với vector. Nếu $AB \ne C$ thì một lần thử bị lừa với xác suất không quá $1/M$; hai lần thử độc lập là $1/M^2$.

Bước 1: nhân ma trận với vector. $O(n^2)$, lấy dư sau mỗi phép cộng; tích hai số dưới $M$ nhỏ hơn $1{,}1 \cdot 10^{18}$, vừa long long.

void mulvec(ll M[N][N], ll *in, ll *out) {
    for (ll i = 0; i < n; i++) {
        ll s = 0;
        for (ll j = 0; j < n; j++) s = (s + M[i][j] * in[j]) % MOD;
        out[i] = s;
    }
}

Bước 2: một lần thử. bv $= Bv$, abv $= A(Bv)$, cv $= Cv$.

for (ll i = 0; i < n; i++) v[i] = rng() % MOD;
mulvec(B, v, bv);
mulvec(A, bv, abv);
mulvec(C, v, cv);

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

const ll N = 1005, MOD = 1e9 + 7;
ll n, A[N][N], B[N][N], C[N][N], v[N], bv[N], abv[N], cv[N];
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

void mulvec(ll M[N][N], ll *in, ll *out) {
    for (ll i = 0; i < n; i++) {
        ll s = 0;
        for (ll j = 0; j < n; j++) s = (s + M[i][j] * in[j]) % MOD;
        out[i] = s;
    }
}

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

    cin >> n;
    for (ll i = 0; i < n; i++)
        for (ll j = 0; j < n; j++) cin >> A[i][j];
    for (ll i = 0; i < n; i++)
        for (ll j = 0; j < n; j++) cin >> B[i][j];
    for (ll i = 0; i < n; i++)
        for (ll j = 0; j < n; j++) cin >> C[i][j];
    for (ll it = 0; it < 2; it++) {
        for (ll i = 0; i < n; i++) v[i] = rng() % MOD;
        mulvec(B, v, bv);
        mulvec(A, bv, abv);
        mulvec(C, v, cv);
        for (ll i = 0; i < n; i++)
            if (abv[i] != cv[i]) {
                cout << "NO\n";
                return 0;
            }
    }
    cout << "YES\n";
}

ChạyInput2
1 0
0 1
2 3
4 5
2 3
4 5

OutputYES

Ước lượng phép tính: hai lần thử, mỗi lần $3n^2 = 3 \cdot 10^6$ phép.

Miller-Rabin

Với $n$ lẻ, viết $n - 1 = 2^s d$ với $d$ lẻ. Nếu $n$ nguyên tố thì với mọi $a$ không chia hết cho $n$, dãy $a^d, a^{2d}, \dots, a^{2^s d}$ theo mô-đun $n$ hoặc bắt đầu bằng 1, hoặc chứa $-1$ trước khi tới 1 (vì $x^2 \equiv 1$ chỉ có nghiệm $\pm 1$ theo mô-đun nguyên tố). Cơ số $a$ vi phạm điều này là nhân chứng cho $n$ là hợp số. Với hợp số, ít nhất $3/4$ số $a$ là nhân chứng, nên vài cơ số ngẫu nhiên đã đủ; với $n < 2^{64}$, 12 cơ số nguyên tố đầu tiên $2, 3, \dots, 37$ đã được kiểm là đủ chắc chắn.

Phép nhân theo mô-đun $n$ tới $2^{63}$ cần __int128. Mỗi cơ số tốn $O(\log n)$ phép nhân.

Ví dụNhận Diện Số Nguyên Tố

Đề. $t \le 500$ số $N \in [2, 2^{63} - 1]$, in YES nếu nguyên tố.

Lời giải. Miller-Rabin với 12 số nguyên tố đầu làm cơ số. Bộ cơ số này đã được kiểm là đúng cho mọi $n < 2^{64}$, nên kết quả chắc chắn đúng chứ không chỉ đúng với xác suất cao.

Bước 1: nhân theo mô-đun lớn. Hai số dưới $2^{63}$ có tích tới $2^{126}$, nên nhân trong __int128.

ll mul(ll a, ll b, ll m) {
    return (__int128)a * b % m;
}

Bước 2: chia thử. Gặp ước nhỏ thì $n$ nguyên tố chỉ khi chính nó là số đó. Bước này cũng loại mọi số chẵn, nên phần sau luôn có $n$ lẻ. Rồi tách $n - 1 = 2^s d$ với $d$ lẻ.

for (ll p : bases)
    if (n % p == 0) return n == p;
ll d = n - 1, s = 0;
while (d % 2 == 0) d /= 2, s++;

Bước 3: kiểm một cơ số. $x = a^d$. Bằng 1 hoặc $n - 1$ thì $a$ không phải nhân chứng. Ngược lại bình phương tối đa $s - 1$ lần; gặp $n - 1$ thì thôi. Không gặp thì $a$ là nhân chứng, $n$ là hợp số.

ll x = power(a, d, n);
if (x == 1 || x == n - 1) continue;
bool comp = 1;
for (ll r = 1; r < s && comp; r++) {
    x = mul(x, x, n);
    if (x == n - 1) comp = 0;
}
if (comp) return 0;

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

ll bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

ll mul(ll a, ll b, ll m) {
    return (__int128)a * b % m;
}

ll power(ll a, ll b, ll m) {
    ll r = 1;
    a %= m;
    while (b) {
        if (b & 1) r = mul(r, a, m);
        a = mul(a, a, m);
        b >>= 1;
    }
    return r;
}

bool isPrime(ll n) {
    if (n < 2) return 0;
    for (ll p : bases)
        if (n % p == 0) return n == p;
    ll d = n - 1, s = 0;
    while (d % 2 == 0) d /= 2, s++;
    for (ll a : bases) {
        ll x = power(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool comp = 1;
        for (ll r = 1; r < s && comp; r++) {
            x = mul(x, x, n);
            if (x == n - 1) comp = 0;
        }
        if (comp) return 0;
    }
    return 1;
}

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

    ll t;
    cin >> t;
    while (t--) {
        ll n;
        cin >> n;
        cout << (isPrime(n) ? "YES" : "NO") << '\n';
    }
}

ChạyInput6
2
29
91
97
221
1000000

OutputYES
YES
NO
YES
NO
NO

Mỗi số $12 \cdot 63$ phép nhân __int128. Số nguyên tố là trường hợp tốn nhất vì phải qua đủ mọi cơ số, nên số cơ số quyết định tốc độ:

Bộ 7 cơ số $2, 325, 9375, 28178, 450775, 9780504, 1795265022$ cũng đã được kiểm là đủ cho mọi số dưới $2^{64}$, và nhanh hơn bộ 12 cơ số khoảng 1,6 lần. Phải kiểm hàng trăm nghìn số lớn thì dùng bộ 7; nhớ lấy $a \bmod n$ và bỏ qua cơ số bằng 0.

Băm trọng số

Muốn biết hai đoạn có cùng đa tập giá trị không, gán mỗi giá trị $v$ một số ngẫu nhiên 64 bit $h_v$ và lấy $H = \sum h_{a_i}$ trên đoạn. Tràn số của unsigned long long chính là cộng theo mô-đun $2^{64}$.

Hai đa tập bằng nhau thì $H$ bằng nhau. Khác nhau thì $H$ bằng nhau với xác suất rất nhỏ, vì hiệu là tổ hợp tuyến tính khác 0 của các $h_v$ độc lập. Tổng tiền tố cho tính $H$ của mọi đoạn trong $O(1)$.

Biến thể:

• XOR trọng số để kiểm “mỗi giá trị xuất hiện chẵn lần”.
• Trọng số theo mô-đun 3 để kiểm số mũ chia hết cho 3.
• Cộng trọng số vào cả hình chữ nhật để đánh dấu vùng.

Ví dụHai Đoạn Cùng Đa Tập

Đề. Hai dãy $A, B$ dài $N \le 2 \cdot 10^5$, giá trị trong $[1, N]$. $Q \le 2 \cdot 10^5$ truy vấn: đoạn $A[l..r]$ có hoán vị được thành $B[L..R]$ không?

Lời giải. Gán cho mỗi giá trị một số ngẫu nhiên 64 bit h[v]. Hai đa tập bằng nhau thì tổng trọng số bằng nhau; khác nhau thì tổng trùng nhau chỉ với xác suất cỡ $2^{-64}$. Tổng trên một đoạn lấy từ tổng tiền tố.

Bước 1: tổng tiền tố trọng số. Kiểu ull nên tràn số chính là cộng theo mô-đun $2^{64}$, không cần lấy dư.

for (ll v = 1; v <= n; v++) h[v] = rng();
for (ll i = 1; i <= n; i++) {
    ll x;
    cin >> x;
    pa[i] = pa[i - 1] + h[x];
}

Bước 2: so. So độ dài và so hiệu hai tổng tiền tố, $O(1)$. So độ dài là cần thiết: hai đoạn khác độ dài không bao giờ hoán vị được thành nhau, dù tổng trọng số có tình cờ bằng nhau.

bool ok = r - l == R - L && pa[r] - pa[l - 1] == pb[R] - pb[L - 1];

#include <bits/stdc++.h>
using namespace std;
#define ll long long
using ull = unsigned long long;

const ll N = 2e5 + 5;
ll n, q;
ull h[N], pa[N], pb[N];
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

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

    cin >> n >> q;
    for (ll v = 1; v <= n; v++) h[v] = rng();
    for (ll i = 1; i <= n; i++) {
        ll x;
        cin >> x;
        pa[i] = pa[i - 1] + h[x];
    }
    for (ll i = 1; i <= n; i++) {
        ll x;
        cin >> x;
        pb[i] = pb[i - 1] + h[x];
    }
    while (q--) {
        ll l, r, L, R;
        cin >> l >> r >> L >> R;
        bool ok = r - l == R - L && pa[r] - pa[l - 1] == pb[R] - pb[L - 1];
        cout << (ok ? "Yes" : "No") << '\n';
    }
}

ChạyInput8 6
1 2 2 3 4 1 5 3
2 1 3 2 4 5 1 3
1 4 1 4
2 5 2 5
5 7 5 7
1 3 6 8
6 8 6 8
1 2 1 3

OutputYes
No
Yes
No
Yes
No

Độ phức tạp: $O(N + Q)$.