Bài 40: Số học nâng cao trong C++ và Python: sàng ước nguyên tố nhỏ nhất, sàng trên đoạn, hàm phi Euler, Euclid mở rộng, định lý thặng dư Trung Hoa

Series: Thuật toán chuyên sâu Độ khó: Khá Bài 40 / Thuật toán chuyên sâu Sàng, phi Euler, Trung Hoa

// Thuật toán chuyên sâu · Bài 40

Số học nâng cao: sàng spf, sàng đoạn, phi Euler và định lý Trung Hoa

Phân tích một triệu số trong nháy mắt, đếm số nguyên tố quanh $10^{12}$, tính lũy thừa với số mũ khổng lồ, giải hệ “chia 3 dư 2, chia 5 dư 3, chia 7 dư 2”: bốn công cụ số học mà đề HSG vòng 2 và HSG quốc gia rất hay dùng.

Lời mở đầu

Ở bài 12 (số học cơ bản) em đã biết liệt kê ước trong $O(\sqrt n)$, thuật toán Euclid, sàng Eratosthenes và phân tích thừa số bằng cách thử chia. Ở bài 26 (modulo và tổ hợp) em đã có lũy thừa nhanh, định lý Fermat nhỏ và Euclid mở rộng để tìm nghịch đảo. Những công cụ đó đủ cho đề cấp trường, nhưng gặp đề vòng 2 thì thường “hụt hơi”: phải phân tích $10^6$ số chứ không phải một số, số nguyên tố cần tìm nằm quanh $10^{12}$ nên không sàng nổi cả mảng, số mũ là $b^c$ với $b, c$ tới $10^9$…

Bài này nâng cấp từng món: sàng ước nguyên tố nhỏ nhất để phân tích mỗi số trong $O(\log)$, sàng trên đoạn $[L, R]$, hàm phi Euler và định lý Euler (tổng quát hoá Fermat cho modulo bất kỳ), phương trình $ax + by = c$ và cuối cùng là định lý thặng dư Trung Hoa, bài toán cổ “Hàn Tín điểm binh” mà máy tính giải trong vài micro giây. Em cần nắm chắc độ phức tạp và tràn số (bài 9) vì phần cuối bài số rất lớn, chỉ sơ ý là tràn long long.

Bài 33Cây, QHĐ trên câyBài 34LCA, nhảy nhị phânBài 35Sắp xếp topo, DAGBài 36TPLT mạnh, khớp, cầuBài 37Sparse Table, LazyBài 38Chia căn, thuật toán MoBài 39Z, Manacher, TrieBài 40Số học nâng caoBài 41Kỹ thuật tìm kiếmBài 42Trò chơi và tổng ônBài 43Ghép cặp, luồngBài 44Small-to-large, Euler tourBài 45Tối ưu QHĐBài 46Hash nâng cao, Suffix ArrayBài 47Heavy-Light Decomposition

Trong bài này bạn sẽ học: sàng ước nguyên tố nhỏ nhất (spf) và phân tích thừa số mỗi số trong $O(\log x)$; công thức số ước, tổng ước từ phân tích; sàng nguyên tố trên đoạn $[L, R]$ với $R$ tới $10^{12}$; hàm phi Euler: công thức, tính cho một số, sàng phi cho cả dãy, định lý Euler và cách hạ số mũ khổng lồ; Euclid mở rộng dạng vòng lặp, giải $ax + by = c$ và họ nghiệm tổng quát; định lý thặng dư Trung Hoa cho 2 rồi $k$ phương trình (kể cả modulo không nguyên tố cùng nhau) và cách chống tràn số. Có khung chạy tay sàng spf và Euclid mở rộng, 3 mô phỏng từng bước và 8 bài tập kèm gợi ý.

Lộ trình bài viết, bấm vào từng phần để nhảy tới:

Phần 1Sàng spf, số ướcPhần 2Sàng trên đoạnPhần 3Hàm phi EulerPhần 4Phương trình ax + by = cPhần 5Định lý Trung HoaBài tập8 bài + gợi ý

Phần 1Sàng ước nguyên tố nhỏ nhất và phân tích nhanh

Ở bài 12, muốn phân tích $x$ ra thừa số nguyên tố ta thử chia cho $2, 3, 4, \ldots$ tới $\sqrt x$. Một số thì nhanh, nhưng nếu đề cho $Q = 10^6$ số, mỗi số tới $10^7$, thì tổng cộng cỡ $10^6 \cdot 3 \cdot 10^3 = 3 \cdot 10^9$ phép chia: quá chậm.

Ý tưởng: sửa sàng Eratosthenes một chút. Thay vì chỉ đánh dấu “là hợp số”, ta ghi lại số nguyên tố đầu tiên gạch nó. Gọi spf[v] (smallest prime factor) là ước nguyên tố nhỏ nhất của $v$:

• Duyệt $i = 2, 3, \ldots, N$. Nếu spf[i] == 0 thì chưa số nguyên tố nào nhỏ hơn chia hết $i$, vậy $i$ là số nguyên tố, gán spf[i] = i.
• Với $i$ nguyên tố, duyệt các bội $j = i^2, i^2 + i, \ldots$; ô nào còn 0 thì gán spf[j] = i. Ô đã có giá trị thì giữ nguyên vì đã có ước nguyên tố nhỏ hơn.

Vì sao bắt đầu từ $i^2$? Bội $k \cdot i$ với $k < i$ có ước nguyên tố $\le k < i$, nên đã được số nhỏ hơn gán rồi. Sàng vẫn tốn $O(N \log \log N)$ như sàng thường.

Phân tích $x$: p = spf[x] là một ước nguyên tố của $x$; chia $x$ cho $p$ tới khi hết, rồi lặp lại với phần còn lại. Mỗi lần chia, $x$ giảm ít nhất một nửa, nên chỉ cần không quá $\log_2 x$ lần chia (với $x \le 10^7$ là tối đa 23 lần).

x	spf[x]	chia	thừa số
360	2	360 → 180 → 90 → 45	$2^3$
45	3	45 → 15 → 5	$3^2$
5	5	5 → 1	$5^1$
1	–	dừng	$360 = 2^3 \cdot 3^2 \cdot 5$

Bảng 1. Phân tích 360 bằng mảng spf: mỗi hàng chỉ cần tra một ô của mảng, không phải thử chia từ 2.

Bấm Tiếp để xem mảng spf được lấp dần: mỗi số nguyên tố $p \le \sqrt N$ “nhận” những bội chưa có chủ, rồi xem số $x$ được phân tích chỉ bằng cách tra bảng. Em thử $x = 64$ (chia 6 lần) và $x = 59$ (số nguyên tố, chỉ một bước):

MÔ PHỎNGSàng ước nguyên tố nhỏ nhất rồi phân tích x trong vài bướcNhập N (10 … 120) và số x ≤ N. Số nhỏ dưới mỗi ô là spf của ô đó.

Nx

ChậmVừaNhanh

số nguyên tố đang xétvừa gán spf / đang chiasố nguyên tốhợp số đã có spf

Số ước và tổng các ước

Có phân tích $x = p_1^{e_1} p_2^{e_2} \cdots p_k^{e_k}$ thì mọi ước của $x$ có dạng $p_1^{f_1} \cdots p_k^{f_k}$ với $0 \le f_i \le e_i$. Mỗi $f_i$ có $e_i + 1$ cách chọn, nên:

$d(x) = (e_1 + 1)(e_2 + 1) \cdots (e_k + 1)$, $\quad \sigma(x) = \prod_{i=1}^{k} (1 + p_i + p_i^2 + \cdots + p_i^{e_i}) = \prod_{i=1}^{k} \dfrac{p_i^{e_i + 1} - 1}{p_i - 1}$.

Với $360 = 2^3 \cdot 3^2 \cdot 5$: $d = 4 \cdot 3 \cdot 2 = 24$ ước, $\sigma = 15 \cdot 13 \cdot 6 = 1170$. Công thức tổng ước có được vì khi khai triển tích $(1 + 2 + 4 + 8)(1 + 3 + 9)(1 + 5)$, mỗi số hạng là đúng một ước của 360.

Xem code C++ / Python: sàng spf, phân tích nhiều số, số ước và tổng ước (bấm để mở)

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

const int N = 10000000;
int spf[N + 1];                 // spf[v] = ước nguyên tố nhỏ nhất của v

int main(){
    ios::sync_with_stdio(false); cin.tie(nullptr);
    for(int i = 2; i <= N; i++){
        if(spf[i] != 0) continue;            // hợp số: đã có spf
        spf[i] = i;                          // i là số nguyên tố
        for(long long j = 1LL * i * i; j <= N; j += i)
            if(spf[j] == 0) spf[j] = i;
    }
    int q;
    cin >> q;
    while(q--){
        int x;
        cin >> x;
        cout << x << " =";
        if(x == 1) cout << " 1";
        long long d = 1, s = 1;
        bool first = true;
        while(x > 1){
            int p = spf[x], e = 0;
            long long pk = 1, sum = 1;       // sum = 1 + p + ... + p^e
            while(x % p == 0){ x /= p; e++; pk *= p; sum += pk; }
            cout << (first ? " " : " * ") << p << "^" << e;
            first = false;
            d *= e + 1;
            s *= sum;
        }
        cout << ", d = " << d << ", sigma = " << s << "\n";
    }
    return 0;
}

import sys
from math import isqrt

def build_spf(n):
    spf = list(range(n + 1))                 # tạm coi spf[v] = v
    r = isqrt(n)
    mark = bytearray(r + 1)                  # sàng nhỏ tìm số nguyên tố <= sqrt(n)
    small = []
    for i in range(2, r + 1):
        if not mark[i]:
            small.append(i)
            mark[i * i::i] = b'\x01' * len(range(i * i, r + 1, i))
    for p in reversed(small):                # p lớn gán trước, p nhỏ ghi đè sau
        spf[p * p::p] = [p] * len(range(p * p, n + 1, p))
    return spf

data = sys.stdin.buffer.read().split()
q = int(data[0])
xs = list(map(int, data[1:q + 1]))
spf = build_spf(max(xs + [2]))
out = []
for x in xs:
    parts, d, s = [], 1, 1
    while x > 1:
        p, e, pk, sm = spf[x], 0, 1, 1
        while x % p == 0:
            x //= p; e += 1; pk *= p; sm += pk
        parts.append('%d^%d' % (p, e))
        d *= e + 1
        s *= sm
    out.append('%d = %s, d = %d, sigma = %d' % (xs[len(out)], ' * '.join(parts) or '1', d, s))
print('\n'.join(out))

Khung chạy tay dưới đây sàng spf cho $N$ nhỏ rồi phân tích một số $x$. Để ý các dòng skip: hợp số được nhận ra ngay vì spf đã khác 0, không cần kiểm tra gì thêm.

C++▶ Chạy tay từng dòng codeSàng spf tới N (ô xanh lá: số nguyên tố), rồi phân tích x bằng cách tra spf[x] liên tục, đồng thời tính số ước d và tổng ước s.

Dữ liệu vào (stdin) – sửa thoải mái30 24

Dữ liệu mẫu:

1#include <bits/stdc++.h>

2using namespace std;

3

4int spf[105];

5

6int main(){

7 int N, x;

8 cin >> N >> x;

9 for(int i = 2; i <= N; i++){

10 if(spf[i] != 0) continue;

11 spf[i] = i;

12 for(int j = i * i; j <= N; j += i)

13 if(spf[j] == 0) spf[j] = i;

14 }

15 long long d = 1, s = 1;

16 cout << x << " =";

17 bool first = true;

18 while(x > 1){

19 int p = spf[x], e = 0;

20 long long pk = 1, sum = 1;

21 while(x % p == 0){ x /= p; e++; pk *= p; sum += pk; }

22 cout << (first ? " " : " * ") << p << "^" << e;

23 first = false;

24 d *= e + 1; s *= sum;

25 }

26 cout << ", d = " << d << ", sigma = " << s;

27 return 0;

28}

ChậmVừaNhanhRất nhanh

Khuôn mẫu sàng spf

for(int i = 2; i <= N; i++) if(spf[i] == 0){
    spf[i] = i;
    for(long long j = 1LL * i * i; j <= N; j += i)
        if(spf[j] == 0) spf[j] = i;
}
// phân tích x
while(x > 1){ int p = spf[x]; int e = 0; while(x % p == 0){ x /= p; e++; } /* dùng (p, e) */ }

Nhớ 1LL * i * i: với $N = 10^7$ thì i * i tính bằng int sẽ tràn khi $i > 46340$. Mảng int spf[10^7] tốn 40 MB, hãy xem giới hạn bộ nhớ của đề trước khi khai báo.

Python và sàng lớn

Vòng for lồng nhau trong Python chỉ chạy nổi tới khoảng $N = 10^6$. Code Python ở trên dùng mẹo gán theo lát cắt spf[p*p::p] = [p] * k (chạy ở tốc độ C) và đi từ p lớn xuống p nhỏ, để số nguyên tố nhỏ ghi đè sau cùng, nên mỗi ô giữ đúng ước nhỏ nhất. Cách này chịu được $N$ cỡ vài triệu. Nếu cần tới $10^7$ trong Python thì nên dùng array hoặc đổi sang C++.

Sàng tuyến tính. Có một biến thể gọi là sàng tuyến tính (linear sieve), mỗi hợp số chỉ bị gán đúng một lần nên chạy đúng $O(N)$, đồng thời cho luôn danh sách số nguyên tố. Trên thực tế sàng spf ở trên đã đủ nhanh: với $N = 10^7$ cả hai đều chạy dưới 0,1 giây trong C++. Em chỉ cần sàng tuyến tính khi muốn tính đồng thời các hàm như $\varphi$, số ước cho mọi số tới $N$.

Phần 2Sàng nguyên tố trên đoạn [L, R]

Đề: đếm (hoặc liệt kê) số nguyên tố trong đoạn $[L, R]$ với $R \le 10^{12}$ nhưng $R - L \le 10^6$. Không thể khai báo mảng dài $10^{12}$; còn kiểm tra từng số bằng thử chia thì tốn $10^6 \cdot 10^6 = 10^{12}$ phép tính.

Nhận xét then chốt: mọi hợp số $v \le R$ đều có một ước nguyên tố $p \le \sqrt v \le \sqrt R = 10^6$. Vậy chỉ cần hai bước:

1. Sàng thường trên $[2, \sqrt R]$ để lấy danh sách các số nguyên tố nhỏ (khoảng 78 nghìn số với $\sqrt R = 10^6$).
2. Khai báo mảng đánh dấu dài $R - L + 1$, ô thứ $k$ ứng với số $L + k$ (dời gốc về $L$). Với mỗi $p$ nhỏ, gạch các bội của $p$ trong đoạn: bắt đầu từ $\max(p^2, \lceil L / p \rceil \cdot p)$, nhảy bước $p$.

Bội đầu tiên $\ge L$ của $p$ là $\lceil L/p \rceil \cdot p$, tính bằng số nguyên: (L + p - 1) / p * p. Lấy max với $p^2$ để không gạch nhầm chính $p$ khi đoạn chứa số nguyên tố nhỏ (ví dụ $L = 2$). Số ô bị gạch là $\sum_p \frac{R - L}{p} \approx (R - L) \log \log R$, cộng thêm $\sqrt R$ cho vòng duyệt $p$: chạy rất nhanh.

Hình 1. Sàng đoạn: sàng thường trên đoạn nhỏ $[2, \sqrt R]$ rồi dùng các số nguyên tố tìm được để gạch bội trong $[L, R]$. Toàn bộ khoảng giữa bị bỏ qua.

Mô phỏng dưới đây chạy đúng hai bước đó với đoạn nhỏ. Em để ý ô đầu tiên bị gạch của mỗi $p$, và thử $L = 1$ để thấy vì sao số 1 phải loại riêng:

MÔ PHỎNGSàng nguyên tố trên đoạn [L, R]Nhập L, R (R ≤ 10⁶, đoạn dài tối đa 60 số). Dòng trên là các số nguyên tố ≤ √R dùng để gạch.

LR

ChậmVừaNhanh

p đang dùngbội của p trong đoạnđã bị gạchsố nguyên tố

Xem code C++ / Python: sàng nguyên tố trên đoạn [L, R] (bấm để mở)

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

typedef long long ll;

int main(){
    ll L, R;
    cin >> L >> R;
    ll lim = sqrtl((long double)R);
    while(lim * lim > R) lim--;
    while((lim + 1) * (lim + 1) <= R) lim++;
    // Bước 1: sàng thường các số nguyên tố <= sqrt(R)
    vector<char> comp(lim + 1, 0);
    vector<ll> primes;
    for(ll i = 2; i <= lim; i++) if(!comp[i]){
        primes.push_back(i);
        for(ll j = i * i; j <= lim; j += i) comp[j] = 1;
    }
    // Bước 2: gạch bội của từng p trên đoạn, chỉ số dời đi L
    vector<char> isP(R - L + 1, 1);
    for(ll p : primes){
        ll start = max(p * p, (L + p - 1) / p * p);
        for(ll j = start; j <= R; j += p) isP[j - L] = 0;
    }
    if(L == 1) isP[0] = 0;                  // 1 không phải số nguyên tố
    int cnt = 0;
    for(char c : isP) cnt += c;
    cout << cnt;
    return 0;
}

from math import isqrt

L, R = map(int, input().split())
lim = isqrt(R)
comp = bytearray(lim + 1)
primes = []
for i in range(2, lim + 1):                 # sàng thường các số <= sqrt(R)
    if not comp[i]:
        primes.append(i)
        comp[i * i::i] = b'\x01' * len(range(i * i, lim + 1, i))
seg = bytearray(b'\x01') * (R - L + 1)      # seg[k] ứng với số L + k
for p in primes:
    start = max(p * p, (L + p - 1) // p * p)
    if start <= R:
        seg[start - L::p] = bytes(len(range(start, R + 1, p)))
if L == 1:
    seg[0] = 0                              # 1 không phải số nguyên tố
print(sum(seg))

Ba lỗi hay gặp khi sàng đoạn

1. Quên trường hợp $L = 1$: số 1 không bị số nào gạch nhưng không phải số nguyên tố.
2. Bắt đầu từ $\lceil L/p \rceil \cdot p$ mà không lấy max với $p^2$: nếu $L \le p$ thì chính $p$ bị gạch.
3. Tràn số: p * p, j, L + p - 1 đều phải là long long. Tính $\sqrt R$ bằng sqrt số thực có thể lệch 1, nên chỉnh lại bằng hai vòng while như trong code.

Mở rộng: phân tích mọi số trong đoạn. Cũng với các số nguyên tố $p \le \sqrt R$, thay vì gạch, ta chia hẳn rest[j - L] cho $p$ tới khi hết (bắt đầu từ $\lceil L/p \rceil \cdot p$, không cần max với $p^2$). Sau khi duyệt hết các $p$, ô nào còn rest > 1 thì phần còn lại đó là một số nguyên tố lớn hơn $\sqrt R$ (không thể có hai thừa số như vậy). Nhờ vậy ta đếm được số ước, tính $\varphi$… của cả triệu số quanh $10^{12}$.

Phần 3Hàm phi Euler và định lý Euler

Hàm phi Euler $\varphi(n)$ là số các số nguyên trong $[1, n]$ nguyên tố cùng nhau với $n$. Ví dụ $\varphi(12) = 4$ vì chỉ có 1, 5, 7, 11. Với $p$ nguyên tố thì $\varphi(p) = p - 1$.

n	1	2	3	4	5	6	7	8	9	10	11	12
$\varphi(n)$	1	1	2	2	4	2	6	4	6	4	10	4

Bảng 2. Mười hai giá trị đầu của hàm phi. Ở các số nguyên tố 2, 3, 5, 7, 11, giá trị là $p - 1$.

Công thức. Nếu $n = p_1^{e_1} \cdots p_k^{e_k}$ thì

$\varphi(n) = n \left(1 - \dfrac{1}{p_1}\right)\left(1 - \dfrac{1}{p_2}\right) \cdots \left(1 - \dfrac{1}{p_k}\right)$.

Hiểu theo kiểu “sàng”: trong $n$ số, bỏ đi $\frac{1}{p_1}$ số chia hết cho $p_1$, trong phần còn lại bỏ tiếp $\frac{1}{p_2}$ số chia hết cho $p_2$… (các ước nguyên tố khác nhau “độc lập” với nhau, chứng minh bằng nguyên lý bù trừ). Ví dụ $\varphi(36) = 36 \cdot \frac12 \cdot \frac23 = 12$. Chú ý công thức chỉ phụ thuộc vào các ước nguyên tố khác nhau, không phụ thuộc số mũ.

Tính cho một số $n \le 10^{12}$: phân tích bằng thử chia tới $\sqrt n$ như bài 12, gặp ước nguyên tố $p$ thì làm res -= res / p (chính là nhân với $1 - \frac1p$ mà không cần số thực, phép chia luôn chia hết). Độ phức tạp $O(\sqrt n)$.

Xem code C++ / Python: tính phi Euler của một số n trong $O(\sqrt n)$ (bấm để mở)

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

long long phi(long long n){
    long long res = n;
    for(long long p = 2; p * p <= n; p++)
        if(n % p == 0){
            while(n % p == 0) n /= p;
            res -= res / p;                 // res = res * (1 - 1/p)
        }
    if(n > 1) res -= res / n;               // còn một ước nguyên tố lớn
    return res;
}

int main(){
    long long n;
    cin >> n;
    cout << phi(n);
    return 0;
}

def phi(n):
    res, p = n, 2
    while p * p <= n:
        if n % p == 0:
            while n % p == 0:
                n //= p
            res -= res // p                 # res = res * (1 - 1/p)
        p += 1
    if n > 1:
        res -= res // n                     # còn một ước nguyên tố lớn
    return res

print(phi(int(input())))

Sàng phi cho mọi số tới $N$. Khởi tạo ph[i] = i. Duyệt $p$ từ 2: nếu ph[p] == p (chưa bị số nguyên tố nào nhỏ hơn đụng tới) thì $p$ là số nguyên tố, với mọi bội $j$ của $p$ làm ph[j] -= ph[j] / p. Mỗi $j$ được nhân đúng một lần với $(1 - \frac1p)$ cho từng ước nguyên tố $p$ của nó. Độ phức tạp $O(N \log \log N)$ như sàng Eratosthenes.

Xem code C++ / Python: sàng phi Euler cho 1 … N (bấm để mở)

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

int ph[1000005];

int main(){
    int n;
    cin >> n;
    for(int i = 0; i <= n; i++) ph[i] = i;
    for(int p = 2; p <= n; p++)
        if(ph[p] == p)                      // chưa bị đụng tới: p nguyên tố
            for(int j = p; j <= n; j += p)
                ph[j] -= ph[j] / p;         // nhân thêm (1 - 1/p)
    for(int i = 1; i <= n; i++) cout << ph[i] << " ";
    return 0;
}

n = int(input())
ph = list(range(n + 1))
for p in range(2, n + 1):
    if ph[p] == p:                          # chưa bị đụng tới: p nguyên tố
        for j in range(p, n + 1, p):
            ph[j] -= ph[j] // p             # nhân thêm (1 - 1/p)
print(*ph[1:])

Định lý Euler: hạ bậc số mũ khổng lồ

Định lý Euler: nếu $\gcd(a, m) = 1$ thì $a^{\varphi(m)} \equiv 1 \pmod m$. Khi $m = p$ nguyên tố, $\varphi(p) = p - 1$ và ta có lại định lý Fermat nhỏ của bài 26. Hệ quả quan trọng nhất cho lập trình: khi $\gcd(a, m) = 1$,

$a^{b} \equiv a^{\,b \bmod \varphi(m)} \pmod m$.

Ví dụ: $3^{100} \bmod 10$. Vì $\gcd(3, 10) = 1$ và $\varphi(10) = 4$, $100 \bmod 4 = 0$ nên $3^{100} \equiv 3^0 = 1$. Kiểm tra: $3^4 = 81 \equiv 1$, đúng là cứ 4 bước lặp lại. Áp dụng hay gặp nhất: tính $a^{b^c} \bmod (10^9 + 7)$. Không thể tính $b^c$ (có tới hàng tỉ chữ số), nhưng chỉ cần $e = b^c \bmod (10^9 + 6)$ bằng lũy thừa nhanh, rồi tính $a^e \bmod (10^9 + 7)$ (bài tập 5).

Thêm một tính chất đẹp, dùng trong bài tập 8: $\sum_{d \mid n} \varphi(d) = n$. Lý do: chia các số $1..n$ theo giá trị $\gcd(i, n) = g$; số các $i$ có $\gcd(i, n) = g$ đúng bằng $\varphi(n/g)$.

Khi $\gcd(a, m) \ne 1$: định lý Euler mở rộng

Công thức $a^b \equiv a^{b \bmod \varphi(m)}$ sai nếu $a$, $m$ không nguyên tố cùng nhau. Ví dụ $2^{4} \bmod 8 = 0$ nhưng $\varphi(8) = 4$ cho $2^{0} = 1$. Dạng đúng cho mọi $a$: nếu $b \ge \log_2 m$ thì $a^b \equiv a^{(b \bmod \varphi(m)) + \varphi(m)} \pmod m$. Phải cộng thêm $\varphi(m)$, và phải biết $b$ có “đủ lớn” hay không, nên khi tính tháp lũy thừa cần giữ thêm cờ “đã vượt ngưỡng”.

Phần 4Euclid mở rộng và phương trình ax + by = c

Ở bài 26 em đã viết Euclid mở rộng bằng đệ quy để tìm $x, y$ sao cho $ax + by = \gcd(a, b)$ rồi lấy nghịch đảo modulo. Bài này ta học thêm dạng vòng lặp (không lo tràn ngăn xếp, dễ chạy tay) và dùng nó để giải trọn vẹn phương trình $ax + by = c$.

Dạng vòng lặp. Ta giữ một “bảng” các hàng $(r, x, y)$ luôn thỏa bất biến $r = a \cdot x + b \cdot y$. Hai hàng đầu là $(a, 1, 0)$ và $(b, 0, 1)$ (hiển nhiên đúng). Mỗi bước, giống Euclid thường, lấy $q = \lfloor r_0 / r_1 \rfloor$ và tạo hàng mới bằng hàng trên trừ $q$ lần hàng dưới: $r_2 = r_0 - q r_1$, $x_2 = x_0 - q x_1$, $y_2 = y_0 - q y_1$. Vì phép trừ giữ nguyên đẳng thức nên hàng mới vẫn đúng bất biến. Khi $r$ về 0, hàng ngay trước cho $g$ và cặp $(x, y)$ cần tìm.

r	q	x	y	kiểm tra: 240x + 46y
240		1	0	240
46	5	0	1	46
10	4	1	−5	240 − 230 = 10
6	1	−4	21	−960 + 966 = 6
4	1	5	−26	1200 − 1196 = 4
2	2	−9	47	−2160 + 2162 = 2
0		23	−120	0

Bảng 3. Euclid mở rộng dạng bảng cho $a = 240$, $b = 46$. Cột $q$ ghi ở hàng chia; hàng tô xanh là hàng cuối có $r \ne 0$: $\gcd = 2 = 240 \cdot (-9) + 46 \cdot 47$.

C++▶ Chạy tay từng dòng codeEuclid mở rộng dạng vòng lặp: ba cột r, x, y cùng biến đổi theo một quy tắc “hàng trên trừ q lần hàng dưới”, nên luôn giữ r = a·x + b·y.

Dữ liệu vào (stdin) – sửa thoải mái240 46

Dữ liệu mẫu:

1#include <bits/stdc++.h>

2using namespace std;

3

4int main(){

5 long long a, b;

6 cin >> a >> b;

7 long long r0 = a, r1 = b;

8 long long x0 = 1, x1 = 0;

9 long long y0 = 0, y1 = 1;

10 while(r1 != 0){

11 long long q = r0 / r1;

12 long long r2 = r0 - q * r1; r0 = r1; r1 = r2;

13 long long x2 = x0 - q * x1; x0 = x1; x1 = x2;

14 long long y2 = y0 - q * y1; y0 = y1; y1 = y2;

15 }

16 cout << r0 << " " << x0 << " " << y0;

17 return 0;

18}

ChậmVừaNhanhRất nhanh

Giải phương trình $ax + by = c$ trong số nguyên

Gọi $g = \gcd(a, b)$. Vế trái luôn chia hết cho $g$, nên có nghiệm khi và chỉ khi $g \mid c$. Khi đó:

• Một nghiệm riêng: Euclid mở rộng cho $a x_g + b y_g = g$. Nhân hai vế với $\frac{c}{g}$: $x_0 = x_g \cdot \frac cg$, $y_0 = y_g \cdot \frac cg$.
• Mọi nghiệm: $x = x_0 + k \cdot \frac bg$, $y = y_0 - k \cdot \frac ag$ với $k \in \mathbb{Z}$. (Tăng $x$ thêm $\frac bg$ làm vế trái tăng $\frac{ab}{g}$, giảm $y$ đi $\frac ag$ làm vế trái giảm đúng $\frac{ab}{g}$.)

Ví dụ $6x + 9y = 21$: $g = 3 \mid 21$. Từ $6 \cdot 2 + 9 \cdot (-1) = 3$ được nghiệm riêng $(14, -7)$, họ nghiệm $x = 14 + 3k$, $y = -7 - 2k$. Nghiệm có $x$ không âm nhỏ nhất là $x = 14 \bmod 3 = 2$, $y = 1$: $12 + 9 = 21$. ✔

Đề hay hỏi: nghiệm có $x \ge 0$ nhỏ nhất, số nghiệm với $x, y \ge 0$ (bài tập 6), nghiệm có $|x| + |y|$ nhỏ nhất… Tất cả đều quy về chọn $k$ trong họ nghiệm tổng quát.

Xem code C++ / Python: giải ax + by = c, in nghiệm có x không âm nhỏ nhất (bấm để mở)

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

typedef long long ll;

// Euclid mở rộng dạng vòng lặp: trả về g, và x, y với a*x + b*y = g
ll extgcd(ll a, ll b, ll &x, ll &y){
    ll x0 = 1, x1 = 0, y0 = 0, y1 = 1;
    while(b != 0){
        ll q = a / b, t;
        t = a - q * b; a = b; b = t;
        t = x0 - q * x1; x0 = x1; x1 = t;
        t = y0 - q * y1; y0 = y1; y1 = t;
    }
    x = x0; y = y0;
    return a;
}

int main(){
    ll a, b, c;
    cin >> a >> b >> c;
    ll x, y;
    ll g = extgcd(a, b, x, y);
    if(c % g != 0){ cout << -1; return 0; }   // vô nghiệm
    ll bg = b / g;                             // bước nhảy của x
    // x0 = x * (c / g), lấy mod bg từng thừa số để không tràn số
    ll xm = ((x % bg + bg) % bg) * ((c / g) % bg) % bg;
    ll ym = (c - a * xm) / b;
    cout << xm << " " << ym;
    return 0;
}

def extgcd(a, b):
    x0, x1, y0, y1 = 1, 0, 0, 1
    while b:
        q = a // b
        a, b = b, a - q * b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return a, x0, y0                        # a*x0 + b*y0 = g

a, b, c = map(int, input().split())
g, x, y = extgcd(a, b)
if c % g != 0:
    print(-1)                               # vô nghiệm
else:
    bg = b // g                             # bước nhảy của x
    xm = x * (c // g) % bg                  # Python: % luôn ra số không âm
    print(xm, (c - a * xm) // b)

Chống tràn số khi tính nghiệm riêng

$x_g$ có thể cỡ $b$ (tới $10^9$), $\frac cg$ tới $10^{18}$: tích tới $10^{27}$, tràn long long. Vì chỉ cần $x$ theo modulo $\frac bg$, hãy lấy modulo từng thừa số trước rồi mới nhân: (x % bg + bg) % bg * ((c / g) % bg) % bg, mỗi thừa số dưới $10^9$ nên tích dưới $10^{18}$. Trong Python không có tràn số, và % luôn trả về số không âm.

Phương trình đồng dư $ax \equiv b \pmod m$ chính là $ax + my = b$, đã giải ở bài tập 7 của bài 26: có nghiệm khi $g = \gcd(a, m) \mid b$, và khi đó có đúng $g$ nghiệm trong $[0, m)$, cách nhau $\frac mg$. Đây cũng là bước cốt lõi của phần tiếp theo.

Phần 5Định lý thặng dư Trung Hoa

Bài toán “Hàn Tín điểm binh”: một đội quân xếp hàng 3 thì dư 2 người, hàng 5 dư 3, hàng 7 dư 2. Hỏi đội có ít nhất bao nhiêu người? Tức là tìm $x$ nhỏ nhất thỏa

$x \equiv 2 \pmod 3$, $\quad x \equiv 3 \pmod 5$, $\quad x \equiv 2 \pmod 7$.

Định lý thặng dư Trung Hoa (CRT): nếu các modulo $m_1, \ldots, m_k$ đôi một nguyên tố cùng nhau thì hệ luôn có nghiệm, và nghiệm là duy nhất theo modulo $M = m_1 m_2 \cdots m_k$. Đáp án ví dụ trên là $x = 23$ (và 128, 233, … cách nhau $105 = 3 \cdot 5 \cdot 7$).

Ghép hai phương trình. Thay vì nhớ công thức tổng, cách dễ cài và dùng được cả khi modulo không nguyên tố cùng nhau là ghép dần từng cặp. Giả sử đã biết $x \equiv A \pmod M$ và cần thêm $x \equiv a \pmod m$:

1. Mọi nghiệm của phương trình đầu có dạng $x = A + M t$. Thay vào phương trình sau: $M t \equiv a - A \pmod m$.
2. Đây là phương trình đồng dư ở phần 4. Gọi $g = \gcd(M, m)$. Nếu $g \nmid (a - A)$ thì hệ vô nghiệm.
3. Ngược lại, Euclid mở rộng cho $M p + m q = g$, suy ra $t \equiv \frac{a - A}{g} \cdot p \pmod{\frac mg}$.
4. Hệ mới: $x \equiv A + M t \pmod{\operatorname{lcm}(M, m)}$, với $\operatorname{lcm}(M, m) = M \cdot \frac mg$.

Bắt đầu từ $A = 0, M = 1$ (mọi số đều thỏa) và ghép lần lượt $k$ phương trình, tổng độ phức tạp $O(k \log M)$.

Ghép	g	Phương trình cho t	t	Hệ sau khi ghép
$x \equiv 2 \pmod 3$	1	$1 \cdot t \equiv 2 \pmod 3$	2	$x \equiv 2 \pmod 3$
$x \equiv 3 \pmod 5$	1	$3 t \equiv 1 \pmod 5$	2	$x \equiv 2 + 3 \cdot 2 = 8 \pmod{15}$
$x \equiv 2 \pmod 7$	1	$15 t \equiv -6 \pmod 7$	1	$x \equiv 8 + 15 = 23 \pmod{105}$

Bảng 4. Ghép từng phương trình cho bài Hàn Tín điểm binh. Ở hàng 3: $15 \equiv 1$ nên $t \equiv -6 \equiv 1 \pmod 7$.

Trong mô phỏng, hàng ô trên cùng là các ứng viên $A + M t$ của hệ hiện tại, con số nhỏ bên dưới là số dư khi chia cho $m$ mới. Ta không thử từng ô mà dùng bảng Euclid mở rộng để nhảy thẳng tới $t$ đúng. Em thử thêm hệ không nguyên tố cùng nhau “1 4 3 6” (có nghiệm 9) và “0 4 1 6” (vô nghiệm, vì một số chẵn không thể chia 6 dư 1):

MÔ PHỎNGĐịnh lý Trung Hoa: ghép từng phương trình bằng Euclid mở rộngNhập các cặp “a m” (tối đa 4 cặp, m ≤ 60). Mỗi lần ghép có bảng Euclid mở rộng để tìm t.

Các cặp a m

ChậmVừaNhanh

giá trị t được chọn / thỏa phương trình mớiứng viên x = A + M·t

Xem code C++ / Python: định lý thặng dư Trung Hoa cho k phương trình, modulo bất kỳ (bấm để mở)

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

typedef long long ll;
typedef __int128 lll;

ll extgcd(ll a, ll b, ll &x, ll &y){
    if(b == 0){ x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = extgcd(b, a % b, x1, y1);
    x = y1; y = x1 - (a / b) * y1;
    return g;
}

int main(){
    int k;
    cin >> k;
    ll A = 0, M = 1;                          // hệ hiện tại: x ≡ A (mod M)
    for(int i = 0; i < k; i++){
        ll a, m;
        cin >> a >> m;
        ll p, q;
        ll g = extgcd(M, m, p, q);            // M*p + m*q = g
        if((a - A) % g != 0){ cout << -1; return 0; }
        ll mg = m / g;
        // t = (a - A)/g * p  (mod m/g), nhân bằng __int128 cho chắc
        ll t = (ll)((lll)((a - A) / g % mg) * (p % mg) % mg);
        if(t < 0) t += mg;
        ll newM = M * mg;                     // BCNN(M, m), đề đảm bảo <= 1e18
        A = (ll)(((lll)A + (lll)M * t) % newM);
        M = newM;
    }
    cout << A;
    return 0;
}

import sys

def extgcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extgcd(b, a % b)
    return g, y1, x1 - (a // b) * y1

data = sys.stdin.read().split()
k = int(data[0])
A, M = 0, 1                                 # hệ hiện tại: x ≡ A (mod M)
ok = True
for i in range(k):
    a, m = int(data[1 + 2 * i]), int(data[2 + 2 * i])
    g, p, q = extgcd(M, m)                  # M*p + m*q = g
    if (a - A) % g != 0:
        ok = False
        break
    mg = m // g
    t = (a - A) // g * p % mg               # số nguyên Python không bao giờ tràn
    A = (A + M * t) % (M * mg)
    M *= mg
print(A if ok else -1)

Tràn số trong định lý Trung Hoa

Đề thường cho $m_i \le 10^9$ và đảm bảo BCNN $\le 10^{18}$. Hai chỗ nguy hiểm: tích $\frac{a - A}{g} \cdot p$ (tới $10^{18} \cdot 10^9$) và $M \cdot t$ (tới $10^{18} \cdot 10^9$ nếu chưa lấy modulo). Cách an toàn: lấy $\frac{a - A}{g} \bmod \frac mg$ trước, nhân bằng __int128, và nhớ $t < \frac mg$ nên $M t < M \cdot \frac mg = $ BCNN mới. Nếu BCNN có thể vượt $10^{18}$ thì chỉ Python (số nguyên không giới hạn) hoặc __int128 mới chứa nổi.

Khi nào nghĩ tới định lý Trung Hoa?

1. Đề cho hệ “chia … dư …”, các chu kỳ lặp lại (đèn nhấp nháy, xe buýt, hành tinh thẳng hàng).
2. Cần tính một đại lượng modulo hợp số $m = p_1^{e_1} \cdots$ mà mỗi modulo $p_i^{e_i}$ thì dễ tính: tính riêng từng phần rồi ghép lại bằng CRT (ví dụ $C(n, k) \bmod 10^6$).
3. Hai phương trình $x \equiv a \pmod m$ với modulo không nguyên tố cùng nhau: đừng vội nghĩ vô nghiệm, chỉ cần kiểm tra $g \mid (a_2 - a_1)$.

Luyện tập trên web

Các bài dưới đây trên blog đều có lời giải chi tiết bằng C++ và Python, xếp từ dễ tới khó. Em hãy tự làm trước rồi mới xem lời giải nhé.

• Chia Quà (CHIAQUA) (Dễ)
• Future Number 1 (FUTURENUM) (Trung bình)
• Phân Số (Trung bình)
• Tích Các Ước (MULDIV) (Khá)
• Số Có 9 Ước Số (Bản Khó) (Khá)
• Cặp Ước Đẹp (Khá)
• Tìm Số Có N Ước (NDIVISORS) (Khá)

Xem toàn bộ bài cùng chủ đề: nhãn Số học.

Bài tậpLuyện số học nâng cao

Bài 1Ước nguyên tố lớn nhấtDễ

Cho $Q$ số nguyên $x_1, \ldots, x_Q$ ($Q \le 2 \cdot 10^5$, $2 \le x_i \le 10^6$). Với mỗi số, in ra ước nguyên tố lớn nhất của nó, mỗi số trên một dòng.

Input

4
12
17
100
999983

Output

3
17
5
999983

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Thử chia tới $\sqrt x$ cho từng số thì tốn $2 \cdot 10^5 \cdot 10^3$ phép chia, vẫn chạy được nhưng sát giờ. Gọn hơn: làm sàng giống spf nhưng ghi đè: với mỗi số nguyên tố $p$ (theo thứ tự tăng dần), gán big[j] = p cho mọi bội $j$ của $p$ (bắt đầu từ chính $p$). Số nguyên tố gán sau là số lớn hơn, nên cuối cùng big[j] là ước nguyên tố lớn nhất. Một số $p$ là nguyên tố khi tới lượt nó mà big[p] vẫn bằng 0. Ví dụ $100 = 2^2 \cdot 5^2$ nên đáp án là 5; 999983 là số nguyên tố.

Xem code tham khảo (C++ / Python)

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

const int N = 1000000;
int big[N + 1];               // big[v] = ước nguyên tố lớn nhất của v

int main(){
    ios::sync_with_stdio(false); cin.tie(nullptr);
    for(int p = 2; p <= N; p++)
        if(big[p] == 0)                       // p nguyên tố
            for(int j = p; j <= N; j += p) big[j] = p;   // p sau ghi đè p trước
    int q;
    cin >> q;
    while(q--){
        int x;
        cin >> x;
        cout << big[x] << "\n";
    }
    return 0;
}

import sys

data = sys.stdin.buffer.read().split()
q = int(data[0])
xs = list(map(int, data[1:q + 1]))
n = max(xs)
big = [0] * (n + 1)
for p in range(2, n + 1):
    if big[p] == 0:                           # p nguyên tố
        big[p::p] = [p] * len(range(p, n + 1, p))   # p sau ghi đè p trước
print('\n'.join(str(big[x]) for x in xs))

Bài 2Số ước và tổng ướcDễ

Cho $Q$ số nguyên dương ($Q \le 2 \cdot 10^5$, mỗi số $\le 10^6$). Với mỗi số $x$, in trên một dòng hai số: số ước $d(x)$ và tổng các ước $\sigma(x)$.

Input

3
1
12
360

Output

1 1
6 28
24 1170

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Đúng phần 1: sàng spf một lần tới $10^6$, rồi phân tích mỗi số trong $O(\log x)$. Với mỗi thừa số $p^e$: nhân $d$ với $e + 1$, nhân $\sigma$ với $1 + p + \cdots + p^e$ (cộng dồn trong lúc chia). $12 = 2^2 \cdot 3$: $d = 3 \cdot 2 = 6$, $\sigma = 7 \cdot 4 = 28$. Với $x \le 10^6$ thì $\sigma(x) < 5 \cdot 10^6$, vẫn nên dùng long long cho quen tay.

Xem code tham khảo (C++ / Python)

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

const int N = 1000000;
int spf[N + 1];

int main(){
    ios::sync_with_stdio(false); cin.tie(nullptr);
    for(int i = 2; i <= N; i++) if(spf[i] == 0){
        spf[i] = i;
        for(long long j = 1LL * i * i; j <= N; j += i)
            if(spf[j] == 0) spf[j] = i;
    }
    int q;
    cin >> q;
    while(q--){
        int x;
        cin >> x;
        long long d = 1, s = 1;
        while(x > 1){
            int p = spf[x], e = 0;
            long long pk = 1, sum = 1;
            while(x % p == 0){ x /= p; e++; pk *= p; sum += pk; }
            d *= e + 1; s *= sum;
        }
        cout << d << " " << s << "\n";
    }
    return 0;
}

import sys
from math import isqrt

data = sys.stdin.buffer.read().split()
q = int(data[0])
xs = list(map(int, data[1:q + 1]))
n = max(xs + [2])
spf = list(range(n + 1))
r = isqrt(n)
small = [p for p in range(2, r + 1) if all(p % k for k in range(2, isqrt(p) + 1))]
for p in reversed(small):                     # p nhỏ ghi đè sau cùng
    spf[p * p::p] = [p] * len(range(p * p, n + 1, p))
out = []
for x in xs:
    d = s = 1
    while x > 1:
        p, e, pk, sm = spf[x], 0, 1, 1
        while x % p == 0:
            x //= p; e += 1; pk *= p; sm += pk
        d *= e + 1; s *= sm
    out.append('%d %d' % (d, s))
print('\n'.join(out))

Bài 3Cặp số nguyên tố sinh đôiTrung bình

Hai số nguyên tố $p$ và $p + 2$ gọi là một cặp sinh đôi. Cho $L, R$ ($1 \le L \le R \le 10^{12}$, $R - L \le 10^6$). Đếm số cặp sinh đôi $(p, p + 2)$ nằm trọn trong đoạn, tức $L \le p$ và $p + 2 \le R$.

Input

1 20

Output

4

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Sàng đoạn $[L, R]$ như phần 2 để có mảng isP, rồi đếm các vị trí $k$ mà isP[k] và isP[k + 2] đều đúng. Nhớ loại số 1 khi $L = 1$. Trong $[1, 20]$ có 4 cặp: (3, 5), (5, 7), (11, 13), (17, 19).

Xem code tham khảo (C++ / Python)

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

typedef long long ll;

int main(){
    ll L, R;
    cin >> L >> R;
    ll lim = sqrtl((long double)R);
    while(lim * lim > R) lim--;
    while((lim + 1) * (lim + 1) <= R) lim++;
    vector<char> comp(lim + 1, 0);
    vector<char> isP(R - L + 1, 1);
    for(ll p = 2; p <= lim; p++) if(!comp[p]){
        for(ll j = p * p; j <= lim; j += p) comp[j] = 1;
        for(ll j = max(p * p, (L + p - 1) / p * p); j <= R; j += p) isP[j - L] = 0;
    }
    if(L == 1) isP[0] = 0;
    int cnt = 0;
    for(ll k = 0; k + 2 <= R - L; k++)
        if(isP[k] && isP[k + 2]) cnt++;
    cout << cnt;
    return 0;
}

from math import isqrt

L, R = map(int, input().split())
lim = isqrt(R)
comp = bytearray(lim + 1)
seg = bytearray(b'\x01') * (R - L + 1)
for p in range(2, lim + 1):
    if not comp[p]:
        comp[p * p::p] = b'\x01' * len(range(p * p, lim + 1, p))
        start = max(p * p, (L + p - 1) // p * p)
        if start <= R:
            seg[start - L::p] = bytes(len(range(start, R + 1, p)))
if L == 1:
    seg[0] = 0
print(sum(1 for k in range(R - L - 1) if seg[k] and seg[k + 2]))

Bài 4Phân số tối giảnTrung bình

Cho $n$ ($1 \le n \le 10^6$). Đếm số phân số tối giản $\frac pq$ với $1 \le p < q \le n$.

Input

5

Output

9

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Với mỗi mẫu $q$, số tử $p \in [1, q)$ để $\frac pq$ tối giản là số $p$ nguyên tố cùng nhau với $q$, đúng bằng $\varphi(q)$ (với $q \ge 2$). Đáp án $= \sum_{q=2}^{n} \varphi(q)$, dùng sàng phi ở phần 3. Với $n = 5$: $\varphi(2) + \varphi(3) + \varphi(4) + \varphi(5) = 1 + 2 + 2 + 4 = 9$. Đáp án cỡ $0{,}3 n^2$ nên phải dùng long long.

Xem code tham khảo (C++ / Python)

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

int ph[1000005];

int main(){
    int n;
    cin >> n;
    for(int i = 0; i <= n; i++) ph[i] = i;
    long long ans = 0;
    for(int p = 2; p <= n; p++){
        if(ph[p] == p)
            for(int j = p; j <= n; j += p) ph[j] -= ph[j] / p;
        ans += ph[p];                 // ph[p] đã hoàn chỉnh khi tới lượt p
    }
    cout << ans;
    return 0;
}

n = int(input())
ph = list(range(n + 1))
for p in range(2, n + 1):
    if ph[p] == p:
        for j in range(p, n + 1, p):
            ph[j] -= ph[j] // p
print(sum(ph[2:]))

Bài 5Lũy thừa tầngTrung bình

Có $T$ truy vấn ($T \le 10^5$), mỗi truy vấn gồm ba số $a, b, c$ ($1 \le a \le 10^9$, $1 \le b, c \le 2 \cdot 10^9$). Tính $a^{(b^c)} \bmod (10^9 + 7)$.

Input

3
3 7 1
2 3 2
5 2 3

Output

2187
512
390625

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Gọi $p = 10^9 + 7$ (nguyên tố). Vì $1 \le a < p$ nên $\gcd(a, p) = 1$, theo định lý Euler (ở đây chính là Fermat) $a^{p-1} \equiv 1$, nên chỉ cần số mũ theo modulo $p - 1$: $e = b^c \bmod (p - 1)$, đáp án $= a^e \bmod p$. Hai lần lũy thừa nhanh. Bẫy: lấy $b^c \bmod p$ (thay vì $p - 1$) là sai. Ví dụ 3: $2^3 = 8$, $5^8 = 390625$.

Xem code tham khảo (C++ / Python)

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

typedef long long ll;
const ll P = 1000000007;

ll power(ll a, ll b, ll m){
    ll res = 1 % m; a %= m;
    while(b > 0){
        if(b & 1) res = res * a % m;
        a = a * a % m;
        b >>= 1;
    }
    return res;
}

int main(){
    ios::sync_with_stdio(false); cin.tie(nullptr);
    int t;
    cin >> t;
    while(t--){
        ll a, b, c;
        cin >> a >> b >> c;
        ll e = power(b, c, P - 1);        // số mũ lấy modulo phi(P) = P - 1
        cout << power(a, e, P) << "\n";
    }
    return 0;
}

import sys

P = 10**9 + 7
data = sys.stdin.buffer.read().split()
t = int(data[0])
out = []
for i in range(t):
    a, b, c = int(data[1 + 3 * i]), int(data[2 + 3 * i]), int(data[3 + 3 * i])
    e = pow(b, c, P - 1)                  # số mũ lấy modulo phi(P) = P - 1
    out.append(pow(a, e, P))
print('\n'.join(map(str, out)))

Bài 6Đếm nghiệm không âmKhá

Cho ba số $a, b, c$ ($1 \le a, b \le 10^9$, $0 \le c \le 10^{18}$). Đếm số cặp số nguyên không âm $(x, y)$ thỏa $ax + by = c$.

Input

2 3 12

Output

3

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Nếu $g = \gcd(a, b) \nmid c$ thì đáp án 0. Ngược lại, như phần 4, tìm $x_{\min}$ là nghiệm $x \ge 0$ nhỏ nhất (lấy modulo $\frac bg$ từng thừa số để khỏi tràn). Nếu $a \cdot x_{\min} > c$ thì $y < 0$, đáp án 0. Nếu không, các nghiệm có $x = x_{\min} + k \frac bg$, mỗi lần tăng $k$ thì $a x$ tăng thêm $\operatorname{lcm}(a, b) = \frac ag b$, và cần $a x \le c$, nên số nghiệm là $\left\lfloor \frac{c - a x_{\min}}{\operatorname{lcm}(a, b)} \right\rfloor + 1$. Ví dụ: (0, 4), (3, 2), (6, 0).

Xem code tham khảo (C++ / Python)

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

typedef long long ll;

ll extgcd(ll a, ll b, ll &x, ll &y){
    ll x0 = 1, x1 = 0, y0 = 0, y1 = 1;
    while(b != 0){
        ll q = a / b, t;
        t = a - q * b; a = b; b = t;
        t = x0 - q * x1; x0 = x1; x1 = t;
        t = y0 - q * y1; y0 = y1; y1 = t;
    }
    x = x0; y = y0;
    return a;
}

int main(){
    ll a, b, c;
    cin >> a >> b >> c;
    ll x, y;
    ll g = extgcd(a, b, x, y);
    if(c % g != 0){ cout << 0; return 0; }
    ll bg = b / g;
    ll xm = ((x % bg + bg) % bg) * ((c / g) % bg) % bg;   // x nhỏ nhất >= 0
    if(a * xm > c){ cout << 0; return 0; }                // y sẽ âm
    ll lcm = a / g * b;                                   // <= 1e18
    cout << (c - a * xm) / lcm + 1;
    return 0;
}

def extgcd(a, b):
    x0, x1, y0, y1 = 1, 0, 0, 1
    while b:
        q = a // b
        a, b = b, a - q * b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return a, x0, y0

a, b, c = map(int, input().split())
g, x, y = extgcd(a, b)
if c % g:
    print(0)
else:
    xm = x * (c // g) % (b // g)          # x nhỏ nhất >= 0
    if a * xm > c:
        print(0)                          # y sẽ âm
    else:
        print((c - a * xm) // (a // g * b) + 1)

Bài 7Đèn nhấp nháyKhá

Có $k$ ngọn đèn ($k \le 10^5$). Đèn thứ $i$ sáng đúng vào các giây $t$ thỏa $t \equiv a_i \pmod{m_i}$ ($1 \le m_i \le 10^9$, $0 \le a_i < m_i$). Cho thời điểm $T$ ($0 \le T \le 10^{18}$), tìm giây $t \ge T$ sớm nhất mà tất cả đèn cùng sáng, hoặc in $-1$ nếu không bao giờ xảy ra. Dữ liệu đảm bảo BCNN của các $m_i$ không vượt quá $10^{18}$.

Dòng đầu gồm $k$ và $T$, $k$ dòng sau mỗi dòng gồm $a_i, m_i$.

Input

3 30
2 3
3 5
2 7

Output

128

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Ghép $k$ phương trình bằng định lý Trung Hoa (phần 5) được $t \equiv A \pmod M$ (hoặc vô nghiệm, in $-1$). Các thời điểm cả dãy cùng sáng là $A, A + M, A + 2M, \ldots$ Nếu $A \ge T$ thì đáp án là $A$; ngược lại là $A + \left\lceil \frac{T - A}{M} \right\rceil M$. Đáp án $< T + M \le 2 \cdot 10^{18}$, vừa trong long long. Ví dụ: hệ Hàn Tín cho $A = 23$, $M = 105$; giây đầu tiên $\ge 30$ là $23 + 105 = 128$.

Xem code tham khảo (C++ / Python)

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

typedef long long ll;
typedef __int128 lll;

ll extgcd(ll a, ll b, ll &x, ll &y){
    if(b == 0){ x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = extgcd(b, a % b, x1, y1);
    x = y1; y = x1 - (a / b) * y1;
    return g;
}

int main(){
    ios::sync_with_stdio(false); cin.tie(nullptr);
    int k;
    ll T;
    cin >> k >> T;
    ll A = 0, M = 1;
    bool ok = true;
    for(int i = 0; i < k; i++){
        ll a, m;
        cin >> a >> m;
        if(!ok) continue;
        ll p, q;
        ll g = extgcd(M, m, p, q);
        if((a - A) % g != 0){ ok = false; continue; }
        ll mg = m / g;
        ll t = (ll)((lll)((a - A) / g % mg) * (p % mg) % mg);
        if(t < 0) t += mg;
        ll newM = M * mg;
        A = (ll)(((lll)A + (lll)M * t) % newM);
        M = newM;
    }
    if(!ok){ cout << -1; return 0; }
    if(A >= T) cout << A;
    else cout << A + (T - A + M - 1) / M * M;
    return 0;
}

import sys

def extgcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extgcd(b, a % b)
    return g, y1, x1 - (a // b) * y1

data = sys.stdin.buffer.read().split()
k, T = int(data[0]), int(data[1])
A, M, ok = 0, 1, True
for i in range(k):
    a, m = int(data[2 + 2 * i]), int(data[3 + 2 * i])
    g, p, q = extgcd(M, m)
    if (a - A) % g:
        ok = False
        break
    mg = m // g
    t = (a - A) // g * p % mg
    A = (A + M * t) % (M * mg)
    M *= mg
if not ok:
    print(-1)
elif A >= T:
    print(A)
else:
    print(A + (T - A + M - 1) // M * M)

Bài 8Tổng ước chung lớn nhấtKhá

Cho $n$ ($1 \le n \le 10^{12}$). Tính $S = \gcd(1, n) + \gcd(2, n) + \cdots + \gcd(n, n)$.

Input

6

Output

15

📝 Nộp bài, chấm tự động trên ClueOJ ↗

💡 Gợi ý

Gom các $i$ theo giá trị $g = \gcd(i, n)$: $\gcd(i, n) = g$ khi và chỉ khi $i = g j$ với $\gcd(j, \frac ng) = 1$, có đúng $\varphi(\frac ng)$ số như vậy. Vậy $S = \sum_{g \mid n} g \cdot \varphi(\frac ng)$. Hàm này nhân tính, với $n = \prod p^e$ ta có công thức gọn: $S = \prod \big((e + 1) p^e - e \, p^{e-1}\big)$. Chỉ cần phân tích $n$ trong $O(\sqrt n)$. Với $n = 6 = 2 \cdot 3$: $(2 \cdot 2 - 1)(2 \cdot 3 - 1) = 3 \cdot 5 = 15$, khớp với $1 + 2 + 3 + 2 + 1 + 6$.

Xem code tham khảo (C++ / Python)

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

typedef long long ll;

int main(){
    ll n;
    cin >> n;
    ll S = 1;
    for(ll p = 2; p * p <= n; p++)
        if(n % p == 0){
            ll e = 0, pe = 1;
            while(n % p == 0){ n /= p; e++; pe *= p; }
            S *= (e + 1) * pe - e * (pe / p);      // (e+1)p^e - e p^(e-1)
        }
    if(n > 1) S *= 2 * n - 1;                       // thừa số nguyên tố cuối, e = 1
    cout << S;
    return 0;
}

n = int(input())
S, p = 1, 2
while p * p <= n:
    if n % p == 0:
        e, pe = 0, 1
        while n % p == 0:
            n //= p; e += 1; pe *= p
        S *= (e + 1) * pe - e * (pe // p)       # (e+1)p^e - e p^(e-1)
    p += 1
if n > 1:
    S *= 2 * n - 1                               # thừa số nguyên tố cuối, e = 1
print(S)

Bài tiếp theo · Bài 41

Kỹ thuật tìm kiếm nâng cao

Chặt nhị phân đáp án dạng cực tiểu hoá giá trị lớn nhất, chặt tam phân cho hàm lồi, gặp nhau ở giữa (meet in the middle) với n = 40 và nhánh cận giúp quay lui nhanh gấp nghìn lần.

Học bài 41 →

← Bài 39: Z-function, Manacher, TrieXem cả series Thuật toán chuyên sâu

Rất vui khi anh em nào đã đọc đến những dòng cuối cùng này. Anh em đọc chắc cũng mỏi mắt luôn chứ 😅 vậy anh em nghĩ người viết ra hết đống ở trên còn thế nào nữa 🤪. Vậy nên đừng tiếc gì cho mình một lời bình luận hoặc chia sẻ bài viết để động viên mình nè. Chúc anh em học tốt và gặt hái nhiều thành công 😇

Bản quyền thuộc về Thầy Nguyễn Việt Nam Sơn. Bài viết, lời giải, hình minh hoạ và mã nguồn trên trang là công sức biên soạn của thầy qua nhiều năm giảng dạy. Vui lòng không sao chép, đăng lại hay sử dụng dưới mọi hình thức khi chưa được sự đồng ý của tác giả. Nếu thấy hữu ích, hãy chia sẻ bằng đường link bài gốc nhé. Cảm ơn bạn!