Stern-Brocot树/Farey序列

What · 是什么

Stern–Brocot 树是枚举全部正的最简分数的二叉搜索树。以三元组 (0/1, 1/1, 1/0) 为根，每个节点 (a/b, p/q, c/d) 的左右子为 (a/b, (a+p)/(b+q), p/q) 与 (p/q, (p+c)/(q+d), c/d)。相邻分数满足 bc-ad=1。应用：在树上二分逼近实数，找分母不超过某上限的最优分数逼近；Farey 序列(分母≤n 的 [0,1] 最简分数)是它的有限前缀；与连分数紧密相关，可 O(log(p+q)) 定位分数。

Why · 为什么

How · 实现

伪代码

# 分数逼近: 在 (a/b, c/d) 区间中位数逼近 x, 分母上限 D
# 朴素: 每次比较中位数与 x
a, b, c, d = 0, 1, 1, 0      # 0/1 与 1/0
while True:
  p, q = a + c, b + d        # 中位分数
  if b + d > D: break        # 分母超限
  if p / q < x:
    a, b = p, q              # 取右半
  elif p / q > x:
    c, d = p, q              # 取左半
  else:
    a, b, c, d = p, q, p, q  # 恰好命中
# 比较 (a/b) 与 (c/d) 谁更接近 x

Python 实现

def find_fraction(x, D):
    """在分母<=D 的最简分数中, 二分逼近实数 x。
    返回 (p, q) 为最佳逼近(比较左右边界距离)。
    """
    a, b, c, d = 0, 1, 1, 0  # 区间 (a/b, c/d)
    while True:
        p, q = a + c, b + d
        if q > D:
            break
        if p < x * q:
            a, b = p, q
        elif p > x * q:
            c, d = p, q
        else:
            a, b = p, q
            c, d = p, q
            break
    # 比较左右边界哪个更接近 x
    lv = abs(a - x * b) * d   # 交叉乘比较 |a/b - x| 与 |c/d - x|
    rv = abs(c - x * d) * b
    if lv <= rv:
        return a, b
    return c, d


def farey_sequence(n):
    """生成第 n 阶 Farey 序列([0,1] 分母<=n 的最简分数升序)。"""
    a, b, c, d = 0, 1, 1, n
    res = [(a, b)]
    while c <= n:
        k = (n + b) // d
        a, c, b, d = c, k * c - a, d, k * d - b
        res.append((a, b))
    return res


def stern_diatomic(n):
    """Stern 双原子序列 a_n (OEIS A002487):
    a0=0,a1=1, a2n=a_n, a2n+1=a_n+a_{n+1}。
    用 fusc 迭代实现。
    """
    if n == 0:
        return 0
    a, b = 1, 0
    while n:
        if n & 1:
            b += a
        else:
            a += b
        n >>= 1
    return b


if __name__ == '__main__':
    # 3/4 = 0.75, 分母<=4 最佳逼近应为 3/4 本身
    assert find_fraction(0.75, 4) == (3, 4)
    # 无理数 sqrt(2)≈1.4142, 分母<=10 的最佳逼近
    import math
    sq2 = math.sqrt(2)
    p, q = find_fraction(sq2, 10)
    # 已知分母<=10 的最佳为 17/12? 分母超界; 手动枚举验证
    best_p, best_q, best_err = 0, 1, float('inf')
    for dd in range(1, 11):
        pp = round(sq2 * dd)
        if abs(pp / dd - sq2) < best_err:
            best_p, best_q, best_err = pp, dd, abs(pp / dd - sq2)
    assert p == best_p and q == best_q, (p, q, best_p, best_q)

    # Farey 序列
    f4 = farey_sequence(4)
    expected = [(0, 1), (1, 4), (1, 3), (1, 2), (2, 3), (3, 4), (1, 1)]
    assert f4 == expected, f4
    assert len(farey_sequence(5)) == 11  # 1+sum_{k<=5} phi(k)=1+10

    # Stern 双原子序列
    sd = [stern_diatomic(i) for i in range(16)]
    assert sd == [0, 1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4], sd

    # 分数逼近: 1/pi 或任意有理数
    assert find_fraction(1 / 3, 20) == (1, 3)
    assert find_fraction(0.1, 50) in ((1, 10), (1, 10))

    print('sqrt(2) best frac (denom<=10):', find_fraction(sq2, 10))
    print('Farey(4):', f4)
    print('Stern diatomic [0..15]:', sd)
    print('All stern-brocot tests passed!')

易错点

• 中位分数 (a+c)/(b+d) 在分母超限时即为最佳逼近的候选, 需与两端点比距离, 别直接返回中位数
• 交叉相乘比较距离: |a/b - x| 与 |c/d - x| 用 abs(a-x*b)*d 与 abs(c-x*d)*b, 别用浮点除法
• 区间初始为 (0/1, 1/0), 1/0 是无穷大记号, 分母为0 不是错误, 求 q 前判断
• Farey 递推 k=(n+b)//d 的下取整必须精确, 别四舍五入
• Stern 双原子序列的递推 a2n=a_n 是基于编号从1开始, 与 Calkin-Wilf 编号对应, 下标从0开始需特判

例题

信源

• OI Wiki✓ https://oi-wiki.org/math/number-theory/stern-brocot/
• OI Wiki✓ https://raw.githubusercontent.com/OI-wiki/OI-wiki/master/docs/math/number-theory/stern-brocot.md
• CP-Algorithms✓ https://cp-algorithms.com/others/stern_brocot_tree_farey_sequences.html
• web✓ https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree