Convert a rational number to a floating-point number exactly

Question

We have two integers, $n$ and $d$. They are coprime (the only positive integer that is a divisor of both of them is $1$). They may be implemented as something that fits in a machine register, or they may be "big integers". Together, $n$ and $d$ represent the numerator and denominator of a rational number.

The goal is to convert the rational number into the corresponding floating-point number in a fixed-width binary floating-point format (for example, 32-bit binary IEEE 754 floating-point, known as float in the C programming language).

As long as using "big floats" (arbitrary-precision floating-point, for example with the MPFR C library) for implementing the algorithm is allowed, the solution is easy:

1. convert both $n$ and $d$ to big floats with a precision that is high enough so both integers would be exactly represented
2. divide the big float numerator with the big float denominator, yielding a big float number
3. convert the above big float number to the required fixed-precision floating-point format

I wonder, however, what would a more from-first-principles solution look like, one that wouldn't use big floats?

Answer 1

Big floats are nothing but big ints with a scaling factor, which you can handle by yourself.

If the two numbers are given as big ints and you are only allowed integer division, multiply the numerator by a power of $2$ large enough that the quotient has at least $24$ (resp. $53$) significant bits. Then normalize the result to single (resp. double)-precision floating-point.

If the numbers are given as integers, you will need to resort to extended precision arithmetic, i.e. 64 bits over 32 or 128 over 64. Append a zero word to the numerator, divide and normalize the result.

Answer 2

I am investigating a floating point number bug in CCL (Clozure Common Lisp), which since Common Lisp is using an old standard by default 32 bit IEEE 753 float known as single-float in most implementation is used.

The bug is demonstrate by the observation that converting a rational number to single-float could be way off:

(float 41107100000541273/100000000000) ;; => 411071.03
(float 41107100000541273/100000000000 0.0D0) ;; => 411071.00000541273D0

This is because CCL used the naive implementation that when neither the numerator or denominator are bignums or if their bit length are both less than single-float-bias that is 126 bits, convert both to single-float and do division is used. This gives:

41107100000541273 converted to single float is 1.0 * 0.57047564 * 2^56 100000000000 is 1.0 * 0.72759575 * 2^37 Division gives: 1.0 * 0.78405577 * 2^19 which in decimal is 4.1107103E+5.

Note that 41107100000541273 is only 56 bits so their implementation is definitely wrong.

So I lookup SBCL's implementation which does not suffer from this inaccuracy issue and turns out the algorithm they used is the following, and despite the algorithm is also intended to work for bignum integers, not even bigfloat or double float is involved:

The only noun ANSI Common Lisp standard function is

(ccl::make-short-float-from-fixnums significand biased-exp sigN)

which should be self explanatory. labels defines local functions and loop just does infinity loop here, which the help of a reference book to CL this algorithm can easily be converted to C.

(defun short-float-ratio (x)
  (let* ((signed-num (numerator x))
         (plusp (plusp signed-num))
         (num (if plusp signed-num (- signed-num)))
         (den (denominator x))
         (digits 24) ; single-float-digits
         (scale 0))
    (declare (fixnum digits scale))
    ;; Strip any trailing zeros from the denominator and move it into the scale
    ;; factor (to minimize the size of the operands.)
    (let ((den-twos (1- (integer-length (logxor den (1- den))))))
      (declare (fixnum den-twos))
      (decf scale den-twos)
      (setq den (ash den (- den-twos))))
    ;; Guess how much we need to scale by from the magnitudes of the numerator
    ;; and denominator. We want one extra bit for a guard bit.
    (let* ((num-len (integer-length num))
           (den-len (integer-length den))
           (delta (- den-len num-len))
           (shift (1+ (the fixnum (+ delta digits))))
           (shifted-num (ash num shift)))
      (declare (fixnum delta shift))
      (decf scale delta)
      (labels ((float-and-scale (bits)
                 (let* ((bits (ash bits -1))
                        (len (integer-length bits)))
                   (cond ((> len digits)
                          (assert (= len (the fixnum (1+ digits))))
                          (scale-float (floatit (ash bits -1)) (1+ scale)))
                         (t
                          (scale-float (floatit bits) scale)))))
               (floatit (bits)
                 (let ((sign (if plusp 1 -1)))
                   (ccl::make-short-float-from-fixnums bits 126 ; single-float-bias
                                      sign))))
        (declare (inline floatit))
        (loop
          (multiple-value-bind (fraction-and-guard rem)
              (truncate shifted-num den)
            (let ((extra (- (integer-length fraction-and-guard) digits)))
              (declare (fixnum extra))
              (cond ((/= extra 1)
                     (assert (> extra 1)))
                    ((oddp fraction-and-guard)
                     (return
                       (if (zerop rem)
                           (float-and-scale
                            (if (zerop (logand fraction-and-guard 2))
                                fraction-and-guard
                                (1+ fraction-and-guard)))
                           (float-and-scale (1+ fraction-and-guard)))))
                    (t
                     (return (float-and-scale fraction-and-guard)))))
            (setq shifted-num (ash shifted-num -1))
            (incf scale)))))))

Unfortunately, TAOCP 2nd Volume and Handbook of floating-point arithmetic does not give useful information on this topic so I cannot give a confirmation on how reliable the method is. I upvoted the problem since this is actually a meaning real world problem and has pitfalls.