* Re: guile-core-20020426 and IEEE 754 arithmetic
@ 2002-05-16 13:06 Nelson H. F. Beebe
2002-05-16 15:17 ` John W. Eaton
0 siblings, 1 reply; 5+ messages in thread
From: Nelson H. F. Beebe @ 2002-05-16 13:06 UTC (permalink / raw)
Cc: beebe, jwe, Marius Vollmer
Thanks to Marius Vollmer <mvo@zagadka.ping.de> and John W. Eaton
<jwe@bevo.che.wisc.edu> for the reports of much improved handling of
floating-point output in the current guile development sources!
John asked about this fragment of my report:
>> ...
>> | (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.
>> ...
Apologies for not making my diff of format.scm clearer: the changed
file increased the buffer size from 200 to 400, fixing the problem
with the display of large floating-point numbers.
>> If (- 0.0) is a valid way to create a signed zero,...
Another way to get a negative zero is (/ -1 Infinity): here's a
summary of what my hoc, and guile and some other Lisp systems, do with
this:
hoc-7.0.x:
hoc> x = 0
hoc> __hex(x)
0x00000000_00000000 == 0 0
hoc> y = -0
hoc> __hex(y)
0x80000000_00000000 == 0 0
hoc> y
-0
hoc> z = 1 / -Inf
hoc> __hex(z)
0x80000000_00000000 == 0 0
bigloo (2.5a)
1:=> -0.0
0.0 !WRONG
1:=> (- 0.0)
0.0 !WRONG
1:=> (/ -1 (/ 1.0 0.0))
0.0 !WRONG
guile-core-20020426 (1.5.6):
guile> -0.0
0.0 !WRONG
guile> (- 0.0)
0.0 !WRONG
guile> (/ -1 (/ 1.0 0.0))
0.0 !WRONG
mzscheme 200alpha12 (PLT Scheme):
> -0.0
-0.0 !OK
> (- 0.0)
-0.0 !OK
> (/ -1 (/ 1.0 0.0))
-0.0 !OK
Scheme Release 7.3.1 Microcode Version 11.146 Runtime 14.166 (MIT Scheme):
1 ]=> -0.0
;Value: 0. !WRONG
1 ]=> (- 0.0)
;Value: 0. !WRONG
1 ]=> (/ -1 (/ 1.0 0.0))
;Value: 0. !WRONG
[By the way, that implementation produces #[+inf] and #[NaN] for
Infinity and NaN.]
emacs-21.2.2:
(- 0.0)
-0.0 !OK
-0.0
-0.0 !OK
(/ -1 (/ 1.0 0.0))
-0.0 !OK
clisp 2.28 (released 2002-03-03):
[1]> (- 0.0)
0.0 !WRONG
[1]> -0.0
0.0 !WRONG
[2]> (/ -1.0e-300 1.0e+300)
*** - floating point underflow !BOTCH
[3] (/ -1 (/ 1.0 0.0))
*** - division by zero !BOTCH
gcl 2.2.2:
>-0.0
0.0 !WRONG
>(- 0.0)
-0.0 !OK
>(/ -1.0e-300 1.0e+300)
-0.0 !OK
>(/ -1 (/ 1.0 0.0))
Error: Zero divisor. !BOTCH
A caution for implementers: when I tested hoc on several different
platforms, I found that on a few, it produced 0 instead of -0: that is
due to the underlying C printf() function's defective handling of
negative zero. I will add code in the next hoc release to handle that
special case.
Similar experiments in with 50+ C, 35_ C++, and 30+ Fortran compilers,
produce mixed results: correct handling of IEEE 754 negative zero
remains elusive in those languages, despite that fact that the
arithmetic system is now more than two decades old, and is in
universal use on 100M+ desktop computer systems.
-------------------------------------------------------------------------------
- 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: beebe@math.utah.edu -
- Department of Mathematics, 110 LCB beebe@acm.org beebe@computer.org -
- 155 S 1400 E RM 233 beebe@ieee.org -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe -
-------------------------------------------------------------------------------
_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: guile-core-20020426 and IEEE 754 arithmetic
2002-05-16 13:06 guile-core-20020426 and IEEE 754 arithmetic Nelson H. F. Beebe
@ 2002-05-16 15:17 ` John W. Eaton
2002-05-22 13:36 ` Marius Vollmer
0 siblings, 1 reply; 5+ messages in thread
From: John W. Eaton @ 2002-05-16 15:17 UTC (permalink / raw)
Cc: guile-devel
On 16-May-2002, Nelson H. F. Beebe <beebe@math.utah.edu> wrote:
| guile-core-20020426 (1.5.6):
| guile> -0.0
| 0.0 !WRONG
| guile> (- 0.0)
| 0.0 !WRONG
| guile> (/ -1 (/ 1.0 0.0))
| 0.0 !WRONG
With the following patch, guile produces
guile> -0.0
0.0
guile> (- 0.0)
-0.0
guile> (/ -1 (/ 1.0 0.0))
-0.0
It looks like guile is dropping the sign on input. I would have tried
to fix that problem too, but I couldn't figure out where the string ->
double conversion takes place.
jwe
ChangeLog:
2002-05-16 John W. Eaton <jwe@bevo.che.wisc.edu>
* configure.in (AC_CHECK_FUNCS): Check for copysign.
libguile/ChangeLog:
2002-05-16 John W. Eaton <jwe@bevo.che.wisc.edu>
* numbers.c (idbl2str): Don't omit sign when printing negative zero.
--- configure.in 2002/05/16 15:04:05 1.1
+++ configure.in 2002/05/16 15:03:54
@@ -505,7 +505,7 @@
AC_CHECK_HEADERS(floatingpoint.h ieeefp.h nan.h)
-AC_CHECK_FUNCS(finite isinf isnan)
+AC_CHECK_FUNCS(finite isinf isnan copysign)
# When testing for the presence of alloca, we need to add alloca.o
# explicitly to LIBOBJS to make sure that it is translated to
--- libguile/numbers.c 2002/05/16 15:01:34 1.1
+++ libguile/numbers.c 2002/05/16 15:02:45
@@ -2076,7 +2076,16 @@
int exp = 0;
if (f == 0.0)
- goto zero; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
+ {
+#ifdef HAVE_COPYSIGN
+ double sgn = copysign (1.0, f);
+
+ if (sgn < 0.0)
+ a[ch++] = '-';
+#endif
+
+ goto zero; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
+ }
if (xisinf (f))
{
_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel
^ permalink raw reply [flat|nested] 5+ messages in thread
* guile-core-20020426 and IEEE 754 arithmetic
@ 2002-05-08 19:02 Nelson H. F. Beebe
2002-05-16 5:12 ` John W. Eaton
0 siblings, 1 reply; 5+ messages in thread
From: Nelson H. F. Beebe @ 2002-05-08 19:02 UTC (permalink / raw)
Cc: beebe, Marius Vollmer
guile developer Marius Vollmer <mvo@zagadka.ping.de> 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: beebe@math.utah.edu -
- Department of Mathematics, 110 LCB beebe@acm.org beebe@computer.org -
- 155 S 1400 E RM 233 beebe@ieee.org -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe -
-------------------------------------------------------------------------------
_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel
^ permalink raw reply [flat|nested] 5+ messages in thread
* guile-core-20020426 and IEEE 754 arithmetic
2002-05-08 19:02 Nelson H. F. Beebe
@ 2002-05-16 5:12 ` John W. Eaton
0 siblings, 0 replies; 5+ messages in thread
From: John W. Eaton @ 2002-05-16 5:12 UTC (permalink / raw)
Cc: guile-devel, Marius Vollmer
On 8-May-2002, Nelson H. F. Beebe <beebe@math.utah.edu> wrote:
| guile developer Marius Vollmer <mvo@zagadka.ping.de> 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
_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel
^ permalink raw reply [flat|nested] 5+ messages in thread
end of thread, other threads:[~2002-05-22 13:36 UTC | newest]
Thread overview: 5+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2002-05-16 13:06 guile-core-20020426 and IEEE 754 arithmetic Nelson H. F. Beebe
2002-05-16 15:17 ` John W. Eaton
2002-05-22 13:36 ` Marius Vollmer
-- strict thread matches above, loose matches on Subject: below --
2002-05-08 19:02 Nelson H. F. Beebe
2002-05-16 5:12 ` John W. Eaton
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).