Version 17 of Notes on continued fractions

Updated 2002-01-16 20:22:42

MS: These are some notes on my playing around with Fraction Math and the reference to the nice notes on Continued Fractions at http://www.inwap.com/pdp10/hbaker/hakmem/cf.html [L1 ].

I have started playing with these things and resisted (for now) looking for other sources. Some of the "errors" I found in the reference may be due to the informal presentation of results there, which might be inaccurate or misunderstood by me: caveat emptor.

DEFINITIONS: a rational approximation p/q to a real number x is "best" iff, for every integer r and s,

  (s <= q) ==> (|x - p/q| <= |x - r/s|)

It is "best on its side" if

  ((s <= q) & (sgn(x-p/q) = sgn(x-r/s))) ==> (|x - p/q| <= |x - r/s|)

i.e, if no other fraction on the same side of x with a lesser denominator comes closer.

NOTATION: let us identify a (positive) real number x with its regular continued fraction representation

  x = {x[0] x[1] x[2] x[3] ...} 

Define the truncation of x after (n+1) terms

  a(x,n) = {x[0] x[1] x[2] x[3] ... x[n]} 

and let

  b(x,n,i) = {x[0] x[1] x[2] x[3] ... x[n-1] i}

be a(x,n) with the last element replaced by 0 < i <= x[n]. Note that

  a(x,n) = b(x,n,x[n])

NOTE: [quotient_rep] in Fraction Math computes the highest-order truncation requiring no integers larger than maxint in the rational representation p/q.


Item 101A (3) clearly says (AFAIU, not all claims reproduced here):

   A - a(x,n) is "best"
   B - b(x,n,i) is "best" if i>1
   C - b(x,n,1) is never "best" if (x[n] != 1)

Let me provide counterexamples to both B and C, thus showing that they are not true.

The example provided in the text actually contains counterexamples to B. Let x = pi = {3 7 15 1 292 ...}

   . b(x,0,2) = 2/1 = {2} is not best (3/1 is better)
   . b(x,1,2) = 7/2 = {3 2} and b(x,1,3) = 10/3 = {3 3} are not best (3/1 is better) 
   . b(x,2,2) = 47/15 = {3 7 2} is not best (22/7 is better)

For another counterexample to B, easier to follow by hand, consider

  x = 0.51 = {0 1 1 24 2}

The number

  b(x,3,4) = 5/9 = {0 1 1 4} = 0.555...

is not best, as 1/2 = {0 1 1} = {0 2} is better.

For a counterexample to C, consider x = 7/10 = {0 1 2 3}; now

  b(x,2,1) = 1/2 = {0 1 1} = {0 2} 

is best.


So, in light of these counterexamples, one proposition I can hope could be true is:

   D - The set of "best on its side" approximations to 0<x<1 coincides with 
       the set of numbers of the form b(x,n,i)

Remark that a simple corollary would be:

   E - if r is a "best" approximation to x, then r = b(x,n,i) for some (n,i)

as a "best" approx has to be "best on its side".

I think (D) can be proved from ITEM 101C in [L2 ] (which may be true AFAIK). I'm working on the details, maybe a restriction (i!=1) will be necessary. After a while I'll stop playing and start reading ...


The example I provided in Fraction Math showing that quotient_rep does not always provide the best approximation is

   quotient_rep 3.1416305 500 --> 355/113

Now, 3.1416305 = {3 7 16 2 ...},and quotient_rep produced

   355/113 = {3 7 16} = a(3.1416305,2)

But the fraction

   377/120 = {3 7 16 1} = b(3.1416305,3,1)

is closer - a new counterexample to C.

Remark that the next truncation involves integers that are too large:

   a(3.1416305,3) = {3 7 16 2} = 732/233

Arjen Markus Some care must be taken, as I discovered myself. I tried to deduce the continued fraction for sqrt(N^2+1)+N (forms such as sqrt(2)+1, sqrt(5)+2, ...):

   (sqrt(N^2+1)+N) * (sqrt(N^2+1)-N = 1

So:

                        1               
   sqrt(N^2+1)+N = -------------
                   sqrt(N^2+1)-N

                             1
                 = ---------------------
                   -2N + (sqrt(N^2+1)+N)

                               1
                 = -------------------------
                                  1
                   -2N + ---------------------
                         -2N + (sqrt(N^2+1)+N)

                 = ...

Now, the approximations are (N=1):

    -0.5, -0.4, ...

How will this ever end up in 2.4141...?

It took me the better part of an evening to realise that I had made a stupid, though understandable, mistake: the continued fraction does not get smaller! Hence my approximations are worthless. They in fact will give 1-sqrt(2).

What is interesting is that these continued fractions have a natural representation in Tcl, namely lists. Try doing that with plain C!

KBK - OK, ok, here's a solution for the continued fraction representation of a quadratic surd. Let

   x = ( sqrt(D) - U ) / V, and D not be a perfect square.

# We need a little auxiliary procedure to get the greatest common divisor of p and q

 proc gcd { p q } {
         while { $q != 0 } {
             foreach { p q } [list $q [expr { $p % $q }]] break
         }
     return $p
 }

# Expand the quadratic surd x = ( -U + sqrt( D ) ) / V as a continued fraction truncating after n partial quotients

 proc cfsurd { U D V {n 42} } {

     # Renormalize so that V divides D-U*U evenly

     set j [gcd $V [expr { $D - $U * $U }]]
     set l [expr { $V / $j }]
     set D [expr { $l * $l * $D }]
     set U [expr { $l * $U }]
     set V [expr { $l * $V }]

     # Start the iteration by finding the integer part of the
     # surd.

     set A [expr { int( ( sqrt( $D ) - $U ) / $V ) }]
     set result [list $A]
     set f [expr { int( sqrt( $D ) ) }]
     set U [expr { $U + $A * $V }]

     # Expand on results by iteration

     while { [llength $result] < $n } {
         set V [expr { ( $D - $U * $U ) / $V }]
         if { $V > 0 } {
             set A [expr { int( ( $f + $U ) / $V ) }]
         } else {
             set A [expr { int( ( $f + 1 + $U ) / $V ) }]
         }
         set U [expr { $A * $V - $U }]
         lappend result $A
     }
     return $result
 }

# As a special case, expand a square root as a continued fraction to n partial quotients.

 proc cfsqrt { x { n 42 } } {
     return [cfsurd 0 $x 1 $n]
 }

# And, if you like, the square root of p/q:

 proc cfsqrtrat { p q { n 42 } } {
     return [cfsurd 0 [expr $p * $q] $q $n]
 }

# Some usage examples:

 puts [cfsqrt 2]       ;# sqrt(2) = 1 + / 2, 2, 2, 2, ... /
 puts [cfsqrt 3]       ;# sqrt(3) = 1 + / 2, 1, 2, 1, ... /
 puts [cfsurd -1 5 2]  ;# (1 + sqrt(5)) / 2 = 1 + / 1, 1, 1, 1, ... /
 puts [cfsqrtrat 8 29] ;# sqrt(8/29) =
                        #   / 1 ; 1, 9, 2, 2, 3, 2, 2, 9, 1, 2, ... /
                        #   repeating with period 10
 puts [cfsurd -1 2 1]  ;# 1 + sqrt(2) (Arjen's example)
                        #       = 2 + / 2, 2, 2, 2, 2, ... /

Some code to deal with continued fractions

a lot of this is ripped off from KBK's code in Fraction Math. Some testing of inputs would also be needed.

   #
   # proc cont2frac
   #
   # Returns the fraction corresponding to the list
   # of denominators in a (truncated) regular continued
   # fraction. There is no check for overflow.
   #

   proc cont2frac {lst} {
       foreach {p q p0 q0} {1 0 0 1} break
       foreach a $lst {
           foreach {p p0} [list [expr {$a*$p+$p0}] $p] break 
           foreach {q q0} [list [expr {$a*$q+$q0}] $q] break 
      }
      list $p $q
   }


   #
   # proc frac2cont
   #
   # Returns list of denominators in the regular continued
   # expansion of frac. Cannot overflow.
   #

   proc frac2cont {frac} {
       foreach {p0 q0} {1 0} break
       foreach {p q} $frac break
       set clist {}

       while {$q} {
           set a [expr {$p/$q}]
           lappend clist $a
           foreach {p q} [list $q [expr {$p - $q * $a}]] break
       }
       set clist
   }

   # 
   # proc num2cont (new version, approximating)
   #
   # Returns a list of two lists:
   #   1. the list of denominators in the best regular continued
   #      fraction approx of num (assuming D above!) with both
   #      numerator and denominator <= maxint.
   #   2. the representation of the above list as a fraction
   #

   proc num2cont {num { maxint 2147483647 } } {
       #
       # Do the first iteration outside the loop: this handles
       # the differences when $num is (<0 | >=0) or (<=1 | >1),
       # always sends a number 0<=x<1 to the loop, and insures that 
       # a *regular* continued fraction is produced in every case:
       # the first element carries the sign, every other element is
       # strictly positive.
       #

       set a [expr {int(floor($num))}]
       set x [expr {$num - $a}]
       set clist [list $a]
       foreach {p q p0 q0} [list $a 1 1 0] break

       #
       # Instead of testing at each step both numerator and denominator,
       # link the parameters now to the coeffs of the larger one, and
       # only test this one.
       #

       if {($a == 0) || ($a == -1)} {
           upvar 0 q0 r0 q r
       } else {
           upvar 0 p0 r0 p r
           if {$a < 0} {
               # Avoid absolute values for amax: set a positive
               # numerator and negative denominator.
               foreach {p q p0 q0} [list [expr {-$a}] -1 1 0] break
           }
       }

       set dist0 0
       set lim [expr {1.0/$maxint}]

       while {$x > $lim} {
           set x [expr {1.0 / $x}]
           set a [expr {int($x)}]
           set amax [expr {double($maxint-$r0)/$r}]

           if {$a > $amax} {
               set a [expr {int($amax)}]
               if {!$a} break
               set dist0 [expr {abs($num - double($p0)/$q0)}]
           }

           foreach {p q p0 q0} \
                   [list [expr {$a * $p + $p0}] [expr {$a * $q + $q0}] $p $q] \
                   break

           if {$dist0} {
               set dist1 [expr {abs($num - double($p)/$q)}]
               if {$dist1 < $dist0} {
                   lappend clist $a
               } else {
                   foreach {p q} [list $p0 $q0] break
               }
               break
           }
           lappend clist $a
           set x [expr {$x - $a}]
       }

       list $clist  [list $p $q]
   }

   # 
   # proc num2cont0 (original version, truncating)
   #
   # Returns a list of two lists:
   #   1. the list of denominators in the longest truncated
   #      regular continued fraction expansion of num with both
   #      numerator and denominator <= maxint.
   #   2. the representation of the above list as a fraction
   #
   # Remark that "quotient_rep $num $maxint]" is equivalent to
   # "lindex [num2cont0 $num $maxint] 1"
   #

    proc num2cont0 {num { maxint 2147483647 } } {
        foreach {p q p0 q0} {1 0 0 1} break
        set clist {}

        while {1} {
           set a [expr {int($num)}]
           set fract [expr {$num - $a}]

           if {(1.0 * $a * $p + $p0 > $maxint) \
                   || (1.0 * $a * $q + $q0 > $maxint)} {
               break
           }
           lappend clist $a
           foreach {p p0} [list [expr {$a * $p + $p0}] $p] break 
           foreach {q q0} [list [expr {$a * $q + $q0}] $q] break 

           if {abs($fract * $maxint) < 1} break
           set num [expr {1.0 / $fract}]
        }
        list $clist [list $p $q]
    }