HAKMEM -- CONTINUED FRACTIONS -- DRAFT, NOT YET PROOFEDBeeler, M., Gosper, R.W., and Schroeppel, R. HAKMEM. MIT AIMemo 239, Feb. 29, 1972.Retyped and converted to html ('Web browser format) byHenry Baker,April, 1995.CONTINUED FRACTIONSPreviousUpNextITEM 97 (Schroeppel): CONTINUED FRACTIONSSimple proofs that certain continued fractions are sqrt(2), sqrt(3),etc. Proof for sqrt(2):X = [1, 2, 2, 2, ...](X-1) (X+1) = [0, 2, 2, 2, ...] * [2, 2, 2, 2, ...] = 1 2X - 1 = 1X = sqrt(2)Proof for sqrt(3):Y = [1, 1, 2, 1, 2, ...](Y + 1) (Y - 1) = [2, 1, 2, 1, 2, ...] * [0, 1, 2, 1, 2, ...] = 2 * [1, 2, 1, 2, 1, ...] * [0, 1, 2, 1, 2, ...] = 2 2Y - 1 = 2Y = sqrt(3)Similar proofs exist for sqrt(5) and sqrt(6); but sqrt(7) is hairy.ITEM 98 (Schroeppel):The continued fraction expansion of the positive minimum of thefactorial function (about 0.46) is[0, 2, 6, 63, 135, 1, 1, 1, 1, 4, 1, 43, ...].ITEM 99 (Schroeppel):The value of a continued fraction with partial quotients increasing inarithmetic progression is I (2/D) A/D[A+D, A+2D, A+3D, ...] = ------------- I (2/D) 1+(A/D)where the I's are Bessel functions. A special case is I (2) 0[1, 2, 3, 4, ...] = ----- . I (2) 1ITEM 100 (Perron): n /===\ ! ! 1 ! ! (1 + --) = ! ! Ak k = 1 1 (A1 + 1)A1 (A2 + 1)A2 (A(n-1) + 1)A(n-1)1 + ----- ---------- ---------- ... ------------------ . A1 - A1+A2+1 - A2+A3+1 - A(n-1)+An+1ITEM 101A (Gosper):On the theory that continued fractions are underused, probably becauseof their unfamiliarity, I offer the following propaganda session onthe relative merits of continued fractions versus other numericalrepresentations. For a good cram course in continued fractions, seeKnuth, Volume 2, page 316 [1st Edition; s.4.5.3, p. 339 in 2ndEdition]. (In what follows, "regular" means that all numerators are1, and any radix can be read in place of decimal.)0) pi is 3. But not really 3, more like 3 + 1/7. But notreally 7, more like 7 + 1/15. But not really 15, ... . So theregular continued fraction for pi is written 3 7 15 1 292 1 1... .1) The continued fractions for rational numbers always come out even,and rather quickly. Thus, the number of inches per meter is exactly100/2.54 or 39 2 1 2 2 1 4. The corresponding decimal fraction is39.3700787... has period 42, making it almost impossible to tell ifthe number is rational. (But if our data are ALL rational, theordered pair 5000/127 is even more concise.)2) Quadratic surds, which are of course inexpressible as rationals, aregenerally unrecognizable in decimal. Their continued fractions, on the otherhand, are periodic. Nth roots of e^2, ratios of Bessel functions, and ratiosof linear functions of these all have regular continued fractions formed byinterleaving one or more arithmetic sequences. These special properties willshow up regardless of number base. You might recognize 5.436563... as 2e, buteven Schroeppel might not notice that 6.1102966796... was 2/34 e - 2---------- 2/3 e - 1until he wrote it as 6 9 15 21 27 33 ... .The familiar transcendental functions of rational arguments also havesimple continued fractions, but these are generally not regular andcannot be reconstructed from numerical values by a simple algorithm,since nonregular representations aren't unique. The point is,however, that numbers like e, pi, sqrt(2), sin .5,sqrt(7) arctan sqrt(7), etc. can be expressed to unlimited precisionby simple programs which produce the terms on demand.3) If we define a rational approximation to be "best" if it comescloser than any other rational with such a small denominator, thencontinued fractions give the complete set of best rationalapproximations to the value which they represent. That is, if youtruncate a (regular) continue fraction at any point, then theresulting rational number is a best approximation. Furthermore, thisremains true if the last term of the approximation is replaced by anysmaller positive integer other than 1. All best approximations can begenerated in this manner, in order of increasing denominators (ornumerators). For example, the approximants to pi = 3 7 15 1292 ... are: 3: 1/1, 2/1, 3/1 7: (4/1), 7/2, 10/3, 13/4, 16/5, 19/6, 22/715: (25/8), 47/15, 69/22, 91/29, 113/36, ... 311/99, 333/106 1: 355/113... ...Note that they are all automatically in lowest terms. the size of adenominator is greater than the product of the terms involved and lessthan the product of the numbers 1 greater than the terms. Theapproximations are low if the number of terms is odd, high if it'seven. (Note that if a 1 ends a continued fraction, it should be addedin to the previous term. Thus, to "round off" a continued fractionafter a certain term, add in the next term iff it is +-1. In theabove, 4/1 and 25/8 correspond to termination with a 1 and are not"best"; 355/113 is "best" because the corresponding term really shouldbe 1.) The error is smaller than 1 over the product of thedenominator squared and the first neglected term, so that the totalnumber of digits (numerator and denominator) is usually slightlysmaller than with equally accurate decimal fractions. 355/113 is goodto 7.5 places instead of 5.5, due to the unusually large term (292)which follows.4) Numerical comparison of continue fractions is slightly harder thanin decimal, but much easier than with rationals -- just invert thedecision as to which is larger whenever the first discrepant terms areeven-numbered. Contrast this with the problem of comparing therationals 113/36 and 355/113.5) Regular continued fractions are in 1 to 1 correspondence with thereal numbers, unlike decimal (.5 = .49999...) or rationals (2/3 = 6/9,sqrt(6) = ?). Even infinity has a continued fraction, namely,the empty one! (Minus and plus infinity are the same in continuedfraction notation.)6) Each representation favors certain operations. Decimal favorsmultiplication by powers of 10. Rationals favor reciprocation, as do continuedfractions. To reciprocate a regular continued fraction, add (or if possible,delete) and initial 0 term. To negate, negate all the terms, optionallyobserving that-a, -b, -c, -d ... = -a-1, 1, b-1, c, d ... .7) The strongest argument for positional (e.g., decimal or floating)representation for non-integers is that arithmetic is easy. Rationalnumber arithmetic often loses because numerators and denominators growso large as to require icky multiprecision. Algorithms for arithmeticon continued fractions seem generally unknown. The next itemsdescribe how to arithmetically combine continued fractions to producenew ones, one term at a time.Unfortunately, the effort required to perform these operationsmanually is several times that for decimal, but the rewards formachine implementation are considerable (which can also be said offloating point). Specifically, these rewards will be seen to be:unlimited significance arithmetic without multiprecisionmultiplication or division, built in error analysis, immorally easycomputation of algebraic functions, no unnecessary computations, nodiscarding of information (as with roundoff and truncation),reversibility of computations, and the terms of the answer start tocome out right away and continue to do so until shut off.ITEM 101B (Gosper): Continued Fraction ArithmeticContinued fractions let us perform numerical calculations a little ata time without ever introducing any error, such as roundoff ortruncation. As if this weren't enough, the calculations provideautomatic error analysis, and obviate most forms of successiveapproximation. This means we can start with an arithmetic expressionlike sqrt(3/pi^2 + e)-----------------------tanh(sqrt(5)) - sin(69)and immediately begin to produce the value as a sequence of continuedfraction terms (or even decimal digits, if we should be soreactionary), limited only by time and storage. If there arequantities in the expression which are known only approximately, thecalculation can provide error bounds on the answer as well as identifythe quantity that limited the significance.All this is possible because each operation (+, /, -, sqrt) in thearithmetic expression requests terms from the continued fractions ofits operands only when necessary, and consequently produces terms ofits own value as soon as possible. Numbers like pi ande and functions like sin and tanh havecontinued fraction terms in simple sequences which can be produced byshort programs. Imprecise quantities can also be programs whichdeliver terms until they run out of confidence, whereupon theyinitiate special action. By then, the last guaranteeable term of theoverall expression will have already been produced.We see then that no calculation is performed unnecessarily, so that,for example, a subexpression which happened to be multiplied by zerowould never be evaluated. Also, an operation detecting a deficiencyin two or more of its operands provides a natural mechanism forallocating multiprocessor resources, should you have some.Here are the algorithms for the elementary arithmetic operations oncontinued fractions.Let x be a continued fractionp0 + q0/(p1 + q1/( ... = p0 + q0/x'where x' is again a continued fraction and the p'sand q's are integers. We shall call a (p q) pair a"term" of the continued fraction for x. Often, only thep's are mentioned, in which case the q's areimplicitly all 1, and x is called a "regular" continued fraction.Instead of a list of p's and q's, let x bea computer subroutine which produces its next p andq each time it is called. Thus on its first usage xwill "output" p0 and q0 and, in effect, changeitself into x'. Similarly, let y be anotherprocedurally represented continued fraction r0 + s0/y'. Ourproblem will be solved if we can write such subroutines forz(x,y) = x+y, x-y, x y, andx/y. When called upon to output a term of z, thesubroutine might in turn call for (or "input") terms from xand y until it is satisfied that the unread portions ofx and y cannot affect the pending term ofz. Then it would output this term and change itself intoz', so that it could produce the next term next time.Unfortunately, when we try to do this, our expressions quicklycomplicate. Let us preempt this complication by computing instead themore general functionz(x,y) = (a x y + b x + c y + d) / (e x y + f x + g y + h)(or(a b c d) / (e f g h)for short)where a through h are integer variables whoseinitial values we are free to choose. Various choices express addition: x + y = (0 1 1 0) / (0 0 0 1), subtraction: x - y = (0 1 -1 0) / (0 0 0 1),multiplication: x y = (1 0 0 0) / (0 0 0 1), division: x / y = (0 1 0 0) / (0 0 1 0).As we shall see, the process of inputting terms of x andy and outputting terms of z will reduce to replacingthe eight integers a through h with linearcombinations of each other.When z inputs a term of x, z becomes a newfunction of x'. To see how this happens, substitute p +q/x' for every occurrence of x in the expression forz(x,y), then multiply numerator and denominator through byx':z(x',y) = (pa+c pb+d qa qb) / (pe+g pf+h qe qf).If x was rational and has run out of terms, it has in effectbecome infinite:z(infinity,y) = (0 0 a b) / (0 0 e f)If instead we input a term of y by substituting r +s/y' for every occurrence of y:z(x,y') = (ra+b sa rc+d sc) / (re+f se rg+h sg).If y runs out of terms:z(x,infinity) = (0 a 0 c) / (0 e 0 g)To output the term (t u), so that z = t + u/z'(i.e., z' = u/(z-t)):z'(x,y) = (ue uf ug uh) / (a-te b-tf c-tg d-th).Thus this basic eight variable form is preserved by all threeoperations, which can be performed in any order since they representindependent substitutions.For simplicity, let us assume that z will output in standardform, that is, every u = 1 (regular) and every output termt >= 1 except perhaps the first. This means thatz' will always exceed 1 and thus 0 <= u/z' <1, so that the integer t = z - u/z' must = [z], thegreatest integer <= z.Since z generally varies with x and y, itshould not output unless [z] is constant for the range ofpossible x and y. We can easily compute the rangeof given the ranges of x and y if we represent eachrange by the endpoints of an interval (in either order), along with abit indicating Inside or Outside. Thus if z isin standard form, we can say that z will always be (Inside 1infinity) (or (Outside -infinity 1)) after the first term. Ifz were to always output its nearest integer instead of isgreatest, then none of the terms after the first would be 1, althoughthey would probably vary in sign. In this case, z would be(Outside -2 2).Now hold y fixed and examine the behavior of z withx. If x is (Inside a b) thenz is (Inside z(a) z(b)) unless thedenominator of z changes sign between a andb (i.e., z has its pole in this interval), whereuponz is (Outside z(a) z(b)). Symmetrically, when x is(Outside a b) then z is (Outsidez(a) z(b)) unless the signs of the denominators ofz(a) and z(b) differ, whereupon z is(Inside z(a) z(b)). This argument still holds withx and y interchanged.Now suppose that with y fixed at one if its endpoints,x constrains z (Inside 1 2), and at y'sother extreme, z(x) is z(y) is (Outside 0 3).Suppose further that at the two extremes of x, z(y)is (Inside 1 3) and (Outside 0 2). Then z(x,y) is (Outside 01), the union of the four ranges. (Outside 0 2) is the widest,indicating that z will probably get more information from aterm of y than a term of x. (Topology hackersshould recognize this Inside-Outside nonsense as ordinary intervals intoroidal space. The clue is that both plus and minus infinity aredenoted by the empty continued fraction.)Due to the basically monotonic behavior of z, we canguarantee that the actual range of z will be the union ofthese four ranges, and that this range will be Inside or Outside someinterval. If it is (Inside z1 z2) and [z1]= [z2], z can output the term t =[z1]. Otherwise, z must input a term fromx or y, whichever was associated with the widest ofthe four ranges of z. (Outside narrowness) is wider than(Outside wideness) is wider than (Inside wideness) is wider than(Inside narrowness).Evaluating z on these endpoints may be facilitated by keepingestimates for the integer variables in floating point.Even if z doesn't produce a term, narrowing the range ofpossible z will still help in computing the range of afunction of z, especially if z gets stuck trying tooutput the last term of a rational number resulting from irrationalx and y. (There is no way to guarantee thatx or y won't eventually deviate, whereuponz would egest a gigantic term.)z can produce its value as decimal digits by multiplying by10 instead of reciprocating, after outputting t = [z]:z'(x,y) = (10(a-te) 10(b-tf) 10(c-tg) 10(d-th)) / (e f g h).Strange to say, it is not serious if z for some reasonoutputs the terms 7 5 1 when it should have produced 6 9. As soon aspermitted, it will simply recant with 0 -1 -5 and continue with thecorrection -1 9. The sequence 7 5 1 0 -1 -5 -1 9 is equivalent to 6 9because b 0 c is the same as b+c. In orderto undo these computations, z violates the condition (Outside-1 1) when it is 0 -1 -5 ... . This condition is obeyed by nearly allconvergent continued fractions after their first term, and itsviolation will very probably cause further retractions among thefunctions dependent upon z.This computation reversal trick is also handy for mechanizing anddenoting imprecise quantities. Instead of 2.997930+-.000003, we have2 1 481 0 2, meaning between 2 1 481 and 2 1 483. Similarly, 137 26 01 replaces 137.0373+-.0006.Successive approximations methods benefit considerably from notrequesting terms until needed. Consider Newton's method for algebraicroots. We expect successive approximations to have about twice asmany correct terms each time. Since the production of these termscannot be aided by reading incorrect terms, the additional correctterms must be produced before the bad ones of the previousapproximation are used. But this means that there is no need to readin the bad ones at all. By feeding back the output terms in place ofthe approximation, we get the correct answer directly! (69% of thecredit for this goes to Schroeppel.)The basic eight variable form exemplified above by z(x,y) isnot the only form preserved by continued fraction term transactions.We need only four variables and a single interval check to compute a x + bz(x) = ------- , c x + dthe homographic function of one argument. On the other hand,z(w,x,y) (linear in all three arguments) requiressixteen variables and a twelve way interval check. Eachof these forms can be solved for x in terms of zetc. to get a function of the same form. This is not true of 2 a x + b x + cz(x) = -------------- , 2 d x + e x + ffor example, even though this form is also preserved. This form isnot guaranteed monotone, thus theoretically invalidating the intervalcheck algorithm, but it hardly every errs. Even if it did, it wouldquickly correct itself anyway. This form is not only more economicalthan z(x,x), it is essential for the success of the Newton'smethod feedback trick, which must know when two variables are reallythe same one.By choosing the eight coefficients a through hproperly, it should be possible to rewrite arithmetic expressions ascompositions of considerably fewer of these forms than one for each +,-, *, and /. The reader is invited to investigate the problem oftrying to find minimal representations. Depending on the metric forminimality, the question can be complicated by allowing higher powersof x and y. If the highest powers of x,y, z, ... in an invariant form are i,j, k, .., then the number of integer variablesrequired for the coefficients (mostly because of all the cross terms)is2 (i+1) (j+1) (k+1) ... .It is awkward in this system to evaluate transcendental functions ofirrational arguments. The problem is that you may need any number ofcontinued fraction (or series) terms which, instead of being numbers,are symbolic functions of x, some infinite continuedfraction. My suggestion is to represent each symbolic term of thefunction by a subroutine which is a function of x and thenext term, with this next term really a dummy until actually calledupon for output, whereupon it replaces itself with a full-fledged termsubroutine which in turn refers to x and a new dummy.Sad to say, the integer variables in these algorithms do not usuallyshrink on outputs as much as they grow on inputs. Fortunately, theoperations for input and output only require (besides addition)multiplication by terms which are almost invariably small. (I havenot seen a term exceed 20776 except in specially constructed numbers.)It is fairly safe, then to declare any function which has gotten(Outside -2^35 2^35) to be infinite, thus terminating its continuedfraction. Better still, note that the term 20776 is equivalent to theterms 20000 0 700 0 70 0 6, i.e., a very large term can be transmittedpiecewise. Although this is just thinly disguised multiprecisionmultiplication, that first piece of the term will probably satisfy itsrecipient for quite some time.In some special cases, the integer variables will becomeperiodic rather than large, especially when all by one of thearguments to a function have terminated. Then, we have the form a x + bz(x) = ------- , c x + dknown as a homographic function. If ad-bc is +-1,then a, b, c, d will eventuallybecome 1, 0, 0, 1, whereupon z will output the terms ofx unmodified. Periodicity will also occur when x isa Hurwitz number, i.e., when the terms of x are thevalues of one or more polynomials evaluated on consecutive integersand then interleaved. Coth 1/69, sqrt(105), and e are Hurwitznumbers whose polynomials are linear or constant. Hurwitzness ispreserved by homographic functions. If one can show that pi isa Hurwitz number, one confirms the long standing conjectures thate*pi, e+pi, e/pi, etc. areall irrational.If z, x, and y are all regular, then itgenerally won't be possible to reduce z by finding a GCD ofa through h which is > 1. However, it has beendetermined empirically that much reduction is often possible in othercases. This reduction is almost always by a divisor of an input oroutput term numerator (or 10 if output is decimal digits) and can befacilitated by keeping certain of the integer variables around modulothese quantities.ITEM 101C (Gosper):Problem: Given an interval, find in it the rational number with thesmallest numerator and denominator.Solution: Express the endpoints as continued fractions. Find thefirst term where they differ and add 1 to the lesser term, unless it'slast. Discard the terms to the right. What's left is the continuedfraction for the "smallest" rational in the interval. (If onefraction terminates but matches the other as far as it goes, append aninfinity and proceed as above.)PreviousUpNext |
|
