I’ve been thinking about the display code from the last posting, and realized it’s still not enough to handle the vagaries of floating-point numbers.  Even if we manually round a number up, there is no guarantee that it won’t be rounded down again as it is coerced to a double-float prior to passing through the format statement.  There’s always the chance that our result will be wrong in the final displayed digit.  So, we’re going to be making some changes to the code to enforce more strongly the use of rational numbers when possible.

The first place this will appear is in the printing code we worked on last time.  We will now enforce the requirement that numbers to be printed are always rationals.  Then, we will need to be able to print a rational to the desired precision with absolute confidence that it will be rendered correctly in all digits.  To this end, we’re going to write our own equivalent of format with the ~,vE option.  This function is presented below.

So, what this does is to extract the numerator and denominator of the absolute value of the number to be printed.  An exponential power of 10 is initialized to zero.

Next, we adjust the numerator or denominator by factors of 10, and record these changes in the exponent, until the numerator is at least as large as the denominator, and strictly less than 10 times the denominator.

At this point, we do grade-school long-hand division.  We compute and record the divisor, subtract the appropriate number from the numerator to get the remainder, then shift the decimal and keep going.

Now, for rounding.  If the remainder is at least 1/2 of the denominator, we want to round up.  We set the carry flag to 1 and then apply that flag to the collection of digits from least to most significant.  The carry flag can eventually become zero, at which point no further digits are adjusted.  If we exit that mapcar with the carry flag set to 1, then we have carried all the way to the most significant digit, which we converted from a 9 to a 0.  So, we’re going to have to put a 1 in front of that, after first dropping the lowest digit, as it’s below the requested precision.  Then, we reverse the list we’ve been keeping, so that the digits are from most to least significant.

Finally, we put together the digits and the exponent to produce a string in scientific notation.

Here is the function:

(defun render-rational-as-sci (rval n-digits)
  "In the end, we can't trust format because of its rounding rules, and the coercion.  Even going to double-float can't guarantee that we won't slip a digit in the last place.  So, here we 'display' a rational number by longhand division."
  (when (= rval 0)
    (return-from render-rational-as-sci
      (format nil "~,vE" n-digits 0.0)))

  (let* ((result '())
         (rv (make-string-output-stream))
         (sign (if (> rval 0) 1 -1))
         (num (numerator (abs rval)))
         (den (denominator (abs rval)))
         (exponent 0))
    
    (do ()
        ((>= num den))
      (setf num (* 10 num))
      (decf exponent))
    (do ()
        ((< num (* 10 den)))
      (setf den (* 10 den))
      (incf exponent))

    (dotimes (i (1+ n-digits))
      (let ((digit (floor (/ num den))))
        (push digit result)
        (decf num (* digit den))
        (setf num (* 10 num))))

    (let ((carry (if (>= (/ num den) (/ 1 2)) 1 0)))
      (setf result (mapcar #'(lambda (x)
                               (setf x (+ carry x))
                               (cond
                                 ((= 10 x)
                                  (setf x 0))
                                 (t
                                  (setf carry 0)))
                               x) result))

      (cond
        ((= carry 1)
         ;; rounded all the way to the beginning. Fix.
         (pop result)
         (setf result (reverse result))
         (incf exponent)
         (push 1 result))
        (t
         (setf result (reverse result)))))

    (format rv "~A~D."
            (if (= sign -1) "-" "")
            (car result))
    (format rv "~{~D~}" (cdr result))
    (format rv "e~A~2,'0D"
            (if (< exponent 0) "-" "")
            (abs exponent))

    (get-output-stream-string rv)))

Here is some sample output:

CL-USER> (render-rational-as-sci (/ 1005 100) 2)
"1.01e01"
CL-USER> (render-rational-as-sci (/ -1005 100) 2)
"-1.01e01"
CL-USER> (render-rational-as-sci (/ 99999 100000) 2)
"1.00e00"
CL-USER> (render-rational-as-sci (/ 99999 100000) 3)
"1.000e00"
CL-USER> (render-rational-as-sci (/ 99999 100000) 4)
"9.9999e-01"
CL-USER> (render-rational-as-sci (/ 99999 100000) 5)
"9.99990e-01"

This function is found in the display.lisp module in the git repository, under the tag v2014-11-06.