tests/Rational_Approximation.test.cpp

View this file on GitHub
Last update: 2025-11-28 10:18:48+07:00
Problem: https://judge.yosupo.jp/problem/rational_approximation

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/rational_approximation"

#include "../misc/macros.h"
#include "../math/SternBrocot.h"

void solve() {
  int n, x, y;
  cin >> n >> x >> y;
  auto f = [&](Frac a) {
    return a.p * y <= a.q * x;
  };
  auto [lo, hi] = SternBrocot::bound(f, n, n);
  if(lo.p * y == lo.q * x) hi = lo;
  cout << lo.p << ' ' << lo.q << ' ' << hi.p << ' ' << hi.q << '\n';
}

signed main() {
  ios::sync_with_stdio(false);
  cin.tie(0);
  int tc = 1;
  cin >> tc;
  while (tc--) solve();
}

#line 1 "tests/Rational_Approximation.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/rational_approximation"

#line 1 "misc/macros.h"
// #pragma GCC optimize("Ofast,unroll-loops")       // unroll long, simple loops
// #pragma GCC target("avx2,fma")                   // vectorizing code
// #pragma GCC target("lzcnt,popcnt,abm,bmi,bmi2")  // for fast bitset operation

#include <bits/extc++.h>
#include <tr2/dynamic_bitset>

using namespace std;
using namespace __gnu_pbds;  // ordered_set, gp_hash_table
// using namespace __gnu_cxx; // rope

// for templates to work
#define all(x) (x).begin(), (x).end()
#define sz(x) (int) (x).size()
#define pb push_back
#define eb emplace_back
using i32 = int32_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;
using ld = long double;
using pii = pair<i32, i32>;
using vi = vector<i32>;

// fast map
const int RANDOM = chrono::high_resolution_clock::now().time_since_epoch().count();
struct chash {  // customize hash function for gp_hash_table
  int operator()(int x) const { return x ^ RANDOM; }
};
gp_hash_table<int, int, chash> table;

/* ordered set
    find_by_order(k): returns an iterator to the k-th element (0-based)
    order_of_key(k): returns the number of elements in the set that are strictly less than k
*/
template <class T>
using ordered_set = tree<T, null_type, less<T>, rb_tree_tag, tree_order_statistics_node_update>;

/*  rope
    rope <int> cur = v.substr(l, r - l + 1);
    v.erase(l, r - l + 1);
    v.insert(v.mutable_begin(), cur);
*/
#line 1 "math/SternBrocot.h"
struct Frac { i64 p, q; };
using Path = vector<pair<char, i64>>;

namespace SternBrocot {
  // 1. Encode: 1/1 -> p/q
  Path encode(i64 p, i64 q) {
    Path res;
    while (p != q) {
      i64 k = (p < q) ? q/p - !(q%p) : p/q - !(p%q);
      res.push_back({p < q ? 'L' : 'R', k});
      if (p < q) q -= k * p; else p -= k * q;
    }
    return res;
  }
  // 2. Decode: Path -> Frac
  Frac decode(const Path& path) {
    i64 lp = 0, lq = 1, rp = 1, rq = 0;
    for (auto [c, k] : path)
      if (c == 'L') rp += k * lp, rq += k * lq;
      else          lp += k * rp, lq += k * rq;
    return {lp + rp, lq + rq};
  }
  // 3. LCA: a/b vs c/d
  Frac lca(i64 a, i64 b, i64 c, i64 d) {
    Path p1 = encode(a, b), p2 = encode(c, d), res;
    for (int i = 0; i < min((int)p1.size(), (int)p2.size()) && p1[i].first == p2[i].first; ++i)
      res.push_back({p1[i].first, min(p1[i].second, p2[i].second)}),
      p1[i].second == p2[i].second ? 0 : (i = p1.size());
    return decode(res);
  }
  // 4. Ancestor: Trả về nút ở độ sâu k trên đường đi từ 1/1 đến p/q
  // k=0 -> 1/1, k=1 -> con trực tiếp của 1/1...
  Frac ancestor(i64 k, i64 p, i64 q) {
    Path path = encode(p, q);
    i64 lp = 0, lq = 1, rp = 1, rq = 0;
    for (auto [c, step] : path) {
      i64 take = min(k, step);
      if (c == 'L') rp += take * lp, rq += take * lq;
      else          lp += take * rp, lq += take * rq;
      k -= take;
      if (k == 0) return {lp + rp, lq + rq};
    }
    return {-1, -1}; // k lớn hơn độ sâu của p/q
  }
  // 5. Range: Tìm khoảng con (L, R) chứa p/q
  pair<Frac, Frac> range(i64 p, i64 q) {
    if (p == 0) return {{0, 1}, {1, 0}};
    i64 lp = 0, lq = 1, rp = 1, rq = 0;
    while (p != q) {
      i64 k = (p < q) ? q/p - !(q%p) : p/q - !(p%q);
      if (p < q) rp += k * lp, rq += k * lq, q -= k * p;
      else       lp += k * rp, lq += k * rq, p -= k * q;
    }
    return {{lp, lq}, {rp, rq}};
  }
  // 6. Bound: Tìm {L, R} sát nhất trên SBT thỏa mãn giới hạn và hàm f
  // f(Frac) -> bool: hàm đơn điệu trên cây (VD: f(x) = x <= target)
  // Trả về {L, R} là 2 phân số kẹp giữa ranh giới T/F của f
  template <class Func>
  pair<Frac, Frac> bound(Func f, i64 maxp, i64 maxq) {
    Frac l{0, 1}, r{1, 0}, m;
    int dir = 1; 
    if (f({1, 1}) == f(l)) l = {1, 1}; else r = {1, 1}, dir = 0;
    while (true) {
      Frac &cur = dir ? l : r, &del = dir ? r : l;
      i64 k = 0;
      for (i64 step = 1; step; ) {
        i64 nk = k + step;
        i64 np = cur.p + nk * del.p, nq = cur.q + nk * del.q;
        if (np <= maxp && nq <= maxq && f({np, nq}) == dir) k = nk, step *= 2;
        else step /= 2;
      }
      cur.p += k * del.p, cur.q += k * del.q;
      m = {l.p + r.p, l.q + r.q};
      if (m.p > maxp || m.q > maxq) break;
      if (f(m)) l = m, dir = 1; else r = m, dir = 0;
    }
    return {l, r};
  }
}
#line 5 "tests/Rational_Approximation.test.cpp"

void solve() {
  int n, x, y;
  cin >> n >> x >> y;
  auto f = [&](Frac a) {
    return a.p * y <= a.q * x;
  };
  auto [lo, hi] = SternBrocot::bound(f, n, n);
  if(lo.p * y == lo.q * x) hi = lo;
  cout << lo.p << ' ' << lo.q << ' ' << hi.p << ' ' << hi.q << '\n';
}

signed main() {
  ios::sync_with_stdio(false);
  cin.tie(0);
  int tc = 1;
  cin >> tc;
  while (tc--) solve();
}