guile-core-20020426 and IEEE 754 arithmetic

[Nelson H. F. Beebe]

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.

Among these expectations are the nonstop computing model from allowing
generation of IEEE 754 NaN (not-a-number, e.g., from 0.0/0.0) and
Infinity (e.g., from 1.0/0.0), support for signed zero, and support
for subnormal numbers (when implemented by the underlying
architecture).  

Ultimately, there should also be access to IEEE 754 rounding control
and floating-point status flags, but those two advanced issues are not
being addressed yet in the material below

The extended hoc program cited below is being developed to provide an
interactive testbed, in a simple language that resembles a small
subset of C, for students to learn about floating-point arithmetic in
general, and IEEE 754 arithmetic in particular, as well as to provide
a research vehicle for trying out ideas proposed by the IEEE 754
committee for suitable programming-language primitives.  More details
can be found at my IEEE 754 test software site:

        http://www.math.utah.edu/~beebe/software/ieee

Marius has already reported that changes in progress in
as-yet-unreleased guile code to address some of the issues that we
have been discussing, so please consider this posting a summary of
problems that are likely to be resolved shortly.

(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

(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

This is almost correct: however, the last number should be infinity,
not 0.10.

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

I would prefer +Inf or +Infinity instead of the peculiar +#.#, since
that is what IEEE 754 recommends, C99 requires, and humans expect.

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

guile evidently uses #.# for NaN and +#.# or -#.# for Infinity (see
below).

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

        (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

(In)equality tests for Inf and NaN are done correctly:

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

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

        (format #t "~a  ~a~%" (= Inf Inf) (not (= Inf Inf)))
        #t  #f
        #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

All four results should be #f.  Compare similar experiments in
extended hoc (http://www.math.utah.edu/pub/hoc/) produce correct
output:

        hoc> (NaN == NaN)
        0
        hoc> (NaN != NaN)
        1
        hoc> (Inf < NaN)
        0
        hoc> (Inf <= NaN)
        0
        hoc> (Inf > NaN)
        0
        hoc> (Inf >= NaN)
        0

Basic arithmetic in guile is consistent with hoc:

        guile> (- NaN NaN)
        #.#
        guile> (+ NaN NaN)
        #.#
        guile> (* NaN NaN)
        #.#
        guile> (/ NaN NaN)
        #.#
        guile> (sqrt NaN)
        #.#

        hoc> (NaN - NaN)
        +QNaN
        hoc> (NaN + NaN)
        +QNaN
        hoc> (NaN / NaN)
        +QNaN
        hoc> (NaN * NaN)
        +QNaN
        hoc> sqrt(NaN)
        +QNaN

        guile> (- Inf Inf)
        #.#
        guile> (+ Inf Inf)
        +#.#
        guile> (- (+ Inf Inf))
        -#.#
        guile> (- (- Inf) Inf)
        -#.#

        hoc> (Inf - Inf)
        +QNaN
        hoc> (Inf + Inf)
        +Inf
        hoc> -(Inf + Inf)
        -Inf
        hoc> (-Inf - Inf)
        -Inf

I tried similar experiments in GNU Common Lisp (2.2.2), Clisp 2.28,
and Emacs (21.2) Lisp.  Unfortunately, all catch the zero divide and
abort the computation of NaN and Inf with the usual straightforward
code:

        (setq Inf (/ 1.0 0.0))
        (setq NaN (/ 0.0 0.0))

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!

I did find out that gcl and emacs lisp can print a signed
floating-point zero, but clisp loses the sign.

(sqrt -0.0) and (sqrt (- 0.0)) correctly produce -0.0 with gcl and
emacs lisp, but both fail with clisp: it incorrectly produces 0.0.

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- Center for Scientific Computing       FAX: +1 801 585 1640, +1 801 581 4148 -
- University of Utah                    Internet e-mail: address@hidden  -
- Department of Mathematics, 110 LCB        address@hidden  address@hidden -
- 155 S 1400 E RM 233                       address@hidden                    -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe  -
-------------------------------------------------------------------------------