[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 [http://www.inwap.com/pdp10/hbaker/hakmem/cf.html]. 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 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] }