First, the Farey sequence $F_n$ represents the set of all completely reduced proper fractions with denominators not exceeding $n$, arranged in ascending order of value. Formally speaking,↵
↵
$$F_n = \left\{\frac{p}{q} \bigg\vert 0 < p < q \le n, \, \gcd(p, q) = 1\right\}$$↵
↵
The asymptotic size of this sequence is $\sum\limits_{i = 1}^{n} \varphi(n) \sim \frac{3}{\pi^2}n^2 + o(n \log n)$.↵
↵
Regarding the Farey sequence, there are two classic conclusions:↵
↵
- If $\frac{a}{b}$ and $\frac{c}{d}$ are two consecutive elements in $F_n$, then $\frac{a + c}{b + d}$ is also a valid Farey fraction.↵
↵
- If $\frac{a}{b}$ and $\frac{c}{d}$ are two consecutive elements in $F_n$, then their next element $\frac{p}{q}$ satisfies↵
$$\begin{cases}p = \left\lfloor\frac{n + b}{d}\right\rfloor c - a \\ q = \left\lfloor\frac{n + b}{d}\right\rfloor d - b\end{cases}$$↵
↵
We will only prove the second theorem above. For the Farey sequence of order $n$, two adjacent terms $\frac{a}{b} < \frac{c}{d}$ satisfy $bc - ad = 1$. Since the adjacent terms $\frac{a}{b}, \, \frac{c}{d}$ are both on the Stern-Brocot tree, the next term $\frac{p}{q}$ satisfies $\frac{c}{d} = \frac{a + p}{b + q}$, which means $(p + a)d = (q + b)c$. Therefore, there exists $k$ such that↵
↵
$$\begin{cases} kc = a + p \\ kd = b + q \end{cases}$$↵
↵
To make $p, \, q$ highly precise (i.e., closely packed), the value of $k$ should be as large as possible. Therefore, concerning the upper bound $n$, $p, \, q$ are defined by the following equations:↵
↵
$$\begin{cases}p = kc - a \le n \\ q = kd - b \le n\end{cases}$$↵
↵
Obviously, from the definition of a proper fraction, we know $q > p$, so the maximum $k = \left\lfloor\frac{n + b}{d}\right\rfloor$, thus↵
↵
$$\begin{cases}p = \left\lfloor\frac{n + b}{d}\right\rfloor c - a \\ q = \left\lfloor\frac{n + b}{d}\right\rfloor d - b\end{cases}$$↵
↵
There are two computational methods to generate this sequence. The first one is based on the Stern-Brocot tree from the first theorem, which is very easy to implement.↵
↵
```cpp↵
int build(int a, int b, int c, int d, int n) {↵
if (b + d > n) return 0;↵
return 1 + build(a, b, a + c, b + d, n) + build(a + c, b + d, c, d, n);↵
}↵
```↵
↵
The other one is based on the recurrence from the second theorem, which is also quite easy to write.↵
↵
```cpp↵
PII nxt_fraction(PII fac1, PII fac2) {↵
auto [a, b] = fac1;↵
auto [c, d] = fac2;↵
return {(n + b) / d * c - a, (n + b) / d * d - b};↵
}↵
```↵
↵
Obviously, they both use an $O(n^2)$ time complexity to compute specific elements of the Farey sequence.↵
↵
To facilitate solving the parent problem, we present the following two problems:↵
↵
> Given $n, \, k$, find the $k$-th term of $F_n$.↵
↵
> Given $n$ and a completely reduced proper fraction $\frac{p}{q} (p \le n)$, determine its rank in $F_n$.↵
↵
Let's consider the first problem first. We find that the first problem can be processed through the following steps:↵
↵
- Binary search for $j$ in $\frac{j}{n}$, and calculate $\text{rank}(\frac{j}{n})$↵
↵
- Let $r = \text{rank}(\frac{j}{n})$. If $r < k$, search upwards; otherwise, search downwards.↵
↵
- Find an interval $\left[\frac{j}{n}, \, \frac{j + 1}{n}\right)$ such that the target fraction falls within this range.↵
↵
- Count the fractions within this interval and find a valid fraction such that its rank equals the given $k$.↵
↵
- Note that there is at most one fraction for each different denominator within this interval, because $\frac{1}{n}$ is the minimum spacing.↵
↵
- For a fraction with a denominator $q$, the numerator of the only possible fraction can only be $\left\lfloor\frac{(j + 1)q - 1}{n}\right\rfloor$.↵
↵
- Find the smallest fraction $\frac{p}{q}$ strictly greater than $\frac{j}{n}$, and use the recurrence from conclusion two until the given fraction reaches rank $k$.↵
↵
We can see that problem one can be transformed into problem two to be implemented. Therefore, we now consider how to accomplish problem two.↵
↵
More trivially, problem two can be reduced to the following problem:↵
↵
> Given a real number $x$, find the number of completely reduced proper fractions $\frac{p}{q} \le x$ with denominators $q \le n$.↵
↵
Let $A_q$ denote the number of completely reduced proper fractions with denominator $q$ that satisfy the above conditions. Obviously, the original problem can be expressed as finding such a set:↵
↵
$$S = \left\{\frac{p}{q} \bigg\vert \frac{p}{q} \le x, \, q \le n, \, \gcd(p, \, q) = 1\right\}$$↵
↵
That is, calculating $\sum\limits_{i = 1}^{n}A_i = |S|$, and we know that↵
↵
$$\left\lfloor x \cdot q \right\rfloor = \sum\limits_{d \mid q}A_d$$↵
↵
This tells us that↵
↵
$$A_q = \left\lfloor x \cdot q \right\rfloor - \sum\limits_{d \mid q, \, d < q}A_d$$↵
↵
We can simply compute this iteratively in ascending order. The complexity of each calculation is $O(n\log{n})$, making the complexity of problem one $O(n\log^2{n})$. Is there a more optimal approach?↵
↵
Considering the essential form of $A_q$, by utilizing[Möbius inversion](https://www.cnblogs.com/YipChipqwq/p/18464696 "莫比乌斯反演"), we can obtain↵
↵
$$A_q = \sum\limits_{i = 1}^{\lfloor x \cdot q \rfloor}[\gcd(i, q) = 1] = \sum\limits_{i = 1}^{\lfloor x \cdot q \rfloor}\sum\limits_{d \mid \gcd(i, q)}\mu(d)$$↵
↵
If we sum up all $A_q$, then we have↵
↵
$$\begin{aligned}↵
|S| &= \sum_{i = 1}^{n}A_i \\↵
&= \sum_{i = 1}^{n}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}[\gcd(i, \, j) = 1] \\↵
&= \sum_{i = 1}^{n}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}\sum_{d \mid \gcd(i, j)}\mu(d) \\↵
&= \sum_{d \mid i, \, d \mid j}\sum_{i = 1}^{n}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}\mu(d) \\↵
&= \sum_{d = 1}^{n}\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}\mu(d) \\↵
&= \sum_{d = 1}^{n}\mu(d)\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\lfloor x \cdot i \rfloor↵
\end{aligned}$$↵
↵
Consider setting $S(n) = \sum\limits_{i = 1}^{n}\lfloor x \cdot i \rfloor$, then what we equivalently need to calculate is↵
↵
$$\sum_{d = 1}^{n}\mu(d)S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)$$↵
↵
By calculating the prefix sums of $\lfloor x \cdot i \rfloor$ and precomputing the Möbius function, we can solve this problem in $O(n)$ time complexity.↵
↵
Thus, in solving Farey sequence problems, we have:↵
↵
- Precomputation: $O(n)$↵
↵
- Finding a fraction given its $\text{rank}$: $O(n\log{n})$↵
↵
- Finding the $\text{rank}$ given a fraction: $O(n)$↵
↵
- Finding the number of elements in $F_n$ smaller than a certain real number: $O(n)$↵
↵
- Rational approximation of positive real numbers: $O(n\log{n})$↵
↵
We provide a template code for reference.↵
↵
```cpp↵
#include<bits/stdc++.h>↵
using namespace std;↵
#define ll long long↵
#define PII pair<ll, ll>↵
const int N = 2e6 + 10;↵
ll n;↵
ll primes[N], st[N], mu[N], cnt;↵
↵
// precalculate Mu↵
void init_mu() {↵
mu[1] = 1;↵
for (int i = 2; i < N; i ++ ) {↵
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;↵
for (int j = 0; i * primes[j] < N; j ++ ) {↵
st[i * primes[j]] = 1;↵
if (i % primes[j] == 0) break;↵
mu[i * primes[j]] = -mu[i];↵
}↵
}↵
}↵
↵
// O(n + logn) calculate next one by one↵
PII nxt_fraction(PII fac) {↵
ll l = 0, r = n;↵
while (l < r) {↵
ll mid = l + r + 1 >> 1;↵
if (mid * fac.second <= n * fac.first) l = mid;↵
else r = mid - 1;↵
}↵
PII res = {r + 1, n};↵
for (int i = 1; i <= n; i ++ ) {↵
PII fac_i = {((r + 1) * i - 1) / n, i};↵
if (fac_i.first * fac.second <= fac_i.second * fac.first) continue;↵
if (fac_i.first * res.second < fac_i.second * res.first) res = fac_i;↵
}↵
ll d = __gcd(res.first, res.second);↵
return {res.first / d, res.second / d};↵
}↵
↵
// O(1) calculate next one by two nearly↵
PII nxt_fraction(PII fac1, PII fac2) {↵
auto [a, b] = fac1;↵
auto [c, d] = fac2;↵
return {(n + b) / d * c - a, (n + b) / d * d - b};↵
}↵
↵
ll fraction_to_rank(PII fac) {↵
static ll A[N];↵
for (int i = 1; i <= n; i ++ ) A[i] = fac.first * i / fac.second + A[i - 1];↵
ll res = 0;↵
for (int i = 1; i <= n; i ++ ) res += mu[i] * A[n / i];↵
return res;↵
}↵
↵
PII rank_to_fraction(ll k) {↵
ll l = 0, r = n;↵
while (l < r) {↵
ll mid = l + r + 1 >> 1;↵
if (fraction_to_rank({mid, n}) <= k) l = mid;↵
else r = mid - 1;↵
}↵
k -= fraction_to_rank({r, n});↵
PII a = {r / __gcd(r, n), n / __gcd(r, n)}, b = nxt_fraction(a);↵
if (!k) return a;↵
while (k -- ) a = nxt_fraction(a, b), swap(a, b);↵
return a;↵
}↵
↵
void solve() {↵
cin >> n;↵
PII now = {1, n};↵
while (now.first != now.second) {↵
ll rk = fraction_to_rank(now);↵
PII fac = rank_to_fraction(rk);↵
cout << rk << " : " << fac.first << " " << fac.second << " | " << now.first << " " << now.second << "\n";↵
now = nxt_fraction(now);↵
}↵
}↵
↵
int main() {↵
ios::sync_with_stdio(false);↵
cin.tie(0), cout.tie(0);↵
init_mu();↵
int T = 1;↵
// cin >> T;↵
while (T -- ) solve();↵
return 0;↵
}↵
```↵
↵
## Expansion↵
↵
Knowing one term, the time complexity of solving for the next term is $O(n + \log{n})$. To find the closest fraction $\frac{c}{d}$ that satisfies $bc - ad = 1$, we can use the Extended Euclidean Algorithm. Since the difference between two adjacent terms is $\frac{1}{bd}$, we only need to maximize $d$.↵
↵
On the other hand, it is easy to see that we rarely encounter situations where $x$ is an irrational number in practical problems. Therefore, we can naturally assume $x = \frac{p}{q}$, giving us↵
↵
$$\sum\limits_{d = 1}^{n}\mu(d)\sum\limits_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\lfloor x \cdot i \rfloor = \sum\limits_{d = 1}^{n}\mu(d)\sum\limits_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor \frac{pi}{q} \right\rfloor$$↵
↵
It is easy to see that the latter is the standard form for our Euclidean-like algorithm. By using number theory block division (square root decomposition), we can achieve $O(\sqrt{n}\log^2{n})$ per query.↵
↵
Furthermore, you can use the Du Sieve to sum the preceding Möbius function, reducing the precomputation time complexity to $O(n^{\frac{2}{3}})$.↵
↵
If you still feel this method is not optimal, you can utilize Dirichlet prefix sums with block processing, achieving an astonishing overall time complexity of $O(n^{\frac{2}{3}} + \sqrt{n}\log^{\frac{3}{2}}{n})$.↵
↵
```cpp↵
#include<bits/stdc++.h>↵
using namespace std;↵
#define ll long long↵
#define PII pair<ll, ll>↵
const int N = 1e5;↵
ll n, k;↵
ll primes[N], st[N], mu[N], cnt;↵
unordered_map<ll, ll> summu;↵
↵
// precalculate Mu↵
void init_mu() {↵
mu[1] = 1;↵
for (int i = 2; i < N; i ++ ) {↵
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;↵
for (int j = 0; i * primes[j] < N; j ++ ) {↵
st[i * primes[j]] = 1;↵
if (i % primes[j] == 0) break;↵
mu[i * primes[j]] = -mu[i];↵
}↵
}↵
for (int i = 2; i < N; i ++ ) mu[i] += mu[i - 1];↵
}↵
↵
ll get_mu(ll n) {↵
if (n < N) return mu[n];↵
if (summu.count(n)) return summu[n];↵
ll res = 1;↵
for (ll l = 2, r; l <= n; l = r + 1) {↵
r = n / (n / l);↵
res -= (r - l + 1) * get_mu(n / l);↵
}↵
return summu[n] = res;↵
}↵
↵
ll f(ll a, ll b, ll c, ll n) {↵
if (!a) return b / c * (n + 1);↵
if (a < c && b < c) {↵
ll m = (a * n + b) / c;↵
if (!m) return 0;↵
return n * m - f(c, c - b - 1, a, m - 1);↵
}↵
return f(a % c, b % c, c, n) + (n + 1) * n / 2 * (a / c) + (n + 1) * (b / c);↵
}↵
↵
// O(logn) calculate next one by one↵
PII nxt_fraction(ll n, PII fac) {↵
ll a = fac.first, b = fac.second, c, d;↵
exgcd(b, a, c, d);↵
d = -d;↵
if (d <= 0) {↵
ll k = (-d + b - 1) / b;↵
c += k * a;↵
d += k * b;↵
}↵
ll k = (n - d) / b;↵
return {c + k * a, d + k * b};↵
}↵
↵
// O(1) calculate next one by two nearly↵
PII nxt_fraction(ll n, PII fac1, PII fac2) {↵
auto [a, b] = fac1;↵
auto [c, d] = fac2;↵
return {(n + b) / d * c - a, (n + b) / d * d - b};↵
}↵
↵
// O(sqrt(n)logn) calculate fraction's rank↵
ll fraction_to_rank(ll n, PII fac) {↵
ll res = 0;↵
for (ll l = 1, r; l <= n; l = r + 1) {↵
r = n / (n / l);↵
res += (get_mu(r) - get_mu(l - 1)) * f(fac.first, 0, fac.second, n / l);↵
}↵
return res;↵
}↵
↵
// O(sqrt(n)log(n)^2) calculate the fraction of k-th rank↵
PII rank_to_fraction(ll n, ll k) {↵
ll l = 0, r = n;↵
while (l < r) {↵
ll mid = l + r + 1 >> 1;↵
if (fraction_to_rank(n, {mid, n}) <= k) l = mid;↵
else r = mid - 1;↵
}↵
k -= fraction_to_rank(n, {r, n});↵
PII a = {r / __gcd(r, n), n / __gcd(r, n)}, b = nxt_fraction(n, a);↵
if (!k) return a;↵
while (k -- ) a = nxt_fraction(n, a, b), swap(a, b);↵
return a;↵
}↵
↵
void solve() {↵
cin >> n;↵
PII now = {1, n};↵
while (now.first != now.second) {↵
ll rk = fraction_to_rank(n, now);↵
PII fac = rank_to_fraction(n, rk);↵
cout << rk << " : " << fac.first << " " << fac.second << " | " << now.first << " " << now.second << "\n";↵
now = nxt_fraction(n, now);↵
}↵
}↵
↵
int main() {↵
ios::sync_with_stdio(false);↵
cin.tie(0), cout.tie(0);↵
init_mu();↵
int T = 1;↵
// cin >> T;↵
while (T -- ) solve();↵
return 0;↵
}↵
```
↵
$$F_n = \left\{\frac{p}{q} \bigg\vert 0 < p < q \le n, \, \gcd(p, q) = 1\right\}$$↵
↵
The asymptotic size of this sequence is $\sum\limits_{i = 1}^{n} \varphi(n) \sim \frac{3}{\pi^2}n^2 + o(n \log n)$.↵
↵
Regarding the Farey sequence, there are two classic conclusions:↵
↵
- If $\frac{a}{b}$ and $\frac{c}{d}$ are two consecutive elements in $F_n$, then $\frac{a + c}{b + d}$ is also a valid Farey fraction.↵
↵
- If $\frac{a}{b}$ and $\frac{c}{d}$ are two consecutive elements in $F_n$, then their next element $\frac{p}{q}$ satisfies↵
$$\begin{cases}p = \left\lfloor\frac{n + b}{d}\right\rfloor c - a \\ q = \left\lfloor\frac{n + b}{d}\right\rfloor d - b\end{cases}$$↵
↵
We will only prove the second theorem above. For the Farey sequence of order $n$, two adjacent terms $\frac{a}{b} < \frac{c}{d}$ satisfy $bc - ad = 1$. Since the adjacent terms $\frac{a}{b}, \, \frac{c}{d}$ are both on the Stern-Brocot tree, the next term $\frac{p}{q}$ satisfies $\frac{c}{d} = \frac{a + p}{b + q}$, which means $(p + a)d = (q + b)c$. Therefore, there exists $k$ such that↵
↵
$$\begin{cases} kc = a + p \\ kd = b + q \end{cases}$$↵
↵
To make $p, \, q$ highly precise (i.e., closely packed), the value of $k$ should be as large as possible. Therefore, concerning the upper bound $n$, $p, \, q$ are defined by the following equations:↵
↵
$$\begin{cases}p = kc - a \le n \\ q = kd - b \le n\end{cases}$$↵
↵
Obviously, from the definition of a proper fraction, we know $q > p$, so the maximum $k = \left\lfloor\frac{n + b}{d}\right\rfloor$, thus↵
↵
$$\begin{cases}p = \left\lfloor\frac{n + b}{d}\right\rfloor c - a \\ q = \left\lfloor\frac{n + b}{d}\right\rfloor d - b\end{cases}$$↵
↵
There are two computational methods to generate this sequence. The first one is based on the Stern-Brocot tree from the first theorem, which is very easy to implement.↵
↵
```cpp↵
int build(int a, int b, int c, int d, int n) {↵
if (b + d > n) return 0;↵
return 1 + build(a, b, a + c, b + d, n) + build(a + c, b + d, c, d, n);↵
}↵
```↵
↵
The other one is based on the recurrence from the second theorem, which is also quite easy to write.↵
↵
```cpp↵
PII nxt_fraction(PII fac1, PII fac2) {↵
auto [a, b] = fac1;↵
auto [c, d] = fac2;↵
return {(n + b) / d * c - a, (n + b) / d * d - b};↵
}↵
```↵
↵
Obviously, they both use an $O(n^2)$ time complexity to compute specific elements of the Farey sequence.↵
↵
To facilitate solving the parent problem, we present the following two problems:↵
↵
> Given $n, \, k$, find the $k$-th term of $F_n$.↵
↵
> Given $n$ and a completely reduced proper fraction $\frac{p}{q} (p \le n)$, determine its rank in $F_n$.↵
↵
Let's consider the first problem first. We find that the first problem can be processed through the following steps:↵
↵
- Binary search for $j$ in $\frac{j}{n}$, and calculate $\text{rank}(\frac{j}{n})$↵
↵
- Let $r = \text{rank}(\frac{j}{n})$. If $r < k$, search upwards; otherwise, search downwards.↵
↵
- Find an interval $\left[\frac{j}{n}, \, \frac{j + 1}{n}\right)$ such that the target fraction falls within this range.↵
↵
- Count the fractions within this interval and find a valid fraction such that its rank equals the given $k$.↵
↵
- Note that there is at most one fraction for each different denominator within this interval, because $\frac{1}{n}$ is the minimum spacing.↵
↵
- For a fraction with a denominator $q$, the numerator of the only possible fraction can only be $\left\lfloor\frac{(j + 1)q - 1}{n}\right\rfloor$.↵
↵
- Find the smallest fraction $\frac{p}{q}$ strictly greater than $\frac{j}{n}$, and use the recurrence from conclusion two until the given fraction reaches rank $k$.↵
↵
We can see that problem one can be transformed into problem two to be implemented. Therefore, we now consider how to accomplish problem two.↵
↵
More trivially, problem two can be reduced to the following problem:↵
↵
> Given a real number $x$, find the number of completely reduced proper fractions $\frac{p}{q} \le x$ with denominators $q \le n$.↵
↵
Let $A_q$ denote the number of completely reduced proper fractions with denominator $q$ that satisfy the above conditions. Obviously, the original problem can be expressed as finding such a set:↵
↵
$$S = \left\{\frac{p}{q} \bigg\vert \frac{p}{q} \le x, \, q \le n, \, \gcd(p, \, q) = 1\right\}$$↵
↵
That is, calculating $\sum\limits_{i = 1}^{n}A_i = |S|$, and we know that↵
↵
$$\left\lfloor x \cdot q \right\rfloor = \sum\limits_{d \mid q}A_d$$↵
↵
This tells us that↵
↵
$$A_q = \left\lfloor x \cdot q \right\rfloor - \sum\limits_{d \mid q, \, d < q}A_d$$↵
↵
We can simply compute this iteratively in ascending order. The complexity of each calculation is $O(n\log{n})$, making the complexity of problem one $O(n\log^2{n})$. Is there a more optimal approach?↵
↵
Considering the essential form of $A_q$, by utilizing
↵
$$A_q = \sum\limits_{i = 1}^{\lfloor x \cdot q \rfloor}[\gcd(i, q) = 1] = \sum\limits_{i = 1}^{\lfloor x \cdot q \rfloor}\sum\limits_{d \mid \gcd(i, q)}\mu(d)$$↵
↵
If we sum up all $A_q$, then we have↵
↵
$$\begin{aligned}↵
|S| &= \sum_{i = 1}^{n}A_i \\↵
&= \sum_{i = 1}^{n}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}[\gcd(i, \, j) = 1] \\↵
&= \sum_{i = 1}^{n}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}\sum_{d \mid \gcd(i, j)}\mu(d) \\↵
&= \sum_{d \mid i, \, d \mid j}\sum_{i = 1}^{n}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}\mu(d) \\↵
&= \sum_{d = 1}^{n}\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j = 1}^{\lfloor x \cdot i \rfloor}\mu(d) \\↵
&= \sum_{d = 1}^{n}\mu(d)\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\lfloor x \cdot i \rfloor↵
\end{aligned}$$↵
↵
Consider setting $S(n) = \sum\limits_{i = 1}^{n}\lfloor x \cdot i \rfloor$, then what we equivalently need to calculate is↵
↵
$$\sum_{d = 1}^{n}\mu(d)S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)$$↵
↵
By calculating the prefix sums of $\lfloor x \cdot i \rfloor$ and precomputing the Möbius function, we can solve this problem in $O(n)$ time complexity.↵
↵
Thus, in solving Farey sequence problems, we have:↵
↵
- Precomputation: $O(n)$↵
↵
- Finding a fraction given its $\text{rank}$: $O(n\log{n})$↵
↵
- Finding the $\text{rank}$ given a fraction: $O(n)$↵
↵
- Finding the number of elements in $F_n$ smaller than a certain real number: $O(n)$↵
↵
- Rational approximation of positive real numbers: $O(n\log{n})$↵
↵
We provide a template code for reference.↵
↵
```cpp↵
#include<bits/stdc++.h>↵
using namespace std;↵
#define ll long long↵
#define PII pair<ll, ll>↵
const int N = 2e6 + 10;↵
ll n;↵
ll primes[N], st[N], mu[N], cnt;↵
↵
// precalculate Mu↵
void init_mu() {↵
mu[1] = 1;↵
for (int i = 2; i < N; i ++ ) {↵
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;↵
for (int j = 0; i * primes[j] < N; j ++ ) {↵
st[i * primes[j]] = 1;↵
if (i % primes[j] == 0) break;↵
mu[i * primes[j]] = -mu[i];↵
}↵
}↵
}↵
↵
// O(n + logn) calculate next one by one↵
PII nxt_fraction(PII fac) {↵
ll l = 0, r = n;↵
while (l < r) {↵
ll mid = l + r + 1 >> 1;↵
if (mid * fac.second <= n * fac.first) l = mid;↵
else r = mid - 1;↵
}↵
PII res = {r + 1, n};↵
for (int i = 1; i <= n; i ++ ) {↵
PII fac_i = {((r + 1) * i - 1) / n, i};↵
if (fac_i.first * fac.second <= fac_i.second * fac.first) continue;↵
if (fac_i.first * res.second < fac_i.second * res.first) res = fac_i;↵
}↵
ll d = __gcd(res.first, res.second);↵
return {res.first / d, res.second / d};↵
}↵
↵
// O(1) calculate next one by two nearly↵
PII nxt_fraction(PII fac1, PII fac2) {↵
auto [a, b] = fac1;↵
auto [c, d] = fac2;↵
return {(n + b) / d * c - a, (n + b) / d * d - b};↵
}↵
↵
ll fraction_to_rank(PII fac) {↵
static ll A[N];↵
for (int i = 1; i <= n; i ++ ) A[i] = fac.first * i / fac.second + A[i - 1];↵
ll res = 0;↵
for (int i = 1; i <= n; i ++ ) res += mu[i] * A[n / i];↵
return res;↵
}↵
↵
PII rank_to_fraction(ll k) {↵
ll l = 0, r = n;↵
while (l < r) {↵
ll mid = l + r + 1 >> 1;↵
if (fraction_to_rank({mid, n}) <= k) l = mid;↵
else r = mid - 1;↵
}↵
k -= fraction_to_rank({r, n});↵
PII a = {r / __gcd(r, n), n / __gcd(r, n)}, b = nxt_fraction(a);↵
if (!k) return a;↵
while (k -- ) a = nxt_fraction(a, b), swap(a, b);↵
return a;↵
}↵
↵
void solve() {↵
cin >> n;↵
PII now = {1, n};↵
while (now.first != now.second) {↵
ll rk = fraction_to_rank(now);↵
PII fac = rank_to_fraction(rk);↵
cout << rk << " : " << fac.first << " " << fac.second << " | " << now.first << " " << now.second << "\n";↵
now = nxt_fraction(now);↵
}↵
}↵
↵
int main() {↵
ios::sync_with_stdio(false);↵
cin.tie(0), cout.tie(0);↵
init_mu();↵
int T = 1;↵
// cin >> T;↵
while (T -- ) solve();↵
return 0;↵
}↵
```↵
↵
## Expansion↵
↵
Knowing one term, the time complexity of solving for the next term is $O(n + \log{n})$. To find the closest fraction $\frac{c}{d}$ that satisfies $bc - ad = 1$, we can use the Extended Euclidean Algorithm. Since the difference between two adjacent terms is $\frac{1}{bd}$, we only need to maximize $d$.↵
↵
On the other hand, it is easy to see that we rarely encounter situations where $x$ is an irrational number in practical problems. Therefore, we can naturally assume $x = \frac{p}{q}$, giving us↵
↵
$$\sum\limits_{d = 1}^{n}\mu(d)\sum\limits_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\lfloor x \cdot i \rfloor = \sum\limits_{d = 1}^{n}\mu(d)\sum\limits_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor \frac{pi}{q} \right\rfloor$$↵
↵
It is easy to see that the latter is the standard form for our Euclidean-like algorithm. By using number theory block division (square root decomposition), we can achieve $O(\sqrt{n}\log^2{n})$ per query.↵
↵
Furthermore, you can use the Du Sieve to sum the preceding Möbius function, reducing the precomputation time complexity to $O(n^{\frac{2}{3}})$.↵
↵
If you still feel this method is not optimal, you can utilize Dirichlet prefix sums with block processing, achieving an astonishing overall time complexity of $O(n^{\frac{2}{3}} + \sqrt{n}\log^{\frac{3}{2}}{n})$.↵
↵
```cpp↵
#include<bits/stdc++.h>↵
using namespace std;↵
#define ll long long↵
#define PII pair<ll, ll>↵
const int N = 1e5;↵
ll n, k;↵
ll primes[N], st[N], mu[N], cnt;↵
unordered_map<ll, ll> summu;↵
↵
// precalculate Mu↵
void init_mu() {↵
mu[1] = 1;↵
for (int i = 2; i < N; i ++ ) {↵
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;↵
for (int j = 0; i * primes[j] < N; j ++ ) {↵
st[i * primes[j]] = 1;↵
if (i % primes[j] == 0) break;↵
mu[i * primes[j]] = -mu[i];↵
}↵
}↵
for (int i = 2; i < N; i ++ ) mu[i] += mu[i - 1];↵
}↵
↵
ll get_mu(ll n) {↵
if (n < N) return mu[n];↵
if (summu.count(n)) return summu[n];↵
ll res = 1;↵
for (ll l = 2, r; l <= n; l = r + 1) {↵
r = n / (n / l);↵
res -= (r - l + 1) * get_mu(n / l);↵
}↵
return summu[n] = res;↵
}↵
↵
ll f(ll a, ll b, ll c, ll n) {↵
if (!a) return b / c * (n + 1);↵
if (a < c && b < c) {↵
ll m = (a * n + b) / c;↵
if (!m) return 0;↵
return n * m - f(c, c - b - 1, a, m - 1);↵
}↵
return f(a % c, b % c, c, n) + (n + 1) * n / 2 * (a / c) + (n + 1) * (b / c);↵
}↵
↵
// O(logn) calculate next one by one↵
PII nxt_fraction(ll n, PII fac) {↵
ll a = fac.first, b = fac.second, c, d;↵
exgcd(b, a, c, d);↵
d = -d;↵
if (d <= 0) {↵
ll k = (-d + b - 1) / b;↵
c += k * a;↵
d += k * b;↵
}↵
ll k = (n - d) / b;↵
return {c + k * a, d + k * b};↵
}↵
↵
// O(1) calculate next one by two nearly↵
PII nxt_fraction(ll n, PII fac1, PII fac2) {↵
auto [a, b] = fac1;↵
auto [c, d] = fac2;↵
return {(n + b) / d * c - a, (n + b) / d * d - b};↵
}↵
↵
// O(sqrt(n)logn) calculate fraction's rank↵
ll fraction_to_rank(ll n, PII fac) {↵
ll res = 0;↵
for (ll l = 1, r; l <= n; l = r + 1) {↵
r = n / (n / l);↵
res += (get_mu(r) - get_mu(l - 1)) * f(fac.first, 0, fac.second, n / l);↵
}↵
return res;↵
}↵
↵
// O(sqrt(n)log(n)^2) calculate the fraction of k-th rank↵
PII rank_to_fraction(ll n, ll k) {↵
ll l = 0, r = n;↵
while (l < r) {↵
ll mid = l + r + 1 >> 1;↵
if (fraction_to_rank(n, {mid, n}) <= k) l = mid;↵
else r = mid - 1;↵
}↵
k -= fraction_to_rank(n, {r, n});↵
PII a = {r / __gcd(r, n), n / __gcd(r, n)}, b = nxt_fraction(n, a);↵
if (!k) return a;↵
while (k -- ) a = nxt_fraction(n, a, b), swap(a, b);↵
return a;↵
}↵
↵
void solve() {↵
cin >> n;↵
PII now = {1, n};↵
while (now.first != now.second) {↵
ll rk = fraction_to_rank(n, now);↵
PII fac = rank_to_fraction(n, rk);↵
cout << rk << " : " << fac.first << " " << fac.second << " | " << now.first << " " << now.second << "\n";↵
now = nxt_fraction(n, now);↵
}↵
}↵
↵
int main() {↵
ios::sync_with_stdio(false);↵
cin.tie(0), cout.tie(0);↵
init_mu();↵
int T = 1;↵
// cin >> T;↵
while (T -- ) solve();↵
return 0;↵
}↵
```
