Next Home Previous

23 and me

Knuth studies continued fractions in the context of Euclid’s algorithm for greatest common divisor. We'll explore the proper fractions with denominator \( 23 \) just for the facts they reveal.

It's easy to work out the regular continued fractions by hand. Here is one example: \[ 7 / 23 = \frac {1} {\frac{23}{7}} = \frac {1} { 3 + \frac{2}{7} } = \frac {1} { 3 + \frac{1}{\frac{7}{2}} } = \frac {1} { 3 + \frac{1}{3 + \frac{1}{2}} } \ \ . \]

Here is the whole table of expansions: \[ \newcommand{\sslash}{\mathbin{/\mkern-6mu/}} \newcommand{\tslash}{\mathbin{/\mkern-5mu/}} \newcommand{\uslash}{\mathbin{/\mkern-4mu/}} \begin{array}{ll} 1 / 23 = \sslash 23 \sslash & 22 / 23 = \sslash 1, 22 \sslash \\ 2 / 23 = \sslash 11, 2 \sslash & 21 / 23 = \sslash 1, 10, 2 \sslash \\ 3 / 23 = \sslash 7, 1, 2 \sslash & 20 / 23 = \sslash 1, 6, 1, 2 \sslash \\ 4 / 23 = \sslash 5, 1, 3 \sslash & 19 / 23 = \sslash 1, 4, 1, 3 \sslash \\ 5 / 23 = \sslash 4, 1, 1, 2 \sslash & 18 / 23 = \sslash 1, 3, 1, 1, 2 \sslash \\ 6 / 23 = \sslash 3, 1, 5 \sslash & 17 / 23 = \sslash 1, 2, 1, 5 \sslash \\ 7 / 23 = \sslash 3, 3, 2 \sslash & 16 / 23 = \sslash 1, 2, 3, 2 \sslash \\ 8 / 23 = \sslash 2, 1, 7 \sslash & 15 / 23 = \sslash 1, 1, 1, 7 \sslash \\ 9 / 23 = \sslash 2, 1, 1, 4 \sslash & 14 / 23 = \sslash 1, 1, 1, 4 \sslash \\ 10 / 23 = \sslash 2, 3, 3 \sslash & 13 / 23 = \sslash 1, 1, 3, 3 \sslash \\ 11 / 23 = \sslash 2, 11 \sslash & 12 / 23 = \sslash 1, 1, 11 \sslash \end{array} \]

Several patterns emerge:

All of the expansions of fractions greater than \( 1 / 2 \) have the form \( \sslash 1, \ldots \sslash \).
Al of the expansions end with a term \( x_n > 1 \). That's obvious because \( \sslash x_1, \ldots, x_n, 1 \sslash = \sslash x_1, \ldots, x_n + 1 \sslash \)
The arrangement of the table sugggests the identity \( 1 - \sslash x_1, \ldots, x_n \sslash = \sslash 1, x_1 - 1, \ldots, x_n \sslash \).
The zero element that arises when the previous formula is applied to expansions for which \( x_1 = 1 \) illustrates a special case of the identity \( \sslash x_1, \ldots, x_k, 0, x_{k+1}, \ldots, x_n \) = \( \sslash x_1, \ldots, x_k + x_{k+1}, \ldots, x_n \sslash \).
The left-hand column shows that for expansions less than \( 1 / 2 \), whenever \( \sslash x_1, \ldots, x_t \sslash \) appears, so does \( \sslash x_t, \ldots, x_1 \sslash \).
The stronger statement of the previous item is that the elements in the left-hand column are generated by all the \( K_t(A_1, \ldots, A_t) \) where \( A_1 > 1 \), \(A_t > 1 \), and \( K_t(A_1, \ldots, A_t) \) is a divisor or, in this case, 23.

Now with a dash of floating point arithmetic

Rational arithmetic is the ideal medium for computing with values like \( 7 / 23 \) and their expansions into regular continued fractions. Every value is a ratio of two integers \( a / b \) and in the best case the integers may be arbitrarily large. But this site has a special focus on floating point arithmetic, so it's interesting to carry out the above computations in the 65-bit double type of the IEEE 754 standard.

When we recompute the regular continued fraction expansions of the proper fractions with denominator \( 23 \) things come out a bit different, as shown in the third column of this table: \[ \small{ \begin{array}{lclcl} 1/23 & = & \tslash 23 \tslash & = & \scriptsize{\uslash 23 \uslash} \\ 2/23 & = & \tslash 11, 2 \tslash & = & \scriptsize{\uslash 11, 2 \uslash} \\ 3/23 & = & \tslash 7, 1, 2 \tslash & \approx & \scriptsize{\uslash 7, 1, 2, 375299968947541, 3, 4, 1, 225179981368524, \ldots \uslash} \\ 4/23 & = & \tslash 5, 1, 3 \tslash & \approx & \scriptsize{\uslash 5, 1, 3, 1125899906842624 \uslash} \\ 5/23 & = & \tslash 4, 1, 1, 2 \tslash & \approx & \scriptsize{\uslash 4, 1, 1, 1, 1, 75059993789508, 3, 1, \ldots \uslash} \\ 6/23 & = & \tslash 3, 1, 5 \tslash & \approx & \scriptsize{\uslash 3, 1, 5, 140737488355328 \uslash} \\ 7/23 & = & \tslash 3, 3, 2 \tslash & \approx & \scriptsize{\uslash 3, 3, 1, 1, 140737488355328 \uslash} \\ 8/23 & = & \tslash 2, 1, 7 \tslash & \approx & \scriptsize{\uslash 2, 1, 7, 281474976710656 \uslash} \\ 9/23 & = & \tslash 2, 1, 1, 4 \tslash & \approx & \scriptsize{\uslash 2, 1, 1, 4, 56294995342131, 4, 1, 12, \ldots \uslash} \\ 10/23 & = & \tslash 2, 3, 3 \tslash & \approx & \scriptsize{\uslash 2, 3, 3, 37529996894754, 7, 1, 1, 8, \ldots \uslash} \\ 11/23 & = & \tslash 2, 11 \tslash & \approx & \scriptsize{\uslash 2, 11, 93824992236885, 3, 21, 13403570319555, 21, 3, \ldots \uslash} \\ 12/23 & = & \tslash 1, 1, 11 \tslash & \approx & \scriptsize{\uslash 1, 1, 11, 93824992236885, 3, 21, 13403570319555, 21, \ldots \uslash} \\ 13/23 & = & \tslash 1, 1, 3, 3 \tslash & \approx & \scriptsize{\uslash 1, 1, 3, 2, 1, 57738456761160, 4, 1, \ldots \uslash} \\ 14/23 & = & \tslash 1, 1, 1, 4 \tslash & \approx & \scriptsize{\uslash 1, 1, 1, 1, 3, 1, 46912496118442, 1, \ldots \uslash} \\ 15/23 & = & \tslash 1, 1, 1, 7 \tslash & \approx & \scriptsize{\uslash 1, 1, 1, 7, 40210710958665, 7, 9, 31274997412295, \ldots \uslash} \\ 16/23 & = & \tslash 1, 2, 3, 2 \tslash & \approx & \scriptsize{\uslash 1, 2, 3, 1, 1, 140737488355328 \uslash} \\ 17/23 & = & \tslash 1, 2, 1, 5 \tslash & \approx & \scriptsize{\uslash 1, 2, 1, 4, 1, 37529996894754, 7, 1, \ldots \uslash} \\ 18/23 & = & \tslash 1, 3, 1, 1, 2 \tslash & \approx & \scriptsize{\uslash 1, 3, 1, 1, 1, 1, 28147497671065, 1, \ldots \uslash} \\ 19/23 & = & \tslash 1, 4, 1, 3 \tslash & \approx & \scriptsize{\uslash 1, 4, 1, 3, 1125899906842624 \uslash} \\ 20/23 & = & \tslash 1, 6, 1, 2 \tslash & \approx & \scriptsize{\uslash 1, 6, 1, 1, 1, 18764998447377, 15, 17, \ldots \uslash} \\ 21/23 & = & \tslash 1, 10, 2 \tslash & \approx & \scriptsize{\uslash 1, 10, 2, 20105355479332, 1, 1, 3, 18, \ldots \uslash} \\ 22/23 & = & \tslash 1, 22 \sslash & \approx & \scriptsize{\uslash 1, 22, 46912496118442, 1, 1, 1, 42, 6701785159777, \ldots \uslash} \end{array} } \]

What's going on? The Python code reflects the computation of \( A_{n+1} \) and \( X_{n+1} \), when \( X_n \neq 0 \):

        new_a = math.floor(1.0 / x)
        new_x = (1.0 / x) - new_a

The computation is doomed from the start. \( 7 / 23 \) rounds up and the difference ripples through the calculation. Let's revisit the earlier calculation, in the presence of rounding. The Greek letters represent tiny positive values. \[ \begin{array}{ll} \frac{7}{23} + \alpha & = \frac {1} { 3 + (\frac{2}{7} - \beta) } = \frac {1} { 3 + \frac{1}{3 + (\frac{1}{2} + \gamma) }} = \frac {1} { 3 + \frac{1}{3 + \frac{1}{1 + (\frac{1}{1} - \delta)}}} \\ & = \frac {1} { 3 + \frac{1}{3 + \frac{1}{1 + \frac{1}{1 + \epsilon}}}} = \frac {1} {3 + \frac{1}{3 + \frac{1}{1 + \frac{1}{1 + \frac{1}{HUGE}}}}} \ \ . \end{array} \] These intermediate values printed by the Python program are annotated with references to the formulas just above.

    new_a = 3  new_x = 0.2857142857142856  2/7 - beta
    new_a = 3  new_x = 0.5000000000000018  1/2 + gamma
    new_a = 1  new_x = 0.9999999999999929  1 - delta
    new_a = 1  new_x = 7.105427357601002e-15  epsilon = 32 ulps of 1
    new_a = 140737488355328  new_x = 0.0  1/epsilon is exact

We can get lucky, as with \( 1 / 23 \), and have the rounding errors cancel and leave us with an exact result. In most cases, the error will ripple down until what ought to be the last step of \( \frac{1}{M} \) turns out instead to be \( \frac{1}{M \pm \delta} \). When it's \( M - \delta \), as happened just above, the continued fraction expansion becomes \( \sslash x_1, \ldots, x_n - 1, 1, HUGE, \ldots \sslash \). In the case of \( M + \delta \) the expansion is correct up to the huge value: \( \sslash x_1, \ldots, x_n, HUGE, \ldots \sslash \). The case \( 8/23 \) exhibits that behavior.

In the step after the \( \delta \) step, we compute \( 1 / \epsilon \). If \( \epsilon \) is \( 2^k \) units in the last place of \( 1 \), then its reciprocal is exact and the computation terminates because new_x == 0. An ulp of \( 1 \) is \( 2^{-52} \). In the case of \( 7/23 \), \( \epsilon \) is 32 ulps of \( 1 \). The huge integer assigned to new_a is exactly \( 2^{47} \). In the case of \( 9/23 \), however, \( \epsilon \) is 80 ulps of \( 1 \), so the computation wanders on.

When HUGE becomes \( \infty \)

On previous page, we computed a bound on the error in a continued fraction expansion. The appearance of the huge values in the computed expansions puts the error down in the ulps of the intendd result.

What happens in the limit, when the huge value jumps to \( \infty \)? Here is the case of \( 7 / 23 \): \[ \small{ \begin{array}{lclcl} 7/23 & = & \tslash 3, 3, 2 \tslash & = & \scriptsize{\uslash 3, 3, 1, 1, \infty \uslash} \end{array} } \ \ . \] The continued fraction expands to \[ \begin{array}{l} \frac{7}{23} & = \frac {1} { 3 + \frac{1}{3 + \frac{1}{ 1 + \frac{1}{1 + \frac{1}{ \infty }}}}} = \frac {1} { 3 + \frac{1}{3 + \frac{1}{ 1 + \frac{1}{1}}}} = \frac {1} { 3 + \frac{1}{3 + \frac{1}{ 2}}} \end{array} \] from which we see that \[ \normalsize{ 7/23 = \sslash 3, 3, 1, 1, \infty \sslash = \sslash 3, 3, 1, 1\sslash = \sslash 3, 3, 2 \sslash } \ \ . \] Although \( \infty \) doesn't actually arise in the calculations on this page, IEEE 754 arithmetic will handle the Python value math.inf when it does appear in a continued fraction evaluation.

The composite case

It's interesting to look at the table of regular continued fraction expansions for a highly composite divisor like \( 24 \). \[ \begin{array}{ll} 1 / 24 = \sslash 24 \sslash & 23 / 24 = \sslash 1, 23 \sslash \\ 2 / 24 = \sslash 12 \sslash & 22 / 24 = \sslash 1, 11 \sslash \\ 3 / 24 = \sslash 8 \sslash & 21 / 24 = \sslash 1, 7 \sslash \\ 4 / 24 = \sslash 6 \sslash & 20 / 24 = \sslash 1, 5 \sslash \\ 5 / 24 = \sslash 4, 1, 4 \sslash & 19 / 24 = \sslash 1, 3, 1, 4 \sslash \\ 6 / 24 = \sslash 4 \sslash & 18 / 24 = \sslash 1, 3 \sslash \\ 7 / 24 = \sslash 3, 2, 3 \sslash & 17 / 24 = \sslash 1, 2, 2, 3 \sslash \\ 8 / 24 = \sslash 3 \sslash & 16 / 24 = \sslash 1, 2 \sslash \\ 9 / 24 = \sslash 2, 1, 2 \sslash & 15 / 24 = \sslash 1, 1, 1, 2 \sslash \\ 10 / 24 = \sslash 2, 2, 2 \sslash & 14 / 24 = \sslash 1, 1, 2, 2 \sslash \\ 11 / 24 = \sslash 2, 5, 2 \sslash & 13 / 24 = \sslash 1, 1, 5, 2 \sslash \\ 12 / 24 = \sslash 2 \sslash & \; \end{array} \] Looking at the various \( \sslash x_1, \ldots, x_t \sslash \) in the form \( K_{t-1}(x_2, \ldots, x_t) / K_t(x_1, \ldots, x_t) \), we see that the denominators generate all of the factors of \( 24 \) greater than one.

And so ends an exploration of regular continued fraction expansions and a look at the effects of rounding on what deserve to be exact calculations.

Next Home Previous