Hướng dẫn cho LQDOJ CUP 2022 - Round 6 - SUSSET
Chép code từ bài hướng dẫn để nộp bài là hành vi có thể dẫn đến khóa tài khoản.
Authors: ,
Subtask \(1\) (\(10\%\) số điểm): \(n \leq 10\)
Tutorial
Vì \(k\) số đều phải phân biệt và các số đều là các số nguyên thuộc tập \(\{1, 2, \ldots, n\}\), nên ta chỉ cần duyệt qua toàn bộ các tập con.
Với mỗi tập con thì ta xây dựng đồ thị và sử dụng thuật toán để tìm đường đi ngắn nhất.
Sau đó giá trị bí ẩn nhỏ nhất và đếm số tập có cùng giá trị đó
Ta dùng duyệt bit hoặc đệ quy để tạo ra các tập hợp \(S\) gồm \(k\) đỉnh thoả \(S_i \in \{1, 2, \ldots, n\}\).
Rồi dựng đồ thị theo yêu cầu đề bài, xét 3 trường hợp
- \(u = v\)
- \(u \neq v\) và \((u, v)\) là cặp bí ẩn
- \(u \neq v\) và \((u, v)\) không phải là cặp bí ẩn
Sau đó ta sử dụng thuật toán Floyd-Warshall để tìm đường đi ngắn nhất giữa tất cả cặp đỉnh.
Sau khi tính được độ bí ẩn của một tập hợp \(S\) bất kì, ta duyệt qua từng tập hợp để tìm giá trị bí ẩn nhỏ nhất, và đếm số tập \(S\) có cùng giá trị đó.
Độ phức tạp: \(\mathcal{O} (q \times 2^n \times k^3)\).
Solution
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
const long long INF = 1e18;
int numQuery;
long long n, k;
long long dist[128][128];
int gcd(int a, int b) {
return b == 0 ? a : gcd(b, a % b);
}
int lcm(int x, int y) {
return x / gcd(x, y) * y;
}
int getBit(int num, int pos) {
return (num >> pos) & 1;
}
bool isSus(int x, int y) {
int v = x * y, a = 0, b = 0;
for (int u = 1; u <= v; ++u) {
if (v % u == 0) {
if (u % 2 != 0) {
++a;
} else {
++b;
}
}
}
return ((gcd(a, b) + lcm(a, b)) % 2) == 1;
}
long long calcSus(const vector<int> &S, int n) {
for (int u = 1; u <= n; ++u) {
for (int v = 1; v <= n; ++v) {
dist[u][v] = INF;
}
}
for (int u : S) {
for (int v : S) {
if (isSus(u, v)) {
dist[u][v] = 1;
}
}
}
for (int u : S) {
dist[u][u] = 0;
}
for (int k : S) {
for (int u : S) {
for (int v : S) {
dist[u][v] = min(dist[u][v], dist[u][k] + dist[k][v]);
}
}
}
long long result = 0;
for (int u : S) {
for (int v : S) {
result = max(result, dist[u][v]);
}
}
return result;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
freopen("SUSSET.inp", "r", stdin);
freopen("SUSSET.out", "w", stdout);
cin >> numQuery;
while (numQuery--) {
cin >> n >> k;
if (k > n) {
cout << -1 << " " << 0 << endl;
continue;
}
assert(n <= 10);
long long result = INF, cnt = 0;
for (int mask = 0; mask < (1 << n); ++mask) {
if (__builtin_popcount(mask) != k) {
continue;
}
vector<int> S;
for (int x = 1; x <= n; ++x) {
if (getBit(mask, x - 1)) {
S.push_back(x);
}
}
long long sus = calcSus(S, n);
if (result > sus) {
result = sus;
cnt = 1;
} else if (result == sus) {
cnt = cnt + 1;
}
}
cout << result % MOD << " " << cnt % MOD << '\n';
}
return 0;
}
Subtask \(2\) (\(10\%\) số điểm): \(n \leq 10^5\) và \(k = 2\)
Tutorial
Ta nhận thấy rằng cặp \((x,y)\) là cặp số bí ẩn khi và chỉ khi \(x\times y\) là một số chính phương. Thật vậy, ta có thể chứng minh như sau:
Vì \(1 \leq x, y \leq n\) nên dễ chứng minh được \(x \times y\) là một số nguyên dương.
Điều kiện đề yêu cầu \(GCD(a,b) + LCM(a,b)\) là số lẻ, mà 2 số đều là số nguyên nên \(GCD(a,b)\) và \(LCM(a,b)\) phải khác tính chắn lẻ.
Mà ta cũng có \(LCM(a,b)\ \vdots\ GCD(a,b)\) nên \(LCM(a,b)\) phải là số nguyên chẵn và \(GCD(a,b)\) là số nguyên lẻ.
Khi \(a, b\) cùng tính chẵn lẻ thì \(LCM(a, b)\) và \(GCD(a, b)\) cũng cùng tính chẵn lẻ, nên suy ra \(a\) và \(b\) khác tính chẵn lẻ và \(a + b\) là một số lẻ.
Với mọi số nguyên dương \(n\) thì với mỗi ước \(x\) của \(n\) ta cũng có một và chỉ một ước tương ứng là \(\frac{n}{x}\), nên số lượng ước của \(n\) không phải số chính phương là một số chẵn (vì chỉ có số chính phương dương \(n = p^2\) mới có một ước bị trùng là \(p = \frac{n}{p}\)).
Ta có \(a + b\) (cũng chính là tổng số lượng ước của \(x \times y\)) là một số lẻ nên \(x \times y\) phải là số chính phương.
Gọi \(x'=p_1^{x_1 \text{ mod } 2} \times p_2^{x_2 \text{ mod } 2} \times p_3^{x_3 \text{ mod } 2} \times \ldots\) với \(p_i\) là các số nguyên tố và \(x_i\) là các số tự nhiên.
Gọi \(y'=p_1^{y_1 \text{ mod } 2} \times p_2^{y_2 \text{ mod } 2} \times p_3^{y_3 \text{ mod } 2} \times \ldots\) với \(p_i\) là các số nguyên tố và \(y_i\) là các số tự nhiên.
Ta chứng minh được \(x \times y\) là số chính phương khi và chỉ khi \(x'=y'\).
Vì \(x\) và \(y\) là 2 số nguyên dương, ta luôn có cách phân tích:
\(\begin{cases}x = p_1^{x_1} \times p_2^{x_2} \times p_3^{x_3} \times \ldots \\ y = p_1^{y_1} \times p_2^{y_2} \times p_3^{y_3} \times \ldots \\ \end{cases}\) \(\ \ \ \ \ \text{}\) với \(p_i\) là các số nguyên tố và \(x_i, y_i\) là các số tự nhiên.
Để \(x \times y\) là số chính phương thì phải có số nguyên \(v = p_1^{v_1} \times p_2^{v_2} \times p_3^{v_3} \times \ldots\) (dãy \(v_i\) là các số tự nhiên) thoả \(v^2 = x \times y\) hay \(p_1^{2v_1} \times p_2^{2v_2} \times p_3^{2v_3} \times \ldots = p_1^{x_1+y_1} \times p_2^{x_2+y_2} \times p_3^{x_3+y_3} \times \ldots\)
Hay ứng với mỗi số nguyên tố \(p_i\) thi \(x_i + y_i = 2 \times v_i\) có nghĩa tổng \(x_i + y_i\) là một số nguyên chẵn (vì \(x_i, y_i, v_i \in \mathbb{N}\)), đồng nghĩa với \(x_i + y_i \equiv 0 \pmod 2\), nên ta có \(x_i\ \text{mod}\ 2 = y_i\ \text{mod}\ 2\) .
Vậy \(x \times y\) là số chính phương khi ứng với mỗi số nguyên tố \(p_i\) ta có \(x_i\ \text{mod}\ 2 = y_i\ \text{mod}\ 2\), mà vế phải cũng chính là định nghĩa của \(x' = y'\), vậy ta có điều phải chứng minh.
Với \(k=2\), ta cần đếm số lượng cặp số \(1\leq x,y \leq n\) nào mà \(x\times y\) là số chính phương hay không. Nếu số lượng này lớn hơn \(0\), thì độ bí ẩn nhỏ nhất sẽ là \(1\), ngược lại là \(10^{18}\)
Nhờ các nhận xét ở trên, ta sẽ có thuật toán như sau: Với mỗi số \(x\) thoả \(1 \leq x \leq n\), ta sẽ tìm \(x'\) tương ứng của nó bằng cách phân tích thừa số nguyên tố sử dụng sàng Eratosthene và đếm số lượng số \(x\) có cùng giá trị \(x'\)
Độ phức tạp: \(\mathcal{O}(n \log n)\)
Solution
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 100005;
const int MOD = 1e9 + 7;
int numQuery;
long long n, k;
int cnt[MAX_N];
int prime[MAX_N];
void eratosthenes(int n) {
for (int x = 0; x <= n; ++x) {
prime[x] = x;
}
for (int x = 2; x <= n; ++x) {
if (prime[x] == x) {
for (int v = 2; v * x <= n; ++v) {
prime[x * v] = x;
}
}
}
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
freopen("SUSSET.inp", "r", stdin);
freopen("SUSSET.out", "w", stdout);
eratosthenes(MAX_N);
cin >> numQuery;
while (numQuery--) {
cin >> n >> k;
if (k > n) {
cout << -1 << " " << 0 << '\n';
continue;
}
assert(k == 2);
long long result = 0;
for (int x = 1; x <= n; ++x) {
cnt[x] = 0;
int x0 = 1;
for (int t = x; t > 1;) {
int p = prime[t], f = 0;
for (; t % p == 0; t /= p) {
++f;
}
if (f % 2 == 1) {
x0 *= p;
}
}
result += cnt[x0]++;
}
cout << 1 << " " << result % MOD << '\n';
}
return 0;
}
Subtask \(3\) (\(20\%\) số điểm): \(n \leq 10^7\) và \(k = 3\)
Tutorial
Với \(n \leq 10^7\), việc phân tích thừa số nguyên tố trong \(log(n)\) sẽ không còn khả thi nữa (TLE :v), thay vào đó ta cần phải làm cách khác.
Mà ta thấy rằng mọi số nguyên dương \(x\) đều có thể phân tích thành dạng \(a^2\times b\) với \(a, b \in \mathbb{N}^*\) và \(a\) đạt max (\(b\) đạt min).
Theo cách phân tích ở trên, khi đặt \(b=(a')^2\times b'\) thì \(a' = 1\) luôn thoả.
Chứng minh
Giả sử tồn tại số \(x = a^2 \times b\) để \(a' > 1\).
Thế \(b = (a')^2\times b'\) vào \(x\) ta có \(x = a^2 \times ((a')^2 \times b') = (a \times a')^2 \times b'\).
Nếu \(a' \neq 1\) thì rõ ràng \((a \times a') > a\) là mâu thuẫn với yêu cầu \(a\) phải đạt max.
Ý tưởng của chúng ta sẽ là như sau, với mỗi \(b\), ta sẽ đánh dấu tất cả các số \(x\) mà phân tích thành dạng \(a^2\times b\) như trên. Khi đó, bất kỳ cặp số nào trong các \(x\) kể trên sẽ có tích là số chính phương.
Lưu ý là giờ thay vì chọn \(k = 2\) số thì ta phải chọn \(k = 3\), mà khi cố định \(b\) mình đã biết tất cả các \(x' = a\) nên khi gọi số cách chọn ứng với mỗi số \(b\) cố định là \(cnt_b\) thì kết quả bài toán là \(\Sigma cnt_b\).
Vì với mối số \(x = b \times a^2\) là biểu diễn duy nhất và \(x, b, a\) đều là các số nguyên trong đoạn \(1..n\) nên mỗi số \(x\) đều chỉ bị đánh dấu, duyệt và kiểm tra trong \(O(1)\) nên độ phức tạp sẽ là \(O(n)\) cho toàn bộ \(n\) số \(x = 1..n\).
Để tính số cách chọn \(3\) số trong một tập \(n\) số ta cần tính hàm tổ hợp \(nc3(n) = \frac{n!}{3!(n-3)!} = \frac{n(n-1)(n-2)}{6}\)
Ta không có phép chia modulo, Nhưng dưới modulo nguyên tố \(p = 10^9 + 7\) ta có \(\frac{1}{x} \equiv x^{-1} \equiv x^{p-2} \pmod p\) (theo định lí Fermat nhỏ)
Nên ta có \(nc3(n) \equiv n(n-1)(n-2) \times 6^{-1} \equiv n(n-1)(n-2) \times 6^{p-2} \equiv n(n-1)(n-2) \times 166666668 \pmod p\)
Độ phức tạp: \(\mathcal{O}(n)\)
Solution
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 1e7 + 17;
const int MOD = 1e9 + 7;
const int inv6 = 166666668;
int numQuery;
long long n, k;
bool isSquareFree[MAX_N];
int f(int x) {
return 1LL * x * (x - 1) % MOD * (x - 2) % MOD * inv6 % MOD;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
freopen("SUSSET.inp", "r", stdin);
freopen("SUSSET.out", "w", stdout);
cin >> numQuery;
while (numQuery--) {
cin >> n >> k;
if (k > n) {
cout << -1 << " " << 0 << endl;
continue;
}
if (k > sqrt(n)) {
cout << 49 << " " << f(n) << endl;
continue;
}
assert(k == 3);
long long result = 0;
for (int x = 1; x <= n; ++x) {
isSquareFree[x] = true;
}
for (int b = 1; b <= n; ++b) {
if (isSquareFree[b]) {
int cnt = 0;
for (int a = 1; 1LL * b * a * a <= n; ++a) {
isSquareFree[b * a * a] = false;
++cnt;
}
result += f(cnt);
}
}
cout << 1 << " " << result % MOD << endl;
}
return 0;
}
Subtask \(4\) (\(15\%\) số điểm): \(n,k \leq 10^5\)
Tutorial
Từ subtask này trở đi, ta sẽ cần phải quan tâm tới định nghĩa của độ bí ẩn.
Với \(k > 1\), ta nhận thấy tập hợp có độ bí ẩn nhỏ nhất sẽ là tập hợp có tất cả cặp phần tử thuộc tập hợp là cặp bí ẩn, hay nói cách khác đồ thị chúng ta có sẽ là đồ thị đầy đủ.
Từ đây, bài toán quy thành đếm số tập hợp \(k\) phần tử phân biệt từ \(1\) tới \(n\) mà tất cả các cặp phần tử đều có tích là số chính phương.
Ở subtask này, ta sẽ áp dụng y hệt thuật toán của subtask \(2\), chỉ khác là thay vì tính số cặp có cùng \(x'\), ta sẽ đếm số tập hợp \(k\) phần tử có cùng \(x'\) bằng công thức tổ hợp.
Ta cần tính \(C_n^k\ \text{mod}\ p\) với số modulo nguyên tố \(p = 10^9 + 7\) đủ nhanh cho \(n, k \leq 10^5\).
Thế thì ta có thể dùng công thức \(C_n^k \equiv \frac{n!}{k! \times (n-k)!} \equiv n! \times k!^{-1} \times (n-k)!^{-1} \pmod p\).
Mà theo định lí Fermat nhỏ ta có \(x^{-1} \equiv x^{p-2} \pmod p\) với modulo nguyên tố \(p\) bất kì.
Và ta cũng có thể tiền xử lí \(fact[i] = n!\ \text{mod}\ p\) bằng công thức \(fact[0] = 1\) và \(fact[n] = n \times fact[n - 1]\ \text{mod}\ p\).
Nên ta dễ tính được \(C_n^k \equiv fact[n] \times fact[k]^{p-2} \times fact[n-k]^{p-2} \pmod p\)
Độ phức tạp: \(\mathcal{O}(n \log (n))\)
Solution
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 100005;
const int MOD = 1e9 + 7;
int numQuery;
long long n, k;
int prime[MAX_N];
int cnt[MAX_N];
long long fct[MAX_N], ifct[MAX_N];
void eratosthenes(int n) {
for (int x = 0; x <= n; ++x) {
prime[x] = x;
}
for (int x = 2; x <= n; ++x) {
if (prime[x] == x) {
for (int v = 2; v * x <= n; ++v) {
prime[x * v] = x;
}
}
}
}
long long binPow(long long a, long long b) {
long long result = 1;
while (b) {
if (b & 1) {
result = result * a % MOD;
}
b >>= 1;
a = a * a % MOD;
}
return result;
}
void initialize() {
fct[0] = ifct[0] = 1;
for (int i = 1; i < MAX_N; i++) {
fct[i] = (fct[i - 1] * i) % MOD;
ifct[i] = binPow(i, MOD - 2);
}
}
long long nCk(int n, int k) {
k = min(k, n - k);
if (k < 0) {
return 0;
}
return (((fct[n] * ifct[k]) % MOD) * ifct[n - k]) % MOD;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
freopen("SUSSET.inp", "r", stdin);
freopen("SUSSET.out", "w", stdout);
initialize();
eratosthenes(MAX_N);
cin >> numQuery;
while (numQuery--) {
cin >> n >> k;
if (k > n) {
cout << -1 << " " << 0 << '\n';
continue;
}
if (k == 1) {
cout << 0 << " " << n % MOD << '\n';
continue;
}
if (k * k > n) {
cout << 49 << " " << nCk(n, k) % MOD << '\n';
continue;
}
assert(n <= 1e5 && k <= 1e5);
long long res = 0;
for (int x = 1; x <= n; ++x) {
cnt[x] = 0;
int x0 = 1;
for (int t = x; t > 1;) {
int p = prime[t], f = 0;
for (; t % p == 0; t /= p) {
++f;
}
if (f % 2 == 1) {
x0 *= p;
}
}
res += nCk(cnt[x0]++, k - 1);
}
cout << 1 << " " << res % MOD << '\n';
}
return 0;
}
Subtask \(5\) (\(15\%\) số điểm): \(n,k \leq 10^7\)
Tutorial
Ý tưởng sẽ tương tự subtask 3, chỉ khác là lần này ta sẽ cần phải khử log trong quá trình tính tổ hợp bằng tiền xử lí.
Ta sẽ tiền xử lý để tính nghịch đảo modulo của \(fact[i] = i!\ \text{mod}\ p\) và \(tcaf[i] = fact[i]^{-1}\ \text{mod}\ p\) trong \(O(n + \log p)\) như sau:
const int MAX_N = 1e7;
...
fact[0] = 1;
for (int i = 1; i <= MAX_N; i++) {
fact[i] = 1LL * fact[i - 1] * i % p;
}
tcaf[MAX_N] = power(fact[MAX_N], p - 2);
for (int i = MAX_N - 1; i >= 0; i--) {
tcaf[i] = 1LL * tcaf[i + 1] * (i + 1) % p;
}
...
Ngoài ra cũng có cách khác làm trong \(O(n)\) với \(p\) nguyên tố, đó là
/// Precalculating factorials under prime modulo
int fact[LIM + 10]; /// fact[n] = n!
int invs[LIM + 10]; /// invs[n] = n^(-1)
int tcaf[LIM + 10]; /// tcaf[n] = (n!)^(-1)
void precal_nck(int n = LIM) {
fact[0] = fact[1] = 1;
invs[0] = invs[1] = 1;
tcaf[0] = tcaf[1] = 1;
for (int i = 2; i <= n; ++i) {
fact[i] = (1LL * fact[i - 1] * i) % MOD;
invs[i] = MOD - 1LL * (MOD / i) * invs[MOD % i] % MOD;
tcaf[i] = (1LL * tcaf[i - 1] * invs[i]) % MOD;
}
}
Ghi chú
Ta có \(p = \left \lfloor \frac{p}{x} \right \rfloor \times x + (p \text{ mod } x)\)
Nên ta có \(-\left \lfloor \frac{p}{x} \right \rfloor \times x \equiv (p \text{ mod } x) \pmod p\)
Nhân cả hai vế cho \(x^{-1} \times (p \text{ mod } x)^{-1}\) ta được \(-\left \lfloor \frac{p}{x} \right \rfloor \times (p\ \text{ mod }\ x)^{-1} \equiv x^{-1} \pmod p\)
Vậy khi đặt hàm \(inv[x] = x^{-1}\ \text{ mod }\ p\) ta có \(inv[x] = (-\left \lfloor \frac{p}{x} \right \rfloor \times inv[p\ \text{ mod }\ x])\ \text{ mod }\ p\)
Và ta có \(tcaf[n] = n!^{-1} = (inv[1] \times inv[2] \times \ldots \times inv[n])\ \text{ mod }\ p\)
Độ phức tạp: \(\mathcal{O}(n)\)
Solution
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 10000005;
const int MOD = 1e9 + 7;
int numQuery = 1;
long long n, k;
bool isSquareFree[MAX_N];
long long fct[MAX_N], ifct[MAX_N];
long long binPow(long long a, long long b) {
long long result = 1;
while (b) {
if (b & 1) {
result = result * a % MOD;
}
b >>= 1;
a = a * a % MOD;
}
return result;
}
void initialize() {
fct[0] = 1;
for (int i = 1; i < MAX_N; i++) {
fct[i] = (fct[i - 1] * i) % MOD;
}
ifct[MAX_N - 1] = binPow(fct[MAX_N - 1], MOD - 2);
for (int i = MAX_N - 2; i >= 0; i--) {
ifct[i] = (ifct[i + 1] * (i + 1)) % MOD;
}
}
long long nCk(int n, int k) {
k = min(k, n - k);
if (k < 0) {
return 0;
}
return (((fct[n] * ifct[k]) % MOD) * ifct[n - k]) % MOD;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
freopen("SUSSET.inp", "r", stdin);
freopen("SUSSET.out", "w", stdout);
initialize();
cin >> numQuery;
while (numQuery--) {
cin >> n >> k;
if (k > n) {
cout << -1 << " " << 0 << '\n';
continue;
}
if (k == 1) {
cout << 0 << " " << n % MOD << '\n';
continue;
}
if (k * k > n) {
cout << 49 << " " << nCk(n, k) % MOD << '\n';
continue;
}
assert(n <= 1e7 && k <= 1e7);
long long res = 0;
for (int x = 1; x <= n; ++x) {
isSquareFree[x] = true;
}
for (int b = 1; b <= n; ++b) {
if (isSquareFree[b]) {
int cnt = 0;
for (int a = 1; 1LL * b * a * a <= n; ++a) {
isSquareFree[b * a * a] = false;
++cnt;
}
res += nCk(cnt, k);
}
}
cout << 1 << " " << res % MOD << '\n';
}
return 0;
}
Subtask \(6\) (\(30\%\) số điểm): Không có ràng buộc gì thêm.
Tutorial
Xét các trường hợp \(k = 1, k > n > 0, n \geq k > \sqrt{n}, \sqrt{n} \geq k \geq 2\)
Trường hợp \(k = 1\) thì giá trị bí ẩn là \(0\)
Trường hợp \(k > n\) thì giá trị bí ẩn là \(-1\)
Trường hợp \(n \geq k > \sqrt{n}\) thì giá trị bí ẩn là \(10^{18}\ \text{mod}\ (10^9 + 7)\)
Trường hợp \(\sqrt{n} \geq k \geq 2\) thì thay vì cố định \(a\) để đếm \(b^2\) với \(a = 1 \ldots n\), chúng ta sẽ làm ngược lại bằng cách cố định \(b^2\) để đếm số \(a\) với \(b = 1 \ldots \sqrt{n}\)
Từ sự khác biệt của subtask \(4, 5\) với subtask \(2, 3\) thì ta biết được bài toán có thể chia làm \(3\) phần riêng biệt là phần tiền xử lí, tính toán tổ hợp và tính toán kết quả.
Nên ta có thể xử từng phần cho tối ưu nhất và lắp ghép cho phù hợp là xong.
Trường hợp \(k = 1\): Kết quả bài toán sẽ là \((0, max(0, n)\ \text{mod}\ p)\) (với số modulo nguyên tố \(p = 10^9 + 7\))
Trường hợp \(k > n > 0\): Kết quả bài toán sẽ là \((-1, 0)\)
Trường hợp \(n \geq k > \sqrt{n}\)
Kết quả bài toán sẽ là \((10^{18}\ \text{mod}\ p, C_n^k\ \text{mod}\ p)\) (với số modulo nguyên tố \(p = 10^9 + 7\))
Vậy đoạn quan trọng là mình cần tính \(\binom{n}{k}\) modulo cho số nguyên tố \(p = 10^9 + 7\).
Dựa theo Lucas Theorem ta có \(\binom{n}{k} = \binom{n\ div\ p}{k\ div\ p} + \binom{n\ \text{mod}\ p}{k\ \text{mod}\ p}\) (vì \(p^2 < 10^{12}\)) với div là phép chia lấy nguyên và mod là phép chia lấy dư.
Mà vì \(n\ div\ p, n\ \text{mod}\ p, k\ div\ p, k\ \text{mod}\ p\) đều nhỏ hơn \(p = 10^9 + 7\), nên ta cần có cách để tính \(x!\) với \(x < p\) đủ nhanh.
Một trong các cách đó chính là dùng một mảng \(FACT[]\) có \(10^3\) phần tử lưu giá trị của \(0e6!, 1e6!, 2e6!, 3e6!, \ldots, 1e9!\) (có thể tính toán từ bên ngoài trước rồi copy lại vô).
Khi đó, với bất cứ giá trị \(x!\) nào, ta có thể tìm số \(FACT[i]\) gần với \(x!\) nhất và tính bằng công thức \(x! = x \times (x - 1)! = x \times (x - 1) \times (x - 2)! = x \times (x - 1) \times (x - 2) \times \ldots \times FACT[i]\).
Cụ thể hơn là trong code chứa \(O\left(\frac{p}{X}\right)\) số dạng \((0X)!, (1X)!, (2X)!, (3X)!, \ldots, \left(\left \lfloor \frac{p}{X} \right \rfloor * X\right)!\) và mỗi truy vấn mất \(O(X)\) thời gian để tiền xử lí.
Trong đó \(p\) là hằng số modulo nguyên tố và \(X\) chính là hằng số khoảng cách giữa các giai thừa đã tiền xử lí.
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
cout << "{";
int fact;
for (int i = 0; i <= int(1e9); i++) {
if (i == 0) {
fact = 1;
} else {
fact = fact * (long long)i % MOD;
}
if (i % int(1e6) == 0) {
cout << fact;
if (i < int(1e9)) {
cout << ", ";
}
}
}
cout << "}";
return 0;
}
long long mulMOD(const long long x, const long long y) { /// Tính x * y % MOD
long long res = x * y;
if (abs(res) >= MOD) {
res %= MOD;
}
if (res < 0) {
res += MOD;
}
return res;
}
long long fix(long long res) { /// Đưa res về [0, MOD) mà không dùng phép nhân chia
while (res >= MOD) {
res -= MOD;
}
while (res < 0) {
res += MOD;
}
return res;
}
/// Tiền xử lí tính toán trong O(n)
long long fact[SQRT_LIM + 10]; /// fact[n] = n!
long long invs[SQRT_LIM + 10]; /// invs[n] = n^(-1)
long long tcaf[SQRT_LIM + 10]; /// tcaf[n] = (n!)^(-1)
void precal_nck(int n = SQRT_LIM) {
fact[0] = fact[1] = 1;
invs[0] = invs[1] = 1;
tcaf[0] = tcaf[1] = 1;
for (int i = 2; i <= n; ++i) {
fact[i] = mulMOD(fact[i - 1], i);
invs[i] = fix(MOD - mulMOD(MOD / i, invs[MOD % i]));
tcaf[i] = mulMOD(tcaf[i - 1], invs[i]);
}
}
/// Truy vấn cho nCk % p khi 0 <= n, k < p
long long nck(long long n, long long k) {
return (k < 0 || k > n) ? 0LL : mulMOD(fact[n], mulMOD(tcaf[k], tcaf[n - k]));
}
/// Tiền xử lí 0e6!, 1e6!, 2e6!, 3e6!, ..., 1e9! (dưới mod p)
const int staticFact[] = /* Dãy mà chương trình trên in ra */;
/// Tính n! trong O(1e6)
int longFactorial(int n) {
long long res = staticFact[n / int(1e6)];
for (long long i = n / int(1e6) * int(1e6) + 1; i <= n; ++i) {
if ((res *= i) >= MOD) {
res %= MOD;
}
}
return res;
}
/// Tính nghịch đảo modulo | O(log x / log log x) -> O(x^(1/3)) for general p | O(log x) for p = 1e9 + 7
int inverse(int x) {
if (x <= 1) {
return 1;
}
return (MOD - MOD / x) * (long long)inverse(MOD % x) % MOD;
}
/// Tính nck % p cho trường hợp 0 <= n, k < p^2
int lnck(long long n, long long k) {
if (n < 0 || k < 0 || n < k) {
return 0;
}
if (n >= MOD || k >= MOD) {
return (long long)lnck(n / MOD, k / MOD) * lnck(n % MOD, k % MOD) % MOD;
}
return longFactorial(n) * (long long)inverse(longFactorial(n - k) * (long long)longFactorial(k) % MOD) % MOD;
}
Trường hợp \(k = 2\)
Ta đặt các định nghĩa sau
Ta gọi \(s(x) = n - 1\) chính là số cách chọn cặp số \(b^2\) với một giá trị \(a\) cố định (mặc dù chúng ta đang cố định \(b^2\) để tính toán).
Ta gọi \(f(x)\) chính là số lượng bộ số \((a, b, x)\) thoả \(1 \leq a < a \leq n\) và \((a, b, x)\) tạo thành dãy cấp số nhân.
Ta gọi \(g(x)\) chính là số lượng bộ số \((a, b, x)\) thoả \(1 \leq a \leq b \leq n\) và \((a, b, x)\) tạo thành dãy cấp số nhân.
Lúc này ta sẽ có kết quả bài toán chính là \(F(n) = \underset{p=1}{\overset{n}{\Large \Sigma}} f(p)\)
Chứng minh F(n) là kết quả bài toán
Với mỗi \(f(x)\) thì \(a, b, x\) tạo thành cấp số nhân nên \(b^2 = a \times x\) (vì $(a, b, x) \Leftrightarrow \((a, a \times z, a \times z^2)\) với hằng số \(z\) là công sai)
Mà khi mình tính \(F(n)\), tức là mình đã cố định số \(x\) với \(x = 1 \rightarrow n\).
Khi này chỉ có \(1\) giá trị không phải hằng, đó là số \(a\), và \(F(n)\) chính là đếm số số \(a\) thoả mãn \(a \times x = b^2\) là một số chính phương.
Tại sao cần phải tính \(F(n)\) với lại cố định \(b^2\) thay vì cố định \(a\) ?
Nếu như bình thường, thì chúng ta sẽ có thuật tuyến tính khi cố định \(a\) như sau
int solve(int n) {
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
long long res = 0;
for (int i = 1, j; i <= n; ++i) {
if (is_squarefree[i]) {
for (j = 1; i * j * j <= n; ++j) {
is_squarefree[i * j * j] = false;
}
res += 1LL * (j - 1) * (j - 2) / 2;
}
}
res %= MOD;
return res;
}
Nếu làm thế thì độ phức tạp sẽ phải là \(\Omega(n)\) vì phải duyệt \(a = 1 \rightarrow n\), tức là ta không thể tính toán nhanh hơn \(O(n)\) được.
Nhưng sau khi ta hoán vị lại \(2\) vòng for, tức là cố định \(b^2\) lại
int solve(int n) {
memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
long long res = 0;
int t = sqrt(n);
while (t * t < n) {
++t;
}
while (t * t > n) {
--t;
}
for (int j = t; j > 1; --j) {
for (int i = 1; i * j * j <= n; ++i) {
if (!used[i * j * j]) {
used[i * j * j] = true;
res += j - 1;
}
}
}
res %= MOD;
return res;
}
Thì ta có độ phức tap của xử lí truy vấn là \(\Omega(\sqrt{n})\) vì chỉ cần phải xét các \(b^2 = \{1, 4, \ldots, \left \lfloor \sqrt{n} \right \rfloor \}\) hay \(O(b) = O(\sqrt{n})\)
Tính toán \(F(n)\) bằng công thức tốt hơn
\(F(N) = \underset{p=1}{\overset{n}{\Large \Sigma}} f(p) = \underset{p=1}{\overset{n}{\Large \Sigma}} g(p) - n = \underset{p=1}{\overset{n}{\Large \Sigma}} \underset{d^2 | p}{\Large \Sigma} \phi(d) - n\)
Gọi hàm Iverson \([x = y] = 1\) khi mệnh đề (\(x = y\)) đúng và \([x = y] = 0\) khi (\(x \neq y\))
Ta có \(F(N) = \underset{p=1}{\overset{n}{\Large \Sigma}} \phi(p) \times \underset{d=1}{\overset{n}{\Large \Sigma}} [u \times p^2 \leq n] = \underset{p=1}{\overset{n}{\Large \Sigma}} \phi(p) \times \left \lfloor \frac{n}{p^2} \right \rfloor\)
Dễ thấy \(\left \lfloor \frac{n}{p^2} \right \rfloor = 0\ \forall\ p > \sqrt{n}\) nên ta có \(F(N) = \underset{p=1}{\overset{\left \lfloor \sqrt{n} \right \rfloor}{\Large \Sigma}} \phi(p) \times \left \lfloor \frac{n}{p^2} \right \rfloor\)
Thuật toán \(O(\sqrt{n})\)
Vì ta có thể tính \(\phi(x)\) với \(x = 1 \rightarrow \left \lfloor \sqrt{n} \right \rfloor\) với thuật sàng trong \(O(\sqrt{n})\) tiền xử lí và \(O(1)\) truy vấn nên \(F(N)\) có thể tính trong \(O(\sqrt n)\)
Để làm thế ta có đoạn code sau:
Độ phức tạp: \(\mathcal{O}(\sqrt{n} \log \log \sqrt{n})\)
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
const int MAX_N = 1e7 + 17;
const int SQRT_MAX_N = 3200;
int euler[SQRT_MAX_N];
void sieve_phi(int n) {
iota(euler, euler + n + 1, 0);
for (int x = 2; x <= n; x++)
if (euler[x] == x) {
for (int j = x; j <= n; j += x) {
euler[j] -= euler[j] / x;
}
}
}
int solve(int n) {
sieve_phi(ceil(sqrt(n) + 1) + 1);
long long res = 0;
for (int p = 2; p * p <= n; ++p) {
res += 1LL * euler[p] * (n / (p * p));
}
res %= MOD;
return res;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
int n;
cin >> n;
cout << solve(n);
return 0;
}
Độ phức tạp: \(\mathcal{O}(\sqrt{n})\)
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
vector<int> lpf;
vector<int> prime;
vector<int> euler;
void linear_sieve_phi(int n) {
lpf.assign(n + 1, 0);
euler.assign(n + 1, 1);
for (int x = 2; x <= n; ++x) {
if (lpf[x] == 0) {
prime.push_back(lpf[x] = x);
euler[x] = x - 1;
}
for (int i = 0; i < prime.size() && x * prime[i] <= n; ++i) {
lpf[x * prime[i]] = prime[i];
if (x % prime[i] == 0) {
euler[x * prime[i]] = euler[x] * prime[i];
break;
}
euler[x * prime[i]] = euler[x] * euler[prime[i]];
}
}
}
int solve(int n) {
linear_sieve_phi(ceil(sqrt(n) + 1) + 1);
long long res = 0;
for (int p = 2; p * p <= n; ++p) {
res += 1LL * euler[p] * (n / (p * p));
}
res %= MOD;
return res;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
int n;
cin >> n;
cout << solve(n);
return 0;
}
Trường hợp \(\sqrt{n} \geq k \geq 2\)
Ta gọi \(s_k(n) = C^{n-1}_{k-1}\) chính là số cách chọn bộ số \(k\) số \(b^2\) với một giá trị \(a\) cố định (mặc dù chúng ta đang cố định \(b^2\) để tính toán).
Ta gọi \(f_k(n)\) chính là số lượng tập số \((a_1, a_2, \dots, a_k, n)\) thoả \(1 \leq a_1 < a_2 < \dots < a_k \leq n\) và \((a_1, a_2, \dots, a_k, n)\) tạo thành dãy cấp số nhân bậc \((k+1)\).
Ta gọi \(g_k(n)\) chính là số lượng tập số \((a_1, a_2, \dots, a_k, n)\) thoả \(1 \leq a_1 \leq a_2 \leq \dots \leq a_k \leq n\) và \((a_1, a_2, \dots, a_k, n)\) tạo thành dãy cấp số nhân bậc \((k+1)\).
Với từng số \(a = 1..n\) ta cố định, ta đã duyệt hết toàn bộ các cặp \(b^2\), nên ta cũng có công thức khi cố định \(b^2\) đó là \(F_k(n) = \overset{n}{\underset{p=1}{\Large \Sigma}} f_k(p)\).
Cách chứng minh tương tự với trường hợp \(k = 2\) ở trên nhưng tổng quát hoá ra cho \(k\) số.
Ta có \(f_k(n) = F_k(n) - F_k(n-1) = s_k(g_k(n)) = \underset{d^2 | n}{\Large \Sigma} (\mu \times s_k)(d) = \underset{d^2 | n}{\Large \Sigma} \underset{p | d}{\Large \Sigma} \mu(p) \times s_k\left(\frac{d}{p}\right)\)
Mà \(F_k(n) = \overset{n}{\underset{p=1}{\Large \Sigma}} f_k(p) = F_k(0) + \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{d = 1}{\Large \Sigma}}(\mu \times s_k)(d) \times \left \lfloor \frac{n}{d^2} \right \rfloor = \overset{\left \lfloor \sqrt{n} \right \rfloor}{\underset{d = 1}{\Large \Sigma}} \left \lfloor \frac{n}{d^2} \right \rfloor \times \left ( \underset{p | d}{\Large \Sigma} \mu(p) \times s_k\left(\frac{d}{p}\right) \right)\)
Để tính \(F(n)\) thì chúng ta có thể dùng tổ hợp bao hàm loại trừ với quy hoạch động trong \(\log^2 n\), hoặc Dirichlet Convolution giữa hàm \(\mu(p)\) với \(s_k\left(\frac{d}{p}\right)\) xuống thành \(\log n\), hoặc là dùng Mobius Inversion xuông thành \(\log \log n\)
Độ phức tạp: \(\mathcal{O}(\sqrt{n} \log n)\)
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 5e6 + 56;
const int SQRT_MAX_N = 2300;
const int MOD = 1e9 + 7;
int fact[SQRT_MAX_N + 10];
int invs[SQRT_MAX_N + 10];
int tcaf[SQRT_MAX_N + 10];
void precal_nck(int n = SQRT_MAX_N) {
fact[0] = fact[1] = 1;
invs[0] = invs[1] = 1;
tcaf[0] = tcaf[1] = 1;
for (int i = 2; i <= n; ++i) {
fact[i] = (1LL * fact[i - 1] * i) % MOD;
invs[i] = MOD - 1LL * (MOD / i) * invs[MOD % i] % MOD;
tcaf[i] = (1LL * tcaf[i - 1] * invs[i]) % MOD;
}
}
int nck(int n, int k) {
k = min(k, n - k);
if (k < 0) {
return 0;
}
long long res = fact[n];
res = res * tcaf[k] % MOD;
res = res * tcaf[n - k] % MOD;
return res;
}
/// Linear Sieve
vector<int> prime; /// prime list = A000040
bool isPrime[SQRT_MAX_N + 10]; /// characteristic function = A010051
int lpf[SQRT_MAX_N + 10]; /// lowest prime factor = A020639
int mu[SQRT_MAX_N + 10]; /// mobius = A008683
void linear_sieve(int n) {
if (n < 1) {
return;
}
/// Extension Sieve || You can add something more
memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
fill_n(mu, n + 1, 1);
/// Main Sieve || Without this, you barely able to achive linear complexity
prime.clear();
prime.reserve(n / log(n - 1));
memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
isPrime[0] = isPrime[1] = false;
for (int x = 2; x <= n; ++x) { /// For each number
if (isPrime[x]) { /// Func[Prime]
mu[x] = -1;
lpf[x] = x;
prime.push_back(x);
}
for (int p : prime) { /// Func[Prime * X] <- Func[Prime]
if (p > lpf[x] || x * p > n) {
break;
}
isPrime[x * p] = 0;
lpf[x * p] = p;
mu[x * p] = (lpf[x] == p) ? 0 : -mu[x];
}
}
}
/// Divisor sieve
vector<int> divisors[SQRT_MAX_N];
void precal_div(int n) { /// O(n log n)
for (int u = n; u >= 1; --u) {
divisors[u].clear();
for (int v = u; v <= n; v += u) {
divisors[v].push_back(u);
}
}
}
long long solve(int n, int k) {
int t = ceil(sqrt(n) + 1) + 1;
linear_sieve(t);
precal_nck(t);
precal_div(t);
long long res = 0;
for (int d = 1; d * d <= n; ++d) {
long long sum = 0;
for (int p : divisors[d]) {
sum += mu[d / p] * nck(p - 1, k - 1);
}
sum %= MOD;
res += sum * (n / (d * d));
}
res %= MOD;
return res;
}
Độ phức tạp: \(\mathcal{O}(\sqrt{n} \log \log \sqrt{n})\)
vector<int> prime;
bool isPrime[SQRT_MAX_N + 10];
int lpf[SQRT_MAX_N + 10];
void linear_sieve(int n) {
if (n < 1) return;
prime.clear();
prime.reserve(n / log(n - 1));
memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
isPrime[0] = isPrime[1] = false;
for (int x = 2; x <= n; ++x) {
if (isPrime[x]) {
lpf[x] = x;
prime.push_back(x);
}
for (int p : prime) {
if (p > lpf[x] || x * p > n) {
break;
}
isPrime[x * p] = 0;
lpf[x * p] = p;
}
}
}
long long res[SQRT_MAX_N + 10];
long long solve(int n, int k) {
int t = ceil(sqrt(n) + 1) + 1;
linear_sieve(t);
precal_nck(t);
memset(res, 0, sizeof(res[0]) * (t + 1));
for (int d = 1; d * d <= n; ++d) {
res[d] = nck(d - 1, k - 1);
}
for (int p : prime) {
for (int d = t / p; d > 0; --d) {
res[d * p] -= res[d];
}
}
long long ans = 0;
for (int d = 1; d * d <= n; ++d) {
ans += res[d] * (n / (d * d));
}
ans %= MOD;
return ans;
}
Độ phức tạp: \(\mathcal{O}\left(\sqrt{\left \lfloor \frac{n}{k} \right \rfloor} \log \log \sqrt{\left \lfloor \frac{n}{k} \right \rfloor} - k \right)\)
vector<int> valid;
int cnt[SQRT_MAX_N];
bool is_squarefree[MAX_N];
long long res[SQRT_MAX_N + 10];
int solve(int n, int k) {
int t = ceil(sqrt(n) + 0.5);
if (k > t) {
return 0;
}
linear_sieve(t);
precal_nck(t);
memset(res, 0, sizeof(res[0]) * (t - k + 1));
for (int d = k; d * d <= n; ++d) {
res[d - k] = nck(d - 1, k - 1);
}
for (int p : prime) {
for (int d = t / p; d >= k; --d) {
res[d * p - k] -= res[d];
}
}
long long ans = 0;
for (int d = k; d <= t; ++d) {
ans += res[d - k] * (n / (d * d));
}
ans %= MOD;
return ans;
}
Độ phức tạp: \(\mathcal{O}\left(\sqrt{\left \lfloor \frac{n}{k} \right \rfloor} \log \log \sqrt{\left \lfloor \frac{n}{k} \right \rfloor} - k + 10^6 \right)\)
Solution
#include <bits/stdc++.h>
using namespace std;
const long long MAX_N = 1e12 + 5;
const int SQRT_MAX_N = 1000005;
const int MOD = 1e9 + 7;
const int staticFact[] = {...};
int numQuery;
vector<pair<long long, long long>> queries;
long long fact[SQRT_MAX_N + 10];
long long invs[SQRT_MAX_N + 10];
long long tcaf[SQRT_MAX_N + 10];
vector<int> prime;
bool isPrime[SQRT_MAX_N + 10];
int lpf[SQRT_MAX_N + 10];
int mu[SQRT_MAX_N + 10];
long long res[SQRT_MAX_N + 10];
long long mulMOD(const long long x, const long long y) {
long long result = x * y % MOD;
if (result < 0) {
result += MOD;
}
return result;
}
long long fix(long long result) {
result %= MOD;
while (result < 0) {
result += MOD;
}
return result;
}
void initialize(int n) {
fact[0] = fact[1] = 1;
invs[0] = invs[1] = 1;
tcaf[0] = tcaf[1] = 1;
for (int i = 2; i <= n; ++i) {
fact[i] = mulMOD(fact[i - 1], i);
invs[i] = fix(MOD - mulMOD(MOD / i, invs[MOD % i]));
tcaf[i] = mulMOD(tcaf[i - 1], invs[i]);
}
}
long long nck(long long n, long long k) {
return (k < 0 || k > n) ? 0LL : mulMOD(fact[n], mulMOD(tcaf[k], tcaf[n - k]));
}
void linear_sieve(int n = SQRT_MAX_N) {
if (n < 1) {
return;
}
memset(lpf, 0, sizeof(lpf[0]) * (n + 1));
fill_n(mu, n + 1, 1);
prime.clear();
prime.reserve(n / log(n - 1) + 1);
memset(isPrime, true, sizeof(isPrime[0]) * (n + 1));
isPrime[0] = isPrime[1] = false;
for (int x = 2; x <= n; ++x) {
if (isPrime[x]) {
mu[x] = -1;
lpf[x] = x;
prime.push_back(x);
}
for (int p : prime) {
if (p > lpf[x] || x * p > n) {
break;
}
isPrime[x * p] = 0;
lpf[x * p] = p;
mu[x * p] = (lpf[x] == p) ? 0 : -mu[x];
}
}
}
int longFactorial(int n) {
long long res = staticFact[n / int(1e6)];
for (long long i = n / int(1e6) * int(1e6) + 1; i <= n; ++i) {
if ((res *= i) >= MOD) {
res %= MOD;
}
}
return res;
}
int inverse(int x) {
if (x <= 1) {
return 1;
}
return (MOD - MOD / x) * (long long)inverse(MOD % x) % MOD;
}
int lnck(long long n, long long k) {
if (n < 0 || k < 0 || n < k) {
return 0;
}
if (n >= MOD || k >= MOD) {
return (long long)lnck(n / MOD, k / MOD) * lnck(n % MOD, k % MOD) % MOD;
}
return longFactorial(n) * (long long)inverse(longFactorial(n - k) * (long long)longFactorial(k) % MOD) % MOD;
}
pair<int, int> solve(long long n, long long k) {
if (k == 1) {
return make_pair(0, n %= MOD);
}
if (k > n) {
return make_pair(-1, 0);
}
int t = sqrt(n);
while (1LL * t * t < n) {
++t;
}
while (1LL * t * t > n) {
--t;
}
if (k > t) {
return make_pair(49, lnck(n, k));
}
memset(res + k, 0, sizeof(res[0]) * (t + 1));
for (int d = 1; 1LL * d * d <= n; ++d) {
res[d] = nck(d - 1, k - 1);
}
for (int p : prime) {
for (int d = t / p; d >= k; --d) {
res[d * p] -= res[d];
}
}
long long ans = 0;
for (int d = 1; 1LL * d * d <= n; ++d) {
ans += mulMOD(res[d], n / (1LL * d * d));
while (ans >= MOD) {
ans -= MOD;
}
}
return make_pair(1, ans %= MOD);
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
freopen("SUSSET.inp", "r", stdin);
freopen("SUSSET.out", "w", stdout);
cin >> numQuery;
queries.resize(numQuery);
long long t = 0;
for (auto &query : queries) {
cin >> query.first >> query.second;
t = max(t, query.first);
}
linear_sieve(sqrt(t));
initialize(sqrt(t));
for (pair<long long, long long> &query : queries) {
pair<int, int> result = solve(query.first, query.second);
cout << result.first << " " << result.second << '\n';
}
return 0;
}
Bình luận