2012/09/30

(exact (/ 3.0)) => ?

Since Gauche's floating point numbers (flonums) have limited precision, you can always calculate the exact number the flonum represents. If a flonum d is finite, it can be represented as follows:

 d = s * m * 2^e
    
 where s : 1 or -1
       m,e : integer

If e >= 0, d has no fractional part. If e < 0, we can have an exact rational number for d as follows:

 (s * m) / 2^(-e)

As of Gauche 0.9.3.3, exact works exactly like this. This is correct, but it has annoying behavior.

;; Gauche 0.9.3.3
gosh> (exact (/ 3.0))
6004799503160661/18014398509481984
gosh> (exact 0.1)
3602879701896397/36028797018963968

Neither 1/3 nor 1/10 can be represented accurately by a flonum; internally we use the closest flonum, but whose exact value is such a gigantic rational number.

But does that make sense? An inexact number is, after all, inexact. We can regard a flonum d as an approximation of any real number within a certain range around d. There are infinite number of rationals in the range. Why don't we pick the simpler rational among them?

(Since I sometimes see the argument that floating point numbers are inexact, I tempted to add a note. A floating point number itself can be either inexact or exact; there's nothing inherently inexact in the floating point number representation. It is the operations on them, in general, that are inexact (when rounding, overflow or underflow occur). If we use "exact floating point number", then the above gigantic rational representation is the correct answer. It is just that Gauche uses floating point representation for inexact numbers, which allows us to choose alternative rational representation.)

We already do similar sort of thing in printing flonums; we choose the shortest decimal representation that can uniquely identify given flonum.

Now we have rationalize, we can use the same algorithm to find out the simplest rational number within the precision of a flonum. We can't use rationalize directly, though, for the error range is more complicated.

In general, any real number r that falls between this range will be represented by a flonum d:

 d - q/2  <   r <  d + q/2   ; m is odd
 d - q/2  <=  r <= d + q/2   ; m is even
 
 where d = s * m * 2^e
   and q = 2^e

That is, q is the least quantum representable with the current exponent.

Since Gauche adopts round-to-even rule, the value that exactly falls between two floating point numbers would belong to the even flonum--- the one whose LSB of mantissa is zero.

This alone makes rationalize short for our purpose, since it always consider a closed range.

There are also two special cases we need to consider. If a flonum is a power of two, the flonums below is more "dense" than usual, so we have smaller lower error range:

 d - q/4  <= r <= d + q/2

Also, if d falls under denormalized range, q is fixed to to the smallest denormalized number.

So we implemented a new procedure real->rational, which can take individual upper and lower error boundary and a flag to indicate whether the range is open or not. Then rationalize can be written on top of it. Our new exact and inexact->exact calls real->rational if the argument is inexact real. Here's what you get in HEAD:

gosh> (exact (/ 3.0))
1/3
gosh> (exact 0.1)
1/10

If you need the exact rational representation of d itself, you can call real->rational with zero error tolerance.

gosh> (real->rational (/. 3) 0 0)
6004799503160661/18014398509481984
gosh> (real->rational 0.1 0 0)
3602879701896397/36028797018963968

Tags: Flonums, 0.9.4, exact, inexact