# Boost :

From: Dave Abrahams (abrahams_at_[hidden])
Date: 2000-11-27 19:26:07

> However, there is an issue in the context of rational<>, namely the
need for
> a maximum denominator. I'm trying to implement a constructor from a
double.
> So where do I get a maximum denominator from? I could require the
user to
> specify it as a constructor argument, or I could use the maximum
value of
> the underlying integer type of rational<>. The former approach
seems to me
> not very user-friendly, whereas the latter doesn't avoid the
problem of
> "over-precise" values.
>
> I believe that mathematically, it is possible to view machine
doubles as
> representing a series of discrete points on the mathematical real
line. The
> problem I am looking to solve is to find the simplest rational (a
> well-defined concept) which is closer (in the absolute mathematical
sense)
> to the supplied double than to any other value representable as a
double. I
> believe this is calculable, and furthermore, I believe that there
exist
> implementations of the algorithm. Unfortunately, I don't know of
existing
> code, nor do I have the expertise to implement the algorithm from
first
> principles.
>
> I'd appreciate Boost members' comments. Should I include Reggie's
algorithm
> (and if so, what should I use as a maximum denominator)? Or should
I wait
> until I can find an implementation of the "ideal" algorithm?

Try a web-search at Google for "rationalize algorithm". That turns up
many promising links. I stole the following algorithm (written in
Python) from similar Java code at . I don't know if it's the ideal
algorithm you're looking for, but it seems to do the job:

---------
def rationalize(f):
if f == 0:
return rat(0,1)
sign = 1
if f < 0:
sign = -1
f = math.fabs(f)

# Compute continued fraction expansion */

# Initialize continued fraction expansion
a0=0
a1=1
b0=1
b1=0

# Compute initial term
ipart=int(math.floor(f))
f=f-ipart
a=a1*ipart+a0
b=b1*ipart+b0
a0=a1
b0=b1
a1=a
b1=b

# Loop for continued fraction expansion
try:
while f!=0:
f=1.0/f
ipart=int(math.floor(f))
f=f-ipart
a=a1*ipart+a0
b=b1*ipart+b0
a0=a1
b0=b1
a1=a
b1=b
print 'approximating:',a1,'/', b1

except: # This catches overflow errors, indicating it's time to
stop
pass

return rat(sign*a1,b1)
---------
>>> rats.rationalize(math.pi)
approximating: 22 / 7
approximating: 333 / 106
approximating: 355 / 113
approximating: 103993 / 33102
approximating: 104348 / 33215
approximating: 208341 / 66317
approximating: 312689 / 99532
approximating: 833719 / 265381
approximating: 1146408 / 364913
approximating: 4272943 / 1360120
approximating: 5419351 / 1725033
approximating: 80143857 / 25510582
approximating: 245850922 / 78256779
approximating: 817696623 / 260280919
rat(817696623, 260280919)
>>> rat(math.pi)
rat(884279719003555L, 281474976710656L)
>>> rationalize(1.0/3)
>>> rats.rationalize(1.0/3)
approximating: 1 / 3
rat(1, 3)
>>> rats.rationalize(1.0/5)
approximating: 1 / 5
rat(1, 5)
>>> rats.rationalize(1.0/213)
approximating: 1 / 213
rat(1, 213)
>>> rats.rationalize(7.0/213)
approximating: 1 / 30
approximating: 2 / 61
approximating: 5 / 152
approximating: 7 / 213
rat(7, 213)
>>> rats.rationalize(3)
rat(3)
>>> rats.rationalize(2.14)
approximating: 15 / 7
approximating: 107 / 50
rat(107, 50)
>>> rats.rationalize(1/0.7)
approximating: 3 / 2
approximating: 10 / 7
rat(10, 7)
>>> rats.rationalize(1/7.0)
approximating: 1 / 7
rat(1, 7)