Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add rational approximation #340

Merged
merged 1 commit into from
Jul 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 66 additions & 0 deletions number/rational_approximation.hpp
Original file line number Diff line number Diff line change
@@ -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}};
}
22 changes: 22 additions & 0 deletions number/rational_approximation.md
Original file line number Diff line number Diff line change
@@ -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)
21 changes: 21 additions & 0 deletions number/test/rational_approximation.test.cpp
Original file line number Diff line number Diff line change
@@ -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';
}
}
Loading