Re: rationals initialization bug

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)