## Approximating Rationals

Monday, September 6, 2010

A few months ago, there was an
article
about approximating numbers with fractions using the `approxRational`

function in Haskell. The documentation for the
function states:

`approxRational`

, applied to two real fractional numbers`x`

and`epsilon`

, returns the simplest rational number within`epsilon`

of`x`

. A rational number`y`

is said to be simpler than another`y'`

if

`abs (numerator y) <= abs (numerator y')`

, and`denominator y <= denominator y'`

.Any real interval contains a unique simplest rational; in particular, note that

`0/1`

is the simplest rational of all.

The source
code
for `approxRational`

is:

```
approxRational :: (RealFrac a) => a -> a -> Rational
approxRational rat eps = simplest (rat-eps) (rat+eps)
where simplest x y | y < x = simplest y x
| x == y = xr
| x > 0 = simplest' n d n' d'
| y < 0 = - simplest' (-n') d' (-n) d
| otherwise = 0 :% 1
where xr = toRational x
n = numerator xr
d = denominator xr
nd' = toRational y
n' = numerator nd'
d' = denominator nd'
simplest' n d n' d' -- assumes 0 < n%d < n'%d'
| r == 0 = q :% 1
| q /= q' = (q+1) :% 1
| otherwise = (q*n''+d'') :% n''
where (q,r) = quotRem n d
(q',r') = quotRem n' d'
nd'' = simplest' d' r' d r
n'' = numerator nd''
d'' = denominator nd''
```

I couldn’t find an equivalent algorithm in Factor, so I decided to translate it (using lexical variables to make the comparison between the languages as direct as possible).

```
USING: combinators kernel locals math math.functions ;
:: (simplest) ( n d n' d' -- val ) ! assumes 0 < n/d < n'/d'
n d /mod :> ( q r )
n' d' /mod :> ( q' r' )
{
{ [ r zero? ] [ q ] }
{ [ q q' = not ] [ q 1 + ] }
[
d' r' d r (simplest) >fraction :> ( n'' d'' )
q n'' * d'' + n'' /
]
} cond ;
:: simplest ( x y -- val )
{
{ [ x y > ] [ y x simplest ] }
{ [ x y = ] [ x ] }
{ [ x 0 > ] [ x y [ >fraction ] bi@ (simplest) ] }
{ [ y 0 < ] [ y x [ neg >fraction ] bi@ (simplest) neg ] }
[ 0 ]
} cond ;
: approximate ( x epsilon -- y )
[ - ] [ + ] 2bi simplest ;
```

We can experiment with the function:

```
IN: scratchpad 33/100 1/10 approximate .
1/3
IN: scratchpad 1/3 1/10 approximate .
1/3
```

We can use the code I
posted
a few days ago to convert floating-point numbers to fractions. That will
be useful for running the experiment from the original blog
post
(to approximate `pi`

with varying degrees of accuracy):

```
IN: scratchpad USING: math.constants math.statistics math.vectors ;
IN: scratchpad 10000 [1..b] [ 0.99 swap ^ float>ratio ] map
IN: scratchpad pi float>ratio '[ _ swap approximate ] map
IN: scratchpad sorted-histogram reverse .
{
{ 3+39854788871587/281474976710656 3447 }
{ 3+16/113 572 }
{ 3+464086704149/3277618523158 433 }
{ 3+1/7 297 }
{ 3+244252/1725033 288 }
{ 3+1470786781/10387451211 248 }
{ 3 194 }
{ 3+11080585/78256779 186 }
{ 3+56667685227/400216280932 176 }
{ 3+449714866/3176117225 172 }
}
```

Looking at the top 10 most common approximations, you can see the first
(and most frequent) result is the exact fractional value of the constant
`pi`

, the fourth result is the common “grade school” estimate `22/7`

,
and the second result is the “best” simple estimate `355/113`

.

*Note: the output is slightly different than the Haskell program. Factor
considers several of these rational numbers to be equal when converted to
floating-point:*

`3+39854788871587/281474976710656`

`3+464086704149/3277618523158`

`3+1470786781/10387451211`

`3+11080585/78256779`

`3+56667685227/400216280932`

The code (and some tests) for this is on my GitHub.