The Lumber Room

"Consign them to dust and damp by way of preserving them"

Not all best rational approximations are the convergents of the continued fraction!

with 15 comments

[For more background on continued fractions and why they are so wonderful at approximations (and wonderful generally) — eventually I may edit this post to mention that. For now I just want to quickly clarify something, which surprisingly many popular expositions of continued fractions seem to mislead by leaving out.]

Any real number can be written as a continued fraction. For instance, the square root of 3 is
\sqrt{3} = 1+\cfrac{1}{1+\cfrac{1}{2+\cfrac{1}{1+\cfrac{1}{2+\dots}}}}
written as \sqrt{3} = [1; 1, 2, 1, 2, \dots],

and π is

\pi = 3+\cfrac{1}{7+\cfrac{1}{15+\cfrac{1}{1+\cfrac{1}{292+\dots}}}}

i.e. \pi = [3; 7, 15, 1, 292, \dots] (more terms below),

and e is e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1,\dots] (proved here).

Rational numbers have finite continued fractions, and irrational numbers have infinite continued fractions, of which only quadratic irrationals (like √3 above) have periodic continued fractions. Non-quadratic irrationals may still have some recognisable patterns in their representation, like e above, or not, like π.

The numbers you get by truncating the continued fraction — taking only the first few terms — are called its convergents. The continued fraction for π is [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2…] (sequence A001203 in OEIS), and the sequence of convergents is 3/1, 22/7, 333/106, 355/113, 103993/33102… (A002485/A002486):

\begin{array}{cll} p_1/q_1 &= [3;] &= 3/1 \\ p_2/q_2 &= [3; 7] &= 22/7 \\  p_3/q_3 &= [3; 7, 15] &= 333/106 \\  p_4/q_4 &= [3; 7, 15, 1] &= 355/113 \\ p_5/q_5 &= [3; 7, 15, 1, 292] &= 103993/33102\end{array}

Each of these convergents is a best rational approximation, in the following sense.

Definition: A fraction p/q is a best rational approximation of a real number x if p/q is closer to x than any other fraction with a smaller denominator. That is, for any other fraction p'/q' with q' \le q and p'/q' \ne p/q, we have \left|\frac{p}{q} -x \right| < \left| \frac{p'}{q'} -x \right|.

Thus among all fractions with denominator at most 7, the fraction 22/7 is closest to π. Among all fractions with denominator at most 113, the fraction 355/113 is closest to π. And so on.

Every convergent is a best rational approximation, but these aren’t all of the best rational approximations! The sequence of best approximations of π is actually 3/1, 13/4, 16/5, 19/6, 22/7, 179/57, 201/64, 223/71, 245/78, 267/85, 289/92, 311/99, 333/106, 355/113, 52163/16604… (A063674/A063673), where only the bold terms are convergents. (Aside: Notice how good 355/113 is: it is the best among all fractions with denominator at most 16603. This is because of the large term 292 in the continued fraction of π.) The others are semi-convergents (intermediate fractions), as explained below.

In general, for any real number, the convergents p_k/q_k of its continued fraction are always best rational approximations, but they aren’t all of the best rational approximations. To get all of them, you also have to take the semi-convergents: fractions of the form (p_k + np_{k+1})/(q_k +nq_{k+1}) for some integer n≥1. Specifically, taking n=a_{k+2} — the next term in the continued fraction expansion — gives p_{k+2}/q_{k+2}, and the ones to take are precisely n from \lceil a_{k+2}/2 \rceil (if it’s better than the previous convergent, else \lfloor a_{k+2}/2 \rfloor + 1) to a_{k+2}. This is also mentioned on Wikipedia.

Thus for π = [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2…], the best rational approximations are:

3/1     = [3]

13/4    = [3; 4]
16/5    = [3; 5]
19/6    = [3; 6]
22/7    = [3; 7]

179/57  = [3; 7, 8]
201/64  = [3; 7, 9]
223/71  = [3; 7, 10]
245/78  = [3; 7, 11]
267/85  = [3; 7, 12]
289/92  = [3; 7, 13]
311/99  = [3; 7, 14]
333/106 = [3; 7, 15]

355/113 = [3; 7, 15, 1]

52163/16604  = [3; 7, 15, 1, 146]
52518/16717  = [3; 7, 15, 1, 147]
… (terms from [3; 7, 15, 1, 148] to [3; 7, 15, 1, 291])...
103993/33102 = [3; 7, 15, 1, 292]

104348/33215 = [3; 7, 15, 1, 292, 1]
...

The above is the natural meaning of “best rational approximation”. This is also known as “best rational approximation of the first kind”. There is a different definition of “best rational approximation” under which the convergents are precisely all the best rational approximations: basically, instead of considering the distance |x - p/q|, we consider |qx - p|. That is: A fraction \frac{a}{b} is a best approximation of the second kind of a real number x if for every fraction \frac{c}{d} with 1 \le d \le b and a/b \neq c/d, we have |dx - c| > |bx - a|
This measure turns out to have a nicer theory, but to omit specifying that this non-obvious definition is the one being used is misleading.

Program

For illustration, here’s a C program (why it’s written in C is not important right now) that given a positive real number, generates its continued fraction, its convergents, and the sequence of best rational approximations. The function find_cf finds the continued fraction (putting the terms in a[] and the convergents in p[] and q[] — excuse the global variables), and the function all_best prints all the best rational approximations.

#include <math.h>
#include <stdio.h>
#include <assert.h>

// number of terms in continued fraction.
// 15 is the max without precision errors for M_PI
#define MAX 15
#define eps 1e-9

long p[MAX], q[MAX], a[MAX], len;
void find_cf(double x) {
  int i;
  //The first two convergents are 0/1 and 1/0
  p[0] = 0; q[0] = 1;
  p[1] = 1; q[1] = 0;
  //The rest of the convergents (and continued fraction)
  for(i=2; i<MAX; ++i) {
    a[i] = lrint(floor(x));
    p[i] = a[i]*p[i-1] + p[i-2];
    q[i] = a[i]*q[i-1] + q[i-2];
    printf("%ld:  %ld/%ld\n", a[i], p[i], q[i]);
    len = i;
    if(fabs(x-a[i])<eps) return;
    x = 1.0/(x - a[i]);
  }
}

void all_best(double x) {
  find_cf(x); printf("\n");
  int i, n; long cp, cq;
  for(i=2; i<len; ++i) {
    //Test n = a[i+1]/2 separately. Enough to test when a[i+1] is even, actually.
    n = a[i+1]/2; cp = n*p[i]+p[i-1]; cq = n*q[i]+q[i-1];
    if(fabs(x-(double)cp/cq) < fabs(x-(double)p[i]/q[i])) printf("%ld/%ld,", cp, cq);
    //And print all the rest, no need to test
    for(n = (a[i+1]+2)/2; n<=a[i+1]; ++n) {
      printf("%ld/%ld,", n*p[i]+p[i-1], n*q[i]+q[i-1]);
    }
  }
}

int main(int argc, char **argv) {
  double x;
  if(argc==1) { x = M_PI; } else { sscanf(argv[1], "%lf", &x); }
  assert(x>0); printf("%.15lf\n\n", x);
  all_best(x); printf("\n");
  return 0;
}

Examples

Here’s the output of this program, in about 0.003 seconds (i.e., it’s truly better than looping through all possible denominators!), line-wrapped for readability:

% ./a.out
3.141592653589793
 
3:  3/1
7:  22/7
15:  333/106
1:  355/113
292:  103993/33102
1:  104348/33215
1:  208341/66317
1:  312689/99532
2:  833719/265381
1:  1146408/364913
3:  4272943/1360120
1:  5419351/1725033
14:  80143857/25510582
 
13/4, 16/5, 19/6, 22/7, 179/57, 201/64, 223/71, 245/78, 267/85, 289/92, 311/99,
333/106, 355/113, 52163/16604, 52518/16717, 52873/16830, 53228/16943, 53583/17056,
53938/17169, 54293/17282, 54648/17395, 55003/17508, 55358/17621, 55713/17734,
56068/17847, 56423/17960, 56778/18073, 57133/18186, 57488/18299, 57843/18412,
58198/18525, 58553/18638, 58908/18751, 59263/18864, 59618/18977, 59973/19090,
60328/19203, 60683/19316, 61038/19429, 61393/19542, 61748/19655, 62103/19768,
62458/19881, 62813/19994, 63168/20107, 63523/20220, 63878/20333, 64233/20446,
64588/20559, 64943/20672, 65298/20785, 65653/20898, 66008/21011, 66363/21124,
66718/21237, 67073/21350, 67428/21463, 67783/21576, 68138/21689, 68493/21802,
68848/21915, 69203/22028, 69558/22141, 69913/22254, 70268/22367, 70623/22480,
70978/22593, 71333/22706, 71688/22819, 72043/22932, 72398/23045, 72753/23158,
73108/23271, 73463/23384, 73818/23497, 74173/23610, 74528/23723, 74883/23836,
75238/23949, 75593/24062, 75948/24175, 76303/24288, 76658/24401, 77013/24514,
77368/24627, 77723/24740, 78078/24853, 78433/24966, 78788/25079, 79143/25192,
79498/25305, 79853/25418, 80208/25531, 80563/25644, 80918/25757, 81273/25870,
81628/25983, 81983/26096, 82338/26209, 82693/26322, 83048/26435, 83403/26548,
83758/26661, 84113/26774, 84468/26887, 84823/27000, 85178/27113, 85533/27226,
85888/27339, 86243/27452, 86598/27565, 86953/27678, 87308/27791, 87663/27904,
88018/28017, 88373/28130, 88728/28243, 89083/28356, 89438/28469, 89793/28582,
90148/28695, 90503/28808, 90858/28921, 91213/29034, 91568/29147, 91923/29260,
92278/29373, 92633/29486, 92988/29599, 93343/29712, 93698/29825, 94053/29938,
94408/30051, 94763/30164, 95118/30277, 95473/30390, 95828/30503, 96183/30616,
96538/30729, 96893/30842, 97248/30955, 97603/31068, 97958/31181, 98313/31294,
98668/31407, 99023/31520, 99378/31633, 99733/31746, 100088/31859, 100443/31972,
100798/32085, 101153/32198, 101508/32311, 101863/32424, 102218/32537, 102573/32650,
102928/32763, 103283/32876, 103638/32989, 103993/33102, 104348/33215, 208341/66317,
312689/99532, 833719/265381, 1146408/364913, 3126535/995207,
4272943/1360120, 5419351/1725033, 42208400/13435351, 47627751/15160384,
53047102/16885417, 58466453/18610450, 63885804/20335483, 69305155/22060516,
74724506/23785549, 80143857/25510582, 

All these terms are correct, though if you increase MAX you start getting errors because of precision. I’m myself impressed with how many terms you get with only 13 convergents. (As you can see, there’s a small bug where it sometimes doesn’t print the very first “n/1″ approximation, or prints it incorrectly — I leave it for you to fix!)

You can try with √2, whose continued fraction is [1; 2, 2, 2, 2…]:

% ./a.out 1.41421356237309504880
1.414213562373095
 
1:  1/1
2:  3/2
2:  7/5
2:  17/12
2:  41/29
2:  99/70
2:  239/169
2:  577/408
2:  1393/985
2:  3363/2378
2:  8119/5741
2:  19601/13860
2:  47321/33461
 
3/2, 4/3, 7/5, 17/12, 24/17, 41/29, 99/70, 140/99, 239/169, 577/408, 816/577, 1393/985, 3363/2378, 4756/3363, 8119/5741, 19601/13860, 47321/33461, 

Or the golden ratio φ = (1+√5)/2 whose continued fraction is [1; 1, 1, 1, …]:

% ./a.out 1.61803398874989484820
1.618033988749895
 
1:  1/1
1:  2/1
1:  3/2
1:  5/3
1:  8/5
1:  13/8
1:  21/13
1:  34/21
1:  55/34
1:  89/55
1:  144/89
1:  233/144
1:  377/233
 
2/1, 3/2, 5/3, 8/5, 13/8, 21/13, 34/21, 55/34, 89/55, 144/89, 233/144, 377/233, 

(See the Fibonacci numbers? Here the convergents are all the approximants.)

Or with rational numbers like 4/3 = [1; 3]:

% ./a.out 1.33333333333333333333
1.333333333333333
 
1:  1/1
3:  4/3
 
3/2, 4/3, 

or 14/11 = [1; 3, 1, 2]:

% ./a.out 1.27272727272727272727
1.272727272727273
 
1:  1/1
3:  4/3
1:  5/4
2:  14/11
 
3/2, 4/3, 5/4, 9/7, 14/11, 

[2011-11-04: Fixed the bug and the description.]

About these ads

Written by S

Mon, 2011-01-10 at 02:28:23 +05:30

Posted in mathematics

15 Responses

Subscribe to comments with RSS.

  1. There’s something wrong with your algorithm for generating best rational approximations. It shows up, for example, in the output of the program when you give it the square root of two. Leaving aside the fact that the first result in the list of best rational approximations (2/1) is wrong, when we come to 10/7 there is an error, because 10/7 is further from sqrt(2) than 7/5 is. I believe that the problem is that you have not implemented the test to see whether a_k/2 is admissible as a result when a_k is divisible by 2. This test is described in the wikipedia article and amounts to comparing two continued fractions.

    I’m glad to see you covering this stuff.
    Keep up the good work. :-)

    [2011-11-04: This is now fixed — S.]

    Nontagonist

    Wed, 2011-02-23 at 13:24:57 +05:30

  2. Oh, and BTW, every 4th result after 10/7 in the best approximations for sqrt(2) is probably wrong too, for the same reason.

    Bye.

    [2011-11-04: This is now fixed. — S.]

    Nontagonist

    Wed, 2011-02-23 at 13:28:48 +05:30

    • Oops, you’re right! Very embarrassed. When writing the post I felt guilty about posting something without completely understanding its proof, and my fears have come true. :p

      Thanks very much for your comments. I’ll understand the theorem properly and update the post; for now I’ve updated it to point out the parts that have mistakes.

      [2011-11-04: This is now fixed. The theorem turned out to be simply to check whether the fraction is better than the previous one. — S.]

      S

      Wed, 2011-02-23 at 21:58:54 +05:30

      • The formulation of rule 3 at on Wikipedia and in Grafic Gems V seems unnecessary complicated.
        You just have to keep track of the differences
        d_k=|p_k/q_k-x| for admissible fractions,
        and scrap those pn/qn, which not give an improvement,
        i.e. for which the latest d_k<=|p_n/q_n-x|.
        Allan

        [2011-11-04: You’re right; thanks! This is now fixed — S.]

        Anonymous

        Tue, 2011-10-25 at 17:19:22 +05:30

        • I think you’re right (for some meaning of “think” as applies at 3am!). That is definitely much simpler.

          But a possible issue with that is the precision to which |p_k/q_k – x| can be computed, if we don’t have x in a convenient form. So I think their goal is to give the rule in terms of just the continued fraction.
          Meanwhile, the linked reference in Wikipedia, Graphics Gems 5 does give a simpler test on page 28 (the next page), if you generate the continued fraction using the a_k’s they use. The reference given is Concrete Mathematics, equations 6.131 and 6.135, which I plan to look at sometime when I’m more awake.

          Thanks for your comment!

          [2011-11-04: This is now fixed. The rule in terms of continued fractions is exactly equivalent to comparing the distance of x from the two, so what I said was stupid. — S.]

          S

          Tue, 2011-10-25 at 17:46:31 +05:30

          • Thank you for pointing this out.
            These precision problems is of course serious with the old compilers from 1980-2000;
            but seem obsolate to day with Mathematica and Maple.
            The condition in GEM V: Allow c_k/2 whenever d_(k-2)*a_(k-1)>d_(k-1)*a_k
            is correct; but the a_i are only explained for rationals.
            The refrences in Concrete Mathematics are only about some basic properties.
            Allan

            Anonymous

            Wed, 2011-10-26 at 19:42:19 +05:30

          • Are there better references?
            Do you know about generalizations to higher dimensions?
            Allan.

            Anonymous

            Wed, 2011-10-26 at 20:11:57 +05:30

            • What do you mean by higher dimensions? If you mean multidimensional continued fractions, I know nothing about them, sorry.

              I’m trying to find better references; do inform me if you find one. Meanwhile, the best I can find (and to which all other references seem to lead) is the classic book Continued Fractions by A. Khinchin. He defines the thing we want:
              Definition: A fraction \frac{a}{b} is a best approximation of the first kind of a real number x if for every fraction \frac{c}{d} with 1 \le d \le b and a/b \neq c/d, we have \left|x - \frac{c}{d}\right| > \left|x - \frac{a}{b}\right|
              And he proves that every best approximation is either a convergent on an intermediate fraction (semiconvergent).

              (A fraction \frac{a}{b} is a best approximation of the second kind if |dx - c| > |bx - a| for all such \frac{c}{d}, and these are precisely the convergents.)

              But surprisingly, he doesn’t characterize which semiconvergents (intermediate fractions) to take.

              Even the book by Rockett and Szüsz skips the issue, just saying “it is difficult to characterize the best approximations of the first kind”!

              It is all too clear that sometimes mathematicians pick whatever is convenient and has a better theory, rather than what is natural. :-)

              S

              Wed, 2011-10-26 at 21:36:45 +05:30

            • I found a better reference! See the paper
              “A more precise rounding algorithm for rational numbers”, by M. Thill (2008)
              which has this in its background section. It cites

              Perron O(1954) Die Lehre von den Kettenbrüchen. Bd. I Elementare Kettenbrüche. (German) 3te Aufl. B. G. Teubner Verlagsgesellschaft, Stuttgart

              I’ll update this comment sometime tomorrow or so, to say what it says.

              [2011-11-05: Posted as a separate comment below.]

              S

              Thu, 2011-10-27 at 01:18:14 +05:30

              • Thank you very much for these references.
                I agree in your viewpoint about Khinchin’s book,
                and have found Perron (1913). It is absollutely the best.
                I may try to simplify the stuff.

                Kind regards
                Allan

                Allan

                Tue, 2011-11-01 at 03:53:12 +05:30

            • From Thill’s paper (after renaming some variables):

              Theorem An intermediate fraction \displaystyle N(i,k) = \frac{i \cdot p_{k-1} + p_{k-2}}{i \cdot q_{k-1}+ q_{k-2}} (necessarily with k ≥ 1) is a best approximation of the first kind of x if and only if

              (i) [8, Sect. II.16.III p. 50] N (i, k) is nearer to x than is \displaystyle \frac{p_{k-1}}{q_{k-1}} which is equivalent to

              (ii) [8, Theorem 2.23 p. 51] either 2i > a_k, or else 2i = a_k as well as
              [a_k , a_{k-1} , . . . , a_2 , a_1 ] > [ a_k , a_{k+1} , a_{k+2} , . . . ].

              The reference [8] is Perron. The form (i) is what you said, and (ii) is what is given on Wikipedia. I have edited Wikipedia to include the simpler version, with a reference to Thill’s paper.

              S

              Sat, 2011-11-05 at 23:53:42 +05:30

              • Thank you for refering from Thill’s paper, which not seems to have free access.
                Is it about a numerical precision problem, or somthing about rational approximation?
                On wiki seems needed that a_k should be even.
                Allan

                Allan

                Sun, 2011-11-06 at 06:09:41 +05:30

                • I have emailed the paper to you. The paper is about something else, but it quotes the above theorem (which is what we want), as background, citing Perron’s book. I don’t have access to the 1954 version, but in the 1913 edition available online, see section 2.16 starting here on page 55. On page 60 is the form (ii) of the rule in terms of continued fractions.

                  Here also we need that a_k is even as on Wikipedia (just now I renamed “n” to “k” in the previous comment), because the test only becomes necessary when 2i =a_k. When a_k is odd, the condition is automatically satisfied and doesn’t need to be tested (though testing it does no harm).

                  S

                  Sun, 2011-11-06 at 09:13:11 +05:30

                  • Thank you very much for the references.
                    I have collected my impressions in the link below,
                    where Rule 3 is treated in Theorem 6:

                    http://lanco.host22.com/index.php?p=1_7

                    On wiki under the headline : Best rational approximations.
                    The “half rule” is that the ….
                    could perhaps be changed to
                    If a_k is even, the “half rule” is that the ….

                    Example: a=11/16=[0,1,2,5]
                    b=[0,1,2,3]=7/10 OK with b-a=1/80=0.0125…
                    c=[0,1,2,5/2]=12/17 not OK, as c-a=5/272=0.0183…,
                    because b is better approximation than c, and with
                    smaller denominator.
                    Thus c is not a a best approximation of first kind
                    Allan

                    Allan

                    Thu, 2011-11-10 at 23:29:44 +05:30

                    • Your PDF looks interesting, thanks.

                      As for Wikipedia, it was already mentioned earlier (the third item in the list) that the rule is to be checked only when ak is even. (When ak is odd, ak/2 is not an integer, so the rule doesn’t make sense for simple continued fractions anyway.) Anyway, I’ve re-edited the Wikipedia article to clarify as you suggested.

                      S

                      Sat, 2011-11-12 at 20:07:00 +05:30


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 78 other followers