algorithm - 在 Farey 序列的第 n 阶中查找第 k 个元素

Question

n 阶 Farey 序列是完全减少的分数序列，介于 0 和 1 之间，当最低项时分母小于或等于 n，按大小递增的顺序排列。详解here .

问题

问题是，给定 n 和 k，其中 n = seq 的顺序，k = 元素索引，我们能否从序列中找到特定元素。例如 (n=5, k =6) 的答案是 1/2。

铅

有许多不是最优的可用解决方案，但我正在寻找一个接近最优的解决方案。讨论了一种这样的算法 here ，我无法理解其中的逻辑，因此无法应用示例。

问题

能否请一些人更详细地解释解决方案，最好是举例说明。

谢谢。

Answer 1

我已经阅读了您链接中提供的方法，以及公认的 C++ 解决方案。让我发布它们以供引用:

编辑说明

Several less-than-optimal solutions exist. Using a priority queue, one can iterate through the fractions (generating them one by one) in O(K log N) time. Using a fancier math relation, this can be reduced to O(K). However, neither of these solution obtains many points, because the number of fractions (and thus K) is quadratic in N.

The “good” solution is based on meta-binary search. To construct this solution, we need the following subroutine: given a fraction A/B (which is not necessarily irreducible), find how many fractions from the Farey sequence are less than this fraction. Suppose we had this subroutine; then the algorithm works as follows:

• Determine a number X such that the answer is between X/N and (X+1)/N; such a number can be determined by binary searching the range 1...N, thus calling the subroutine O(log N) times.
  • Make a list of all fractions A/B in the range X/N...(X+1)/N. For any given B, there is at most one A in this range, and it can be determined trivially in O(1).
  • Determine the appropriate order statistic in this list (doing this in O(N log N) by sorting is good enough).

It remains to show how we can construct the desired subroutine. We will show how it can be implemented in O(N log N), thus giving a O(N log^2 N) algorithm overall. Let us denote by C[j] the number of irreducible fractions i/j which are less than X/N. The algorithm is based on the following observation: C[j] = floor(X*B/N) – Sum(C[D], where D divides j). A direct implementation, which tests whether any D is a divisor, yields a quadratic algorithm. A better approach, inspired by Eratosthene’s sieve, is the following: at step j, we know C[j], and we subtract it from all multiples of j. The running time of the subroutine becomes O(N log N).

相关代码

#include <cassert>

#include <algorithm>
#include <fstream>
#include <iostream>
#include <vector>
using namespace std;

const int kMaxN = 2e5;
typedef int int32;
typedef long long int64_x;
// #define int __int128_t
// #define int64 __int128_t

typedef long long int64;


int64 count_less(int a, int n) {
    vector<int> counter(n + 1, 0);
    for (int i = 2; i <= n; i += 1) {
        counter[i] = min(1LL * (i - 1), 1LL * i * a / n);
    }

    int64 result = 0;
    for (int i = 2; i <= n; i += 1) {
        for (int j = 2 * i; j <= n; j += i) {
            counter[j] -= counter[i];
        }
        result += counter[i];
    }

    return result;
}

int32 main() {
//    ifstream cin("farey.in");
//    ofstream cout("farey.out");
    int64_x n, k; cin >> n >> k;
    assert(1 <= n);
    assert(n <= kMaxN);
    assert(1 <= k);
    assert(k <= count_less(n, n));

    int up = 0;
    for (int p = 29; p >= 0; p -= 1) {
        if ((1 << p) + up > n) 
            continue;

        if (count_less((1 << p) + up, n) < k) {
            up += (1 << p);
        }
    }

    k -= count_less(up, n);

    vector<pair<int, int>> elements;
    for (int i = 1; i <= n; i += 1) {
        int b = i;
        // find a such that up/n < a / b and a / b <= (up+1) / n
        int a = 1LL * (up + 1) * b / n;
        if (1LL * up * b < 1LL * a * n) {
        } else {
            continue;
        }

        if (1LL * a * n <= 1LL * (up + 1) * b) {
        } else {
            continue;
        }

        if (__gcd(a, b) != 1) {
            continue;
        }

        elements.push_back({a, b});
    }

    sort(elements.begin(), elements.end(), 
            [](const pair<int, int>& lhs, const pair<int, int>& rhs) -> bool {
                return 1LL * lhs.first * rhs.second < 1LL * rhs.first * lhs.second; 
            });

    cout << (int64_x)elements[k - 1].first << ' ' << (int64_x)elements[k - 1].second << '\n';
    return 0;
}

基本方法论

以上编辑解释导致以下简化版本。让我从一个例子开始。

比方说，我们想找到 N = 5 的 Farey 序列的第 7 个元素。

1. 我们首先编写一个子例程，如解释中所述，它为我们提供“k”值(在给定分数之前存在多少 Farey 序列约简分数 - 给定数字可能会或可能不会减少)<

所以，按照你的 F5 序列:

k = 0,  0/1
k = 1,  1/5
k = 2,  1/4
k = 3,  1/3
k = 4,  2/5
k = 5,  1/2
k = 6,  3/5
k = 7,  2/3
k = 8,  3/4
k = 9,  4/5
k = 10, 1/1

如果我们可以找到一个函数来计算 Farey 序列中先前减少分数的计数，我们可以执行以下操作:

int64 k_count_2 = count_less(2, 5); // result = 4
int64 k_count_3 = count_less(3, 5); // result = 6
int64 k_count_4 = count_less(4, 5); // result = 9

这个函数写在公认的解决方案中。它使用社论最后一段中解释的确切方法。

如您所见，count_less()函数生成相同的 k我们手写列表中的值。

我们使用该函数知道 k = 4、6、9 的约化分数的值。 k = 7 呢？如社论中所述，我们将列出 X/N 和 (X+1)/N 范围内的所有约化分数，此处 X = 3 和 N = 5。

使用已接受解决方案(靠近底部)中的函数，我们列出并排序减少的分数。

之后我们将重新排列我们的 k 值，以适应我们的新数组:

k = -,  0/1
k = -,  1/5
k = -,  1/4
k = -,  1/3
k = -,  2/5
k = -,  1/2
k = -,  3/5 <-|
k = 0,  2/3   | We list and sort the possible reduced fractions 
k = 1,  3/4   | in between these numbers
k = -,  4/5 <-|
k = -, 1/1

(这就是为什么有这段代码:k -= count_less(up, n);，它基本上重新映射了 k 个值)

(而且我们在索引时还减去一个，即:cout << (int64_x)elements[k - 1].first << ' ' << (int64_x)elements[k - 1].second << '\n';。这只是为了基本上调用生成数组中的正确位置。)

因此，对于我们重新映射的新 k 值，对于 N = 5 和 k = 7(原始 k)，我们的结果是 2/3。

(我们在新 map 中选择值 k = 0)

如果你编译并运行接受的解决方案，它会给你这个:

Input:  5 7 (Enter)
Output: 2 3

我相信这是编辑和接受的解决方案的基本要点。