Over the past few posts we’ve discussed the realization that we should be using rational numbers for most operations.  To that end, we wrote code to render a rational as a printable string in any one of the three output modes of the calculator.

We have rewritten some of the stack/memory code so that only rational numbers to 10 digits of precision are stored there.  To truncate the rationals at 10 digits, they are parsed as printable 10-digit scientific-notation strings, then the strings are parsed back into exact rational equivalents.

As mentioned earlier, rational numbers are fine for many purposes, but if a function that requires floating point numbers, like sqrt or log, receives a rational number, that number is promoted to a single-precision float.  Single-precision floats generally have less than 10 digits of precision internally, so this is a lossy operation that must not be permitted.  To avoid this, the stack and memory operations can optionally convert the number they return to a double-precision float before handing it to the caller.  The calling environment must explicitly declare whether it will be requesting doubles or rationals.  The parsing code in key-structs.lisp can either set all its stores and recalls to rationals, or all to doubles; the syntax does not permit the mixing of those two modes.  Therefore, if a rational is requested, and a single-precision or double-precision float is later pushed or stored, an error will be raised, one that is not caught by the handlers in the key forms, as this is a programming bug that must be fixed.  In this way, we catch all operations that convert rationals to single-precision floats.

We’ve also written a function to retrieve the parsed form associated with a defined key, for debugging perusal.  Reading over the forms produced for the keys we’ve defined so far, some parsing errors were found and the code was corrected.

Here are the two new conditions used by the code that protects against single-precision float errors:

(define-condition invalid-float-arrived (error)
  ((val         :initarg value
                :reader get-val))
  (:documentation "A floating-point number was pushed when rationals were promised.  This is a coding bug and should be fixed.")
  (:report (lambda (c s)
             (format s "The float value ~A was encountered."
                     (get-val c)))))

(define-condition single-precision-float (error)
  ((val         :initarg value
                :reader get-val))
  (:documentation "A single-precision floating-point number was seen.  This is below the promised precision of the calculator, and is a programming bug.")
  (:report (lambda (c s)
             (format s "The single-precision float value ~A was encountered."
                     (get-val c)))))

Here is the function that is used by stack and memory store operations to validate the passed value:

(defun fix-input-val (val rflag)
  (when (eq rflag :DOUBLE-FLOAT)
    (when (eq (type-of val) 'single-float)
      (error (make-condition 'single-precision-float
                             :value val)))
         (setf val (rational val)))

  (unless (rationalp val)
    (error (make-condition 'invalid-float-arrived
                           :value val)))

  (round-to-ultimate-precision val))

Here is the function used by stack and memory retrieve operations to convert to double-precision floats if desired:

(defun fix-output-val (val rflag)
  (if (eq rflag :DOUBLE-FLOAT)
      (coerce val 'double-float)
      val))

The current code is in the git repository under the tag v2014-11-10.

The HP-67 calculator has some specialized behaviour when displaying numbers in fixed mode.  Numbers less than 1 always have a leading zero before the decimal, leaving 9 digits of precision for the remainder of the value.  If asked to display 0.001 to two digits of precision in fixed mode, the calculator will display the number in scientific notation instead, as that value rounds to a visual representation of zero.  If asked to display 1234567890.123 to two digits of precision, it will display only 1234567890, as the display has only ten digits.

So, now we need functions to display a rational number in either scientific or fixed mode.  The scientific mode is fairly simple, as we can use the long division code we wrote earlier to produce a sequence of digits, and we know the exponent on the total.  Fixed mode is more complicated, because of the automatic conversion to scientific notation.

A number that is too wide to fit in the 10 digits of the display must be at least 10^10-(1/2).  With the rounding rules, that number will round up to an 11 digit number.  A number that is too small to fit in the number of digits of precision requested must be smaller than 10^-d / 2.  With the rounding rules, that number will round up to 1 in the last digit displayed after the decimal.  In between these values, the magnitude of the number (the power of 10 represented by the most significant non-zero digit) and the precision after the decimal point will, together, determine how many digits total must be displayed.

The temptation is to find the base-10 logarithm of the rational number, take the floor function on that, and you have the magnitude of the rational.  However, this is not safe.  Rational numbers are promoted to single-precision floats, and cannot generally be exactly represented.  Sometimes, the error introduced in the conversion from rational to float can cause a number to cross to a different magnitude, as seen in this example:

CL-USER> (log (/ 1000001 10) 10.d0)
5.000000276521525d0
CL-USER> (log (/ 10000001 10) 10.d0)
6.000000083320534d0
CL-USER> (log (/ 100000001 10) 10.d0)
6.999999890119543d0

So, for safety, after the magnitude is determined with the logarithm, it is verified using purely rational arithmetic, and adjusted if necessary.

The next thing of interest is rounding numbers to lower precision.  It is easier to write the code to produce a certain fixed precision, then successively round the number to the desired precision later.  However, successive rounding has a major issue in the case of the number 0.44445.  If you were to round this successively, you might first get 0.4445, then 0.445, then 0.45, then 0.5.  This sequence is clearly incorrect.  When performing successive rounding, when the digit you’re eliminating is 5, you have to know how you got that 5 digit, whether it was from rounding down, or rounding up.  If it came from rounding down or was exact, then the 5 causes a rounding up of the next-higher-precision number.  If it came from rounding up, then the 5 causes no rounding up.  The long division function has been modified to return that information in a third entry in the returned list, and a new function that understands this format, round-long-division-result, has been written to make use of it.

With these changes, we now are ready to display any rational number in either fixed or scientific notation, in a way that exactly matches the HP-67’s behaviour.

The current code is checked into the git repository under the tag v2014-11-08.

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.

The temptation is to say that formatting output is simple, just use a format statement with ~F or ~E and the appropriate flags.  In fact, things are more complicated than that.  The difficulty is not just because format is permitted by the standard to round 0.5 up or down at its discretion, while the calculator always rounds 0.5 up for positive numbers and down for negative numbers.  There are more awkward problems than that ahead.

The HP-67 calculator, using BCD arithmetic, always had an internal representation that exactly matched the maximum output precision.  Every number that could be displayed could be exactly represented, and there were no left over representations.  Such is not the case with many modern floating-point platforms, such as the familiar IEEE-754 representation used in many modern computers.  Certain decimal values cannot be exactly represented as floats or double-floats, and this leads to some difficulties when trying to emulate the behaviour of the BCD calculator.  I’ve mentioned before the futility of just “throwing more bits at the problem”, this rarely solves the issue, only hides it in more subtle ways.

Here is an example of how things can go wrong.  If you want to represent a small floating-point number to, say, two decimal places, you might be tempted to scale the number up, round it off, then divide it back down.  Watch what happens when I divide two exact powers of ten in double-precision arithmetic:

CL-USER> (/ 1.d99 1.d97)
99.99999999999999d0

This is not helpful.  In fact, a lot of our manipulations of numbers for display are going to have to be in string form.  We round numbers off like grade-schoolers, looking at the digital representation and tweaking it appropriately, everything being passed around as strings, not as numbers.

So, we now have a new module in the tree, display.lisp.  This allows us to display a passed number, either a double-precision float or a rational, in one of three modes.  FIXED displays the number in fixed-point mode, if it can be represented on the 10-digit display with the desired number of digits after the decimal point.  If a number cannot be so displayed, it will be displayed in scientific notation.  An example of such non-displayable numbers might be 0.001 with 2 digits of precision, which would erroneously display as zero in fixed mode, or 100000000000, which has too many digits to display on the screen.

SCIENTIFIC displays in the familiar scientific notation.  If negative, a ‘-‘ is displayed.  Then comes the mantissa which consists of a non-zero digit followed by a decimal point and 0 or more further digits.  After this is either a space, or a minus sign, depending on whether the exponent is positive or negative.  Finally, a 2-digit exponent.

ENGINEERING is much like SCIENTIFIC, but if the exponent is not a multiple of 3, the next higher multiple of 3 is chosen and the mantissa is adjusted to compensate.  The mantissa will always, then, be at least 1 and less than 1000.

Here are the interesting parts of the file.  First, we need to know whether a fixed-mode display has rounded a number to look like zero.  This function scans a string and returns non-nil if the string contains at least one non-zero digit:

(defun string-contains-non-zero-digit (string)
  (dotimes (i (length string))
    (let ((one-char (char string i)))
      (when (and (digit-char-p one-char)
                 (char/= one-char #\0))
        (return-from string-contains-non-zero-digit t))))
  nil)

Next, we’re going to need to manipulate the components of a scientific-notation string, so we have a function that returns a list of the sign of the number, the mantissa, the sign of the exponent, and the exponent:

(defun break-down-sci-notation (string)
  (let* ((negative (char= (char string 0) #\-))
         (epos (position-if #'(lambda (x)
                                (or (char= x #\e)
                                    (char= x #\d))) string))
         (neg-expt (char= (char string (1+ epos)) #\-))
         (mantissa (subseq string
                           (if negative 1 0)
                           epos))
         (expt (subseq string
                       (if neg-expt
                           (+ 2 epos)
                           (1+ epos)))))
    (when (char= (char expt 0) #\+)
      (setf expt (subseq expt 1)))
    (list (if negative "-" " ")
          mantissa
          (if neg-expt "-" " ")
          expt)))

Our engineering notation code is going to have to be able to shift the decimal point up to two digits to the right, padding with zeroes if there aren’t enough characters after the decimal.  It has to be able to handle a bad case that can appear sometimes.  Normally we expect the format statement with ~E to return a mantissa at least one and strictly less than 10.  However, here is what happens sometimes on SBCL v1.1.14:

CL-USER> (format nil "~,8,2E" 1.0d-6)
"10.00000000d-07"

So, the function to shift decimal points has to notice when the point starts in the wrong place, and shift one digit less, while adjusting the exponent appropriately.  That is the what d-pos does in this code:

(defun shift-char-to-right (string start-pos n-shift
                            &key (padding #\0))
  "Moves the character at start-pos n-shift to the right"
  (let ((workspace (copy-seq string))
        (moved (char string start-pos))
        (pad-len (- (+ 1 start-pos n-shift) (length string))))

    (when (> pad-len 0)
      (setf workspace
            (concatenate 'string
                         workspace
                         (make-sequence 'string
                                        pad-len
                                        :initial-element padding))))
    (dotimes (i n-shift)
      (setf (char workspace (+ i start-pos))
            (char workspace (+ i 1 start-pos)))
      (setf (char workspace (+ i 1 start-pos)) moved))
    workspace))

Here, now, is the code to print numbers in fixed mode:

(defun format-for-printing-fix (val digits-after-decimal
                                &key readable)

  (when (= val 0)
    (return-from format-for-printing-fix
      (format nil "~,vF" digits-after-decimal 0.0d0)))
  
  (let* ((negmult (if (< val 0) -1.0d0 1.0d0))
         (scaleup (expt 10.0d0 digits-after-decimal))
         (magnitude (abs val))
         (rounded (* negmult
                     (floor (+ 0.50000000004d0
                               (* magnitude scaleup)))))
         (first-try (format nil "~,v,vF"
                            digits-after-decimal
                            (- digits-after-decimal)
                            rounded))
         (max-width (+ 1 *digits-in-display*
                       (if (< val 0) 1 0))))

    (let ((overrun (- (length first-try) max-width)))
      (cond
        ((and (> overrun 0)
              (<= overrun digits-after-decimal))
         (format-for-printing-fix val
                                  (- digits-after-decimal
                                     overrun)
                                  :readable readable))
        ((> overrun 0)
         (format-for-printing-sci val digits-after-decimal
                                  :readable readable))
        ((and (/= val 0)
              (not (string-contains-non-zero-digit first-try)))
         (format-for-printing-sci val digits-after-decimal
                                  :readable readable))
        (t
         first-try)))))

The code for scientific mode:

(defun format-for-printing-sci (val digits-after-decimal
                                &key readable)
  (when (= 0 val)
    (return-from format-for-printing-sci
      (if readable
          "0.0d0"
          (format nil "~,vE" digits-after-decimal 0.0d0))))
  
  (let* ((magnitude (abs val))
         (first-try (format nil "~A~,v,2E"
                            (if (< val 0) "-" "")
                            digits-after-decimal
                            magnitude))
         formatted)

    (setf first-try (round-sci-notation-to-digits first-try
                                                  digits-after-decimal))

    (unless readable
      (destructuring-bind (sign mantissa e-sign exponent)
          (break-down-sci-notation first-try)

        (setf formatted
              (format nil "~A~vA~A~A"
                      sign
                      (1+ *digits-in-display*)
                      mantissa
                      e-sign
                      exponent))))

    (if readable
        (values first-try first-try)
        (values formatted first-try))))

The code for engineering mode:

(defun format-for-printing-eng (val digits-after-decimal
                                &key readable)
  (multiple-value-bind (junk parsed)
      (format-for-printing-sci val digits-after-decimal
                               :readable readable)
    (declare (ignore junk))
    (when readable
      (return-from format-for-printing-eng parsed))

    (destructuring-bind (sign mantissa e-sign exponent)
        (break-down-sci-notation parsed)

      (let* ((e-num (read-from-string exponent))
             (man-len (length mantissa))
             (shift-num (mod e-num 3)))
        (when (string= e-sign "-")
          (setf shift-num (mod (- 3 shift-num) 3)))
        (when (and (= man-len 3) (= shift-num 2))
          (setf mantissa (format nil "~A0" mantissa)))

        (dotimes (i shift-num)
          (psetf (char mantissa (1+ i)) (char mantissa (+ 2 i))
                 (char mantissa (+ 2 i)) #\.))

        (when (string= e-sign "-")
          (setf e-num (* -1 e-num)))
        (decf e-num shift-num)

        (format nil "~A~vA~A~2,'0D"
                sign
                (1+ *digits-in-display*)
                mantissa
                e-sign
                (abs e-num))))))

This module, and a few supporting changes, are all available in the git repository with the tag v2014-11-04.

At this point, the indirection code was getting unreasonable.  The case-insensitive label we use for indirection, “(i)”,  was starting to show up in too many places.  There’s no reason that the logic for indirection as applied to memory and flags can’t sit entirely in the memory and flag code.  So, this was pushed back into that module.  For simplicity, the I-register was moved from a special value in the structure to just another memory register, one indexed by the label “(i)”.  A new condition was defined for indirection operations that are attempted with the I-register out of its valid domain.  The HP-67 calculator required that the value in the I-register be from 0 to 25, inclusive, for store operations, and 0 to 3, inclusive, for flag operations.  Operations with invalid I-register will now signal a condition that will cause the calculator to enter an error state.

The memory code now looks like this:

(defun canonicalize-memory-name (stack mem-name)
  (when (integerp mem-name)
    (setf mem-name (format nil "~D" mem-name)))
  (assert (stringp mem-name))
  (cond
    ((string-equal mem-name "(i)")
     (multiple-value-bind (junk int-val str-val)
         (get-i-register stack)
       (declare (ignore junk))
       (cond
         ((and (not *unlimited-indirection*)
               (or (< int-val 0) (> int-val 25)))
          (error (make-condition 'i-register-range-error
                                 :value int-val
                                 :min-allowed 0
                                 :max-allowed 25)))
         ((= int-val 25)
          "(i)")
         ((> int-val 19)
          (subseq "ABCDE"
                  (- int-val 20)
                  (- int-val 19)))
         (t
          str-val))))
    (t
     mem-name)))


(defun store-memory-by-name (stack name val)
  "Does no indirection, just stores under the name."
  (setf (stack-memory stack)
        (delete-duplicates
         (push (cons name val)
               (stack-memory stack))
         :key 'car
         :test 'string=
         :from-end t))
  val)

(defun recall-memory-by-name (stack name)
  "Does no indirection, just recalls from the name."
  (let ((record (assoc name
                       (stack-memory stack)
                       :test 'string=)))
    (if record
        (cdr record)
        0)))


(defun store-memory (stack name val)
  (setf name (canonicalize-memory-name stack name))
  (store-memory-by-name stack name val))


(defun recall-memory (stack name)
  (setf name (canonicalize-memory-name stack name))
  (recall-memory-by-name stack name))

(defun set-i-register (stack value)
  (store-memory-by-name stack "(i)" value))

;; Returns 3 values.  The unmodified value of I, the greatest-integer
;; value, and a string holding the greatest-integer value
(defun get-i-register (stack)
  (let ((rval (recall-memory-by-name stack "(i)")))
    (values
     rval
     (floor rval)
     (format nil "~D" (floor rval)))))

The flag code looks like this:

(defun canonicalize-flag-name (stack flag-name)
  (when (integerp flag-name)
    (setf flag-name (format nil "~D" flag-name)))
  (assert (stringp flag-name))
  (cond
    ((string-equal flag-name "(i)")
     (multiple-value-bind (junk int-val str-val)
         (get-i-register stack)
       (declare (ignore junk))
       (cond
         ((and (not *unlimited-indirection*)
               (or (< int-val 0) (> int-val 3)))
          (error (make-condition 'i-register-range-error
                                 :value int-val
                                 :min-allowed 0
                                 :max-allowed 3)))
         (t
          str-val))))
    (t
     flag-name)))
  

(defun set-flag-by-name (stack name &key clear)
  (let ((record (assoc name (stack-flags stack)
                       :test 'string=)))
    (cond
      (record
       (setf (cdr record) (not clear)))
      (t
       (setf (stack-flags stack)
             (push (cons name (not clear))
                   (stack-flags stack)))))))

(defun get-flag-by-name (stack name)
  (let* ((record (assoc name (stack-flags stack)
                       :test 'string=))
         (rval (cdr record)))
    (when (or (string= name "2")
              (string= name "3"))
      (set-flag-by-name stack name :clear t))
    rval))


(defun set-flag-fcn (stack name &key clear)
  (setf name (canonicalize-flag-name stack name))
  (set-flag-by-name stack name :clear clear))

(defun clear-flag-fcn (stack name)
  (set-flag-fcn stack name :clear t))

(defun get-flag-fcn (stack name)
  (setf name (canonicalize-flag-name stack name))
  (get-flag-by-name stack name))

Both are found in “stack.lisp”.

Several more key operations were coded in “calc1.lisp”.  The statistical operations, some flow-control operations, and a key that affects the way data is presented on the screen of the calculator.

The current code is in the git repository, under the tag v2014-11-02.