Finding irreducible fractions between two given fractions in ascending order of denominators and then numerators

Question

Time limit per test: 2 seconds
Memory limit per test: 512 megabytes

You are given two fractions a/b < c/d and a positive number N. Consider all irreducible fractions e/f such that 0 < e, f ≤ N and a/b < e/f < c/d. Let s be a sequence of these fractions sorted in ascending order of denominators and then numerators (fraction e1/f1 precedes e2/f2 if either f1 < f2 or f1 = f2 and e1 < e2). You should print first n terms of the sequence s or the whole sequence s if it consists of fewer than n terms.

Input
The first line of each test contains 6 integers a, b, c, d, N, n (0 ≤ a ≤ 10^18, 1 ≤ b, c, d, N ≤ 10^18, 1 ≤ n ≤ 200 000, a/b < c/d).

Output
First, print how many terms of sequence s you will output. And then output these terms in the right order.

Examples

• Input:

  0 1 1 1 5 10
  

  Output:

  9
  1 2
  1 3
  2 3
  1 4
  3 4
  1 5
  2 5
  3 5
  4 5
  

• Input:

  55 34 68 42 90 1
  

  Output:

  1
  89 55
  

• Input:

  49 33 45 30 50 239
  

  Output:

  0
  

So far, I've only managed to write a solution that iterates over all the denominators from 1 to N, and for each denominator iterates over all the numerators from a*f/b to c*f/d, adding all found irreducible fractions to the answer.

Here is my code:

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

Unfortunately, this solution is too slow. I wonder how to solve this problem more efficiently.

Answer 1

This is a very interesting question and I really enjoyed thinking about it, and I don't know why it received so many downvotes. Anyways, the below is my rough sketch of the solution. I will update my answer later if I add some refinement to it.

---

As suggested by @user3357359 answer, we need a more effective way to generate coprimes. A common technique (that is used in diophantine equation solvers) is Farey Sequence.

Definition: a Farey Sequence of order n between a/b and c/d is the ascending (by value) sequence of irreducible fractions p/q such that a/b <= p/q <= c/d, where q < n.

Properties of Farey Sequence is presented in this paper. So, we have several questions at hand:

• How to effectively generate Farey Sequence of order n?
  Knowing the first element a0/b0 = a/b and the second element a1/b1, one can generate a2/b2 using the following algorithm (with complexity O(1)):

  k = int((n + b0) / b1)
  a2 = k * a1 - a0
  b2 = k * b1 - b0
  

• How to get the second element of Farey Sequence of order n?
  We know that the denominator is somewhere in the range 1..n. We also know that a/b and c/d are neighbors in Farey Sequence if and only if b*c - a*d = 1. So by itering the value of d through 1..n we can find the smallest fraction, which will the the next one in the sequence. Complexity O(n).

• How do we decide which order of Farey Sequence we should generate?
  Blindly generating of order N, when N is in magnitude of 10^18 is stupid. Moreover, you won't have enough memory for that. We only need to generate a Farey Sequence of some order k, such that the length of it is larger than n which is bounded by 200000. This is the hardest part of this algorithm, and right now the tools in number theory only allows us to estimate: |F_k| = 3*k^2 / pi^2. So we have:

  k = ceil(sqrt(n * pi^2)) + C
  

  In page 11 of this paper you can also find the approximation error of this formula, so that you accomodate more room for error, thus the C. Note that, for each a/b and c/d, the length of F_k would be different.

To sum up, the pseudocode for the algorithm is:

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

In the worst case scenario when your estimation in step 1 gave less elements than n, you'll need to generate again, using a bit higher order k'. Which will only increase the execution time of your algorithm by a constant. So, the overall complexity is O(n log n) in average, where n < 200000.

---

Remarks: the bottleneck of this algorithm is estimation of order k. This might be avoided completely by using not Farey Sequence, but the Stern-Brocot Tree. The tree generates exactly in the order that you need, but I doubt that in an arbitrary setting with a/b != 0/1 and c/d != 1/1 it will iterate through all the fractions in a good order.