Brun's algorithm

Brun's algorithm is an extension of continued fractions to more than one number. It can be used to find EDOs that contain good approximations to given just ratios.

Algorithm

Say we would like to find EDOs which have good approximations of 5/4, 3/2, and 2/1. First, write the corresponding cents values, largest to smallest, in a vector [math]\displaystyle{ A_{0} }[/math],

\begin{displaymjx} A_{0} = \begin{bmatrix} 1200.00 \\ 701.96 \\ 386.31 \end{bmatrix} \end{displaymjx}

Now write [math]\displaystyle{ XA = A_{0} }[/math], where to begin with [math]\displaystyle{ X }[/math] is the identity matrix and [math]\displaystyle{ A = A_{0} }[/math],

\begin{displaymjx} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1200.00 \\ 701.96 \\ 386.31 \end{bmatrix} \\ = \begin{bmatrix} 1200.00 \\ 701.96 \\ 386.31 \end{bmatrix} \end{displaymjx}

Each iteration of the algorithm transforms [math]\displaystyle{ X }[/math] and [math]\displaystyle{ A }[/math] such that

• We preserve [math]\displaystyle{ X A = A_{0} }[/math]
• The elements of [math]\displaystyle{ A }[/math] get smaller
• The first column of [math]\displaystyle{ X }[/math] gives a better EDO approximation of [math]\displaystyle{ A_{0} }[/math]

To do one iteration of the algorithm, firstly subtract the second element from the first element in the cents vector [math]\displaystyle{ A }[/math] and add the first column to the second column in the matrix [math]\displaystyle{ X }[/math], giving

\begin{displaymjx} \begin{bmatrix} 1 & 0 + 1 & 0 \\ 0 & 1 + 0 & 0 \\ 0 & 0 + 0 & 1 \end{bmatrix} \begin{bmatrix} 1200.00 - 701.96 \\ 701.96 \\ 386.31 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 498.04 \\ 701.96 \\ 386.31 \end{bmatrix} \end{displaymjx}

Secondly, sort the cents vector [math]\displaystyle{ A }[/math] and rearrange the columns of [math]\displaystyle{ X }[/math] in the corresponding order, giving

\begin{displaymjx} \begin{bmatrix} 1 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 701.96 \\ 498.04 \\ 386.31 \end{bmatrix} \end{displaymjx}

Repeating this process gives:

[math]\displaystyle{ X }[/math]	[math]\displaystyle{ A }[/math]	Ratios
[math]\displaystyle{ \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}1200.00 \\ 701.96 \\ 386.31\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}2/1 \\ 3/2 \\ 5/4\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}1 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}701.96 \\ 498.04 \\ 386.31\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}3/2 \\ 4/3 \\ 5/4\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}2 & 0 & 1 \\ 1 & 0 & 1 \\ 0 & 1 & 0\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}498.04 \\ 386.31 \\ 203.91\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}4/3 \\ 5/4 \\ 9/8\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}2 & 1 & 2 \\ 1 & 1 & 1 \\ 1 & 0 & 0\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}386.31 \\ 203.91 \\ 111.73\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}5/4 \\ 9/8 \\ 16/15\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}3 & 2 & 2 \\ 2 & 1 & 1 \\ 1 & 1 & 0\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}203.91 \\ 182.40 \\ 111.73\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}9/8 \\ 10/9 \\ 16/15\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}5 & 2 & 3 \\ 3 & 1 & 2 \\ 2 & 0 & 1\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}182.40 \\ 111.73 \\ 21.51\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}10/9 \\ 16/15 \\ 81/80\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}7 & 5 & 3 \\ 4 & 3 & 2 \\ 2 & 2 & 1\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}111.73 \\ 70.67 \\ 21.51\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}16/15 \\ 25/24 \\ 81/80\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}12 & 7 & 3 \\ 7 & 4 & 2 \\ 4 & 2 & 1\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}70.67 \\ 41.06 \\ 21.51\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}25/24 \\ 128/125 \\ 81/80\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}19 & 12 & 3 \\ 11 & 7 & 2 \\ 6 & 4 & 1\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}41.06 \\ 29.61 \\ 21.51\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}128/125 \\ 3125/3072 \\ 81/80\end{bmatrix} }[/math]
[math]\displaystyle{ \begin{bmatrix}31 & 3 & 19 \\ 18 & 2 & 11 \\ 10 & 1 & 6\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}29.61 \\ 21.51 \\ 11.45\end{bmatrix} }[/math]	[math]\displaystyle{ \begin{bmatrix}3125/3072 \\ 81/80 \\ 393216/390625\end{bmatrix} }[/math]

Here each row shows one iteration of the algorithm. The Ratios column shows the ratios corresponding to the cents values in [math]\displaystyle{ A }[/math] (from dividing ratios each time we subtract cents values).

The first column of each matrix [math]\displaystyle{ X }[/math] gives an EDO approximation to [math]\displaystyle{ A_{0} }[/math]. For example, the matrix with first column [math]\displaystyle{ [31, 18, 10] }[/math] gives

\begin{displaymjx} \begin{aligned} \frac{1200}{31} \begin{bmatrix} 31 \\ 18 \\ 10 \end{bmatrix} &= \begin{bmatrix} 1200.00 \\ 696.77 \\ 387.10 \end{bmatrix} \\ &\approx \begin{bmatrix} 1200.00 \\ 701.96 \\ 386.31 \end{bmatrix} \end{aligned} \end{displaymjx}

Code

The following Python code implements Brun's algorithm.

import numpy as np

target_ratios = [5/4, 3/2, 2/1]

A = 1200 * np.log2(sorted(target_ratios, reverse=True))
X = np.identity(len(A), dtype=np.int64)

while X.max() <= 100:
    print(f"{X}\n{A}\n")
    A[0] -= A[1]
    X[:, 1] += X[:, 0]
    i = np.argsort(-A)
    A = A[i]
    X = X[:, i]

Other applications

Erv Wilson used the intervals generated by Brun's algorithm to construct scales.

Further reading

• Nicolas Chevallier, Best simultaneous Diophantine approximations and multidimensional continued fraction expansions. 2013.

External links

• bruntosaur - Online Brun's algorithm calculator