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/233Arjen 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 } {
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} {
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]
}
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.
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.
[ Category Mathematics ] [L3 ] [L4 ] [L5 ] [L6 ] [L7 ] [L8 ] [L9 ] [L10 ] [L11 ] [L12 ] [L13 ] [L14 ] [L15 ] [L16 ] [L17 ] [L18 ] [L19 ] [L20 ] [L21 ] [L22 ] [L23 ] [L24 ] [L25 ] [L26 ] [L27 ] [L28 ] [L29 ] [L30 ] [L31 ] [L32 ] [L33 ] [L34 ] [L35 ] [L36 ] [L37 ] [L38 ] [L39 ] [L40 ] [L41 ] [L42 ] [L43 ] [L44 ] [L45 ] [L46 ] [L47 ] [L48 ] [L49 ] [L50 ] [L51 ] [L52 ] [L53 ] [L54 ] [L55 ] [L56 ] [L57 ] [L58 ] [L59 ] [L60 ] [L61 ] [L62 ] [L63 ] [L64 ] [L65 ] [L66 ] [L67 ] [L68 ] [L69 ] [L70 ] [L71 ] [L72 ] [L73 ] [L74 ] [L75 ] [L76 ] [L77 ] [L78 ] [L79 ] [L80 ] [L81 ] [L82 ] [L83 ] [L84 ] [L85 ] [L86 ] [L87 ] [L88 ] [L89 ] [L90 ] [L91 ] [L92 ] [L93 ] [L94 ] [L95 ] [L96 ] [L97 ] [L98 ] [L99 ] [L100 ] [L101 ] [L102 ] [L103 ] [L104 ] [L105 ] [L106 ] [L107 ] [L108 ] [L109 ] [L110 ] [L111 ] [L112 ] [L113 ]
[L114 ] [L115 ] [L116 ] [L117 ] [L118 ] [L119 ] [L120 ] [L121 ] [L122 ] [L123 ] [L124 ] [L125 ] [L126 ] [L127 ] [L128 ] [L129 ] [L130 ] [L131 ] [L132 ] [L133 ] [L134 ] [L135 ] [L136 ] [L137 ] [L138 ] [L139 ] [L140 ] [L141 ] [L142 ] [L143 ] [L144 ] [L145 ] [L146 ] [L147 ] [L148 ] [L149 ] [L150 ] [L151 ] [L152 ] [L153 ] [L154 ] [L155 ] [L156 ] [L157 ] [L158 ] [L159 ] [L160 ] [L161 ] [L162 ] [L163 ] [L164 ] [L165 ] [L166 ] [L167 ] [L168 ] [L169 ] [L170 ] [L171 ] [L172 ] [L173 ] [L174 ] [L175 ] [L176 ] [L177 ] [L178 ] [L179 ] [L180 ] [L181 ] [L182 ] [L183 ] [L184 ] [L185 ] [L186 ] [L187 ] [L188 ] [L189 ] [L190 ] [L191 ] [L192 ] [L193 ] [L194 ] [L195 ] [L196 ] [L197 ] [L198 ] [L199 ] [L200 ] [L201 ] [L202 ] [L203 ] [L204 ] [L205 ] [L206 ] [L207 ] [L208 ] [L209 ] [L210 ] [L211 ] [L212 ] [L213 ] [L214 ] [L215 ] [L216 ] [L217 ] [L218 ] [L219 ] [L220 ] [L221 ] [L222 ] [L223 ] [L224 ] [L225 ] [L226 ] [L227 ] [L228 ] [L229 ] [L230 ] [L231 ] [L232 ] [L233 ] [L234 ] [L235 ] [L236 ] [L237 ] [L238 ] [L239 ] [L240 ] [L241 ] [L242 ] [L243 ] [L244 ] [L245 ] [L246 ] [L247 ] [L248 ] [L249 ] [L250 ] [L251 ] [L252 ] [L253 ] [L254 ] [L255 ] [L256 ] [L257 ] [L258 ] [L259 ] [L260 ] [L261 ] [L262 ] [L263 ] [L264 ] [L265 ] [L266 ] [L267 ] [L268 ] [L269 ] [L270 ] [L271 ] [L272 ] [L273 ] [L274 ] [L275 ] [L276 ] [L277 ] [L278 ] [L279 ] [L280 ] [L281 ] [L282 ] [L283 ] [L284 ] [L285 ] [L286 ] [L287 ] [L288 ] [L289 ] [L290 ] [L291 ] [L292 ] [L293 ] [L294 ] [L295 ] [L296 ] [L297 ] [L298 ] [L299 ] [L300 ] [L301 ] [L302 ] [L303 ] [L304 ] [L305 ] [L306 ] [L307 ] [L308 ] [L309 ] [L310 ] [L311 ] [L312 ] [L313 ] [L314 ] [L315 ] [L316 ] [L317 ] [L318 ] [L319 ] [L320 ] [L321 ] [L322 ] [L323 ] [L324 ] [L325 ] [L326 ] [L327 ] [L328 ] [L329 ] [L330 ] [L331 ] [L332 ] [L333 ] [L334 ] [L335 ] [L336 ] [L337 ] [L338 ] [L339 ] [L340 ] [L341 ] [L342 ] [L343 ] [L344 ] [L345 ] [L346 ] [L347 ] [L348 ] [L349 ] [L350 ] [L351 ] [L352 ] [L353 ] [L354 ] [L355 ] [L356 ] [L357 ] [L358 ] [L359 ] [L360 ] [L361 ] [L362 ] [L363 ] [L364 ] [L365 ] [L366 ] [L367 ] [L368 ] [L369 ] [L370 ] [L371 ] [L372 ] [L373 ] [L374 ] [L375 ] [L376 ] [L377 ] [L378 ] [L379 ] [L380 ] [L381 ] [L382 ] [L383 ] [L384 ] [L385 ] [L386 ] [L387 ] [L388 ]
[L389 ] [L390 ] [L391 ] [L392 ] [L393 ] [L394 ] [L395 ] [L396 ] [L397 ] [L398 ] [L399 ] [L400 ] [L401 ] [L402 ] [L403 ] [L404 ] [L405 ] [L406 ] [L407 ] [L408 ] [L409 ] [L410 ] [L411 ] [L412 ] [L413 ] [L414 ] [L415 ] [L416 ] [L417 ] [L418 ] [L419 ] [L420 ] [L421 ] [L422 ] [L423 ] [L424 ] [L425 ] [L426 ] [L427 ] [L428 ] [L429 ] [L430 ] [L431 ] [L432 ] [L433 ] [L434 ] [L435 ] [L436 ] [L437 ] [L438 ] [L439 ] [L440 ] [L441 ] [L442 ] [L443 ] [L444 ] [L445 ] [L446 ] [L447 ] [L448 ] [L449 ] [L450 ] [L451 ] [L452 ] [L453 ] [L454 ] [L455 ] [L456 ] [L457 ] [L458 ] [L459 ] [L460 ] [L461 ] [L462 ] [L463 ] [L464 ] [L465 ] [L466 ] [L467 ] [L468 ] [L469 ] [L470 ] [L471 ] [L472 ] [L473 ] [L474 ] [L475 ] [L476 ] [L477 ] [L478 ] [L479 ] [L480 ] [L481 ] [L482 ] [L483 ] [L484 ] [L485 ] [L486 ] [L487 ] [L488 ] [L489 ] [L490 ] [L491 ] [L492 ] [L493 ] [L494 ] [L495 ] [L496 ] [L497 ] [L498 ] [L499 ] [L500 ] [L501 ] [L502 ] [L503 ] [L504 ] [L505 ] [L506 ]