## Notes on continued fractions

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 [1 ].

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 x x x ...} `

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

`  a(x,n) = {x x x x ... x[n]} `

and let

`  b(x,n,i) = {x x x x ... 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 [2 ] (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 } {
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.

• 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.

 Category Mathematics