[Tutor] Stern-Brocot tree: more fun math
Pijus Virketis
virketis@fas.harvard.edu
Tue, 8 Jan 2002 16:40:08 -0500
This is a multi-part message in MIME format.
------=_NextPart_000_00F6_01C19863.225ED8B0
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
Dear all,
inspired by the thread about Continued Fractions and armed with the =
reference to "Concrete Mathematics", I've written a few functions to =
play around with another interesting mathematical construct, a =
Stern-Brocot tree. For an introduction to SB trees, see =
http://www.cut-the-knot.com/ctk/SB_tree.html. It can be show (by me, =
even :)) that this tree is a representation of the Rational numbers =
(i.e. a number system). Its nodes are only reduced fractions, and they =
are ordered, so it is easy to make algorithms that traverse the tree. It =
also suggests an interesting notation, because any fraction can be =
unambiguously identified by listing how many left and right turns are =
needed from the first node to get to it. So, 5/7 is L, R, R, L. So, =
here's the code:
# define the L and R matrices
L =3D [[1,1],[0,1]]
R =3D [[1,0],[1,1]]
# some helper functions
def median(S):
"""
Returns the median fraction from a State matrix
"""
return [S[1][0]+S[1][1],S[0][0]+S[0][1]]
=20
def matrix_mult(A,B):
"""
Multiplies two 2x2 matrixes
"""
return ([[A[0][0]*B[0][0]+A[0][1]*B[1][0], =
A[0][0]*B[0][1]+A[0][1]*B[1][1]],
[A[1][0]*B[0][0]+A[1][1]*B[1][0], =
A[1][0]*B[0][1]+A[1][1]*B[1][1]]])
=20
def frac2path(m,n):
"""
Returns the LR address of reduced fraction m/n
"""
path =3D [] =20
S =3D [[1,0],[0,1]] # we start =
at the Identity node
frac =3D float(m)/float(n) # this is =
the fraction we wish to end up at
while frac !=3D float(median(S)[0])/float(median(S)[1]): # while not =
at the correct node
if frac < float(median(S)[0])/float(median(S)[1]): # either turn =
left
path.append("L") # and output =
"L"
S =3D matrix_mult(S,L) # update =
the State Matrix
else: # or turn =
right
path.append("R") # output "R"=20
S =3D matrix_mult(S,R) # update =
the State Matrix
return path
def path2frac(path_list):
"""
Returns the fraction corresponding to given ["L","R"] address
"""
S =3D [[1,0],[0,1]] # we start with Identity matrix
for turn in path_list: # for each turn in the path list
if turn =3D=3D "L": # if we turn left
S =3D matrix_mult(S,L) # multiply the State matrix by the =
Left matrix =20
else: # if we turn right =20
S =3D matrix_mult(S,R) # multiply the State matrix by the =
Right matrix
return median(S) # return the resulting fraction =
(median of State matrix)
=20
We can use this to estimate irrational numbers to arbitrary degree of =
precision. Just as Kirby did in his post a month ago, I am going to take =
a short at root 3. From highschool I remember that its something like =
1.73 ... Let's plug 173/100 into our function that gives back the path =
to any fraction:
>>> sternbrocot.frac2path(173,100)
['R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'L', 'R', 'L']
Is there a pattern here? It seems that it might be R, L, R, R, L, R, R, =
L, ... Irrationals always make patterns in SB representation. So, let's =
make a longer list and try it to see if our conjecture was good:
>>> root_list =3D ['R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', =
'R', 'L']
>>> sternbrocot.path2frac(root_list)
[121635, 70226]
>>> import math
>>> math.sqrt(3) - 121635./70226
1.7560397580496101e-010
So, our estimation is pretty close to the one obtained from math.sqrt(). =
Not bad. We can try it a few more steps further:
>>> root_list =3D ['R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', =
'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', =
'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L']
>>> sternbrocot.path2frac(root_list)
[17083879020L, 9863382151L]
>>> math.sqrt(3) - 17083879020./9863382151
0.0
That's a pretty good estimate of the square root of three. The next =
thing to do is to rewrite all of this as a generator, just to be trendy. =
:)
-PV=20
------=_NextPart_000_00F6_01C19863.225ED8B0
Content-Type: text/html;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=3DContent-Type content=3D"text/html; =
charset=3Diso-8859-1">
<META content=3D"MSHTML 6.00.2712.300" name=3DGENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=3D#ffffff>
<DIV><FONT face=3DArial size=3D2>Dear all,</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>inspired by the thread about Continued =
Fractions=20
and armed with the reference to "Concrete Mathematics", I've written a =
few=20
functions to play around with another interesting mathematical =
construct, a=20
Stern-Brocot tree. For an introduction to SB trees, see <A=20
href=3D"http://www.cut-the-knot.com/ctk/SB_tree.html">http://www.cut-the-=
knot.com/ctk/SB_tree.html</A>.=20
It can be show (by me, even :)) that this tree is a representation of =
the=20
Rational numbers (i.e. a number system). Its nodes are only reduced =
fractions, and they are ordered, so it is easy to make algorithms that =
traverse=20
the tree. It also suggests an interesting notation, because any fraction =
can be=20
unambiguously identified by listing how many left and right turns are =
needed=20
from the first node to get to it. So, 5/7 is L, R, R, L. So, here's the=20
code:</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV> <FONT face=3DArial size=3D2># define the L and R =
matrices<BR>L =3D=20
[[1,1],[0,1]]<BR>R =3D [[1,0],[1,1]]</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2># some helper functions<BR>def=20
median(S):<BR> """<BR> Returns the =
median=20
fraction from a State matrix<BR> =
"""<BR> =20
return [S[1][0]+S[1][1],S[0][0]+S[0][1]]<BR> <BR>def=20
matrix_mult(A,B):<BR> """<BR> =
Multiplies two=20
2x2 matrixes<BR> """<BR> return=20
([[A[0][0]*B[0][0]+A[0][1]*B[1][0],=20
A[0][0]*B[0][1]+A[0][1]*B[1][1]],<BR> =
=20
[A[1][0]*B[0][0]+A[1][1]*B[1][0],=20
A[1][0]*B[0][1]+A[1][1]*B[1][1]]])<BR> <BR>def=20
frac2path(m,n):<BR> """<BR> Returns =
the LR=20
address of reduced fraction m/n<BR> =
"""<BR> =20
path =3D [] =20
<BR> S =3D=20
[[1,0],[0,1]] =
&=
nbsp; &n=
bsp; =20
# we start at the Identity node<BR> frac =3D=20
float(m)/float(n) &n=
bsp; &nb=
sp; =20
# this is the fraction we wish to end up at<BR> while =
frac !=3D=20
float(median(S)[0])/float(median(S)[1]): # while not at the correct=20
node<BR> if frac <=20
float(median(S)[0])/float(median(S)[1]): # either turn=20
left<BR>  =
;=20
path.append("L") &nb=
sp; &nbs=
p; =20
# and output=20
"L"<BR> =
S =3D=20
matrix_mult(S,L) &nb=
sp; &nbs=
p; =20
# update the State Matrix<BR> =20
else: &n=
bsp; &nb=
sp; &nbs=
p; =20
# or turn=20
right<BR> &nbs=
p;=20
path.append("R") &nb=
sp; &nbs=
p; =20
# output "R"=20
<BR> S =
=3D=20
matrix_mult(S,R) &nb=
sp; &nbs=
p; =20
# update the State Matrix<BR> return path</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>def =
path2frac(path_list):<BR> =20
"""<BR> Returns the fraction corresponding to given =
["L","R"]=20
address<BR> """<BR> S =3D=20
[[1,0],[0,1]] =
=20
# we start with Identity matrix<BR> for turn in=20
path_list: # for =
each turn=20
in the path list<BR> if turn =
=3D=3D=20
"L": &nb=
sp; #=20
if we turn=20
left<BR>  =
; S =3D=20
matrix_mult(S,L) # multiply the State matrix by the =
Left=20
matrix =20
<BR> =20
else: &n=
bsp; =20
# if we turn right =20
<BR> S =
=3D=20
matrix_mult(S,R) # multiply the State matrix by the =
Right=20
matrix<BR> return=20
median(S) &nbs=
p; =20
# return the resulting fraction (median of State =
matrix)<BR> =20
<BR>We can use this to estimate irrational numbers to arbitrary degree =
of=20
precision. Just as Kirby did in his post a month ago, I am going to take =
a short=20
at root 3. From highschool I remember that its something like 1.73 ... =
Let's=20
plug 173/100 into our function that gives back the path to any=20
fraction:</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>>>>=20
sternbrocot.frac2path(173,100)<BR>['R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L', 'L',=20
'R', 'L']</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>Is there a pattern here? It seems that =
it might be=20
R, L, R, R, L, R, R, L, ... Irrationals always make patterns in SB=20
representation. So, let's make a longer list and try it to see if our =
conjecture=20
was good:</DIV>
<DIV><BR>>>> root_list =3D ['R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L', 'R',=20
'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', =
'R', 'R',=20
'L']<BR>>>> sternbrocot.path2frac(root_list)<BR>[121635,=20
70226]<BR>>>> import math<BR>>>> math.sqrt(3) -=20
121635./70226<BR>1.7560397580496101e-010</DIV>
<DIV> </DIV>
<DIV>So, our estimation is pretty close to the one obtained from =
math.sqrt().=20
Not bad. We can try it a few more steps further:</DIV>
<DIV><BR>>>> root_list =3D ['R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L', 'R',=20
'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', =
'R', 'R',=20
'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', =
'R', 'L',=20
'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'L', 'R', 'R', =
'L']<BR>>>>=20
sternbrocot.path2frac(root_list)<BR>[17083879020L, =
9863382151L]<BR>>>>=20
math.sqrt(3) - 17083879020./9863382151<BR>0.0</DIV>
<DIV> </DIV>
<DIV>That's a pretty good estimate of the square root of =
three. The=20
next thing to do is to rewrite all of this as a generator, just to be =
trendy.=20
:)</DIV>
<DIV> </DIV>
<DIV>-PV </FONT></DIV></BODY></HTML>
------=_NextPart_000_00F6_01C19863.225ED8B0--