[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Rationalize



I thought I'd put my two cents since I had something to do with this stuff.
First, unless my memory is massively wrong, I was the one whoe first
implemented rationalize in Macsyma.  I believe it was Bill Gosper that
requested it.  Bill correct me if I'm wrong.

Second, while I have your attention let me explain the algorithm since it is
really quite simple.  All positive real numbers can be represented as a
continued fraction, viz.

alpha = a0 + 1/(a1 + 1/(a2 + ... 1/an) ... )

If alpha is not a rational number then the continued fraction does not
terminate.  Initial segments of the continued fraction lead to the best
rational approximations to alpha.  (This is proven in any book on elementary
number theory.)

The algorithm in rationalize consists of two parts that are melded together.
First, one computes the sequence of partial quotients: a0, a1, a2, ...
and then use them to compute rational approximations to alpha.  

The ai are easy to calculate, a0 is the floor (truncate) of alpha,
a1 is the floor of 1/(alpha-a0), etc.  Let alpha[0] = alpha, a0 =
truncate(alpha[0]) and alpha[i] = 1/(alpha(i-1) - ai) then
ai = truncate(alpha(i)).

Computing the rational approximation is a little more complicated.
Let p[0]/q[0] = a0/1, p[1]/q[1] = a0 + 1/a1, etc. (the rational
approximations to alpha).  Also let p[-1] = 0, q[-1] = 1.
Then p[i] = ai p[i-1] + p[i-2]
     q[i] = ai q[i-1] + q[i-2].

It is also useful to know that each pair of p[i] and q[i] is relatively
prime.  

The code given by Skef, intertwines the two phases of this algorithm to
avoid consing or using the stack.  Remember it was originally implemented
for Macsyma which was terribly short on space.