Algebraic Numbers

Any algebraic number has a fractal fraction expansion. It appears that these are suboptimal for computation since the number of numerators and denominators doubles each iteration.

Linear recurrence relations can be constructed for cube roots and more that allow continued-fraction-like expansions. [TODO: See if they can efficiently produce continued fraction expansions.]

With floating-point numbers, using Newton’s method gives us \(n\) bits of any algebraic number in the same time it takes to multiply \(n\)-bit numbers. Does Newton’s method also help calculations with continued fractions?

Newton’s Method

Here’s a sure-fire, obvious method for calculating the continued fraction expansion of a positive real algebraic number \(\alpha\) with minimal polynomial \(c^n x^n + ... + c_0\) over the rationals.

We maintain convergents for \(\alpha, ..., \alpha^{n-1}\). Using Newton’s method, find \(a_1 = \lfloor \alpha \rfloor\), and also our first convergent for \(\alpha\). Use this to compute convergents for the other powers of \(\alpha\). Then \(a_2 = \lfloor \frac {1}{\alpha - a_1}\rfloor\), which we rewrite as \(\lfloor d_n \alpha^{n-1} + ... + d_0 \rfloor\) via abstract algebra (extended GCD algorithm for polynomials).

Substitute our convergents for the powers of \(\alpha\) and provisionally assign \(a_2\) to the floor of this value. Correct it by trying \(a_1 + 1/a_2\) in the minimal polynomial. I haven’t worked out the details: perhaps we should use Newton’s method again; in some cases, we can keep increasing \(a_2\) until the sign changes. Update all convergents, and repeat.

Almost certainly better is using Newton’s method to generate a sequence of rationals whose closeness to \(\alpha\) is known, then compute the continued fraction expansion from this.

Square Roots

If we specialize to polynomials of degree 2, we have more concrete results.

Consider the function:

\[ f(y) = \frac{a y + b}{c y - a} \]

and we want to find \(\alpha\) such that \(\alpha = f(\alpha)\). Then starting with some guess \(\alpha_0\) and successively computing

\[ \alpha_{i+1} = \frac{\alpha_i + f(\alpha_i)}{2} \]

is equivalent to Newton’s method on the quadratic:

\[ c y^2 - 2 a y - b = 0 \]

Let us follow Gosper and see how these statements apply to finding \(\sqrt{\frac{17}{10}}\). This is equivalent to finding a positive fixed point of

\[ f(y) = \frac{17}{10 y} \]

We initialize a table for computing the continued fraction expansion:

17/0

0/10

but we run into a bootstrapping problem: we cannot proceed until we have part of the continued fraction sequence.

To remedy this, we find \(a_1\) via Newton’s method: guess arbitrarily (e.g. 3) and put the guess into the function (i.e. 17/30). If their integer parts differ (they do) take the integer part of their average (i.e. 1) and repeat (we get 17/10, whose integer part agrees, so we have found \(a_1\)).

This starts off our table, as well as our output.

\[ [1;...] \]

1

17/0

0/10

17/10

Since we know the first output is 1, we subtract and invert preemptively. In fact, we could have done this even before applying the recurrence relation (if we choose this route, we subtract 1 from 17/0 by computing (17 - denominator x 1)/ denominator, i.e. \((17-0\times 1)/0 = 17/0\).)

\[ [1;...] \]

1

17/0

0/10/-10

17/10/7

Again we need more of the expansion to continue, and again we use Newton’s method: we want the integer part of the positive solution of

\[ y = \frac{10 y + 10}{7 y - 10} \]

We must guess at least \(10/7\) for the positive root. Say we pick \(2\). Then the right hand side floors to give \(4\), so we take the average, \(3\), and substituting this in the right hand side shows this is correct guess for the integer part of the root. We can now continue, applying the recurrence relation as well as subtracting 3 and inverting:

\[ [1;3, ...] \]

1

3

17/0

0/10/-10

17/10/7/-11

40/11/7

We repeat this process: 3, 2 are the next two correct guesses:

\[ [1;3, 3, 2, ...] \]

1

3

3

2

17/0

0/10/-10

17/10/7/-11

40/11/7/-10

40/10/10/-10

27/10/7

We’re at the same state we were just after our first output, so the expansion is \([1;3, 3, 2, 3, 3, 2, ...]\).

Did you notice…?

  • \(f(y)\) is its own inverse.

  • The homographic function represented by the last two convergents is always its own inverse.

  • Every quadratic considered has roots of opposing signs.

  • The determinant of each Mobius function is \(D = -a^2-b c\); we could have avoided most of the guesswork by finding \(d = \lfloor\sqrt{-D}\rfloor\) and observing

\[ \lfloor y \rfloor = \left\lfloor \frac{a+\sqrt{-D}}{c} \right\rfloor = \left\lfloor \frac{a+d}{c} \right\rfloor\]

All but the last of these facts hold when we consider:

Square Roots of Continued Fractions

Using bihomographic functions we can take square roots of continued fractions, not just rationals. As before these results apply to quadratic equations, not just square roots.

We illustrate with an example: we compute \(\coth 1/2\) using

\[ \frac{(\coth t)^2 + 1}{2 \coth t} = \coth 2t \]

From matrix calculations we know \(\coth 1 = [1;3, 5, ...]\).

The above quadratic has two positive solutions. In general, our initial guesses for Newton’s method will decide where we end up. Once we’ve chosen a root to pursue, the follow-up quadratics always have exactly one positive root greater than 1, because continued fraction expansions are unique.

We wish to find a solution \(\alpha\) of

\[ y = \frac{x y - 1}{y - x} \]

where \(x = \coth 1\). We initialize the three-dimensional table as follows:

1

3

5

7

-1/0

0/-1

0/1

1/0

Observe the left two convergents correspond to a self-inverse homographic function, as do the right two convergents.

We have a similar situation as above. We know the expansion of \(x\), but without that of \(\alpha\) we cannot proceed far. Yet the expansion of \(\alpha\) is exactly what we are seeking.

Thus we open with Newton’s method. Or at least, we attempt to. The right convergents instruct us to find a the integer part of the fixed point of

\[ y= - y \]

which is 0, whilst the left convergents want the same of

\[ y= -1/y \]

which is certainly not 0. We must move along horizontally to obtain consistent equations. In fact, we must do so for a few steps:

1

3

5

7

-1/0

0/-1

-1/-1

-3/-4

-16/-21

0/1

1/0

1/1

4/3

21/16

Now the integer parts of the fixed points of the two functions agree, namely at \(y = 2\). Thus we start the \(y\)-axis with \(2\), output \(2\) and subtract 2 and invert before continuing:

\[ [2;...] \]

1

3

5

7

-1/0

0/-1

-1/-1

-3/-4

-16/-2

0/1

1/0

1/1

4/3/-2

21/16/-11

2

5/2/1

26/11/4

The rightmost two convergents lead to a function with a fixed point between 6 and 7, which disagrees with the function implied by the other two convergents. We must move horizontally again.

\[ [2;...] \]

1

3

5

7

-1/0

0/-1

-1/-1

-3/-4

-16/-2

0/1

1/0

1/1

4/3/-2

21/16/-11

115/-79

2

5/2/1

26/11/4

79/29

Now the two functions both have fixed points between 6 and 7, so we output 6, before subtracting 6 and inverting. And so on.

It is similar to square rooting a rational: we use Newton’s method to find the next term, except that there are two functions at a time to consider, and they must agree on the integer part of the fixed point; if they differ, we must feed in more terms of the input continued fraction to refine the functions.

[TODO: Check this.] Increasing the number of dimensions gives a similar algorithm for quadratics with more than one continued fraction coefficient.

Accuracy of integers in Newton’s method

We were fortunate in the above example that we could find the integer parts of fixed points of functions by evaluating them exclusively at integers. However, this fails in general. Consider:

\[ f(y) = \frac{2y + 3}{y - 2} \]

We have \(f(4) = 11/2 = 5 + 1/2\) and \(f(5) = 13/3 = 4 + 1/3\). It gets worse: we could have \(\lfloor f(8)\rfloor = 6, \lfloor f(6) \rfloor = 7, \lfloor f(7) \rfloor = 8\). In general we need more than integer precision to find the integer part of the answer.

An alternative is to look for sign changes in the corresponding quadratic \(y^2 - 4y - 3\) via a binary search.


Ben Lynn blynn@cs.stanford.edu 💡