Stern–Brocot Tree についてのメモ

この記事の目標

SBT についてざっくり理解して、次の機能を実装できるようになります。

Library Checker - Stern–Brocot Tree
ABC333 G - Nearest Fraction のような問題を SBT 上の探索で解く

おことわり

新規性は特になく、既存の記事を自分用にまとめた程度の内容になっています。
気持ち寄りの説明が多めで、議論がやや雑な部分があるかもしれません。

Stern–Brocot Tree について

定義

Stern–Brocot Tree (SBT) は、次のように定められる、（無限に続く）完全二分木です。

各頂点は非負整数の $4$ つ組 $(p, q, r, s)$ を持ち、正の有理数 $\dfrac{p+r}{q+s}$ が対応する。*1
根は $(0, 1, 1, 0)$ である（対応する有理数は $\dfrac{1}{1}$ ）。
$(p, q, r, s)$ の左の子は $(p, q, p+r, q+s)$ 、右の子は $(p+r, q+s, r, s)$ である。

性質

Stern–Brocot Tree には次の重要な性質があります。

二分探索木である。つまり、 $(左の子孫の値) \lt (親の値) \lt (右の子孫の値)$ である。（cf. 加比の理）
根から右の子に $a_0$ 回、左の子に $a_1$ 回、 $\ldots$ 、( $k-1$ が奇数 ? 右の子 : 左の子) に $a_{k-1}$ 回進んだ頂点が対応する有理数は、連分数 $a _ 0 + \dfrac{1}{a _ 1 + \dfrac{1}{\ddots + \dfrac{1}{a _ {k-1} + 1}}}$ として表される。
正の有理数がすべて現れる。（2. からわかる。すべての正の有理数は連分数展開を持つので）
逆に、各頂点に対応する正の有理数は相異なる。（1. からわかる。あるいは、正の有理数の連分数展開は $a_i \geq 1 \: (i \geq 1)$ に限定すれば一意なので、と思ってもよい）
根から始めて、有理数 $\dfrac{p}{q}$ に対応する頂点に到達するまでのパスを考える。「右の子/左の子に $d$ 回進む」ことをひとまとまりの操作とみなしたときの操作回数（パスの連長圧縮の長さ、2. での $k$ の値 *2）は $O(\log(p+q))$ である。（2. からわかる。連分数展開は実質的に互除法なので）
頂点 $(p, q, r, s)$ の子孫の集合は、開区間 $\left(\dfrac{p}{q}, \dfrac{r}{s} \right)$ になる（ $\dfrac{1}{0}$ は $\infty$ とみなす）。（気持ちとしては、左の子に $d$ 回進むと $\dfrac{p + (dp + r)}{q + (dq + s)}$ 、右の子に $d$ 回進むと $\dfrac{(p + dr) + r}{(q + ds) + s}$ で、 $d \to \infty$ とするとそれぞれ $\dfrac{p}{q}$ 、 $\dfrac{r}{s}$ ）

Library Checker - Stern–Brocot Tree

ENCODE_PATH

有理数から SBT 上のパス（の連長圧縮）を得る問題です。性質 2. を利用して連分数展開すればよいです。（性質 5. より、出力は多くなりません）

DECODE_PATH

SBT 上のパス（の連長圧縮）から有理数を得る問題ですが、素直に根から辿ると有理数だけでなく頂点の $4$ 整数 $(p, q, r, s)$ を得ることができます（そう実装すると後々楽です）。（これも性質 5. より、入力は多くなりません）

LCA

$2$ つの有理数の LCA（に対応する有理数）を求める問題です。 $2$ つの有理数を ENCODE_PATH したのち、根から辿っていってどこまで一致するか見ればよいです。性質 5. より間に合います。

ANCESTOR

有理数 $\dfrac{p}{q}$ の祖先であって深さ $d$ の頂点（に対応する有理数）を求める問題です。これも $\dfrac{p}{q}$ を ENCODE_PATH したのち、根から深さ $d$ まで辿っていけば性質 5. より間に合います。

RANGE

有理数が与えられて、その子孫の下限と上限を求める問題です。性質 6. そのものです。ENCODE_PATH してから DECODE_PATH すると $(p, q, r, s)$ が得られるので、それが答えです。

ABC333 G - Nearest Fraction

SBT 上を探索する話です。

より一般に、次の形で考えましょう。（この定式化は [2], [3] をかなり参考にしています）

関数 $f \colon \mathbb{R} _ {\gt 0} \to \{ \mathrm{true}, \mathrm{false} \}$ は単調性がある。すなわち、ある $\alpha \in \mathbb{R} _ {\gt 0}$ が存在して、 $x \leq \alpha \implies f(x) = \mathrm{true}, x \gt \alpha \implies f(x) = \mathrm{false}$ といった形になっている（等号・不等号や true/false は入れ替わってもよい）。このとき、境界 $\alpha$ を分母・分子が $n$ 以下の有理数で近似したい。 $\dfrac{p}{q} \leq \alpha \lt \dfrac{r}{s}$ *3 を満たす $0 \leq p,q,r,s \leq n$ のうち、 $\alpha - \dfrac{p}{q}$ および $\dfrac{r}{s} - \alpha$ が最小になるものを求めよ。

性質 1. を利用して、二分探索木を探索することを考えます。性質 6. より、探索の途中で頂点 $(p, q, r, s)$ にいるとき、 $(p, q, r, s)$ は $\dfrac{p}{q} \leq \alpha \lt \dfrac{r}{s}$ を満たしています。よって、素直に二分探索木を探索して、 $n$ を超えたら打ち切る、とすると $O(n)$ で解けます。性質 5. を利用して、「その方向にどれだけ進むか」を二分探索 *4 すれば速くなります（log ふたつに見えて実は $O(\log n)$ です。証明は [2]）。

実装で自分が苦労した細かい部分を解説します。まず、頂点 $(p, q, r, s)$ に到達したとき、そこからはもう進まずに $(p, q, r, s)$ を答えとする条件は $p + r > n$ または $q + s > n$ です。 $(p + r, q + s)$ は、左および右に $1$ 個進んだときの $(r, s)$ および $(p, q)$ だからです。

同様に、進む深さ $d$ を二分探索するときの $d$ の上限 $\mathrm{ng}$ は、（左に進むのであれば） $dp + r > n$ または $dq + r > n$ を満たす最小の $d$ 、というようにして決めることができます *5。

こういったことは分母・分子が $n$ 以下という制限が強い場面でなければ（例：[2] のように、真の値が有理数で、かつ分母・分子がある値以下であることが示せる場合）気にしなくてもあまり問題はありませんが、今回の問題のような場面では出力が $n$ を超えてしまうバグの原因になるので、注意して実装する必要があります。

実装例

クリックして展開

#include <iostream>
#include <vector>
#include <functional>
#include <cassert>
using namespace std;

// stern-brocot tree
namespace sbt
{
  // (p, q, r, s) から (isleft ? 左 : 右) に d 個進んだ頂点を返す
  template<class T>
  tuple<T, T, T, T> child(T d, T p, T q, T r, T s, bool isleft)
  {
    if (isleft)
      r += d * p, s += d * q;
    else
      p += d * r, q += d * s;
    return make_tuple(p, q, r, s);
  }

  // (p, q, r, s) の親を返す
  // (p, q, r, s) が根なら (0, 0, 0, 0) を返す
  template<class T>
  tuple<T, T, T, T> parent(T p, T q, T r, T s)
  {
    if (p == 0 && q == 1 && r == 1 && s == 0)
      return make_tuple(0, 0, 0, 0);
    if (p < r || q < s)
      r -= p, s -= q;
    else
      p -= r, q -= s;
    return make_tuple(p, q, r, s);
  }

  // 1/1 から p/q へのパスが右に a_1 回、左に a_2 回、… としたときの a を返す
  // (a_1, …, a_k + 1) は p/q の連分数展開になっている
  template<class T>
  vector<T> encode_path(T p, T q)
  {
    vector<T> a;
    if (p < q)
    {
      a.emplace_back(0);
      swap(p, q);
    }
    while (p != 1)
    {
      a.emplace_back(p / q);
      p %= q;
      swap(p, q);
    }
    if (!a.empty())
    {
      if (a.back() == 1)
        a.pop_back();
      else
        a.back()--;
    }
    return a;
  }

  // 1/1 から右に a_1 回、左に a_2 回、… 移動した頂点 (p, q, r, s) を返す
  // (a_1, …, a_k + 1) は (p+r)/(q+s) の連分数展開になっている
  template<class T>
  tuple<T, T, T, T> decode_path(const vector<T> &a)
  {
    T p = 0, q = 1, r = 1, s = 0;
    for (int i = 0; i < (int)a.size(); i++)
      tie(p, q, r, s) = child(a[i], p, q, r, s, i & 1);
    return make_tuple(p, q, r, s);
  }

  // p/q と r/s の LCA (f, g, h, k) を返す
  template<class T>
  tuple<T, T, T, T> lca(T p, T q, T r, T s)
  {
    vector<T> a = encode_path(p, q), b = encode_path(r, s);
    int n = min(a.size(), b.size());
    T f = 0, g = 1, h = 1, k = 0;
    for (int i = 0; i < n; i++)
    {
      T c = min(a[i], b[i]);
      tie(f, g, h, k) = child(c, f, g, h, k, i & 1);
      if (a[i] != b[i])
        break;
    }
    return make_tuple(f, g, h, k);
  }

  // p/q の祖先であって深さが d のノード (f, g, h, k) を返す
  // 存在しなければ (0, 0, 0, 0)
  template<class T>
  tuple<T, T, T, T> ancestor(T d, T p, T q)
  {
    vector<T> a = encode_path(p, q);
    T f = 0, g = 1, h = 1, k = 0;
    for (int i = 0; i < (int)a.size(); i++)
    {
      T c = min(d, a[i]);
      tie(f, g, h, k) = child(c, f, g, h, k, i & 1);
      d -= c;
      if (d == 0)
        break;
    }
    return d == 0 ? make_tuple(f, g, h, k) : make_tuple(0, 0, 0, 0);
  }

  // p/q の子孫の下限 f/g と上限 h/k を返す
  // p/q に対応する頂点 (f, g, h, k) でもある
  template<class T>
  tuple<T, T, T, T> range(T p, T q) { return decode_path(encode_path(p, q)); }

  // judge: [0, ∞] → {true, false}
  // judge は単調性を持つ、すなわち、ある実数 α を境に true と false が切り替わる
  // 分母・分子が max_value 以下の有理数のうち、α に最も近いもの (下側: p/q, 上側: r/s) を返す
  template<class T>
  tuple<T, T, T, T> search(const function<bool(T, T)> &judge, const T &max_value)
  {
    T p = 0, q = 1, r = 1, s = 0;
    const bool judge_01 = judge(0, 1), judge_10 = judge(1, 0);
    assert(judge_01 != judge_10);
    while (p + r <= max_value && q + s <= max_value)
    {
      const bool judge_now = judge(p + r, q + s);
      const bool isleft = judge_now ^ judge_01;
      T maxd;
      if (isleft)
        maxd = p == 0 ? (max_value - s) / q : min((max_value - r) / p, (max_value - s) / q);
      else
        maxd = s == 0 ? (max_value - p) / r : min((max_value - q) / s, (max_value - p) / r);
      T dl = 1, dr = 2;
      while (dr <= maxd)
      {
        auto [np, nq, nr, ns] = child(dr, p, q, r, s, isleft);
        if (judge(np + nr, nq + ns) != judge_now)
          break;
        dl = dr, dr = min(2 * dr, maxd + 1);
      }
      while (dr - dl > 1)
      {
        T dm = dl + (dr - dl) / 2;
        auto [np, nq, nr, ns] = child(dm, p, q, r, s, isleft);
        if (judge(np + nr, nq + ns) == judge_now)
          dl = dm;
        else
          dr = dm;
      }
      tie(p, q, r, s) = child(dl, p, q, r, s, isleft);
    }
    return make_tuple(p, q, r, s);
  }
}

// 使用例 (https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1208)
int main()
{
  while (true)
  {
    int p, n;
    cin >> p >> n;
    if (p == 0 && n == 0)
      break;
    const function<bool(int, int)> judge = [&](int x, int y) -> bool
    { return y * y * p < x * x; };
    auto [u, v, x, y] = sbt::search(judge, n);
    cout << x << "/" << y << " " << u << "/" << v << endl;
  }
}