1
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ In 2., is it always possible to find an exact representation of the quotient? $\endgroup$
    – greybeard
    May 6, 2023 at 4:21
  • $\begingroup$ @greybeard No, but I'm assuming the result would be correctly rounded, as with MPFR. $\endgroup$ May 6, 2023 at 11:15
  • $\begingroup$ (Resisting the temptation to go grab Knuth's TAoCP volume 2 to look it up...) $\endgroup$
    – vonbrand
    May 6, 2023 at 18:17

2 Answers 2

3
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ The way I'm reading your answer, the algorithm requires performing multiple integer divisions until one of them produces a quotient with the required number of digits. Is there perhaps a way to solve this with a constant number of integer divisions? $\endgroup$ May 7, 2023 at 4:17
  • $\begingroup$ You say "multiply the numerator by a power of $2$ large enough that the quotient has at least $24$ significant bits. then normalize". But if I get more than $24$ significant bits, it's not clear to me how to proceed with getting rid of the extra bits. It also seems like a double-rounding problem, the integer division included the first rounding, getting rid of extra bits after that would be a second rounding. $\endgroup$ May 7, 2023 at 10:47
  • $\begingroup$ I can build an IEEE representation using ldexp, that's not a problem. $\endgroup$ May 7, 2023 at 10:50
  • $\begingroup$ @user2373145: good, solved then. Though I don't see the connection to "from-first-principles" anymore. $\endgroup$
    – user16034
    May 7, 2023 at 10:50
  • 2
    $\begingroup$ @user2373145 "how to proceed with getting rid of the extra bits" After normalization, all that is needed for correct rounding with any IEEE-74 rounding mode is one additional bit (round bit) and all subsequent bits accumulated (ORed) into a sticky bit. It is not necessary to compute the quotient to an infinite number of bits; once there are enough quotient bits + a round bit simply look at the remainder: if it's zero, the sticky bit is zero, otherwise the sticky bit is 1. $\endgroup$
    – njuffa
    Sep 3, 2023 at 2:06
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ It has been thirty years since I last looked at LISP code, but isn't this just using binary longhand division with an extra quotient bit generated for rounding and looking at the remainder for sticky-bit information? $\endgroup$
    – njuffa
    Oct 5, 2023 at 4:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.