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 numbersx
andepsilon
, returns the simplest rational number withinepsilon
ofx
. A rational numbery
is said to be simpler than anothery'
if
abs (numerator y) <= abs (numerator y')
, anddenominator 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