My question concerns the numerical accuracy of a continued fraction expansion. A typical algorithm for computing a continued fraction can be written in Python as :
x0 = sqrt(2)
N = 40
a = [0]*N
u = [0]*N
x = x0
for k in range(N):
a[k] = int(x//1)
u[k] = x % 1 # Often replaced with x - a[k] ???
x = 1/u[k]
print(a)
For $x=\sqrt{2}$, this produces
[1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 3, 3, 1, 3, 1, 1, 2, 1809, 1, 2, 5, 2, 2, 1, 2, 1]
Or for $x = \pi$, I get :
[3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3, 3, 23, 1, 1, 7, 4, 35, 1, 1, 1, 2, 3, 3, 3, 3, 1, 1, 14, 6, 4, 5, 1, 7, 1, 5, 1]
One observation is that round-off error has clearly crept into the calculation. The expansion for $\sqrt{2}$ should be repeating 2s, and the approximation for $\pi$ is accurate up to the first appearance of 14 (12 terms after the leading $a_0=3$ term).
What I also found is that while round off error does creep in, the number of correct terms in the expansion is exactly what is needed to get convergents that agree with $x$ up to machine precision. So for example, convergent 20 of $\sqrt{2}$ is given by
54608393 / 38613965
This convergent approximates $\sqrt{2}$ to $4.44 \times 10^{-16}$. Convergent 21 agrees with $\sqrt{2}$ exactly (in finite precision arithmetic). Interestingly, the expansion has exactly 20 correct terms (i.e. 2s) after the leading $a_0 = 1$.
Similarly, for $\pi$, the convergent 12 is given by
80143857 / 25510582
which agrees with $\pi$ to within an error of $4.44 \times 10^{-16}$. Convergent 13 agrees with $\pi$ exactly (in finite precision arithmetic).
In both cases, the numerators for the convergents agree with the OEIS. See A001333 and A002485.
The above observations led me to these questions :
-
Is there a standard stopping criteria that can be used to determine when the expansion has converged to a desired tolerance?
-
Is it always the case that one will have enough correct terms in the expansion to approximate the desired number to machine precision?
-
Is it possible to detect whether a continued fraction is periodic? Or, if one knows it will be periodic (i.e. $x$ is a root of a quadratic with integer coefficients), is it possible to get the periodic sequence exactly?
It has also occurred to me that nobody would think of computing continued fraction expansions using finite precision arithmetic! I would like to do this problem as an exercise in a beginning computational math course, and was hoping not to go into variable precision arithmetic (an area which is really outside of my wheel house).
Best Answer
Here is a method for $\sqrt n$ that is likely to be what Fermat used in hand computations.
$$ \sqrt { 5} = 2 + \frac{ \sqrt {5} - 2 }{ 1 } $$ $$ \frac{ 1 }{ \sqrt {5} - 2 } = \frac{ \sqrt {5} + 2 }{1 } = 4 + \frac{ \sqrt {5} - 2 }{1 } $$
Simple continued fraction tableau:
$$ \begin{array}{cccccccc} & & 2 & & 4 & & 4 & \\ \\ \frac{ 0 }{ 1 } & \frac{ 1 }{ 0 } & & \frac{ 2 }{ 1 } & & \frac{ 9 }{ 4 } \\ \\ & 1 & & -1 & & 1 \end{array} $$
$$ \begin{array}{cccc} \frac{ 1 }{ 0 } & 1^2 - 5 \cdot 0^2 = 1 & \mbox{digit} & 2 \\ \frac{ 2 }{ 1 } & 2^2 - 5 \cdot 1^2 = -1 & \mbox{digit} & 4 \\ \frac{ 9 }{ 4 } & 9^2 - 5 \cdot 4^2 = 1 & \mbox{digit} & 4 \\ \end{array} $$
$$\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$
$$ \sqrt { 13} = 3 + \frac{ \sqrt {13} - 3 }{ 1 } $$ $$ \frac{ 1 }{ \sqrt {13} - 3 } = \frac{ \sqrt {13} + 3 }{4 } = 1 + \frac{ \sqrt {13} - 1 }{4 } $$ $$ \frac{ 4 }{ \sqrt {13} - 1 } = \frac{ \sqrt {13} + 1 }{3 } = 1 + \frac{ \sqrt {13} - 2 }{3 } $$ $$ \frac{ 3 }{ \sqrt {13} - 2 } = \frac{ \sqrt {13} + 2 }{3 } = 1 + \frac{ \sqrt {13} - 1 }{3 } $$ $$ \frac{ 3 }{ \sqrt {13} - 1 } = \frac{ \sqrt {13} + 1 }{4 } = 1 + \frac{ \sqrt {13} - 3 }{4 } $$ $$ \frac{ 4 }{ \sqrt {13} - 3 } = \frac{ \sqrt {13} + 3 }{1 } = 6 + \frac{ \sqrt {13} - 3 }{1 } $$
Simple continued fraction tableau:
$$ \begin{array}{cccccccccccccccccccccccc} & & 3 & & 1 & & 1 & & 1 & & 1 & & 6 & & 1 & & 1 & & 1 & & 1 & & 6 & \\ \\ \frac{ 0 }{ 1 } & \frac{ 1 }{ 0 } & & \frac{ 3 }{ 1 } & & \frac{ 4 }{ 1 } & & \frac{ 7 }{ 2 } & & \frac{ 11 }{ 3 } & & \frac{ 18 }{ 5 } & & \frac{ 119 }{ 33 } & & \frac{ 137 }{ 38 } & & \frac{ 256 }{ 71 } & & \frac{ 393 }{ 109 } & & \frac{ 649 }{ 180 } \\ \\ & 1 & & -4 & & 3 & & -3 & & 4 & & -1 & & 4 & & -3 & & 3 & & -4 & & 1 \end{array} $$
$$ \begin{array}{cccc} \frac{ 1 }{ 0 } & 1^2 - 13 \cdot 0^2 = 1 & \mbox{digit} & 3 \\ \frac{ 3 }{ 1 } & 3^2 - 13 \cdot 1^2 = -4 & \mbox{digit} & 1 \\ \frac{ 4 }{ 1 } & 4^2 - 13 \cdot 1^2 = 3 & \mbox{digit} & 1 \\ \frac{ 7 }{ 2 } & 7^2 - 13 \cdot 2^2 = -3 & \mbox{digit} & 1 \\ \frac{ 11 }{ 3 } & 11^2 - 13 \cdot 3^2 = 4 & \mbox{digit} & 1 \\ \frac{ 18 }{ 5 } & 18^2 - 13 \cdot 5^2 = -1 & \mbox{digit} & 6 \\ \frac{ 119 }{ 33 } & 119^2 - 13 \cdot 33^2 = 4 & \mbox{digit} & 1 \\ \frac{ 137 }{ 38 } & 137^2 - 13 \cdot 38^2 = -3 & \mbox{digit} & 1 \\ \frac{ 256 }{ 71 } & 256^2 - 13 \cdot 71^2 = 3 & \mbox{digit} & 1 \\ \frac{ 393 }{ 109 } & 393^2 - 13 \cdot 109^2 = -4 & \mbox{digit} & 1 \\ \frac{ 649 }{ 180 } & 649^2 - 13 \cdot 180^2 = 1 & \mbox{digit} & 6 \\ \end{array} $$
$$\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc\bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$