Organization: Penn State University
Date: Wed, 12 Aug 1992 15:00:39 EDT
From: H. D. Knoble <HDK@psuvm.psu.edu>
Newsgroups: comp.lang.fortran
Subject: Re: SUMMARY: precision of real numbers

Two excellent (less formal than ACM, but well qualified) practical papers
dealing with precision of real numbers are published in the SHARE 75
PROCEEDINGS, Volume 1, 14 August 1990, are:

1) John R. Ehrman, "Floating Point Numbers and Arithmetic, a Tutorial".
   This gives algebra of the Real Field of numbers that often do not hold in
   the Floating Point approximations; then the IBM 370 and IEEE Floating Point
   models are compared and many practical examples and consequences given
   along with some very interesting questions about how to make the IEEE
   model better.

2) Herman D. Knoble, "Reality of REAL Arithmetic's Really Wrong Results".
   This gives several examples and analyses of some catastrophic cases.
   Lively teaching examples. Results shown for both IBM 370 and IEEE
   Floating Point Models. Show these to the decision makers who believe
   everything on paper with holes coming down the sides!

Example 1 and 3 in the second paper follow here. They may be useful
if for no other benchmark then one in awareness that program results should
be thoroughly checked before standing behind or believing any of them.

Example 1: For an apparently simple dollars and cents example whose result
           algebraically is 1.0, the result computed using the IBM Floating
           Floating Point model with the suggested input values of A=100.,
           X=0.01 is not 1.0 but rather -39.06... ; using the IBM PC IEEE
           Model the result is: 9.76... .  Please note both results are
           "really wrong."  This is better than close or subtly wrong!
           Try reversing the suggested input values for A and X and look
           at the output, P; the result is P=0.999999...  and 1.00000...
           respectively. Also try A=100. and X=0.0078125; on the IBM 370
           model the result P=0.999999... ; using the IBM PC IEEE model
           the result is P=0.000000... . This program is also useful as a
           teaching tool to wake Computer Science College Seniors up just
           before they graduate! Ask what one would do to prevent and to
           detect problems like this as part of that exercise.

C==== Example 1: Dollars and Cents numeric breakdown; Try A=100., X=0.01.
C     NOTE: P is algebraically 1.0 for all A and X<>0.
C  Penn State University Center for Academic Computing
C  H. D. Knoble, August 1979.
       DO 1 I=1,99999
         WRITE(6,*) 'Please enter A and X (or press Enter to Quit):'
         WRITE(6,*) '(Try values 100., 0.01, then 0.01, 100.)'
         READ(5,*,END=99) A,X
         P=((A+X)**2 - A**2 - 2.*A*X)/X**2
         WRITE(6,100) P,A,X
100      FORMAT(' P=',F14.7,' A=',F14.7,' X=',F14.7)
1      CONTINUE
99     WRITE(6,*) '** Example 1 Ended **'
       STOP
       END

Example 3. Results may be right (or wrong) seemingly by chance, even in the
           same program. Using the IBM 370 Floating Point Model, the result,
           "Swartz's Inequality holds" is output for X=Y=1 to X=Y=2901 by 1.
           Then it "holds" and "fails" rather sporadically (but in a pattern)
           from X=Y=2902. to 8190. by 1. Then it "holds" again from X=Y=8191
           to X=Y=12743 by 1. It then appears to fail sporadically (but with
           a changing pattern) above X=Y=12743. Using Double Precision only
           changes where the "holding" and "failing" begin and end. Using
           the IEEE Floating Point Model on the IBM PC, the result is "holds"
           indefinitely, even well beyond where the limits where the
           associative law of algebra breaks down for the arithmetic of this
           model. So seemingly correct results can be obtained by a very
           shaky program by chance alone; and floating point arithmetic
           number systems are not closed with respect to addition and
           multiplication. That is, there are X's and Y's in every
           computer's floating point system such that X+Y, X-Y, X*Y,
           and/or X/Y do not belong to this system. For example X=90010.
           and Y=20.03125 (0.03125 is 1/32nd) both belong to the IBM 370
           single precision floating point system; but neither their
           sum, difference, product, or ratio do. Banks and (S & L's!)
           do run into this problem trying to divide money into units
           that don't exist; ever wonder in whose favor the results are
           adjusted? Swartz's Inequality failing to hold shows that
           floating point numbers do not obey both the associative and
           distributive laws of algebra. That is, for all computers there
           there are many X, Y, and Z's belonging to the system such that
           (X+Y)+Z .NE. X+(Y+Z) and (X*Y)*Z .NE. X*(Y*Z) and (X+Y)*Z .NE.
           X*Z+Y*Z. This too is a good example for students of numerical
           computations to study; have them count and come up with a
           simple ... notation to enumerate the floating point numbers
           for a specific platform and model and then approximate what
           proportion of that set obey some simple rules of algebra.
C
C==== Example 3, Swartz's Inequality.
C  Example shows while it is true mathematically doesn't necessarily
C  hold in straight-forward Fortran.
C  Penn State University Center for Academic Computing
C  H. D. Knoble, August 1979.
      REAL X,Y, FROM,TO
      WRITE(6,*)
     * 'Please give 2 loop limits to run Swartz''s Inequality thru:'
      WRITE(6,*) '(Try limits like: 1., 10.  2900., 2910.  8190., 8210.'
      READ(5,*,END=99) FROM,TO
      DO 1 X=FROM,TO,1.
        Y=X
        IF (2.*(X**2+Y**2) .GE. (X+Y)**2) THEN
          WRITE(6,100) X,Y
100       FORMAT(' Swartz''s Inequality holds for X=',F6.1,' Y=',F6.1)
        ELSE
          WRITE(6,200) X,Y
200       FORMAT(' Swartz''s Inequality fails for X=',F6.1,' Y=',F6.1)
        ENDIF
1     CONTINUE
99    WRITE(6,*) '** Example 3 Ended **'
      STOP
      END
---------------------------------------------------------End of examples