/[dtapublic]/to_be_filed/uculib01/doc/manual/c_rla1/c_rla1.tex
ViewVC logotype

Contents of /to_be_filed/uculib01/doc/manual/c_rla1/c_rla1.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 30 - (show annotations) (download) (as text)
Sat Oct 8 07:22:17 2016 UTC (7 years, 5 months ago) by dashley
File MIME type: application/x-tex
File size: 53824 byte(s)
Initial commit.
1 %$Header: /home/dashley/cvsrep/uculib01/uculib01/doc/manual/c_rla1/c_rla1.tex,v 1.11 2010/01/28 21:18:33 dashley Exp $
2
3 \chapter{Rational Linear Approximation}
4
5 \label{crla1}
6
7 \section{Introduction}
8 %Section Tag: INT
9 \label{crla1:sint0}
10
11 Inexpensive microcontrollers CPUs often possess integer multiplication
12 and/or division machine instructions.
13 These instructions are always faster and more compact than corresponding
14 integer multiplication and division subroutines written without the
15 machine instructions.
16
17 In the case a processor possesses both integer multiplication and integer division
18 machine instructions, it is economical to approximate an integer-valued
19 function $f(x)$ of a non-negative real constant $r_I$ and an
20 non-negative integer $x$ defined
21 by
22
23 \begin{equation}
24 \label{eq:crla1:sint0:01}
25 f(x) = \lfloor r_I x \rfloor
26 \end{equation}
27
28 \noindent{}by choosing a non-negative integer $h$ and
29 a positive integer $k$ so as to place
30 the rational number $r_A$ close to $r_I$\@.
31 (\ref{eq:crla1:sint0:01}) can then be
32 approximated by
33
34 \begin{equation}
35 \label{eq:crla1:sint0:02}
36 g(x) = \lfloor r_A x \rfloor
37 = \left\lfloor \frac{hx}{k} \right\rfloor.
38 \end{equation}
39
40 Note that (\ref{eq:crla1:sint0:02}) can be evaluated directly
41 by a processor with both integer multiplication and integer division
42 instructions, and that the
43 \emph{floor($\cdot$)} function (as only non-negative $h$, $k$, $x$ are
44 considered) approximates the behavior of most processor integer division
45 instructions, which provide a remainder separately from the quotient.
46
47 In order to control (or ``place'')
48 the approximation error,
49 it is also possible and common to add a non-negative integer $z$ to the numerator
50 of the approximation to yield
51
52 \begin{equation}
53 \label{eq:crla1:sint0:03}
54 g(x) =
55 \left\lfloor \frac{hx + z}{k} \right\rfloor
56 =
57 \left\lfloor r_Ax + \frac{z}{k} \right\rfloor.
58 \end{equation}
59
60 In the event that the processor has an integer multiplication instruction but
61 no integer division instruction, it is common to choose
62
63 \begin{equation}
64 \label{eq:crla1:sint0:04}
65 r_A = \frac{h}{2^q}
66 \end{equation}
67
68 \noindent{}so that the division can be implemented by a bit shifting
69 operation.
70
71 In the event that the processor has an integer division instruction but
72 no integer multiplication instruction (rare), $r_A$ can be chosen as
73
74 \begin{equation}
75 \label{eq:crla1:sint0:05}
76 r_A = \frac{2^p}{k}
77 \end{equation}
78
79 \noindent{}so that the multiplication can be implemented by a bit-shifting
80 operation.
81
82 This chapter addresses the implementation of approximations of
83 the form of (\ref{eq:crla1:sint0:02}). Choosing a rational number
84 $r_A = h/k$ as close as possible to an arbitrary $r_I$ requires
85 results from number theory (a branch of mathematics); and these
86 results are presented here.
87
88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91
92 \section{Nomenclature}
93 %Section Tag: NOM
94 \label{crla1:snom0}
95
96 We assume that the domain of the approximation is
97
98 \begin{equation}
99 \label{eq:crla1:snom0:01}
100 x \in [0, x_{MAX}] .
101 \end{equation}
102
103 \noindent{}$x_{MAX}$ may in some applications be smaller than the maximum value that
104 can be represented in
105 the data operand to the processor's integer multiplication instruction; but
106 the most usual case is that $x_{MAX} = 2^w-1$, where $w$ is the number of bits
107 in the data operand accepted by the processor's multiplication instruction.
108
109 We also assume that $h$ (of Equation \ref{eq:crla1:sint0:02} and related
110 equations) is constrained by
111
112 \begin{equation}
113 \label{eq:crla1:snom0:02}
114 h \in [0, h_{MAX}] .
115 \end{equation}
116
117 \noindent{}The most typical constraint on $h$ is due to the size of operand
118 that an integer multiplication instruction will accept; but there may be other reasons
119 for constraining $h$ (for example, to allow an arbitrary choice of $z$
120 of Equation \ref{eq:crla1:sint0:03} without
121 requiring test and branch logic).
122
123 Finally, we assume that $k$ is also constrained
124
125 \begin{equation}
126 \label{eq:crla1:snom0:03}
127 k \in [1, k_{MAX}] ,
128 \end{equation}
129
130 \noindent{}again with the most typical reason being the size of divisor that
131 a machine integer division instruction will accept.
132
133 We don't constrain $z$ except to assume that $z \geq 0$
134 (details of possible overflow of
135 $hx+z$ are left to the programmer).
136
137
138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141
142 \section[Choosing $r_A = h/k \approx r_I$]
143 {Choosing \mbox{\boldmath $r_A = h/k \approx r_I$}}
144
145 %Section Tag: lcr0
146 \label{crla1:slcr0}
147
148 This section presents the details of how to choose $r_A = h/k$ so that
149 it is as close as possible to $r_I$. Proofs are generally omitted, as they
150 would rely on material not presented for reasons of brevity.
151
152
153 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
156
157 \subsection{The Farey Series}
158
159 %Subsection Tag: fry0
160 \label{crla1:slcr0:sfry0}
161
162 The \emph{Farey}\footnote{Named after geologist John Farey. whose letter about
163 these series was published in the \emph{Philosophical Magazine} in 1816.} \emph{series
164 of order $N$},\index{Farey series} denoted $F_{N}$,\index{F@$F_N$}
165 is the ordered set of all irreducible
166 rational numbers $h/k$ in the interval
167 [0,1]
168 with denominator $k\leq N$.
169 As examples, the Farey series of
170 orders 1 through 7, $F_1$ through $F_7$, are shown
171 in (\ref{eq:crla1:slcr0:sfry0:eq0001a}) through (\ref{eq:crla1:slcr0:sfry0:eq0001g}).
172
173 \begin{equation}
174 \label{eq:crla1:slcr0:sfry0:eq0001a}
175 F_1 = \left\{ {\frac{0}{1},\frac{1}{1}} \right\}
176 \end{equation}
177
178 \begin{equation}
179 \label{eq:crla1:slcr0:sfry0:eq0001b}
180 F_2 = \left\{ {\frac{0}{1},\frac{1}{2},\frac{1}{1}} \right\}
181 \end{equation}
182
183 \begin{equation}
184 \label{eq:crla1:slcr0:sfry0:eq0001c}
185 F_3 = \left\{ {\frac{0}{1},\frac{1}{3},\frac{1}{2},
186 \frac{2}{3},\frac{1}{1}} \right\}
187 \end{equation}
188
189 \begin{equation}
190 \label{eq:crla1:slcr0:sfry0:eq0001d}
191 F_4 = \left\{ {\frac{0}{1},\frac{1}{4},
192 \frac{1}{3},\frac{1}{2},
193 \frac{2}{3},\frac{3}{4},
194 \frac{1}{1}} \right\}
195 \end{equation}
196
197 \begin{equation}
198 \label{eq:crla1:slcr0:sfry0:eq0001e}
199 F_5 = \left\{ {\frac{0}{1},\frac{1}{5},\frac{1}{4},
200 \frac{1}{3},\frac{2}{5},\frac{1}{2},
201 \frac{3}{5},\frac{2}{3},\frac{3}{4},
202 \frac{4}{5},\frac{1}{1}} \right\}
203 \end{equation}
204
205 \begin{equation}
206 \label{eq:crla1:slcr0:sfry0:eq0001f}
207 F_6 = \left\{ {\frac{0}{1},\frac{1}{6},\frac{1}{5},
208 \frac{1}{4},
209 \frac{1}{3},\frac{2}{5},\frac{1}{2},
210 \frac{3}{5},\frac{2}{3},
211 \frac{3}{4},
212 \frac{4}{5},
213 \frac{5}{6},\frac{1}{1}} \right\}
214 \end{equation}
215
216
217 \begin{equation}
218 \label{eq:crla1:slcr0:sfry0:eq0001g}
219 F_7 = \left\{ {\frac{0}{1},\frac{1}{7},\frac{1}{6},\frac{1}{5},
220 \frac{1}{4},\frac{2}{7},
221 \frac{1}{3},\frac{2}{5},\frac{3}{7},\frac{1}{2},
222 \frac{4}{7},\frac{3}{5},\frac{2}{3},
223 \frac{5}{7},\frac{3}{4},
224 \frac{4}{5},
225 \frac{5}{6},\frac{6}{7},\frac{1}{1} } \right\}
226 \end{equation}
227
228 The distribution of Farey rational numbers in
229 [0,1] is repeated
230 in any
231 $[i,i+1]$, $i\in \vworkintset$; so that the distribution of
232 Farey rationals in [0,1] supplies complete
233 information about the distribution in
234 all of $\vworkrealset$. We
235 occasionally abuse the proper nomenclature by referring
236 to sequential rational numbers outside the
237 interval [0,1] as Farey terms or as part of
238 $F_N$, which, technically, they are not.
239 All of the results presented in
240 this chapter, can be shown to apply
241 everywhere in $\vworkrealsetnonneg$, so this abuse
242 is not harmful.
243
244 It can be proved that if $h/k$ is irreducible, then
245 $(h+ik)/k$ is also irreducible (i.e. the analogous terms
246 in $[i, i+1]$ corresponding to the Farey terms in
247 $[0,1]$ are also irreducible).
248
249 Recursive formulas do exist for generating
250 successive terms of the Farey series.
251 Given two successive terms of a Farey series of
252 order $N$, $h_{j-2}/k_{j-2}$ and $h_{j-1}/k_{j-1}$
253 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq01})
254 and
255 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq02})
256 can be applied to generate the next term,
257 $h_{j}/k_{j}$\@. In applying
258 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq01})
259 and
260 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq02}), the
261 two terms used as input must be irreducible.
262
263 \begin{equation}
264 \label{eq:crla1:slcr0:sfry0:thm:01:eq01}
265 h_{j} = \left\lfloor {\frac{{k_{j-2}
266 + N}}{{k_{j - 1} }}} \right\rfloor h_{j - 1} - h_{j-2}
267 \end{equation}
268
269 \begin{equation}
270 \label{eq:crla1:slcr0:sfry0:thm:01:eq02}
271 k_{j} = \left\lfloor {\frac{{k_{j-2} + N}}{{k_{j
272 - 1} }}} \right\rfloor k_{j - 1} - k_{j-2}
273 \end{equation}
274
275 Similarly, given two successive terms of a Farey series of
276 order $N$, $h_{j+1}/k_{j+1}$ and $h_{j+2}/k_{j+2}$
277 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq03})
278 and
279 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq04})
280 can be applied to generate the preceding term,
281 $h_{j}/k_{j}$\@. Again, in applying
282 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq03})
283 and
284 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq04}), the
285 two terms used as input must be irreducible.
286
287 \begin{equation}
288 \label{eq:crla1:slcr0:sfry0:thm:01:eq03}
289 h_j = \left\lfloor {\frac{{k_{j + 2} + N}}{{k_{j + 1} }}}
290 \right\rfloor h_{j + 1} - h_{j + 2}
291 \end{equation}
292
293 \begin{equation}
294 \label{eq:crla1:slcr0:sfry0:thm:01:eq04}
295 k_j = \left\lfloor {\frac{{k_{j + 2} + N}}{{k_{j + 1} }}}
296 \right\rfloor k_{j + 1} - k_{j + 2}
297 \end{equation}
298
299 It might appear on the surface that (\ref{eq:crla1:slcr0:sfry0:thm:01:eq01})
300 through (\ref{eq:crla1:slcr0:sfry0:thm:01:eq04}) provide a viable
301 method for finding best rational approximations (i.e.
302 start with $0/1$ and $1/N$ and generate successive terms until
303 $r_I$ is bracketed), but in fact they do not.
304 The number of terms in the Farey series is approximately quadratic with
305 respect to the order of the series, and so generating terms starting at
306 an integer boundary is $O(N^2)$ and doesn't scale up well
307 finding best rational approximations for processors that can
308 operate on 32- and 64-bit integers.
309
310 The Farey series has a convenient and intuitive graphical interpretation
311 involving the integer lattice\index{integer lattice}\index{lattice}%
312 \index{Farey series!integer lattice interpretation}
313 (see Fig. \ref{fig:crla1:slcr0:sfry0:00},
314 which illustrates this interpretation, but with $h$
315 also restricted).
316 [By integer lattice, we mean the $\vworkrealset{}^2$ plane
317 with each point $(x,y)$, $x,y \in \vworkintset$, marked.]
318 In such an interpretation, each rational number $h/k$ corresponds to
319 a point $(k,h)$ which is $h$ units above and $k$ units
320 to the right of the origin.
321
322 \begin{figure}
323 \centering
324 \includegraphics[width=4.6in]{c_rla1/farey01a.eps}
325 \caption{Graphical Interpretation Of Rational Numbers
326 $h/k$ That Can Be Formed With $h \leq h_{MAX}=3$, $k \leq k_{MAX}=5$}
327 \label{fig:crla1:slcr0:sfry0:00}
328 \end{figure}
329
330 From the graphical interpretation suggested by Fig. \ref{fig:crla1:slcr0:sfry0:00},
331 the following properties are intuitively clear.
332
333 \begin{itemize}
334 \item The angle of a ray drawn from the origin to the point
335 $(k,h)$ corresponding to the rational number $h/k$ is
336 $\theta = tan^{-1} \; h/k$.
337
338 \item Any integer lattice point on a line from
339 the origin drawn at the angle $\theta$
340 has the value $h/k = tan \; \theta$. All points corresponding
341 to rational numbers with the same value will be on such a line,
342 and thus form an equivalence class.
343
344 \item A rational number $h/k$ is irreducible if and only if its corresponding
345 point $(k,h)$ is ``directly'' visible from the origin with
346 no intervening points.
347
348 \item The Farey series of order $N$, $F_N$, can be
349 formed graphically by starting with the
350 set of integer lattice points
351 $(k,h): \; h \in \vworkintsetnonneg \wedge 1 \leq k \leq N$,
352 then sweeping
353 a line extended from the origin, starting with
354 angle $\theta = 0$, through
355 $0 \leq \theta < \pi{}/2$, and recording
356 in order each point directly visible from
357 the origin.\footnote{Note that Fig. \ref{fig:crla1:slcr0:sfry0:00},
358 because it illustrates the case when $h$ is constrained
359 as well, does not show integer lattice points for
360 $h > h_{MAX}$. In principle, if the integer lattice shown
361 in Fig. \ref{fig:crla1:slcr0:sfry0:00} were extended indefinitely
362 ``upward'', every positive irreducible rational number with
363 $k \leq k_{MAX} = 5$ could be found graphically.}
364 \end{itemize}
365
366 Fig. \ref{fig:crla1:slcr0:sfry0:01} illustrates the graphical construction method
367 of $F_5$. Note that only integer lattice points which are directly
368 visible from the origin (with no intervening points) are selected.
369 (Fig. \ref{fig:crla1:slcr0:sfry0:01}, like Fig. \ref{fig:crla1:slcr0:sfry0:00},
370 shows the case of constrained $h$---the integer lattice should be
371 continued ``upward'' to construct the Farey series.)
372
373 \begin{figure}
374 \centering
375 \includegraphics[width=4.6in]{c_rla1/farey01b.eps}
376 \caption{Graphical Interpretation Of Irreducible Rational Numbers
377 $h/k$ That Can Be Formed With $h \leq h_{MAX}=3$, $k \leq k_{MAX}=5$}
378 \label{fig:crla1:slcr0:sfry0:01}
379 \end{figure}
380
381 Note that Figures \ref{fig:crla1:slcr0:sfry0:00}
382 and \ref{fig:crla1:slcr0:sfry0:01} depict the case when
383 \emph{both} $h$ and $k$ are constrained (whereas the Farey series
384 constrains only $k$).
385
386 To give a compact notation, we denote the set of ascending irreducible
387 rational numbers that can be graphically formed
388 from Figures \ref{fig:crla1:slcr0:sfry0:00}
389 and \ref{fig:crla1:slcr0:sfry0:01} as
390 $F_{k_{MAX}, h_{MAX}}$. Using this notation, the graphical construction method
391 depicted in Figure \ref{fig:crla1:slcr0:sfry0:01} identifies $F_{5,3}$.
392
393 The ``corner point'' in Figures \ref{fig:crla1:slcr0:sfry0:00}
394 and \ref{fig:crla1:slcr0:sfry0:01} plays a special role:
395
396 \begin{itemize}
397 \item From 0/1 up through the corner point $h_{MAX}/k_{MAX}$,
398 the terms are the terms of the Farey series of order
399 $k_{MAX}$.
400 \item From $h_{MAX}/k_{MAX}$ up through $h_{MAX}/1$, the terms
401 are the reverse-ordered reciprocals of the terms of the
402 Farey series of order $h_{MAX}$.\footnote{This can be verified
403 by transposing the $h$ and $k$ axes of the figures.}
404 \end{itemize}
405
406 As an example, $F_{5,3}$ identified graphically
407 is Figure \ref{fig:crla1:slcr0:sfry0:01} is
408
409 \begin{equation}
410 \label{eq:crla1:slcr0:sfry0:10}
411 F_5 = \left\{ {\frac{0}{1},\frac{1}{5},\frac{1}{4},
412 \frac{1}{3},\frac{2}{5},\frac{1}{2},
413 \frac{3}{5},\frac{2}{3},\frac{3}{4},
414 \frac{1}{1},\frac{3}{2}, \frac{2}{1},
415 \frac{3}{1} } \right\} .
416 \end{equation}
417
418 It can be seen by comparing (\ref{eq:crla1:slcr0:sfry0:10})
419 with (\ref{eq:crla1:slcr0:sfry0:eq0001c}) and
420 (\ref{eq:crla1:slcr0:sfry0:eq0001e}) that the first seven terms of
421 (\ref{eq:crla1:slcr0:sfry0:10}) come from $F_5$ and the remaining
422 six terms are the reverse-ordered reciprocals of $F_3$ (although $F_3$
423 must be extended from Eq. \ref{eq:crla1:slcr0:sfry0:eq0001c}
424 into [1,2] to include 4/3 and 3/2).
425
426 The symmetry of Figures \ref{fig:crla1:slcr0:sfry0:00} and
427 \ref{fig:crla1:slcr0:sfry0:01} with respect to the corner point
428 is important, because it means that finding best
429 rational approximations for
430 $r_I < h_{MAX}/k_{MAX}$ is the same problem as for
431 $r_I > h_{MAX}/k_{MAX}$. In the case of
432 $r_I < h_{MAX}/k_{MAX}$, the Farey series of order
433 $k_{MAX}$ is used, and in the case of
434 $r_I > h_{MAX}/k_{MAX}$ the Farey series of order
435 $h_{MAX}$ is used (but the reciprocals of the terms are
436 used, the order of the terms is reversed, and $1/r_I$ is used). Thus, if we
437 know how to bracket $r_I < h_{MAX}/k_{MAX}$
438 in $F_{k_{MAX}}$, we can approach
439 the problem of $r_I > h_{MAX}/k_{MAX}$
440 through the inherent symmetry.
441
442
443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
445 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
446
447 \subsection{The Continued Fraction Algorithm}
448
449 %Subsection Tag: cfr0
450 \label{crla1:slcr0:scfr0}
451
452 A \emph{finite simple continued fraction} is a fraction of the form
453
454 \begin{equation}
455 \label{eq:crla1:slcr0:scfr0:00}
456 a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2
457 + \cfrac{1}{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ldots + \cfrac{1}{a_n}}}}
458 =
459 [a_0; a_1, a_2, \ldots , a_n] ,
460 \end{equation}
461
462 \noindent{}where $a_0 \in \vworkintsetnonneg$ and
463 $a_i \in \vworkintsetpos$, $i > 0$. Each integer
464 $a_i$ is called an \index{continued fraction!element}\emph{element} or
465 \index{continued fraction!partial quotient}\emph{partial quotient}
466 of the continued fraction.
467 To ensure a unique representation, we require, except in the case of
468 the continued fraction representation of an integer,
469 that the final element $a_n$ not be equal
470 to 1.
471
472 Continued fractions are quite unwieldly to write and typeset,
473 and so a continued fraction in the form of (\ref{eq:crla1:slcr0:scfr0:00})
474 is written as $[a_0; a_1, a_2, \ldots , a_n]$. Note that the
475 separator between $a_0$ and $a_1$ is a semicolon (`;'), and that all other
476 separators are commas (`,'). In some works, commas are used exclusively; and in
477 other works, the first element is $a_1$ rather than $a_0$. Throughout this
478 work, the notational conventions illustrated in (\ref{eq:crla1:slcr0:scfr0:00}) are
479 followed.
480
481 Continued fractions can be either finite or infinite.
482
483 A finite continued fraction consists of a finite number of elements
484 $[a_0; a_1, a_2, \ldots , a_n]$. It can be proved that
485 every rational number corresponds to a unique
486 finite continued fraction\footnote{So long as the $a_n \neq 1$ convention
487 described earlier is followed.}, and that
488 every finite continued fraction corresponds to a rational number.
489
490 An infinite continued fraction consists of an infinite number
491 of elements $[a_0; a_1, a_2, \ldots]$. Because every rational number
492 corresponds to a finite continued fraction, all irrational numbers have
493 infinite continued fraction representations.
494
495 In engineering work (and due to the general prevalence of computers
496 and calculators), any $r_I$ to be approximated has a known approximate
497 numerical value. Even quantities that are known to be irrational (such
498 as $\pi$ or $\sqrt{2}$) have a numerical value known to a large number
499 of significant digits. For this reason, only the numerical procedure for
500 obtaining the continued fraction representation of a rational number is
501 presented here (the symbolic procedure is not discussed). Numerical values
502 are always rational (for example, 3.1415 is 31,415/10,000).
503
504 \index{continued fraction!convergent}
505 The \emph{kth order convergent} of a continued fraction
506 $[a_0; a_1, \ldots{}, a_n]$ is the irreducible rational number
507 corresponding to $[a_0; a_1, \ldots{}, a_k]$, $k \leq n$.
508 In other words, the $k$th order convergent is the irreducible rational number
509 corresponding to the first $k+1$ partial quotients of a
510 continued fraction.\footnote{``$k+1$'' because the notational
511 numbering
512 for partial quotients starts at 0 rather than 1.}
513
514 An $n$th order continued fraction $[a_0; a_1, \ldots{}, a_n]$
515 has $n+1$ convergents, $[a_0]$,
516 $[a_0; a_1]$, \ldots{}, and $[a_0; a_1, \ldots{}, a_n]$.
517 We denote the $k$th order convergent as $s_k$, with numerator
518 $p_k$ and denominator $q_k$.
519
520 Without proof, we present the following algorithm, Algorithm
521 \ref{alg:crla1:slcr0:scfr0:akgenalg}, for
522 determining the continued fraction representation (i.e. the partial
523 quotients) as well as the convergents of a non-negative
524 rational number $a/b$.
525
526 \begin{vworkalgorithmstatementpar}{Continued Fraction Representation and
527 Convergents of
528 A Rational Number \mbox{\boldmath $a/b$}}
529 \label{alg:crla1:slcr0:scfr0:akgenalg}
530 \begin{alglvl0}
531 \item $k:=-1$.
532 \item $divisor_{-1} := a$.
533 \item $remainder_{-1} := b$.
534
535 \item Repeat
536
537 \begin{alglvl1}
538 \item $k := k + 1$.
539 \item $dividend_k := divisor_{k-1}$.
540 \item $divisor_k := remainder_{k-1}$.
541 \item $a_k := dividend_k \; div \; divisor_k$.
542 \item $remainder_k := dividend_k \; mod \; divisor_k$.
543 \item If $k=0$, $p_0 = a_0$; else if $k=1$, $p_1 = a_0 a_1 + 1$;
544 else $p_i = a_i p_{i-1} + p_{i-2}$.
545 \item If $k=0$, $q_0 = 1$; else if $k=1$, $q_1 = a_1$;
546 else $q_i = a_i q_{i-1} + q_{i-2}$.
547 \end{alglvl1}
548
549 \item Until ($remainder_k = 0$).
550 \end{alglvl0}
551 \textbf{\emph{Note:}} The final $s_k = p_k / q_k$ is the irreducible
552 form of $a/b$. For brevity, this is not proved here.
553 \end{vworkalgorithmstatementpar}
554 %\vworkalgorithmfooter{}
555 \begin{vworkexamplestatement}
556 \label{ex:crla1:slcr0:scfr0:01}
557 Find the continued fraction partial quotients and convergents of
558 $67/29$.
559 \end{vworkexamplestatement}
560 \begin{vworkexampleparsection}{Solution} Table
561 \ref{tbl:crla1:slcr0:scfr0:01} shows the application of
562 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg} to find the
563 continued fraction partial quotients and convergents of $67/29$. From
564 Table \ref{tbl:crla1:slcr0:scfr0:01}, the continued fraction
565 representation of $67/29$ is $[2;3,4,2]$.
566
567 \begin{table}
568 \caption{Continued Fraction Partial Quotients and Convergents of $67/29$ (Example \ref{ex:crla1:slcr0:scfr0:01})}
569 \label{tbl:crla1:slcr0:scfr0:01}
570 \begin{center}
571 \begin{tabular}{|c|c|c|c|c|c|c|}
572 \hline
573 \small{Index} & \small{$dividend_k$} & \small{$divisor_k$} & \small{$a_k$} & \small{$remainder_k$} & \small{$p_k$} & \small{$q_k$} \\
574 \small{($k$)} & & & & & & \\
575 \hline
576 \hline
577 \small{-1} & \small{N/A} & \small{67} & \small{N/A} & \small{29} & \small{N/A} & \small{N/A} \\
578 \hline
579 \small{0} & \small{67} & \small{29} & \small{2} & \small{9} & \small{2} & \small{1} \\
580 \hline
581 \small{1} & \small{29} & \small{9} & \small{3} & \small{2} & \small{7} & \small{3} \\
582 \hline
583 \small{2} & \small{9} & \small{2} & \small{4} & \small{1} & \small{30} & \small{13} \\
584 \hline
585 \small{3} & \small{2} & \small{1} & \small{2} & \small{0} & \small{67} & \small{29} \\
586 \hline
587 \end{tabular}
588 \end{center}
589 \end{table}
590 \end{vworkexampleparsection}
591 \vworkexamplefooter{}
592
593 Finally, we present without proof a theorem that indicates how to bracket
594 a rational number $a/b$ that is not in $F_{k_{MAX}}$ with its two neighbors
595 in $F_{k_{MAX}}$.
596
597 \begin{vworktheoremstatementpar}{Enclosing Neighbors Of \mbox{\boldmath $x \notin F_N$}
598 In \mbox{\boldmath $F_N$}}
599 \label{thm:crla1:slcr0:scfr0:cfenclosingneighbors}
600 For a non-negative rational
601 number $a/b$\footnote{It is not required that $a/b$ be irreducible in
602 order to apply this theorem.
603 For brevity, many properties of convergents were omitted; it is provable that
604 the convergents formed by Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}
605 will be identical for $a/b$ and $ia/ib$.} not in
606 $F_N$ which has a
607 continued fraction representation
608 $[a_0;a_1,a_2,\ldots{} ,a_n]$, the
609 highest-order convergent $s_k = p_k/q_k$ with $q_k \leq N$ is one
610 neighbor\footnote{By neighbors in $F_N$ we mean the rational numbers
611 in $F_N$ immediately to the left and immediately to the right
612 of $a/b$.}
613 to $a/b$ in $F_N$, and the other neighbor in
614 $F_N$ is\footnote{Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors}
615 is a somewhat stronger statement about best approximations
616 than Khinchin makes in \cite{bibref:b:KhinchinClassic}, Theorem 15.
617 We were not able to locate
618 this theorem or a proof in print,
619 but this theorem is understood within the number theory community.
620 It appears on the Web
621 page of David Eppstein \cite{bibref:i:davideppstein} in the form of a
622 `C'-language computer program,
623 \texttt{http://www.ics.uci.edu/\~{}{}eppstein/numth/frap.c}.
624 Although
625 Dr. Eppstein phrases the solution in terms of modifying
626 a partial quotient, his approach is equivalent to
627 (\ref{eq:crla1:slcr0:scfr0:thm:cfenclosingneighbors:01}).}
628
629 \begin{equation}
630 \label{eq:crla1:slcr0:scfr0:thm:cfenclosingneighbors:01}
631 \frac{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}} \right\rfloor}
632 p_k + p_{k - 1} }}{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}}
633 \right\rfloor} q_k + q_{k - 1} }}.
634 \end{equation}
635 \end{vworktheoremstatementpar}
636 \begin{vworktheoremproof}
637 Omitted, as it relies on material not presented for brevity.
638 \end{vworktheoremproof}
639 \vworktheoremfooter{}
640
641 Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors}
642 can also be applied to find the Farey neighbors of an $a/b$ already
643 in $F_{k_{MAX}}$. If Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}
644 is applied to $a/b$, (\ref{eq:crla1:slcr0:scfr0:thm:cfenclosingneighbors:01})
645 will provide one Farey neighbor, and
646 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq01}) through
647 (\ref{eq:crla1:slcr0:sfry0:thm:01:eq04}) can be used to provide the other
648 Farey neighbor. (Again, for brevity, the mathematical basis for this
649 is not presented.)
650
651 Many constants $r_I$ to be approximated are engineering constants based
652 on measurements or arbitrary conventions, and so are known or accepted to
653 only a finite number of significant digits. Such constants are always
654 rational, and
655 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}
656 and
657 Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors}
658 can be applied with no special consideration.
659
660 Some constants, however, are irrational. The question naturally arises
661 of how to be sure that one is using enough decimal digits
662 in applying
663 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}
664 and
665 Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors}.
666 The easiest approach to apply in practice\footnote{\emph{In practice}
667 because some theoretical results may be possible as far as how
668 many significant digits are always adequate or as far as other
669 criteria, but the approach taken here is the easiest practical one.}
670 is to confine the quantity of interest by an inequality and to be
671 sure that the results are the same at both boundaries
672 of the inequality.
673
674 For example, $\pi$ is a transcendental constant, so has a non-terminating
675 decimal representation. $\pi$ to several digits (truncated at the
676 end) is 3.1415926535. It follows that
677
678 \begin{equation}
679 \label{eq:crla1:slcr0:scfr0:30}
680 3.1415926535 < \pi < 3.1415926536
681 \end{equation}
682
683 For compact notation, we denote the left limit as $r_{LEFT}$ and the
684 right limit by $r_{RIGHT}$. We also make the observation that in some
685 applications, the interval is closed rather than open.\footnote{This depends
686 on how much is known about $r_I$---for example, we know that $\pi$ is irrational and
687 can't be equal to any rational number, but we
688 may not know this about other $r_I$ of interest.} In the more
689 general case, $r_I$ is confined by:
690
691 \begin{equation}
692 \label{eq:crla1:slcr0:scfr0:30b}
693 r_{LEFT} \leq r_I < r_{RIGHT} .
694 \end{equation}
695
696 With the the $r_I$ of interest confined as suggested in
697 (\ref{eq:crla1:slcr0:scfr0:30b}), there are two easy approaches
698 to decide if the Farey neighbors of $r_I$ can be determined with
699 the information available.
700
701 \begin{enumerate}
702 \item \emph{Easier:} Locate the Farey neighbors $h_L/k_L$ and $h_R/k_R$ of
703 $r_{LEFT}$, then numerically determine whether
704 $h_L/k_L \leq r_{RIGHT} \leq h_R/k_R$.
705 \item \emph{Harder:} Determine whether $r_{LEFT}$ and $r_{RIGHT}$ have the
706 same convergents up through $p_k/q_k$ with $q_k \leq N$. If so,
707 $r_{LEFT}$ and $r_{RIGHT}$ have the same Farey neighbors.
708 \end{enumerate}
709
710 \noindent{}If
711 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}
712 and
713 Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors} are
714 applied with 31415926535/10000000000 and 31415926536/10000000000
715 separately and yield the same rational numbers $h_L/k_L$ and $h_R/k_R$ as
716 left and right neighbors, then
717
718 \begin{equation}
719 \label{eq:crla1:slcr0:scfr0:31}
720 \frac{h_L}{k_L} < 3.1415926535 < \pi < 3.1415926536 < \frac{h_R}{k_R}
721 \end{equation}
722
723 \noindent{}and it is thus confirmed that $h_L/k_L$ and $h_R/k_R$ are the left
724 and right neighbors of $\pi$ in the Farey series of interest.
725
726 Several examples follow which illustrate the technique and various special
727 cases.
728
729 \begin{vworkexamplestatement}
730 \label{ex:crla1:slcr0:scfr0:10}
731 Find the best rational approximation to $1/e$ in $F_{65535}$.
732 \end{vworkexamplestatement}
733 \begin{vworkexampleparsection}{Solution} Note that $e$ is irrational, implying
734 that $1/e$ is also irrational, thus a bounding technique might best be used
735 to ensure finding the correct Farey neighbors. Using an ordinary scientific pocket
736 calculator, the displayed value of $e^{-1}$ is approximately 0.367879441171.
737 Allowing for some possible imprecision in the last digit\footnote{The guess at
738 how accurate a calculator is likely to be is subjective.}, it is fairly safe to assume that
739
740 \begin{equation}
741 \label{eq:ex:crla1:slcr0:scfr0:10:01}
742 \frac{367879441170}{1000000000000} < \frac{1}{e} < \frac{367879441172}{1000000000000} .
743 \end{equation}
744
745 Table \ref{tbl:crla1:slcr0:scfr0:10a} shows the application of
746 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg} to find the
747 continued fraction partial quotients of the left inequality
748 limit, and Table Table \ref{tbl:crla1:slcr0:scfr0:10b} shows the
749 calculation of the convergents. (The tables are separated due to typesetting
750 limitations---the partial quotients and convergents would normally be tabulated
751 together.)
752
753 \begin{table}
754 \caption{Continued Fraction Partial Quotients of $367,879,441,170/1,000,000,000,000$ (Example \ref{ex:crla1:slcr0:scfr0:10})}
755 \label{tbl:crla1:slcr0:scfr0:10a}
756 \begin{center}
757 \begin{tabular}{|c|c|c|c|c|}
758 \hline
759 \small{Index} & \small{$dividend_k$} & \small{$divisor_k$} & \small{$a_k$} & \small{$remainder_k$} \\
760 \small{($k$)} & & & & \\
761 \hline
762 \hline
763 \small{-1} & \small{N/A} & \small{367,879,441,170} & \small{N/A} & \small{1,000,000,000,000} \\
764 \hline
765 \small{0} & \small{367,879,441,170} & \small{1,000,000,000,000} & \small{0} & \small{367,879,441,170} \\
766 \hline
767 \small{1} & \small{1,000,000,000,000} & \small{367,879,441,170} & \small{2} & \small{264,241,117,660} \\
768 \hline
769 \small{2} & \small{367,879,441,170} & \small{264,241,117,660} & \small{1} & \small{103,638,323,510} \\
770 \hline
771 \small{3} & \small{264,241,117,660} & \small{103,638,323,510} & \small{2} & \small{56,964,470,640} \\
772 \hline
773 \small{4} & \small{103,638,323,510} & \small{56,964,470,640} & \small{1} & \small{46,673,852,870} \\
774 \hline
775 \small{5} & \small{56,964,470,640} & \small{46,673,852,870} & \small{1} & \small{10,290,617,770} \\
776 \hline
777 \small{6} & \small{46,673,852,870} & \small{10,290,617,770} & \small{4} & \small{5,511,381,790} \\
778 \hline
779 \small{7} & \small{10,290,617,770} & \small{5,511,381,790} & \small{1} & \small{4,779,235,980} \\
780 \hline
781 \small{8} & \small{5,511,381,790} & \small{4,779,235,980} & \small{1} & \small{732,145,810} \\
782 \hline
783 \small{9} & \small{4,779,235,980} & \small{732,145,810} & \small{6} & \small{386,361,120} \\
784 \hline
785 \small{10} & \small{732,145,810} & \small{386,361,120} & \small{1} & \small{345,784,690} \\
786 \hline
787 \small{11} & \small{386,361,120} & \small{345,784,690} & \small{1} & \small{40,576,430} \\
788 \hline
789 \small{12} & \small{345,784,690} & \small{40,576,430} & \small{8} & \small{21,173,250} \\
790 \hline
791 \small{13} & \small{40,576,430} & \small{21,173,250} & \small{1} & \small{19,403,180} \\
792 \hline
793 \small{14} & \small{21,173,250} & \small{19,403,180} & \small{1} & \small{1,770,070} \\
794 \hline
795 \small{15} & \small{19,403,180} & \small{1,770,070} & \small{10} & \small{1,702,480} \\
796 \hline
797 \small{16} & \small{1,770,070} & \small{1,702,480} & \small{1} & \small{67,590} \\
798 \hline
799 \small{17} & \small{1,702,480} & \small{67,590} & \small{25} & \small{12,730} \\
800 \hline
801 \small{18} & \small{67,590} & \small{12,730} & \small{5} & \small{3,940} \\
802 \hline
803 \small{19} & \small{12,730} & \small{3,940} & \small{3} & \small{910} \\
804 \hline
805 \small{20} & \small{3,940} & \small{910} & \small{4} & \small{300} \\
806 \hline
807 \small{21} & \small{910} & \small{300} & \small{3} & \small{10} \\
808 \hline
809 \small{22} & \small{300} & \small{10} & \small{30} & \small{0} \\
810 \hline
811 \end{tabular}
812 \end{center}
813 \end{table}
814
815 \begin{table}
816 \caption{Continued Fraction Convergents of $367,879,441,170/1,000,000,000,000$ (Example \ref{ex:crla1:slcr0:scfr0:10})}
817 \label{tbl:crla1:slcr0:scfr0:10b}
818 \begin{center}
819 \begin{tabular}{|c|c|c|c|}
820 \hline
821 \small{Index} & \small{$a_k$} & \small{$p_k$} & \small{$q_k$} \\
822 \small{($k$)} & & & \\
823 \hline
824 \hline
825 \small{-1} & \small{N/A} & \small{N/A} & \small{N/A} \\
826 \hline
827 \small{0} & \small{0} & \small{0} & \small{1} \\
828 \hline
829 \small{1} & \small{2} & \small{1} & \small{2} \\
830 \hline
831 \small{2} & \small{1} & \small{1} & \small{3} \\
832 \hline
833 \small{3} & \small{2} & \small{3} & \small{8} \\
834 \hline
835 \small{4} & \small{1} & \small{4} & \small{11} \\
836 \hline
837 \small{5} & \small{1} & \small{7} & \small{19} \\
838 \hline
839 \small{6} & \small{4} & \small{32} & \small{87} \\
840 \hline
841 \small{7} & \small{1} & \small{39} & \small{106} \\
842 \hline
843 \small{8} & \small{1} & \small{71} & \small{193} \\
844 \hline
845 \small{9} & \small{6} & \small{465} & \small{1,264} \\
846 \hline
847 \small{10} & \small{1} & \small{536} & \small{1,457} \\
848 \hline
849 \small{11} & \small{1} & \small{1,001} & \small{2,721} \\
850 \hline
851 \small{12} & \small{8} & \small{8,544} & \small{23,225} \\
852 \hline
853 \small{13} & \small{1} & \small{9,545} & \small{25,946} \\
854 \hline
855 \small{14} & \small{1} & \small{18,089} & \small{49,171} \\
856 \hline
857 \small{15} & \small{10} & \small{190,435} & \small{517,656} \\
858 \hline
859 \small{16} & \small{1} & \small{208,524} & \small{566,827} \\
860 \hline
861 \small{17} & \small{25} & \small{5,403,535} & \small{14,688,331} \\
862 \hline
863 \small{18} & \small{5} & \small{27,226,199} & \small{74,008,482} \\
864 \hline
865 \small{19} & \small{3} & \small{87,082,132} & \small{236,713,777} \\
866 \hline
867 \small{20} & \small{4} & \small{375,554,727} & \small{1,020,863,590} \\
868 \hline
869 \small{21} & \small{3} & \small{1,213,746,313} & \small{3,299,304,547} \\
870 \hline
871 \small{22} & \small{30} & \small{36,787,944,117} & \small{100,000,000,000} \\
872 \hline
873 \end{tabular}
874 \end{center}
875 \end{table}
876
877 Note that since we are finding best rational approximations
878 in $F_{65535}$, Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors}
879 requires only that we carry the tabular procedure of
880 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg} out until
881 $q_k \geq 65535$ ($k=15$ in this example).
882
883 By
884 Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors}
885 and
886 Table \ref{tbl:crla1:slcr0:scfr0:10b}, one Farey neighbor
887 to the left inequality limit is $p_{14}/q_{14} = 18,089/49,171$.
888 The other Farey neighbor is given by
889 (\ref{eq:crla1:slcr0:scfr0:thm:cfenclosingneighbors:01}), with the calculation
890 detailed below.
891
892 \begin{equation}
893 \label{eq:ex:crla1:slcr0:scfr0:10:50a}
894 \frac{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}} \right\rfloor}
895 p_k + p_{k - 1} }}{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}}
896 \right\rfloor} q_k + q_{k - 1} }}.
897 \end{equation}
898
899 \begin{equation}
900 \label{eq:ex:crla1:slcr0:scfr0:10:50b}
901 = \frac{{\displaystyle{\left\lfloor {\frac{{65535 - q_{13} }}{{q_{14} }}} \right\rfloor}
902 p_{14} + p_{13} }}{{\displaystyle{\left\lfloor {\frac{{65535 - q_{13} }}{{q_{14} }}}
903 \right\rfloor} q_{14} + q_{13} }}.
904 \end{equation}
905
906 \begin{equation}
907 \label{eq:ex:crla1:slcr0:scfr0:10:50c}
908 = \frac{{\displaystyle{\left\lfloor {\frac{{65535 - 25946 }}{{49171 }}} \right\rfloor}
909 18089 + 9545 }}{{\displaystyle{\left\lfloor {\frac{{65535 - 25946 }}{{49171 }}}
910 \right\rfloor} 49171 + 25946 }}.
911 \end{equation}
912
913 \begin{equation}
914 \label{eq:ex:crla1:slcr0:scfr0:10:50d}
915 = \frac{9545}{25946}
916 \end{equation}
917
918 It can be verified by cross-multiplication that $9545/25946 > 18089/49171$, therefore
919
920 \begin{equation}
921 \label{eq:ex:crla1:slcr0:scfr0:10:50e}
922 \frac{18089}{49171} < 0.367879441170 < \frac{9545}{25946} .
923 \end{equation}
924
925 A similar tabulation procedure carried out with
926 the right inequality limit (0.367879441172), not reproduced here for brevity,
927 verifies that it has the same convergents up through $s_{15}$ as the left
928 inequality limit. Therefore it has the same Farey neighbors in $F_{65535}$
929 as the left limit, and it is established that
930
931 \begin{equation}
932 \label{eq:ex:crla1:slcr0:scfr0:10:50f}
933 \frac{18089}{49171} < 0.367879441170 < \frac{1}{e} < 0.367879441172 < \frac{9545}{25946} .
934 \end{equation}
935
936 It can be verified numerically that 18089/49171 and 9545/25946 are both \emph{very}
937 close approximations to $1/e$ (with differences on the order of a couple parts per
938 \emph{billion}).
939 \end{vworkexampleparsection}
940 %\vworkexamplefooter{}
941
942 \begin{vworkexamplestatement}
943 \label{ex:crla1:slcr0:scfr0:11}
944 Find the enclosing neighbors to $8/43$ in $F_{65535}$.
945 \end{vworkexamplestatement}
946 \begin{vworkexampleparsection}{Solution} As hinted earlier,
947 since $8/43 \in F_{65535}$, we can simply treat 8/43 as a number very close to 8/43
948 so that the convergent after $s_k = 8/43$ has a larger denominator than 65535 (the details
949 of why and how this works are omitted for brevity).
950
951 \begin{table}
952 \caption{Continued Fraction Partial Quotients and Convergents of $8/43$ (Example \ref{ex:crla1:slcr0:scfr0:11})}
953 \label{tbl:crla1:slcr0:scfr0:11a}
954 \begin{center}
955 \begin{tabular}{|c|c|c|c|c|c|c|}
956 \hline
957 \small{Index} & \small{$dividend_k$} & \small{$divisor_k$} & \small{$a_k$} & \small{$remainder_k$} & \small{$p_k$} & \small{$q_k$} \\
958 \small{($k$)} & & & & & & \\
959 \hline
960 \hline
961 \small{-1} & \small{N/A} & \small{8} & \small{N/A} & \small{43} & \small{N/A} & \small{N/A} \\
962 \hline
963 \small{0} & \small{8} & \small{43} & \small{0} & \small{8} & \small{0} & \small{1} \\
964 \hline
965 \small{1} & \small{43} & \small{8} & \small{5} & \small{3} & \small{1} & \small{5} \\
966 \hline
967 \small{2} & \small{8} & \small{3} & \small{2} & \small{2} & \small{2} & \small{11} \\
968 \hline
969 \small{3} & \small{3} & \small{2} & \small{1} & \small{1} & \small{3} & \small{16} \\
970 \hline
971 \small{4} & \small{2} & \small{1} & \small{2} & \small{0} & \small{8} & \small{43} \\
972 \hline
973 \end{tabular}
974 \end{center}
975 \end{table}
976
977 Table \ref{tbl:crla1:slcr0:scfr0:11a} shows the
978 application of
979 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}
980 to determine the partial quotients and convergents of 8/43.
981
982 (\ref{eq:crla1:slcr0:scfr0:thm:cfenclosingneighbors:01}) can then be applied
983 to find a Farey neighbor of 8/43, with the calculation
984 detailed below.\footnote{Whether the left or right Farey neighbor is found depends on
985 whether the final convergent has $k$ even or $k$ odd. The explanation is beyond the scope here.}
986
987 \begin{equation}
988 \label{eq:ex:crla1:slcr0:scfr0:11:50a}
989 \frac{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}} \right\rfloor}
990 p_k + p_{k - 1} }}{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}}
991 \right\rfloor} q_k + q_{k - 1} }}.
992 \end{equation}
993
994 \begin{equation}
995 \label{eq:ex:crla1:slcr0:scfr0:11:50b}
996 = \frac{{\displaystyle{\left\lfloor {\frac{{65535 - q_{3} }}{{q_{4} }}} \right\rfloor}
997 p_{4} + p_{3} }}{{\displaystyle{\left\lfloor {\frac{{65535 - q_{3} }}{{q_{4} }}}
998 \right\rfloor} q_{4} + q_{3} }}.
999 \end{equation}
1000
1001 \begin{equation}
1002 \label{eq:ex:crla1:slcr0:scfr0:11:50c}
1003 = \frac{{\displaystyle{\left\lfloor {\frac{{65535 - 16 }}{{43 }}} \right\rfloor}
1004 8 + 3 }}{{\displaystyle{\left\lfloor {\frac{{65535 - 16 }}{{43 }}}
1005 \right\rfloor} 43 + 16 }}.
1006 \end{equation}
1007
1008 \begin{equation}
1009 \label{eq:ex:crla1:slcr0:scfr0:11:50d}
1010 = \frac{12187}{65505}
1011 \end{equation}
1012
1013 It can be verified by cross-multiplication that (\ref{eq:ex:crla1:slcr0:scfr0:11:50d})
1014 is the right Farey neighbor of 8/43. As two consecutive
1015 terms in $F_{65535}$ are now known, (\ref{eq:crla1:slcr0:sfry0:thm:01:eq03})
1016 and (\ref{eq:crla1:slcr0:sfry0:thm:01:eq04}) can be applied to find the left Farey
1017 neighbor, as shown below.
1018
1019 \begin{equation}
1020 \label{eq:ex:crla1:slcr0:scfr0:11:51}
1021 h_j = \left\lfloor {\frac{{k_{j + 2} + N}}{{k_{j + 1} }}}
1022 \right\rfloor h_{j + 1} - h_{j + 2}
1023 \end{equation}
1024
1025 \begin{equation}
1026 \label{eq:ex:crla1:slcr0:scfr0:11:52}
1027 = \left\lfloor {\frac{{65505 + 65535}}{{43 }}}
1028 \right\rfloor 8 - 12187 = 12189
1029 \end{equation}
1030
1031 \begin{equation}
1032 \label{eq:ex:crla1:slcr0:scfr0:11:55}
1033 k_j = \left\lfloor {\frac{{k_{j + 2} + N}}{{k_{j + 1} }}}
1034 \right\rfloor k_{j + 1} - k_{j + 2}
1035 \end{equation}
1036
1037 \begin{equation}
1038 \label{eq:ex:crla1:slcr0:scfr0:11:56}
1039 = \left\lfloor {\frac{{65505 + 65535}}{{43 }}}
1040 \right\rfloor 43 - 65505 = 65516
1041 \end{equation}
1042
1043 Thus, 12189/65516 is the left neighbor to 8/43 in $F_{65535}$.
1044 \end{vworkexampleparsection}
1045 %\vworkexamplefooter{}
1046
1047 \begin{vworkexamplestatement}
1048 \label{ex:crla1:slcr0:scfr0:12}
1049 Create assembly-language code to form an approximation
1050 to multiplication by $\pi$ that
1051 can be implemented using a \texttt{MUL} instruction followed by
1052 a \texttt{DIV} instruction on the Freescale CPU08 core. Assume that:
1053 \begin{itemize}
1054 \item The input is an 8-bit unsigned integer, passed in the accumulator.
1055 \item The output is an 8-bit unsigned integer, returned in the accumulator.
1056 \item The code is ``in-line'' (i.e. not a subroutine), and the code is free to modify
1057 any processor registers or flags.
1058 \item Input data that is too large should cause the output to be clipped at 255 (the
1059 maximum value for an unsigned byte).
1060 \end{itemize}
1061 \end{vworkexamplestatement}
1062 \begin{vworkexampleparsection}{Solution} $\pi$ is transcendental, and based
1063 on the first several digits, can be bounded by
1064
1065 \begin{equation}
1066 \label{eq:ex:crla1:slcr0:scfr0:12:01}
1067 3.1415926535 < \pi < 3.1415926536 .
1068 \end{equation}
1069
1070 We first note that, due to the charactertistics of the CPU08,
1071 $h_{MAX} = 255$ and $k_{MAX} = 255$. We are thus operating in
1072 $F_{255, 255}$, with a corner point of $h/k = 255/255 = 1$.
1073 $\pi > 1$, so we are operating along the top side of Figures
1074 \ref{fig:crla1:slcr0:sfry0:00} and \ref{fig:crla1:slcr0:sfry0:01}
1075 (pages \pageref{fig:crla1:slcr0:sfry0:00} and \pageref{fig:crla1:slcr0:sfry0:01}).
1076 In this region above the corner point, the numerator rather than the denominator
1077 is constrained.
1078
1079 As discussed earlier, to obtain best rational approximations above the corner point,
1080 we may transpose the problem to that of obtaining best rational approximations
1081 to $1/\pi$ in $F_{h_{MAX}}$ (rather than $F_{k_{MAX}}$---it is coincidental in this
1082 case that $h_{MAX} = k_{MAX}$).
1083
1084 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg} requires only a rational number as a starting
1085 point, so it is most expedient to invert the terms of
1086 (\ref{eq:ex:crla1:slcr0:scfr0:12:01}) to yield
1087
1088 \begin{equation}
1089 \label{eq:ex:crla1:slcr0:scfr0:12:02}
1090 \frac{10000000000}{31415926536}
1091 <
1092 \frac{1}{\pi}
1093 <
1094 \frac{10000000000}{31415926535} .
1095 \end{equation}
1096
1097 \begin{table}
1098 \caption{Continued Fraction Partial Quotients and Convergents of $10000000000/31415926536$ (Example \ref{ex:crla1:slcr0:scfr0:12})}
1099 \label{tbl:ex:crla1:slcr0:scfr0:12a}
1100 \begin{center}
1101 \begin{tabular}{|c|c|c|c|c|c|c|}
1102 \hline
1103 \small{$k$} & \small{$dividend_k$} & \small{$divisor_k$} & \small{$a_k$} & \small{$remainder_k$} & \small{$p_k$} & \small{$q_k$} \\
1104 \hline
1105 \hline
1106 \small{-1} & \small{N/A} & \small{10,000,000,000} & \small{N/A} & \small{31,415,926,536} & \small{N/A} & \small{N/A} \\
1107 \hline
1108 \small{0} & \small{10,000,000,000} & \small{31,415,926,536} & \small{0} & \small{10,000,000,000} & \small{0} & \small{1} \\
1109 \hline
1110 \small{1} & \small{31,415,926,536} & \small{10,000,000,000} & \small{3} & \small{1,415,926,536} & \small{1} & \small{3} \\
1111 \hline
1112 \small{2} & \small{10,000,000,000} & \small{1,415,926,536} & \small{7} & \small{88,514,248} & \small{7} & \small{22} \\
1113 \hline
1114 \small{3} & \small{1,415,926,536} & \small{88,514,248} & \small{15} & \small{88,212,816} & \small{106} & \small{333} \\
1115 \hline
1116 \multicolumn{7}{|c|}{\small{It isn't necessary to carry Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg}}} \\
1117 \multicolumn{7}{|c|}{\small{any further, as it has been established that $s_2 = 7/22$ is the last}} \\
1118 \multicolumn{7}{|c|}{\small{convergent with $q_k \leq 255$.}} \\
1119 \hline
1120 \end{tabular}
1121 \end{center}
1122 \end{table}
1123
1124 Table \ref{tbl:ex:crla1:slcr0:scfr0:12a} shows the application of
1125 Algorithm \ref{alg:crla1:slcr0:scfr0:akgenalg} to obtain the partial
1126 quotients and convergents of $10000000000/31415926536$. Note that it is not
1127 necessary to carry out the algorithm any further than establishing the
1128 highest-order convergent with $q_k \leq h_{MAX} = 255$.
1129
1130 By Theorem \ref{thm:crla1:slcr0:scfr0:cfenclosingneighbors} 7/22 is either
1131 the left or right Farey neighbor to 10000000000/31415926536. Numerically,
1132 it can be verified that 7/22 is the left Farey neighbor.
1133
1134 The right Farey neighbor is given by
1135 (\ref{eq:crla1:slcr0:scfr0:thm:cfenclosingneighbors:01}), with the calculation
1136 detailed below.
1137
1138 \begin{equation}
1139 \label{eq:ex:crla1:slcr0:scfr0:12:10a}
1140 \frac{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}} \right\rfloor}
1141 p_k + p_{k - 1} }}{{\displaystyle{\left\lfloor {\frac{{N - q_{k - 1} }}{{q_k }}}
1142 \right\rfloor} q_k + q_{k - 1} }}
1143 \end{equation}
1144
1145 \begin{equation}
1146 \label{eq:ex:crla1:slcr0:scfr0:12:10b}
1147 = \frac{{\displaystyle{\left\lfloor {\frac{{255 - q_{1} }}{{q_{2} }}} \right\rfloor}
1148 p_{2} + p_{1} }}{{\displaystyle{\left\lfloor {\frac{{255 - q_{1} }}{{q_{2} }}}
1149 \right\rfloor} q_{2} + q_{1} }}
1150 \end{equation}
1151
1152 \begin{equation}
1153 \label{eq:ex:crla1:slcr0:scfr0:12:10c}
1154 = \frac{{\displaystyle{\left\lfloor {\frac{{255 - 3 }}{{22 }}} \right\rfloor}
1155 7 + 1 }}{{\displaystyle{\left\lfloor {\frac{{255 - 3 }}{{22 }}}
1156 \right\rfloor} 22 + 3 }}
1157 \end{equation}
1158
1159 \begin{equation}
1160 \label{eq:ex:crla1:slcr0:scfr0:12:10d}
1161 = \frac{78}{245}
1162 \end{equation}
1163
1164 We have shown that
1165
1166 \begin{equation}
1167 \label{eq:ex:crla1:slcr0:scfr0:12:20}
1168 \frac{7}{22}
1169 <
1170 \frac{10000000000}{31415926536}
1171 <
1172 \frac{78}{245} .
1173 \end{equation}
1174
1175 However, we can't bound $1/\pi$ without checking whether the right
1176 limit of
1177 (\ref{eq:ex:crla1:slcr0:scfr0:12:02}) falls between
1178 7/22 and 78/245. It is easy to verify numerically that this is the case.
1179 (If this were not the case, Eq. \ref{eq:ex:crla1:slcr0:scfr0:12:02} should
1180 be rephrased using more digits of $\pi$.)
1181
1182 It is then known that:
1183
1184 \begin{equation}
1185 \label{eq:ex:crla1:slcr0:scfr0:12:21}
1186 \frac{7}{22}
1187 <
1188 \frac{10000000000}{31415926536}
1189 <
1190 \frac{1}{\pi}
1191 <
1192 \frac{10000000000}{31415926536}
1193 <
1194 \frac{78}{245} .
1195 \end{equation}
1196
1197 In order to convert (\ref{eq:ex:crla1:slcr0:scfr0:12:21})
1198 to a form involving constrained $h_{MAX}$ and $\pi$, we must
1199 invert the terms and reverse the order of the inequality.
1200
1201 \begin{equation}
1202 \label{eq:ex:crla1:slcr0:scfr0:12:22}
1203 \frac{245}{78}
1204 <
1205 \frac{31415926536}{10000000000}
1206 <
1207 \pi
1208 <
1209 \frac{31415926535}{10000000000}
1210 <
1211 \frac{22}{7}
1212 \end{equation}
1213
1214 It has thus been shown that, subject to the constraints
1215 $h \leq 255$ and $k \leq 255$, the surrounding rational approximations
1216 to $\pi$ are 245/78 and 22/7. Since 245/78 is the better approximation, that
1217 is used for the assembly-language code (Figure \ref{fig:ex:crla1:slcr0:scfr0:12:10}).
1218
1219 \begin{figure}
1220 \begin{verbatim}
1221 ;Assume input argument in accumulator.
1222 cmp #82
1223 bhs toobig ;Input argument is 82 or larger. Must
1224 ;clip or will overflow MUL and/or DIV.
1225 ldx #245 ;Set up to multiply A by 245.
1226 mul ;X:A now contains MSB:LSB multiplication
1227 ;result.
1228 pshx ;Push/pull best way to get X into H
1229 pulh ;to set up for division.
1230 ldx #78 ;78 is divisor.
1231 div ;Do the division. Result guaranteed to
1232 ;be in range as divisor could be no
1233 ;larger than 19,845 before division.
1234 bra theend ;Branch around clip.
1235 toobig: lda #255 ;Load "clip" value.
1236 theend:
1237 ;Output result now in accumulator.
1238 \end{verbatim}
1239 \caption{Freescale CPU08 Code to Approximate $\pi$ (Example \ref{ex:crla1:slcr0:scfr0:12})}
1240 \label{fig:ex:crla1:slcr0:scfr0:12:10}
1241 \end{figure}
1242
1243 In order to ``clip'' so that any input arguments too large result in an output of 255,
1244 it is necessary to determine which input arguments will be problematic. This is done
1245 in the inequality below, and the result appears in the assembly-language code of
1246 (Figure \ref{fig:ex:crla1:slcr0:scfr0:12:10}).
1247
1248 \begin{equation}
1249 \label{eq:ex:crla1:slcr0:scfr0:12:23}
1250 \left({\left\lfloor{\frac{245 x}{78}}\right\rfloor
1251 \geq 256}\right)
1252 \longrightarrow \left({x \geq 82}\right)
1253 \end{equation}
1254 \end{vworkexampleparsection}
1255 \vworkexamplefooter{}
1256
1257
1258 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1260 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1261
1262 \section[Choosing $r_A = h/2^q \approx r_I$]
1263 {Choosing \mbox{\boldmath $r_A = h/2^q \approx r_I$}}
1264
1265 %Section Tag: lcr0
1266 \label{crla1:slcr1}
1267
1268
1269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1272
1273 \section[Choosing $r_A = 2^p/k \approx r_I$]
1274 {Choosing \mbox{\boldmath $r_A = 2^p/k \approx r_I$}}
1275
1276 %Section Tag: lcr0
1277 \label{crla1:slcr2}
1278
1279
1280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1281 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1282 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1283
1284 \section{End-to-End Approximation Error}
1285
1286 %Section Tag: ete0
1287 \label{crla1:sete2}
1288
1289 %
1290 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1291 %
1292 %\noindent\begin{figure}[!b]
1293 %\noindent\rule[-0.25in]{\textwidth}{1pt}
1294 %\begin{tiny}
1295 %\begin{verbatim}
1296 %$RCSfile: c_rla1.tex,v $
1297 %$Source: /home/dashley/cvsrep/uculib01/uculib01/doc/manual/c_rla1/c_rla1.tex,v $
1298 %$Revision: 1.11 $
1299 %$Author: dashley $
1300 %$Date: 2010/01/28 21:18:33 $
1301 %\end{verbatim}
1302 %\end{tiny}
1303 %\noindent\rule[0.25in]{\textwidth}{1pt}
1304 %\end{figure}
1305 %
1306 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1307 %% $Log: c_rla1.tex,v $
1308 %% Revision 1.11 2010/01/28 21:18:33 dashley
1309 %% a)Chapter start quotes removed.
1310 %% b)Aesthetic comment line added at the bottom of most files.
1311 %%
1312 %% Revision 1.10 2007/10/01 14:20:01 dtashley
1313 %% Example completed.
1314 %%
1315 %% Revision 1.9 2007/10/01 02:02:49 dtashley
1316 %% Edits.
1317 %%
1318 %% Revision 1.8 2007/09/30 21:59:51 dtashley
1319 %% Edits.
1320 %%
1321 %% Revision 1.7 2007/09/29 04:58:42 dtashley
1322 %% Edits.
1323 %%
1324 %% Revision 1.6 2007/09/29 03:15:56 dtashley
1325 %% Edits.
1326 %%
1327 %% Revision 1.5 2007/09/28 19:59:55 dtashley
1328 %% Edits.
1329 %%
1330 %% Revision 1.4 2007/09/28 04:59:49 dtashley
1331 %% Edits.
1332 %%
1333 %% Revision 1.3 2007/09/27 22:54:33 dtashley
1334 %% Edits.
1335 %%
1336 %% Revision 1.2 2007/09/27 21:44:22 dtashley
1337 %% Edits.
1338 %%
1339 %% Revision 1.1 2007/09/27 15:23:31 dtashley
1340 %% Initial checkin.
1341 %%
1342 %%End of $RCSfile: c_rla1.tex,v $.
1343 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1344

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25