Algorithm (Double Binary Search)

1. Compute $I = \inf\{x\in\mathbb{R}: f(x) > K\}$, where $f(x)$ is the number of all possible fractions generated by $A$ that are less than $x$. It’s easy to see that $f$ is a monotonically non-decresing function on $\mathbb{R}$.
2. The desired answer is $\sup_{(p, q)\in A \times A}\{p / q: p / q < I\}$
3. Note that (1) could be solved using binary search because $f$ is monotone.
4. Note that (2) could be solved using binary search because $A$ is sorted. It could also be solved using a two pointer approach, although it is not straightforward to see at first glance. We settle for binary search.

Time Complexity

1. A single evaluation of $f$ would cost $O(N\log N)$. The search range for (1) is $(0, 1]$. In order to reach the terminating condition for binary search, at least a total of $\frac{\log \epsilon^{-1}}{\log 2}$ searchs are required, where $\epsilon$ is the seach precision which we take to be $10^{-9}$. In total, step(1) is $O\left\(\frac{\log \epsilon^{-1}}{\log 2}\cdot N\log N\right\)$.
2. Step (2) would require $O(N\log N)$ since we need to traverse through $A$ and find for each $A[i]$ the desired sub-supremum which we accumulate over max-monoid to get the final answer.

Memory

1. $O(1)$ since we didn’t allocate any new space that is proportional to the input size.
2. Recursive calls are not counted since they can be optimized into loops by the compiler on most occassions.

Code

template <class F>
struct recursive {
  F f;
  template <class... Ts>
  decltype(auto) operator()(Ts&&... ts)  const { return f(std::ref(*this), std::forward<Ts>(ts)...); }

  template <class... Ts>
  decltype(auto) operator()(Ts&&... ts)  { return f(std::ref(*this), std::forward<Ts>(ts)...); }
};

template <class F> recursive(F) -> recursive<F>;
auto const rec = [](auto f){ return recursive{std::move(f)}; };

class Solution {
 public:
  vector<int> kthSmallestPrimeFraction(vector<int>& A, int K) {
    const int n = size(A);
    const double eps = 1e-9;

    auto outer_binary_search = [&](double lo, double hi, int target, auto f) {
      optional<double> result;
      auto search = rec([&](auto&& search, double lo, double hi) -> void {
        const double mid = lo + (hi - lo) / 2;
        if (std::abs(lo - hi) < eps) return;
        else if (f(mid) <= target) search(mid, hi);
        else if (f(mid) > target) (result = mid,  search(lo, mid));
      });
      return (search(lo, hi), result);
    };

    auto inner_binary_search = [&](int lo, int hi, double target, auto g) {
      optional<int> result;
      auto search = rec([&](auto&& search, int lo, int hi) -> void {
        const int mid = lo + (hi - lo) / 2;
        if (lo > hi) return;
        else if (g(mid) < target) (result = mid, search(mid + 1, hi));
        else if (g(mid) >= target) search(lo, mid - 1);
      });
      return (search(lo, hi), result);
    };

    auto eval_frac = [](const double num, const double denom) { return num / denom; };

    auto f = [&](double x, int acc = 0) {
      for (int i = 0; i < n; ++i) 
        acc += inner_binary_search(0, n - 1, x, [&](auto j) { return eval_frac(A[j], A[i]); }).value_or(-1) + 1;
      return acc;
    };

    auto solution = [&](vector<int> acc = {INT_MIN, 1}) -> vector<int> {
      const auto upper_limit = outer_binary_search(0, 1+eps, K, f);
      for (int i = 0; i < n; ++i) {
        if (const auto match = inner_binary_search(0, n - 1, upper_limit.value() - eps,
                                                   [&](auto x) { return eval_frac(A[x], A[i]); });
            match and eval_frac(acc[0], acc[1]) < eval_frac(A[match.value()], A[i]))
          acc = {A[match.value()], A[i]};
      }
      return acc;
    }();

    return solution;
  }
};