浅谈 Farey 数列

YipChip · 2025-05-29 19:16:03 · 题解

首先，Farey 数列 F_n 表示分母不超过 n 的所有既约真分数按大小顺序排列的集合，形式化来说：

F_n = \left\{\frac{p}{q} \bigg\vert 0 < p < q \le n, \, \gcd(p, q) = 1\right\}

这个数列的渐进大小为 \sum\limits_{i = 1}^{n} \varphi(n) \sim \frac{3}{\pi^2}n^2 + o(n \log n)。

关于 Farey 数列，有两个经典结论：

若 \frac{a}{b} 和 \frac{c}{d} 是 F_n 中的两个连续元素，那么 \frac{a + c}{b + d} 也是一个合法的 Faray 数。
若 \frac{a}{b} 和 \frac{c}{d} 是 F_n 中的两个连续元素，那么它们的下一个元素 \frac{p}{q} 满足：
\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}

我们只对上述定理二进行证明，对 n 阶 Farey 数列，相邻两项 \frac{a}{b} < \frac{c}{d} 满足 bc - ad = 1，根据相邻两项 \frac{a}{b}, \, \frac{c}{d} 均在 Stern-Brocot 树上，我们有下一项 \frac{p}{q} 满足 \frac{c}{d} = \frac{a + p}{b + q}，也就是 (p + a)d = (q + b)c，因此存在 k 使得：

\begin{cases} kc = a + p \\ kd = b + q \end{cases}

为了使 p, \, q 精度高，那么 k 的值应该尽可能大，所以对上界 n 来说，p, \, q 由下式定义：

\begin{cases}p = kc - a \le n \\ q = kd - b \le n\end{cases}

显然由真分数定义知 q > p，所以最大 k = \left\lfloor\frac{n + b}{d}\right\rfloor，因此

\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}

生成这个数列有两种计算方式，一是种基于第一个定理的 Stern-Brocot 树，这个很好实现。

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);
}

另一种是基于第二个定理的递推，也比较好写。

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};
}

很显然，他们都是 O(n^2) 的时间复杂度进行 Faray 数列特定元素的计算。

为了方便我们解决母问题，我们给出这样两个问题：

给定 n, \, k，求 F_n 的第 k 项。

给定 n 和既约真分数 \frac{p}{q}(p \le n)，判断其在 F_n 中的排名。

我们先考虑第一个问题，我们发现第一个问题可以通过如下步骤进行：

二分 \frac{j}{n} 中的 j，计算 \text{rank}(\frac{j}{n})。
- 记 r = \text{rank}(\frac{j}{n})，如果 r < k，向上搜索，否则向下搜索。
找到一个区间 \left[\frac{j}{n}, \, \frac{j + 1}{n}\right) 使得目标分数在这个范围以内。
统计这个区间内的分数，找到合法的分数使得其排名为所给 k。
- 注意到这个区间内不同分母的分数至多一个，因为 \frac{1}{n} 是最小间距了。
- 对于一个分母为 q 的分数，唯一可能的分数的分母只能是 \left\lfloor\frac{(j + 1)q - 1}{n}\right\rfloor。
找到严格大于 \frac{j}{n} 的最小分数 \frac{p}{q}，利用结论二进行递推，直到给定分数为 k 排名。

我们可以看出，问题一可以转化为问题二实现，因此我们现在来考虑问题二如何完成。

更平凡的，问题二可以规约为如下问题：

给定实数 x，求分母 q \le n 的既约真分数 \frac{p}{q} \le x 的个数。

记 A_q 表示，以 q 为分母满足上述条件的既约真分数的个数，显然的，可以表示原问题为求这样一个集合：

S = \left\{\frac{p}{q} \bigg\vert \frac{p}{q} \le x, \, q \le n, \, \gcd(p, \, q) = 1\right\}

也就是求 \sum\limits_{i = 1}^{n}A_i = |S|，而我们知道：

\left\lfloor x \cdot q \right\rfloor = \sum\limits_{d \mid q}A_d

这告诉我们：

A_q = \left\lfloor x \cdot q \right\rfloor - \sum\limits_{d \mid q, \, d < q}A_d

按顺序从小到大递推即可，每次计算复杂度 O(n\log{n})，那么问题一的复杂度就是 O(n\log^2{n})，有没有更优的做法呢？

考虑 A_q 的本质形态，利用莫比乌斯反演，我们可以得到：

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)

如果我们对所有 A_q 求和，那么就有：

|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}

考虑设 S(n) = \sum\limits_{i = 1}^{n}\lfloor x \cdot i \rfloor，那么我们相当于要求的东西也就是：

\sum_{d = 1}^{n}\mu(d)S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)

对 \lfloor x \cdot i \rfloor 求前缀和，预处理莫比乌斯函数，我们可以在 O(n) 时间复杂度解决这个问题。

如此，我们在解决 Farey 数列问题上我们有：

预处理：O(n)
给 \text{rank} 求分数：O(n\log{n})
给分数求 \text{rank}：O(n)
求 F_n 中比某实数小的数的个数：O(n)
正实数的有理逼近：O(n\log{n})

我们给出一份模版代码供参考。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define PII pair<ll, ll>
const int N = 4e4 + 10;
ll n, k;
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(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(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, k;
    PII fac = rank_to_fraction(rk);
    cout << fac.first << " " << fac.second;
}

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

已知一项，求解下一项我们的时间复杂度是 O(n + \log{n}) 的，我们为了求解满足 bc - ad = 1 的最接近分数 \frac{c}{d} 可以利用扩展欧几里得算法，因为相邻两项的差为 \frac{1}{bd}，我们只需要 d 取到最大值即可。

另一方面，容易发现，我们实际问题中很难遇到 x 是无理数的情况，那么我们很容易想到，不妨设 x = \frac{p}{q}，有：

\sum_{d = 1}^{n}\mu(d)\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\lfloor x \cdot i \rfloor = \sum_{d = 1}^{n}\mu(d)\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor \frac{pi}{q} \right\rfloor

容易看出，后者是我们类欧几里得算法的标准形式，利用数论分块，可以做到 O(\sqrt{n}\log^2{n}) 每个询问。

更进一步，你可以使用杜教筛对前面的莫比乌斯函数求和，预处理时间复杂度降为 O(n^{\frac{2}{3}})。

如果你仍然对此法感觉不优，你可以利用狄利克雷前缀和分块处理，总体时间复杂度可以达到惊人的 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 = 1200;
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(n + logn) calculate next one by one
PII nxt_fraction(ll n, 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(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 >> k;
    PII fac = rank_to_fraction(n, k);
    cout << fac.first << " " << fac.second;
}

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;
}