Stern-Brocot Trees and their Applications

• Download source code - 6.1 KB

Introduction

Sometimes, in some corners of math, you find things completely different from what you were looking for. This is what happened to a 19^th century French clock-maker, Achilles Brocot who was looking for the best way to make geared wheels for his clocks and found the best rational approximation for irrational numbers. Why would this be of interest to you? Well, despite having floating point units in most processors these days, working with integer numbers sometimes has its advantages. In these cases, it's useful to know that you can approximate π as 355/113 accurate to 7 decimal places. You can get this particular piece of information from Wikipedia, but what if your algorithm requires something like ln(2) or φ.

You can find these answers using the tree structure discovered by Messrs Stern and Brocot. The first few levels of this tree are shown in the image below:

Mathematical Background

If you are interested only in the programming details, you can skip this part and go directly to the applications section.

Rational numbers are those numbers that can be expressed as a fraction p/q. As there are many fractions that represent the same rational number, we will choose only the fraction where P and q don't have any common factors. For instance, 2/3 and 4/6 represent the same number, but we will choose only the first one. If we arrange the integer numbers on two axes, we can imagine rational numbers as the slope of the line from origin to the p/q point. The image below shows the representation of "2/3".

We need to define now the concept of mediant: the mediant of two rational numbers m/n and m'/n' is the rational number (m+m')/(n+n'). Below are represented the numbers 1/2 and 1/3 and their mediant 2/5.

Looking at those numbers as a vector addition, it can be seen that the mediant of two numbers lies somewhere between them. In math terms, if m/n < m'/n' then m/n < (m+m')/(n+n') < m'/n';

We can now start building the Stern-Brocot tree. We start with the numbers 0/1 and 1/0. 1/0 is not technically a fraction but it's useful to consider it as a representation for "infinity".

On each successive level of the tree, we add the mediants of the numbers on previous level. On level 2, we place the mediant of those numbers (0+1)/(1+0) = 1/1. On the next level, we place the medaints of each of those numbers: (0+1)/(1+1) = 1/2 and (1+1)/(1+0) = 2/1.

Next image shows how successive levels of the tree expand to cover the grid of rational numbers.

Grayed out points are those where p and q have common factors.

We've seen how rational numbers can be seen as lines going through one of the points of the grid of integers. Irrational numbers can then be seen as lines that never quite hit one of the grid points. Finding an approximation translates into finding a grid point that is close enough.

Application of Stern-Brocot trees for Rational Approximations

We've seen that our tree has on any level the mediants of the values on the previous level and that the value of a mediant is between the values of its ancestors. In such a tree, searching for a good enough approximation for the number N with tolerance ε can be done with the following algorithm:

1. Let current node T be the root of the tree.
2. If |N-T| < ε exit and the answer is T.
3. If N < T, replace T with the left child of T, otherwise replace T with the right child of T.
4. Proceed to step 2.

This algorithm, or better said the properties of Stern-Brocot tree, guarantees that the approximation p/q found is "best" in the sense that any other p'/q' approximation with the same precision will have p'>p and q'>q.

Usage

The sba program computes the best rational approximation with a given precision. The command line is:

sba <number> <precision>

Here is an example calculating the approximation for π with 7 decimal places:

Implementation

Nothing fancy here! Each tree node is represented by a sbnode structure that contains the p and q values, left and right children and a pointer to the parent node.

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;
  }
}

Constructor and destructor functions take care of initializing the node and cleaning up at end.

The program implements the algorithm described before. However, it doesn't create the tree beforehand. Instead, child nodes are created only when needed.

To show prime factorizations, we create a table of prime numbers using Eratosthenes sieve. This is not the fastest way to search for prime factors but for the small numbers we are dealing with, it is perfectly adequate.

Epilogue - The Gears of Time

After reading this story, you might wonder why a clockmaker was interested in these trees. If two gears, one with p teeth and the other with q teeth are connected, their rotation ratio (gear ratio) is p/q.

Gears can be combined like in the following example of a double gear:

In this case, if the gear ratio of the first couple (red/blue) is p1/q1 and the ratio of the second couple (yellow/green) is p2/q2, the combined gear ratio (red/green) is p1/q1 * p2/q2.

As an example of the problems clockmakers were having in the past, an 18th century clockmaker, Charles-Étienne-Louis Camus, presented the following problem: “To find the number of the teeth…of the wheels and pinions of a machine, which being moved by a pinion, placed on the hour wheel, shall cause a wheel to make a revolution in a mean year, supposed to consist of 365 days, 5 hours, 49 minutes”. If we work everything to minutes, the problem requires designing a gear train with the ratio 720/525949 (the hour wheel makes a complete rotation in 12 hours = 720 minutes and the mean tropical year has 525949 minutes according to Camus).

Building gears with 525949 teeth is hard to imagine and making compound gears in this case is hindered by the fact that 525949 is a prime number. Camus had to find a combination that was close to the required ratio but with numbers that factor nicely.

Our little sba utility provides an answer:

~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

This number can be decomposed as 4/25 * 7/69 * 7/83. It took Camus 20 pages of arduous work to get to the same result.

New Developments - Stern-Brocot tree and the Antikythera Mechanism

Discovered in 1901, the Antikythera mechanism is a Greek astronomical calculator dating to second or first century B.C. Since its discovery, it has challenged many researchers to find the internal workings of this technical marvel.

Recently, a team from UCL (University College of London) has published a complete reconstruction of the mechanism.

In order find the gearing factors the authors propose a process called Parmenides Proposition:

We have developed a new theory about how the Venus and Saturn periods were discovered and apply this to restore the missing planetary periods. A dialogue of Plato (fifth-fourth century BC) was named after the philosopher Parmenides of Elea (sixth-fifth century BC). This describes Parmenides Proposition:

In approximating θ, suppose rationals, p/q and r/s, satisfy p/q < θ < r/s. Then (p + r)/(q + s) is a new estimate between p/q and r/s:

If it is an underestimate, it is a better underestimate than p/q. If it is an overestimate, it is a better overestimate than r/s.

This description matches exactly the definition of the mediant and the search algorithm in Stern-Brocot tree.

Using the Parmenides Proposition, the authors describe the process of finding the gearing for the planet Venus. It must be a number close to 720/1151 with a convenient prime decomposition. Below is a fragment of the output from the sba utility:

The best such number is 289/462.

Continuing, the authors establish different gearing ratios for all planets and create a complete reconstruction of the mechanism. Below is an exploded view of the whole set of gears:

I found it fascinating to trace the structure of the Stern-Brocot tree from Plato's dialogues, to Greek astronomers, 18^th-century clock makers and contemporary floating-point approximations.

References

• [1] Graham, R. L., Knuth, D. E.,, Patashnik, O. (1994). Concrete Mathematics: A Foundation for Computer Science 2nd Edition. Reading: Addison-Wesley.
• [2] Wikipedia page for Stern-Brocot trees (https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree)
• [3] Wikipedia page for approximations of π (https://en.wikipedia.org/wiki/Approximations_of_%CF%80)
• [4] Wikipedia page for gear train (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

History

• 19^th January, 2021 - Initial version
• 22^st March, 2021 - Added code for prime factorization; added section on Antikythera mechanism