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,
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
We will only prove the second theorem above. For the Farey sequence of order $$$n$$$, two adjacent terms $$$\frac{a}{b} \lt \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
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:
Obviously, from the definition of a proper fraction, we know $$$q \gt p$$$, so the maximum $$$k = \left\lfloor\frac{n + b}{d}\right\rfloor$$$, thus
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.
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.
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 \lt 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:
That is, calculating $$$\sum\limits_{i = 1}^{n}A_i = |S|$$$, and we know that
This tells us that
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, we can obtain
If we sum up all $$$A_q$$$, then we have
Consider setting $$$S(n) = \sum\limits_{i = 1}^{n}\lfloor x \cdot i \rfloor$$$, then what we equivalently need to calculate is
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.
#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
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})$$$.
#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;
}
