guile-core-20020426 and IEEE 754 arithmetic

[John W. Eaton]

On  8-May-2002, Nelson H. F. Beebe <address@hidden> wrote:

| guile developer Marius Vollmer <address@hidden> requested that I
| repost a summary of some discussions that we've been having offline
| about build problems of guile-core-20020426 on multiple UNIX
| platforms, and particularly, on its floating-point arithmetic behavior
| and the expectations of how IEEE 754 arithmetic should work on the
| part of the IEEE 754 designers, and the numerical analysis community.

I think guile-core.unstable-20020515 (and presumably the current
sources from CVS) are better.

Here are the results of the examples for the cases you mention as
having problems:

| (1) Output of small numbers is incorrect:
| 
|       (use-modules (ice-9 format))
|       (let ((ratio (- 1.0 (expt 2 -53))))
|         (while (> ratio (/ ratio 2))
|           (format #t "~25,15g~%" ratio)
|           (set! ratio (* ratio 2)))
|         (format #t "~25,15g~%" ratio))
| 
| At Marius' suggestion, this was mostly fixed by this patch:
| 
|       % diff /usr/local/share/guile/site/ice-9/format.scm*
|       1445c1445
|       < (define format:fn-max 400)            ; max. number of number digits
|       ---
|       > (define format:fn-max 200)            ; max. number of number digits

It's not clear from the diff which is the old and new.  If it is
supposed to be 400, then this change appears to be in the newer
sources.

| (2) The little test loop
| 
|       (use-modules (ice-9 format))
|       (let ((ratio (- 1.0 (expt 2 -53))))
|         (while (> ratio (/ ratio 2))
|           (format #t "~25,15g~%" ratio)
|           (set! ratio (* ratio 2)))
|         (format #t "~25,15g~%" ratio))
| 
| now produces:
| 
|       ...
|          8.988465674311730E+307
|          1.797693134862350E+308
|           0.100000000000000
|       #t

With the newer sources:

        ...
           8.988465674311730E+307
           1.797693134862350E+308
             +inf.0              
        #t

| The ~g and ~a formats are not being handled consistently:
| 
|       (format #t "~a  ~a~%" 1.79E+308 (* 2 1.79E+308))
|       1.79000000000003e308  +#.#
|       #t
| 
|       (format #t "~f  ~f~%" 1.79E+308 (* 2 1.79E+308))
|       
179000000000003000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0
  0.1
|       #t
| 
|       (format #t "~25,17g  ~25,17g~%" 1.79E+308 (* 2 1.79E+308))
|        1.79000000000003000E+308    0.10000000000000000
|       #t

With the newer sources:

        (format #t "~a  ~a~%" 1.79E+308 (* 2 1.79E+308))
        1.79000000000003e308  +inf.0
        #t

        (format #t "~f  ~f~%" 1.79E+308 (* 2 1.79E+308))
        
179000000000003000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.0
  +inf.0
        #t

        (format #t "~25,17g  ~25,17g~%" 1.79E+308 (* 2 1.79E+308))
         1.79000000000003000E+308     +inf.0                
        #t

So now we have infinity.  In this context, it might be better to have
the string "Infinity", but I think +inf.0 is clear and has the
advantage of being something that guile can now read as a number.

| Here is another test that shows that ~f and ~g format items are wrong
| for NaN, which must be distinct from Inf, and also from any finite
| value:
| 
|       (format #t "~a  ~a~%"  0.0 (/ 0.0 0.0))
|       0.0  #.#
|       #t
| 
|       (format #t "~f  ~f~%"  0.0 (/ 0.0 0.0))
|       0.0  0.0
|       #t
| 
|       (format #t "~25,17g  ~25,17g~%" 0.0 (/ 0.0 0.0))
|         0.00000000000000000        0.00000000000000000
|       #t

With the newer sources:

        (format #t "~a  ~a~%"  0.0 (/ 0.0 0.0))
        0.0  +nan.0
        #t

        (format #t "~f  ~f~%"  0.0 (/ 0.0 0.0))
        0.0  +nan.0
        #t

        (format #t "~25,17g  ~25,17g~%" 0.0 (/ 0.0 0.0))
          0.00000000000000000         +nan.0                
        #t

| Signed zero is not handled properly by any of ~a, ~f, or ~g:

If (- 0.0) is a valid way to create a signed zero, then this is still
not correct:

        guile> (format #t "~a  ~a~%"  0.0 (- 0.0))
        0.0  0.0
        #t

        (format #t "~25,17g  ~25,17g~%"  0.0 (- 0.0))
          0.00000000000000000        0.00000000000000000    
        #t

        (format #t "~f  ~f~%" 0.0 (- 0.0))
        0.0  0.0
        #t

| Comparisons are suspect however:
| 
|       (format #t "~a  ~a~%" (< Inf NaN) (<= Inf NaN))
|       #f  #t
|       #t
| 
|       (format #t "~a  ~a~%" (> Inf NaN) (>= Inf NaN))
|       #f  #t
|       #t

With the newer sources:

        (define NaN (/ 0.0 0.0))                     
        (define Inf (/ 1.0 0.0))                                         

        (format #t "~a  ~a~%" (< Inf NaN) (<= Inf NaN))
        #f  #f

        (format #t "~a  ~a~%" (> Inf NaN) (>= Inf NaN))
        #f  #f

BTW, we now also have

  nan
  inf
  inf?
  nan?

procedures for creating and testing inf and nan values.

| Attempts to generate Inf by multiplication produced numeric-overflow
| aborts.  Evidently, neither implementation of Common Lisp, nor of
| Emacs Lisp, has been done with an understanding of IEEE 754
| floating-point arithmetic.  Let's make guile/Scheme do it better!

Some people are trying!  :-)

jwe

-- 
www.octave.org        | Unfortunately we were hopelessly optimistic in 1954
www.che.wisc.edu/~jwe | about the problems of debugging FORTRAN programs.
                      |                                       -- J. Backus