Stern-Brocot 树及其应用

• 下载源代码 - 6.1 KB

引言

有时，在数学的某些角落，你会发现与你所寻找的东西完全不同的事物。19 世纪的法国钟表匠 Achilles Brocot 在寻找制作钟表齿轮的最佳方法时，就发生了这样的事情，他找到了有理数的最优有理数逼近。为什么这会引起你的兴趣呢？嗯，尽管如今大多数处理器都有浮点单元，但处理整数有时也有其优势。在这些情况下，了解你可以将 π 近似为 355/113，精确到小数点后 7 位，这是很有用的。你可以从维基百科上获得这些信息，但如果你的算法需要 ln(2) 或 φ 呢？

你可以使用 Stern 和 Brocot 先生发现的树状结构来找到这些答案。树的前几层在下图所示

数学背景

如果你只对编程细节感兴趣，可以跳过这部分，直接到应用部分。

有理数是可以表示为 p/q 分数的数字。由于存在许多表示同一有理数的约数，我们将只选择 p 和 q 没有公因数的约数。例如，2/3 和 4/6 表示相同的数字，但我们将只选择前者。如果我们把整数放在两个轴上，我们可以将有理数想象成从原点到 p/q 点的直线的斜率。下图显示了“2/3”的表示。

我们现在需要定义中项的概念：两个有理数 m/n 和 m'/n' 的中项是 (m+m')/(n+n')。下面表示了数字 1/2 和 1/3 以及它们的中项 2/5。

将这些数字视为向量加法，可以看到两个数字的中项介于它们之间。用数学术语来说，如果 m/n < m'/n'，则 m/n < (m+m')/(n+n') < m'/n'；

现在我们可以开始构建 Stern-Brocot 树。我们从数字 0/1 和 1/0 开始。1/0 在技术上不是一个分数，但将其视为“无穷大”的表示很有用。

在树的每一层，我们都添加上一层的数字的中项。在第二层，我们放置这些数字的中项 (0+1)/(1+0) = 1/1。在下一层，我们放置这些数字的每对中项：(0+1)/(1+1) = 1/2 和 (1+1)/(1+0) = 2/1。

下一张图片展示了树的连续层如何扩展以覆盖整数网格。

灰色点是 p 和 q 有公因数的地方。

我们已经看到有理数如何被看作是通过整数网格的一个点延伸出去的直线。无理数则可以看作是永远无法精确击中一个整数点的直线。寻找一个近似值就转化为寻找一个足够接近的整数点。

Stern-Brocot 树在有理数逼近中的应用

我们已经看到，我们的树的任何一层都包含前一层值的中间项，并且中间项的值介于其祖先的值之间。在这种树中，使用以下算法可以搜索数字 N 的一个足够好的近似值，容差为 ε。

1. 令当前节点 T 为树的根。
2. 如果 |N-T| < ε，则退出，答案为 T。
3. 如果 N < T，将 T 替换为 T 的左子节点，否则将 T 替换为 T 的右子节点。
4. 转到步骤 2。

该算法，或者更准确地说，Stern-Brocot 树的性质，保证找到的近似值 p/q 是“最佳”的，因为任何其他具有相同精度的 p'/q' 近似值都将具有 p'>p 和 q'>q。

用法

sba 程序计算具有给定精度的最佳有理数逼近。命令行是

sba <number> <precision>

以下是计算 π 的 7 位小数近似值的示例

实现

这里没什么花哨的！每个树节点都由一个 `sbnode` 结构表示，该结构包含 p 和 q 值、左右子节点以及指向父节点的指针。

struct sbnode {
  int p, q;
  struct sbnode *left, *right, *parent;
  sbnode (sbnode *parent_) 
  {
    p = q = 0; 
    left = right = 0; 
    parent = parent_;
  }
  ~sbnode () 
  {
    if (left)
      delete left;
    if (right)
      delete right;
  }
}

构造函数和析构函数负责初始化节点和清理。

该程序实现了前面描述的算法。但是，它不预先创建树。相反，子节点仅在需要时创建。

为了显示素数分解，我们使用 Eratosthenes 筛法创建素数表。这不是寻找素数因子的最快方法，但对于我们处理的小数字来说，它完全足够了。

结语 - 时间之轮

读完这个故事后，你可能会想，为什么一个钟表匠会对这些树感兴趣。如果两个齿轮，一个有 p 齿，另一个有 q 齿，连接在一起，它们的转动比（齿轮比）就是 p/q。

齿轮可以像下面的双齿轮示例一样组合

在这种情况下，如果第一对（红/蓝）的齿轮比是 `p1/q1`，第二对（黄/绿）的齿轮比是 `p2/q2`，那么组合齿轮比（红/绿）是 `p1/q1 * p2/q2`。

作为过去钟表匠遇到的问题的一个例子，一位 18 世纪的钟表匠 Charles-Étienne-Louis Camus 提出了以下问题：“确定一个机构的齿轮和齿轮的齿数…当一个齿轮被放置在小时轮上的齿轮驱动时，将使一个齿轮在一个平均年内转动一圈，平均年被假定为 365 天 5 小时 49 分钟”。如果我们把所有东西都换算成分钟，这个问题就需要设计一个齿轮系，其比率为 720/525949（根据 Camus 的说法，小时轮在 12 小时 = 720 分钟内完成一次旋转，而平均回归年有 525949 分钟）。

制造有 525949 个齿的齿轮是难以想象的，而且在这种情况下制造复合齿轮会受到 525949 是一个素数的阻碍。Camus 不得不找到一个接近所需比率但数字易于分解的组合。

我们的小 `sba` 工具提供了一个答案

~sba -f 0.0013689540240594 10
Finding approximation of 0.0013689540 with 10 decimals
....
Current approximation 97/70857 = 0.0013690 (err=-3.49e-10)
97 is prime
Factors of 70857 = 3^2*7873

Current approximation 130/94963 = 0.0013690 (err=-2.00e-10)
Factors of 130 = 2*5*13
Factors of 94963 = 11*89*97

Current approximation 163/119069 = 0.0013690 (err=-1.12e-10)
163 is prime
119069 is prime

Found fraction 196/143175 = 0.00136895408
Error= -5.31e-11
Factors of 196 = 2^2*7^2
Factors of 143175 = 3*5^2*23*83

这个数字可以分解为 `4/25 * 7/69 * 7/83`。Camus 用 20 页艰苦的工作才得到相同的结果。

新进展 - Stern-Brocot 树与安提凯希拉机械

安提凯希拉机械于 1901 年被发现，是一个希腊天文计算器，可追溯到公元前二世纪或一世纪。自发现以来，它一直让许多研究人员着迷，试图找出这个技术奇迹的内部运作。

最近，伦敦大学学院（UCL）的一个团队在《自然》杂志上发表了一篇关于该机械完整重建的文章。

为了找到齿轮比，作者提出了一种称为 Parmenides 命题的工艺。

我们开发了一种关于金星和土星周期是如何被发现的新理论，并将其应用于恢复丢失的行星周期。柏拉图（公元前五至四世纪）的一段对话以哲学家 Elea 的 Parmenides（公元前六至五世纪）命名。这描述了Parmenides 命题。

在近似 θ 时，假设有理数 p/q 和 r/s 满足 p/q < θ < r/s。那么 (p + r)/(q + s) 是一个介于 p/q 和 r/s 之间的新估计。

如果它是低估，那么它比 p/q 更好。如果它是高估，那么它比 r/s 更好。

这个描述与 Stern-Brocot 树中中间项的定义和搜索算法完全匹配。

使用Parmenides 命题，作者描述了寻找金星齿轮比的过程。它必须是接近 720/1151 的数字，并且具有方便的素数分解。下面是 `sba` 工具输出的一个片段。

最好的数字是 289/462。

接着，作者为所有行星建立了不同的齿轮比，并创建了该机械的完整重建。下面是整个齿轮组的爆炸图：

我发现从柏拉图的对话、希腊天文学家、18 世纪的钟表匠和当代的浮点近似追踪 Stern-Brocot 树的结构非常迷人。

参考文献

• [1] Graham, R. L., Knuth, D. E.,, Patashnik, O. (1994). Concrete Mathematics: A Foundation for Computer Science 2nd Edition. Reading: Addison-Wesley。
• [2] 维基百科 Stern-Brocot 树页面 (https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree)
• [3] 维基百科 π 的近似值页面 (https://en.wikipedia.org/wiki/Approximations_of_%CF%80)
• [4] 维基百科齿轮列页面 (https://en.wikipedia.org/wiki/Gear_train)
• [5] Brian Hayes, "On the Teeth of Wheels" in American Scientist, Vol. 88, No. 4, July-August 2000, pages 296-300。
• [6] Freeth, T., Higgon, D., Dacanalis, A. et al. A Model of the Cosmos in the ancient Greek Antikythera Mechanism. Sci Rep 11, 5821 (2021)。https://doi.org/10.1038/s41598-021-84310-w

历史

• 2021 年 1 月 19 日 - 初始版本
• 2021 年 3 月 22 日 - 添加了素数分解代码；添加了安提凯希拉机械部分