Click here to Skip to main content
15,867,453 members
Articles / General Programming / Algorithms

Stern-Brocot Trees and their Applications

Rate me:
Please Sign up or sign in to vote.
4.96/5 (25 votes)
1 May 2021MIT7 min read 21K   291   24   19
A data structure for finding best rational approximations
Rational approximations are useful when trying to avoid floating point rounding errors. This article describes a method for finding the 'best' such approximation.

Introduction

Sometimes, in some corners of math, you find things completely different from what you were looking for. This is what happened to a 19th 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:

tree

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".

Rational numbers

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.

Mediant construction

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.

Tree levels

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:

PI computation

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.

C++
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

Gears

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:

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:
Results for Venus

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: Exploded model of Cosmos gearing. From right to left: b1, mean Sun, Nodes, Mercury, Venus;true Sun and superior planets gearing; CP and shared gears; Ring Display; Dragon Hand; Moon position and phasemechanism

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

References

History

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

License

This article, along with any associated source code and files, is licensed under The MIT License


Written By
Canada Canada
Mircea is the embodiment of OOP: Old, Opinionated Programmer. With more years of experience than he likes to admit, he is always opened to new things, but too bruised to follow any passing fad.

Lately, he hangs around here, hoping that some of the things he learned can be useful to others.

Comments and Discussions

 
QuestionThis is a great article Pin
HoshiKata25-Jun-21 3:41
HoshiKata25-Jun-21 3:41 
AnswerRe: This is a great article Pin
Mircea Neacsu25-Jun-21 5:27
Mircea Neacsu25-Jun-21 5:27 
GeneralMy vote of 5 Pin
pierrecoach3-May-21 21:15
professionalpierrecoach3-May-21 21:15 
GeneralMy vote of 5 Pin
Espen Harlinn1-May-21 6:31
professionalEspen Harlinn1-May-21 6:31 
GeneralRe: My vote of 5 Pin
Mircea Neacsu1-May-21 7:08
Mircea Neacsu1-May-21 7:08 
GeneralMy vote of 5 Pin
DEK4665625-Mar-21 5:25
DEK4665625-Mar-21 5:25 
GeneralRe: My vote of 5 Pin
Mircea Neacsu25-Mar-21 5:58
Mircea Neacsu25-Mar-21 5:58 
GeneralMuch entertaining while going though the nice article Pin
Sanit Sharma23-Mar-21 18:24
professionalSanit Sharma23-Mar-21 18:24 
GeneralRe: Much entertaining while going though the nice article Pin
Mircea Neacsu23-Mar-21 23:58
Mircea Neacsu23-Mar-21 23:58 
GeneralMy vote of 5 Pin
Member 332128223-Mar-21 11:53
Member 332128223-Mar-21 11:53 
GeneralRe: My vote of 5 Pin
Mircea Neacsu23-Mar-21 12:19
Mircea Neacsu23-Mar-21 12:19 
GeneralMy Vote of 5+ Pin
AndyHo23-Mar-21 5:19
professionalAndyHo23-Mar-21 5:19 
GeneralRe: My Vote of 5+ Pin
Mircea Neacsu23-Mar-21 5:26
Mircea Neacsu23-Mar-21 5:26 
GeneralMy vote of 5 Pin
BillWoodruff23-Mar-21 4:43
professionalBillWoodruff23-Mar-21 4:43 
GeneralRe: My vote of 5 Pin
Mircea Neacsu23-Mar-21 5:09
Mircea Neacsu23-Mar-21 5:09 
GeneralMy vote of 5 Pin
Jan Heckman25-Jan-21 3:07
professionalJan Heckman25-Jan-21 3:07 
AnswerRe: My vote of 5 Pin
Mircea Neacsu25-Jan-21 3:54
Mircea Neacsu25-Jan-21 3:54 
GeneralMy vote of 5 Pin
Greg Utas19-Jan-21 7:16
professionalGreg Utas19-Jan-21 7:16 
GeneralRe: My vote of 5 Pin
Mircea Neacsu19-Jan-21 7:19
Mircea Neacsu19-Jan-21 7:19 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.