# 2012/09/25

## Way to rationalize

I confess. I claimed Gauche conformed R5RS, but it didn't.
In the early stage of development, I somehow thought `rationalize`

was an optional procedure and
I didn't bother to implement it.
But it wasn't optional; it was a library procedure. Ugh.

So one day I decided to write one, and here's my little journey.

What `(rationalize x e)`

returns
is the simplest rational number *p/q*, where
|*x* - *p/q*| is within the specified error bound *e*.

A quick search finds a method
to use Stern-Brocot tree.
Suppose we have *a/b* and *c/d* where *a/b < x < c/d*.
We calculate the mediant *M = (a+c)/(b+d)*, then we have the following cases:

- If
*M*is close enough to*x*, it's the answer. - If
*M*is above*x*, shrink the search range to*a/b < x < M*and recurse. - If
*M*is below*x*, shrink the search range to*M < x < c/d*and recurse.

Reasonable. This is the straight implementation (for simplicity,
I made it always return an exact rational, and also omitted
handling exceptional cases like `+inf.0`

).

Looks good, yeah?

gosh> (use math.const) #<undef> gosh> (rationalize1 pi 1/100) 22/7 gosh> (rationalize1 pi 1/1000) 201/64 gosh> (rationalize1 pi 1/100000) 355/113 gosh> (rationalize1 pi 1/10000000) 75948/24175 gosh> (rationalize1 (/ 3.0) 1e-17) 1/3

Well, not quite. This doesn't work well when the error bound is
small and *x* falls just outside of error bound from
either *a/b* or *c/d*. Suppose *x* is 1/1000000 and error
bound is 1/2000000. We start from 0/1 and 1/1, but the mediant is
always greater than *x* so only higher bound is moving towards
*x*. The denominator of the lower bound is fixed to 1, so the
mediant's denominator increases only 1 at every iteration.
It'll take 1000000 steps to get the answer.

So I searched a bit more. Someone described an algorithm about calculating (regular) continued fraction up to the point the error gets small enough. Calculating continued fraction is easy; it's basically the same as Euclidean algorithm.

gosh> (continued-fraction (exact pi)) (3 7 15 1 292 1 1 1 2 1 3 1 14 3 3 2 1 3 3 7 2 1 1 3 2 42 2) gosh> (continued-fraction (exact e)) (2 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12 1 1 11 1 1 1 11 5 1 1 2 1 4 2 1 1 9 17 3)

That is, π = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + ...)))) and so on.

(Note1: Since we have limited precision in the floating number constant pi and e, the above expansions don't match the real expansion of pi and e after first dozen items or so.)

(Note2: Each rational number has two regular
continued fraction representation, for you can always
decrement the last item and append '1' to the tail; e.g.
3/8 can be [0; 2, 1, 2] and [0; 2, 1, 1, 1].
Our `continued-fraction`

procedure always returns
the shorter one.)

Suppose we have a continued fraction [*a0*; *a1*, *a2*, ...], and
its calculated value up to *n*-th term is *p_n/q_n*, then
it is known that the value up to *n+1*-th term *p_{n+1}/q_{n+1}*
is *(a_{n+1}*p_{n} + p_{n-1})/(a_{n+1}*q_{n} + q_{n-1})*,
where *p_0/q_0* = *a0* and *p_1/q_1* = *a0 + 1/a1*.
So we can keep two previous approximated value to calcuate
the next approximation, instead of calculating the continued fraction
every time we add the term.

Here's an implementation. Again, we always return an exact result
for the simplicity.
Note that `continued-fraction`

returns a lazy sequence,
so we don't calculate the continued fraction more than necessary.

(If an gets empty, Qn must be exactly the same as x, so we never need to recurse further.)

Is this good enough?
Unfortunately, this converges too quickly and possibly
misses simpler (but less closer) approximation. In the following example,
`rationalize1`

finds the desired approximation (1/500)
but `rationalize2`

returns closer but less simple rational (1/999).

gosh> (rationalize1 2/1999 1/1000) 1/500 gosh> (rationalize2 2/1999 1/1000) 1/999

The Wikipedia entry of continued fraction has some
good explanation.
Basically, if we have a truncated continued fraction
[a0; a1, a2, ..., an-1, an], then
we have to check [a0; a1, a2, ... , an-1, ceiling(an/2)],
[a0; a1, a2, ... , an-1, ceiling(an/2)+1], ... etc.
This would make convergence slow again, when we have very large *an*.

The solution is in the next headings of the same Wikipedia article. We can find the simplest rational between two numbers, by comparing continued fractions of the two. Let's try some examples.

gosh> (continued-fraction 31415/10000) (3 7 14 1 8 2) gosh> (continued-fraction 31416/10000) (3 7 16 11)

To find a simplest rational between 31415/10000 and 31416/10000, we look for the first item that differ between the two sequences---in this case, it's the third item, 14 vs 16. We replace the item with (+ (min 14 16) 1), which is 15, and discard the rest; that is, [3; 7, 15] is the answer.

gosh> (+ 3 (/ (+ 7 (/ 15)))) 333/106

Given *x* and *e*, the answer of rationalize is
the simplest rational between *x-e* and *x+e*, hence
we get this implementation:

The actual code in Gauche includes handling of inexactness, infinites and NaNs.

(Note: The above implementation finds the answer in
the *closed* range between *x-e* and *x+e*,
which is ok for `rationalize`

, but if you want
to search in the *open* range, the `refine`

part
will be a bit more complicated. Besides, the above
implementation doesn't handle the case when e >= 1.)

Tags: r5rs, rationalize

Post a comment