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:

( scratchpad ) 33/100 1/10 approximate . 1/3 ( 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):

( scratchpad ) USING: math.constants math.statistics math.vectors ; ( scratchpad ) 10000 [1,b] [ 0.99 swap ^ float>ratio ] map ( scratchpad ) pi float>ratio '[ _ swap approximate ] map ( 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.

## No comments:

Post a Comment