From 38e58192a3aff4794fb488f8b640d21b083cb8ce Mon Sep 17 00:00:00 2001
From: hitonanode <32937551+hitonanode@users.noreply.github.com>
Date: Sun, 14 Jul 2024 00:23:24 +0900
Subject: [PATCH] add rational approximation

---
 number/rational_approximation.hpp           | 66 +++++++++++++++++++++
 number/rational_approximation.md            | 22 +++++++
 number/test/rational_approximation.test.cpp | 21 +++++++
 3 files changed, 109 insertions(+)
 create mode 100644 number/rational_approximation.hpp
 create mode 100644 number/rational_approximation.md
 create mode 100644 number/test/rational_approximation.test.cpp

diff --git a/number/rational_approximation.hpp b/number/rational_approximation.hpp
new file mode 100644
index 00000000..bc1102fe
--- /dev/null
+++ b/number/rational_approximation.hpp
@@ -0,0 +1,66 @@
+#pragma once
+
+#include <cassert>
+#include <utility>
+#include <vector>
+
+// Rational approximation based on Stern–Brocot tree
+// Input: N > 0, num >= 0, den >= 0 (num > 0 or den > 0)
+// return {{p, q}, {r, s}} where p/q <= num/den <= r/s. Strictly,
+// - p/q = max(n / d | n / d <= num / den, 0 <= n <= N, 1 <= d <= N)
+// - r/s = min(n / d | num / den <= n / d, 0 <= n <= N, 1 <= d <= N)
+template <class Int>
+std::pair<std::pair<Int, Int>, std::pair<Int, Int>> rational_approximation(Int N, Int num, Int den) {
+    assert(N >= 1);
+    assert(num >= 0);
+    assert(den >= 0);
+    assert(num > 0 or den > 0);
+
+    if (num == Int(0)) return {{Int(0), Int(1)}, {Int(0), Int(1)}}; // 0
+    if (den == Int(0)) return {{Int(1), Int(0)}, {Int(1), Int(0)}}; // inf
+
+    // p/q <= num/den <= r/s
+    Int p = 0, q = 1, r = 1, s = 0;
+
+    bool is_left = false;
+    while (true) {
+        Int max_steps = num / den;
+
+        if (is_left) {
+            // r + steps * p <= N
+            // s + steps * q <= N
+
+            if (p > Int(0)) max_steps = std::min(max_steps, (N - r) / p);
+            max_steps = std::min(max_steps, (N - s) / q);
+
+            r += max_steps * p;
+            s += max_steps * q;
+        } else {
+            // p + steps * r <= N
+            // q + steps * s <= N
+
+            max_steps = std::min(max_steps, (N - p) / r);
+            if (s > Int(0)) max_steps = std::min(max_steps, (N - q) / s);
+
+            p += max_steps * r;
+            q += max_steps * s;
+        }
+
+        if (is_left and !max_steps) break;
+
+        num -= max_steps * den;
+
+        if (num == Int(0)) {
+            if (is_left) {
+                return {{r, s}, {r, s}};
+            } else {
+                return {{p, q}, {p, q}};
+            }
+        }
+
+        std::swap(num, den);
+        is_left = !is_left;
+    }
+
+    return {{p, q}, {r, s}};
+}
diff --git a/number/rational_approximation.md b/number/rational_approximation.md
new file mode 100644
index 00000000..89f646b8
--- /dev/null
+++ b/number/rational_approximation.md
@@ -0,0 +1,22 @@
+---
+title: Rational approximation （有理数近似）
+documentation_of: ./rational_approximation.hpp
+---
+
+与えられた非負の有理数 `num / den` に対して，分子・分母がともに $N$ 以下の分数で両側から近似する．
+
+## 使用方法
+
+```cpp
+int N, num, den;
+const auto [l, r] = rational_approximation<int>(N, num, den);
+
+// lnum / lden <= num / den <= rnum / rden
+// max(lnum, lden, rnum, rden) <= N
+const auto [lnum, lden] = l;
+const auto [rnum, rden] = r;
+```
+
+## 問題例
+
+- [Library Checker: Rational Approximation](https://judge.yosupo.jp/problem/rational_approximation)
diff --git a/number/test/rational_approximation.test.cpp b/number/test/rational_approximation.test.cpp
new file mode 100644
index 00000000..fe875630
--- /dev/null
+++ b/number/test/rational_approximation.test.cpp
@@ -0,0 +1,21 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/rational_approximation"
+#include "../rational_approximation.hpp"
+
+#include <cassert>
+#include <iostream>
+using namespace std;
+
+int main() {
+    cin.tie(nullptr), ios::sync_with_stdio(false);
+
+    int T;
+    cin >> T;
+    while (T--) {
+        int N, x, y;
+        cin >> N >> x >> y;
+        const auto [l, r] = rational_approximation<int>(N, x, y);
+        const auto [lnum, lden] = l;
+        const auto [rnum, rden] = r;
+        cout << lnum << ' ' << lden << ' ' << rnum << ' ' << rden << '\n';
+    }
+}