# Re: Factor

Factor: the language, the theory, and the practice.

## 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

1. `abs (numerator y) <= abs (numerator y')`, and
2. `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
{
{ 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.