按分母和分子的升序查找两个给定分数之间的不可约分数答案

【问题标题】：Finding irreducible fractions between two given fractions in ascending order of denominators and then numerators按分母和分子的升序查找两个给定分数之间的不可约分数
【发布时间】：2020-04-29 00:41:28
【问题描述】：

每次测试的时间限制：2 秒
每次测试的内存限制：512 MB

给你两个分数a/b c/d 和一个正数N。考虑所有不可约分数e/f 使得0 < e, f ≤ N 和 a/b e/f c/d。让s 是这些分数的序列按分母升序排序，然后是分子（分数如果f1 < f2 或f1 = f2 and e1 < e2，e1/f1 在e2/f2 之前）。您应该首先打印序列s 的n 项或整个序列s，如果它包含少于n 项。

输入
每个测试的第一行包含6个整数a、b、c、d、N、n（0 ≤ a ≤ 10^18、1 ≤ b, c, d, N ≤ 10^18、1 ≤ n ≤ 200 000、@987653）。

输出
首先，打印您将输出多少个序列s。然后以正确的顺序输出这些术语。

示例
输入：
0 1 1 1 5 10
输出：
9
1 2
1 3
2 3
1 4
3 4
1 5
2 5
3 5
4 5
输入：
55 34 68 42 90 1
输出：
1
89 55
输入：
49 33 45 30 50 239
输出：
0

到目前为止，我只设法编写了一个解决方案，它遍历从 1 到 N 的所有分母，并且对于每个分母，遍历从 a*f/b 到 c*f/d 的所有分子，添加所有找到答案的不可约分数。

这是我的代码：

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

long long a, b, c, d, N, n;
vector<pair<long long, long long>> result;

long long gcd(long long a, long long b) {
    while (b) {
        a %= b;
        swap(a, b);
    }
    return a;
}

void computeResult() {
    for (long long f = 1; f <= N; f++) {
        long long eMax = c*f / d;
        if (c*f % d != 0) eMax++;
        eMax = min(eMax, N);
        for (long long e = a*f / b + 1; e < eMax; e++) {
            if (gcd(e, f) == 1) {
                result.push_back(make_pair(e, f));
                if (result.size() == n)
                    return;
            }
        }
    }
}

int main() {
    cin >> a >> b >> c >> d >> N >> n;  
    computeResult();
    cout << result.size() << endl;
    for (pair<long long, long long> fraction : result)
        cout << fraction.first << " " << fraction.second << endl;
}

很遗憾，这个解决方案太慢了。我想知道如何更有效地解决这个问题。

【问题讨论】：

我认为 SO 不是这个问题的正确平台。如果您想查看您的代码/算法（为了提高效率等），Code Review 是适合它的平台。请将您的问题移到那里。
@AjayDabas 我不想审查我的代码/算法。我知道它很慢。我正在寻找另一种更有效的解决方案。
能否给出这个问题的根源？可以在线测试解决方案吗？
@Evg，很遗憾，报名人数有限。
在 main 末尾迭代时，make 为 & fraction 以避免无用的副本。也可以在 main 的开头尝试为 result 预分配空间。哦，不要使用endl

标签： c++ algorithm performance optimization fractions

【解决方案1】：

这是一个非常有趣的问题，我非常喜欢思考它，但我不知道为什么它会收到如此多的反对票。无论如何，下面是我对解决方案的粗略草图。如果我添加一些改进，我稍后会更新我的答案。

正如@user3357359 answer 所建议的，我们需要一种更有效的方法来生成互质数。一种常用技术（在丢番图方程求解器中使用）是Farey Sequence。

定义：n 介于 a/b 和 c/d 之间的 Farey 序列是不可约分数 p/q 的升序（按值）序列，使得 a/b <= p/q <= c/d，其中 @ 987654331@.

Farey 序列的属性在this paper 中提供。所以，我们手头有几个问题：

如何有效生成法雷序列n？
知道第一个元素a0/b0 = a/b和第二个元素a1/b1，可以使用以下算法生成a2/b2（复杂度为O(1)）：
```
k = int((n + b0) / b1)
a2 = k * a1 - a0
b2 = k * b1 - b0
```
如何得到法雷序列的第二个元素n？
我们知道分母在1..n 范围内。我们还知道a/b 和c/d 是Farey 数列中的邻居当且仅当b*c - a*d = 1。因此，通过迭代d 到1..n 的值，我们可以找到最小的分数，它将是序列中的下一个分数。复杂性O(n)。
我们如何决定我们应该生成哪个顺序的法雷序列？
当N 的大小为10^18 时，盲目地生成订单N 是愚蠢的。此外，您将没有足够的内存。我们只需要生成某个阶k 的Farey 序列，使得它的长度大于n 以200000 为界。这是该算法中最难的部分，目前数论中的工具只允许我们估计：|F_k| = 3*k^2 / pi^2。所以我们有：
```
k = ceil(sqrt(n * pi^2)) + C
```
在page 11 of this paper你也可以找到这个公式的近似误差，让你有更多的误差空间，所以C。请注意，对于每个a/b 和c/d，F_k 的长度会有所不同。

综上所述，算法的伪代码为：

1. Estimate the order k of the Farey Sequence to generate, such that |F_k| >= n
2. Calculate the second element of F_k. O(k), where k << n and 1 < n < 200000
3. Generate the whole Farey Sequence. O(n), where 1 < n < 200000
4. Sort by the requirements. O(n log n), where 1 < n < 200000

在最坏的情况下，当您在 step 1 中的估计给出的元素少于 n 时，您需要使用更高阶的 k' 再次生成。这只会使算法的执行时间增加一个常数。因此，总体复杂度平均为O(n log n)，其中n < 200000。

备注：该算法的瓶颈在于k的阶数估计。这可以通过使用Stern-Brocot Tree 而不是法雷序列来完全避免。树完全按照您需要的顺序生成，但我怀疑在a/b != 0/1 和c/d != 1/1 的任意设置中，它会以良好的顺序遍历所有分数。

【讨论】：

刚刚意识到如何更多地指定算法并使序列顺序具有确定性，总体复杂度O((log n + log N) * n)。明天会更新我的答案。
Stern-Brocot 树也不错。它应该可以工作，因为它生成与法雷序列相同的值，只是顺序不同。

【解决方案2】：

您应该遍历从 1 到 N 的所有分母，但您不需要遍历给定范围内的所有分子。

对于任何不可约分数分子/分母，分子和分母都是coprimes。因此，您需要有效地生成 denominator 的所有互质数，其中 co-prime/denominator 位于给定范围内。

例如，对于分母 = 8，互质数是 1,3,5,7。 [0,1] 范围内的不可约分数是 1/8、3/8、5/8、7/8。您可以通过向它们添加一个整数来扩展这个不可约分数的范围：1+1/8、1+3/8、1+5/8、1+7/8 => 9/8、11/8、13/ 8, 15/8是以8为分母的不可约分数，范围为[0+1,1+1] = > [1,2]

用于生成 x 的所有互质数的伪代码，使得 a/b

prime_factors_of_x = get_prime_factors(x)
is_coprime = boolean array of size x
set every element of is_coprime to True
// next, set to False all elements indexed by multiples of a prime factor of x
for every factor in prime_factors_of_x:
   // next loop can be a vectorial operation in some languages
   for i in 0..(x/factor):
      is_coprime[i*factor] = False

all_coprime_list = [] // empty list
min_coprime = floor(a/b * x)+1
max_coprime = floor(c/d * x)-1
for i in min_coprime..max_coprime:
   if is_coprime[i mod x]:
      all_coprimes_list.append(i)

这是总体思路。

【讨论】：

我如何有效地生成所有互质数/分母位于给定范围内的分母互质数？
@Bananon，我希望我的伪代码有帮助。如果可以的话，我会用 C++ 编写它。我是一个蟒蛇人。
感谢这个优化想法。我会赞成你的回答。但我想指出，只有当max_coprime - min_coprime 相对较大时，您的算法才会比为min_coprime..max_coprime 中的每个i 检查gcd(x, i) 更有效。所以我建议根据x 和max_coprime - min_coprime 的值动态选择我们的算法之一。不幸的是，这种优化仍然不足以通过至少a, b, c, d, N ≤ 10^9, n ≤ 1000 组测试。
在我的 python 测试中，一旦你可以将内部循环作为矢量操作执行，这个算法就会更快。它变得几乎是线性的。在 python/numpy 中，我们编写 is_coprime[::factor] = [False]*len(is_coprime[::factor]) 。在第一个解决方案中，在嵌套循环中调用了 GCD，这似乎是二次复杂度。无论如何，为了进一步改进它，我会缓存一些计算结果。

【解决方案3】：

除了另一个问题，你会注意到如果N=10^9它排序是有效的，但是如果N=10^18那么你的代码在乘以f时会溢出，一切都会崩溃。你可以做long long eMax = (long long)(((long double)c) * f / d); （同样适用于下限）来减少这个问题，但它会减慢它的速度。

【讨论】：