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.