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 } { lassign [list $q [expr { $p % $q }]] p q } 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 } { break ;# D was really the square of an integer } elseif { $V > 0 } { set A [expr { ( $f + $U ) / $V }] } else { set A [expr { ( $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} { lassign {1 0 0 1} p q p0 q0 foreach a $lst { lassign [list [expr {$a*$p+$p0}] $p] p p0 lassign [list [expr {$a*$q+$q0}] $q] q q0 } list $p $q } # # proc frac2cont # # Returns list of denominators in the regular continued # expansion of frac. Cannot overflow. # proc frac2cont {frac} { lassign {1 0} p0 q0 lassign $frac p q set clist {} while {$q} { set a [expr {$p/$q}] lappend clist $a lassign [list $q [expr {$p - $q * $a}]] p q } 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] lassign [list $a 1 1 0] p q p0 q0 # # 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. lassign [list [expr {-$a}] -1 1 0] p q p0 q0 } } 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)}] } lassign [list [expr {$a * $p + $p0}] [expr {$a * $q + $q0}] $p $q] p q p0 q0 if {$dist0} { set dist1 [expr {abs($num - double($p)/$q)}] if {$dist1 < $dist0} { lappend clist $a } else { lassign [list $p0 $q0] p q } 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 } } { lassign {1 0 0 1} p q p0 q0 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 lassign [list [expr {$a * $p + $p0}] $p] p p0 lassign [list [expr {$a * $q + $q0}] $q] q q0 if {abs($fract * $maxint) < 1} { break } set num [expr {1.0 / $fract}] } list $clist [list $p $q] }

AM (19 march 2004) Here is a small script to calculate the value of a continued fraction:

# contfrac.tcl -- # Compute continued fractions like: # a = 1 + 1 # ------------------ # 2 + 1 # ------------- # 3 + 1 # -------- # 4 + 1 # --- # 5 + ... # # With the script below this can be calculated (approximated) by: # set a [// 1 2 3 4 5 ..] # // -- # Compute the numerical value of a continued fraction # # Arguments: # values List of coefficients # Result: # Numerical value # proc // {args} { if { [llength $args] == 1 } { set values [lindex $args 0] } else { set values $args } # # First reverse the list # set rvalues {} foreach v [lrange $values 1 end] { set rvalues [concat $v $rvalues] } # # Then do the computation # set result 0.0 foreach r $rvalues { set result [expr {1.0/($r+$result)}] } expr {[lindex $values 0] + $result} } # main -- # Simple tests # catch {console show} puts "1/1: [// 1]" puts "1/(1+1/2)=2/3: [// {1 2}]" puts "1/(1+1/(2+1/3))=7/10: [// {1 2 3}]" puts "Alternative call: \[// 1 2 3\]: [// 1 2 3]" puts "Successive approximations to sqrt(2):" set values 1 for {set i 0} {$i < 20} {incr i} { lappend values 2 puts "[// $values] / sqrt(2) = [expr {[// $values]/sqrt(2.0)}] ([llength $values] terms)" } puts "Successive approximations to e = exp(1):" set values {2} for {set i 2} {$i < 20} {incr i 2} { lappend values 1 $i 1 puts "[expr {([// $values])}] / exp(1) = [expr {[// $values]/exp(1.0)}] \ ([llength $values] terms)" } puts "Successive approximations:" set values {} for {set i 1} {$i < 10} {incr i} { lappend values $i puts "[set r [// $values]] - $values" }

GWM Here are my recursive evaluators for Continued Fraction series - it is quite short since terms are automatically evaluated in reverse order.

proc contfrac { ai } { ;# CFs are a0+1/(a1+1/(a2+1/(a3...) # ai is a list; called recursively set s [lindex $ai 0] ;# a0 if [llength $ai]>1 { set s [expr $s+1.0/[contfrac [lrange $ai 1 end] ] ] } return $s } # evaluation of continued fractions- show sum of each approximation term proc successive {label ai} { ;# show successive approximations for {set i 1} { $i < [llength $ai] } { incr i } { puts "$label approximation $i [contfrac [lrange $ai 0 $i] ]" } return [contfrac $ai] } proc findContFrac {value nterms} { ;# find the nterms approximation series for real value #see http://en.wikipedia.org/wiki/Continued_fraction#Calculating_continued_fraction_representations list ai lappend ai [expr int($value)] ;# first part is the integer part of the value if {[expr $nterms > 0] && [expr fmod($value,1.0)] } {;# next term in approx is 1/ fractional part of the value append ai " [findContFrac [expr 1.0/fmod($value,1.0)] $nterms-1 ]" } return $ai } # sample fractions series are at http://www.research.att.com/~njas/sequences/Sindx_Con.html#confC puts "root(2) [contfrac {1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2}]" puts "root(3) [contfrac {1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 }]" puts "Pi [contfrac {3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2}]" puts "e [contfrac {2 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12 1 1}]" # show successive approximations from the truncated series puts "golden ratio [successive \"grat\" {1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 } ]" puts "Pi [successive Pi {3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2}]" # test the findContFrac by feeding into contfrac: puts "Sum [contfrac [findContFrac [expr sqrt(2)] 12] ]" # find continued fractions for each square root up to 100 for {set i 2} { $i<100} {incr i} { puts "sq root($i) = [findContFrac [expr sqrt($i)] 12]" }

Lars H: Sorry to be negative, GWM, but those procs are *terribly* bad.

- The expr and if expressions aren't braced.
- Both procedures are recursive, for a problem where there are perfectly good iterative solutions of the same length. (In the case of
**findContFrac**the iterative variant is even much more obvious.) - In
**findContFrac**you forget to evaluate an expression in the recursive call and there are exprs inside expressions. - Again in
**findContFrac**, you seem to think list declares a variable as a list and mixes append and lappend to build up a list. - Both
**contfrac**and**findContFrac**are asymptotically quadratic, for a problem that can be done in linear time.

Note that I don't claim that the procs don't work; I only say that they're doing things in a *very* far from optimal way.