公式

G - Nearest Fraction 解説 by MMNMM


\(r\) を既約分数として表したときの分母 \(D\) が \(N\) 以下の場合、明らかです。 以下では \(N\lt D\) とします。


正整数 \(X\) について、Stern–Brocot 木の頂点のうち

  • 値が区間 \([0,1]\) に含まれる
  • 既約分数 \(p/q\) として表したとき、\(q\leq X\)

を満たす頂点のみをすべて残した木を \(T _ X\) とします。 これは、どのような \(X\) に対しても二分探索木の条件を満たしています。

\(T _ N\) のノードのうち値が \(r\) を超えない最大のもの \(x\) と、値が \(r\) 以上である最小のもの \(y\) について考えます。

例えば、\(N=5,r=0.45\) としたときの \(T _ N\) および \(x,y\) は以下の図のようになります。

\(x,y\) は \(T _ D\) における根 \(\dfrac01\) から \(r\) までのパス上のノードのいずれかと等しいです。

\(r\) の連分数展開を \([0;a _ 1,a _ 2,\ldots,a _ K,a _ {K+1}=1]\) としたとき、根から \(r\) までのパス上のノードの値は \([0;a _ 1,\ldots,a _ k,n]\ (0\leq k\leq K,1\leq n\leq a _ {k+1})\) の形で書ける既約分数です(証明は省略します)。
また、同じ \(k\) に対する \([0;a _ 1,\ldots,a _ k,n]\ (1\leq n\leq a _ {k+1})\) と \(r\) との大小関係は \(k\) の偶奇のみにより、\(n\) が大きいほうが \(r\) に近いです。 \(k\) の偶奇が同じ場合、\(k\) が大きいほうが \(r\) に近いです。 よって、\(x,y\) は \([0;a _ 1,\ldots,a _ k,n]\) を既約分数で表したときの分母が \(N\) を超えないような辞書順最大の \((k,n)\ (0\leq k\leq K,1\leq n\leq a _ {k+1})\) について、\([0;a _ 1,\ldots,a _ k,n]\) と \([0;a _ 1,\ldots,a _ k]\) です。

よって、\(r\) の連分数展開がわかれば \(x,y\) を求めることができます。

求めるべき答えは \(x,y\) のいずれか(の分子と分母からなる整数の \(2\) つ組)です。 \(|r-x|\) と \(|r-y|\) の大小関係を決定することで答えを求めることができます。
この部分の判定を \(64\operatorname{bit}\) 整数や倍精度浮動小数点数・拡張倍精度浮動小数点数のみを用いて正しく行うのは面倒あるいは難しいことに注意してください(実装例では \(64\operatorname{bit}\) 整数を用いて判定を行っています )。 より大きい適切な桁数もしくは適切な精度の数値型を用いるほうが実装が簡単になります。

実装例は以下のようになります。 この実装では、\(r\) の連分数展開を行いながら連分数展開を途中で打ち切ったものの分子・分母を計算しています。

\(r\) を既約分数で表さなくても、\(r\) の連分数展開が「連分数展開を打ち切ったものの分母が \(N\) より大きくなる」より前に終了するかで \(D\leq N\) かを判定することもできます。

#include <utility>
#include <vector>
#include <iostream>
#include <ranges>
#include <atcoder/modint>
#include <atcoder/math>

// -ad+bc = 1 であるような既約分数 x = a/b, y = c/d と有理数 r = p/q (x < r < y)について、
// p/q が a/b と c/d のどちらに近いかを判定する
//
// p/q - a/b < c/d - a/b = 1 / bd より
// (p/q - a/b) * qbd = pbd - qad が (0, q) の範囲の整数になることを利用し、適当な mod で評価したのち復元する
template<class Fraction>
bool compare(const Fraction &r, const Fraction &x, const Fraction &y) {
    const auto&[p, q]{r};
    const auto&[a, b]{x};
    const auto&[c, d]{y};

    using modint1 = atcoder::static_modint<1000000007>;
    using modint2 = atcoder::static_modint<1000000009>;

    std::vector<long long> reminder, modulo;
    reminder.emplace_back((modint1{p} * b * d - modint1{q} * a * d).val());
    modulo.emplace_back(modint1::mod());

    reminder.emplace_back((modint2{p} * b * d - modint2{q} * a * d).val());
    modulo.emplace_back(modint2::mod());

    // |r-x| <= |r-y| iff (p/q - a/b) * qbd <= q/2
    return atcoder::crt(reminder, modulo).first * 2 <= r.second;
}

int main() {
    using namespace std;
    using fraction = pair<unsigned long, unsigned long>;

    string r_str;
    unsigned long N;
    cin >> r_str >> N;

    // r を有理数に
    fraction r{0, 1};
    auto&&[p, q]{r};

    for (const auto digit : r_str | views::drop(2)) {
        (p *= 10) += digit - '0';
        q *= 10;
    }

    // 既約分数に直す
    {
        const auto g{gcd(p, q)};
        p /= g;
        q /= g;
    }

    // すでに分母・分子が N 以下ならそれが答え
    if (q <= N) {
        cout << p << " " << q << endl;
        return 0;
    }

    // x < r < y を満たす最大の x, 最小の y を求める
    const auto&[x, y]{[](const unsigned long N, fraction frac) -> pair<fraction, fraction> {
        auto&&[p, q]{frac};

        pair<unsigned long, unsigned long> numerator{1, 0}, denominator{0, 1};
        auto&&[n1, n2]{numerator};
        auto&&[d1, d2]{denominator};

        bool depth_parity{};
        // 連分数展開
        while (q) {
            const auto quo{p / q};
            const auto max_q{d1 ? (N - d2) / d1 : (N - n2) / n1};

            // 途中で N を越えるとき、x, y がわかる
            if (quo >= max_q)
                if (depth_parity)
                    return {{n1, d1}, {n1 * max_q + n2, d1 * max_q + d2}};
                else
                    return {{n1 * max_q + n2, d1 * max_q + d2}, {n1, d1}};

            numerator = make_pair(n1 * quo + n2, n1);
            denominator = make_pair(d1 * quo + d2, d1);
            frac = make_pair(q, p % q);

            depth_parity ^= 1;
        }
        return {};
    }(N, r)};

    if (compare(r, x, y)) // r <= (x+y)/2 なら
        // x を出力
        cout << x.first << " " << x.second << endl;
    else // そうでなければ
        // y を出力
        cout << y.first << " " << y.second << endl;

    return 0;
}

投稿日時:
最終更新: