[comp.lang.fortran] IEEE 754 vs Fortran arithmetic

burley@world.std.com (James C Burley) (10/24/90)

In article <9010230628.AA22160@admin.ogi.edu> John Roberts <robertsj@admin.ogi.edu> writes:

   Does anyone know of any references that compare or contrast IEEE 754
   floating point arithmetic with other standards, such as Fortran?

   [... I know of no reason that an IEEE implementation of F77 would be
   nonconforming. -John]

I don't know any references, but I do know we ran into this problem
implementing Fortran on a machine using an IEEE 754 math chip:

      REAL R(...)
      DATA R/0.5,1.5,2.5,3.5,.../
      DO I=1,...
         PRINT *,NINT(R)
      END DO
      END

Fortran specifies that the following values must be output:

1, 2, 3, 4,...

However, the IEEE 754 defines nearest-integer so that using its function
instead of Fortran's definition of NINT produces:

0, 2, 2, 4,...

This is because Fortran literally specifies that NINT(X) is INT(X+0.5)
(assuming for convenience here that X>=0), while IEEE 754 specifies that
nearest-integer is to round to the nearest EVEN integer.  (Perhaps it allows
alternatives -- my memory is vague here -- but those alternative have to do
with rounding to 0, infinity, and so on, and somehow I think they didn't come
into play with round-to-nearest-integer as a function and in any case there
didn't seem to be any rounding mode equivalent to Fortran's definition of
NINT.)

Also, Fortran specifically prohibits zero from being negative (or being
significantly negative -- for example, I think, given that X is 0.0 and Y is
-0.0, IEEE 754 specifies that (Y.LT.X) whereas Fortran specifies that
(Y.EQ.X)).  I might be wrong about the specifics, but there is some issue
here even if I can't remember just what it is!  (Love these definitive
pronouncements, don't y'all?!  :-)

James Craig Burley, Software Craftsperson    burley@world.std.com
[754 says that -0 and +0 are equal, but the point about NINT is good, one
can't implement NINT with a single instruction unless the rounding mode is
set to chop. -John]
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

henry@zoo.toronto.edu (Henry Spencer) (10/24/90)

> [our moderator writes]
>... I know of no reason that an IEEE implementation of F77 would be
>nonconforming. ...

I can think of at least one:  F77 flatly denies the existence of -0,
while IEEE demands it.  (One of Dr. Kahan's favorite examples in his
talks is an algorithm which, when implemented straightforwardly, does
the right thing if -0 is implemented properly and screws up bizarrely
if not, so yes, it does matter.)
-- 
Henry Spencer at U of Toronto Zoology, henry@zoo.toronto.edu   utzoo!henry
[Given that +0 = -0, it's not clear to me that the existence of -0 breaks
anything.  Keep in mind that F77 is a permissive standard, extensions are
permitted so long as conforming programs do the right thing. -John]
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

tim@ksr.com (Tim Peters) (10/25/90)

> [Given that +0 = -0, it's not clear to me that the existence of -0
> breaks anything.  Keep in mind that F77 is a permissive standard,
> extensions are permitted so long as conforming programs do the right
> thing. -John]

John, the main problem here is that F77 specifically forbids formatted
output from printing "-0" (pg 13-9, lines 3-5), and while 754/854 never
actually say a -0 has to print as "-0", everyone believes they do
<grin>.  I don't think -0 causes any other problems, except for the
inevitable curse of coming to rely on non-portable semantics.

But there are other problems.  E.g., the meaning of

	X.EQ.Y

is *defined* by F77 to be the value of

       ((X-Y).EQ.0)

(pg 6-10, lines 4-14), and Fortran's usual "must respect parens" rules
prevent a conforming processor from transforming the subtraction into a
"mathematically equivalent" comparison operation (although every vendor
with a comparison instruction cheats here ...).  If X and Y are, e.g.,
both +infinity, following the F77 std here would at best signal an
invalid operation exception and yield a nonsense result, while 754/854
insist that the outcome be an exception-free .TRUE..

Perhaps the most widespread "problem" is that 754 vendors tend to like
to evaluate entire expressions in an extended precision and cut back the
result to "storage precision" only at the end.  This violates my reading
of both the F77 and 754 stds.

None of the Fortran problems will be repaired by F90.

while-the-marriage-may-not-have-been-made-in-hell-that's-where-the-
   honeymoon-is-taking-place-ly y'rs  - tim
   

Tim Peters   Kendall Square Research Corp
tim@ksr.com,  ksr!tim@harvard.harvard.edu
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

dik@cwi.nl (Dik T. Winter) (10/25/90)

In article <BURLEY.90Oct24025053@world.std.com> burley@world.std.com (James C Burley) writes:

About Fortran NINT and IEEE NINT.  That is correct, but it only tells us
that Fortran NINT does not map on IEEE NINT, so if you implement Fortran
NINT you might do it exactly as the Fortran standard says using IEEE
primitives.

Further:
 > Also, Fortran specifically prohibits zero from being negative (or being
 > significantly negative -- for example, I think, given that X is 0.0 and Y is
 > -0.0, IEEE 754 specifies that (Y.LT.X) whereas Fortran specifies that
 > (Y.EQ.X)).  I might be wrong about the specifics, but there is some issue
 > here even if I can't remember just what it is!

In the Ada Numerics Working Groups in Europe and in the US/Canada there is
also discussion going on about -0.0 and +0.0 (Ada does not distinguish them
like Fortran).  However, in a full implementation of IEEE including the
recommendations the only way to distinguish -0.0 and +0.0 is by either
calculating 1/X or by use of the copysign function.  Both ways do not
conflict with the languages (in Fortran and Ada calculating 1/X for X=zero
is undefined or defined to trap, depending on parameters, the Fortran
copysign function is undefined for a zero sign argument).

As Henry Spencer notes, Kahan has good arguments for -0.0 and +0.0 (not too
surprising as he was one of the leading forces behind IEEE 754).  On the
other hand in his use -0.0 stands for negative infinitiseminals and +0.0
for zero and positive infinitiseminals.  In my opinion (and I am not the
only one) it would have been better to have a single zero and two signed
infinitiseminals.  (And as an aside, there have been machines that had only
representations for signed infinitiseminals but not for zero, still not such
a very bad idea.)  But that is taking us quite a bit from compilers.

What is relevant with respect to compilers (or the run-time system writers,
but must we distinguish?) is the following.  And this is why the point
came up in the Ada Numerics Working Groups.  We are working on the
standardization of elementary mathematical functions in Ada.  The current
status is that the basic package of functions like SINE, COSINE etc. is
very near to standardization.  One of the functions included is:
	ARCTAN(Y, X)
which returns the arctangent of Y/X.  (Fortran users will recognize the
ATAN2 function.)  The specification tells us that the result is the range
from -PI to +PI (approximately).  The problem is, what is the result of
ARCTAN(Y, zero).  Does it depend on the sign of zero?  Offhand I do not
know what the Fortran standard tells us.

Summary:  IEEE 754 does not constrain languages that refuse to distinguish
the sign of zero, but it gets hairy if some mathematical functions are
considered.  And of course, if you want all possibilities Kahan wants
the language needs to know about inifinities and NaNs.
--
dik t. winter, cwi, amsterdam, nederland
dik@cwi.nl
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

wsb@eng.Sun.COM (Walt Brainerd) (10/25/90)

In article <BURLEY.90Oct24025053@world.std.com>, burley@world.std.com (James C Burley) writes:
> 
> I don't know any references, but I do know we ran into this problem
> implementing Fortran on a machine using an IEEE 754 math chip:
> 
>       REAL R(...)
>       DATA R/0.5,1.5,2.5,3.5,.../
>       DO I=1,...
>          PRINT *,NINT(R)
>       END DO
>       END
> 
> Fortran specifies that the following values must be output:
> 
> 1, 2, 3, 4,...

We have had this discussion before some here, but to be a bit nit-picking,
"Fortran" (i.e., the standard) does not specify such things, as it
does not even specify what + means.  It certainly does encourage
such things and a vendor must be prepared to answer to the customer,
but not worry about strict standard conformace in this area.
> 
> However, the IEEE 754 defines nearest-integer so that using its function
> instead of Fortran's definition of NINT produces:
> 
> 0, 2, 2, 4,...
> 
> Also, Fortran specifically prohibits zero from being negative (or being
> significantly negative

The appendix (not a legal part of the standard) says:

"A processor must not consider a negative zero to be different
from a positive zero."

I would take this as a SUGGESTION to make 0 and I-I compare true,
however the result of the subtraction is represented.

The standard (13.5.9) does say that a "negative signed zero"
must not be produced when doing numeric output into a formatted record.
So, from the point of view of the standard, printing 4-4 as 7 is OK,
but printing it as -0 (with say an I2 format) is not!
It's just a matter of who you run to if something doesn't work
the way you think it should.
--
Walt Brainerd        Sun Microsystems, Inc.
wsb@eng.sun.com      MS MTV 5-40
                     Mountain View, CA 94043
                     415/336-5991
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

eggert@twinsun.com (Paul Eggert) (10/25/90)

>[Given that +0 = -0, it's not clear to me that the existence of -0 breaks
>anything.  Keep in mind that F77 is a permissive standard....  -John]

Not permissive enough, unfortunately: IEEE's -0.0 breaks F77's input/output.
ANSI X3.9-1978 says that -0.0 and +0.0 must print the same way.  IEEE 754-1985
says that -0.0 and +0.0 are distinct, and that printing any number and reading
it back in must yield the original number if the proper precision is used.
These provisions contradict each other.  An implementation cannot conform to
both standards without unnatural contortions.  I suspect most F77
implementations simply ignore the IEEE requirement here.
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

wsb@eng.Sun.COM (Walt Brainerd) (10/26/90)

In article <9010242205.AA04208@lunch.ksr.com>, tim@ksr.com (Tim Peters) writes:
> 
> But there are other problems.  E.g., the meaning of
> 
> 	X.EQ.Y
> 
> is *defined* by F77 to be the value of
> 
>        ((X-Y).EQ.0)
> 
> (pg 6-10, lines 4-14), and Fortran's usual "must respect parens" rules
> prevent a conforming processor from transforming the subtraction into a
> "mathematically equivalent" comparison operation (although every vendor
> with a comparison instruction cheats here ...).

The referenced lines begin:

"If the two arithmetic expressions are of different types, ..."
                                           ^^^^^^^^^^^^^^^
so this does not apply if X and Y both are real, for example.
This text is there to explain how to do type conversion when
X and Y are different types.  Thus, I don't think compilers that
use a relational operator are "cheating".

> Perhaps the most widespread "problem" is that 754 vendors tend to like
> to evaluate entire expressions in an extended precision and cut back the
> result to "storage precision" only at the end.  This violates my reading
> of both the F77 and 754 stds.

Since there is nothing in the F77 standard that indicates how much precision
should be used for anything, such an evaluation cannot violate F77. (IMHO)
Yes, there are some words that indicate when you should convert from one
precision to another, but there is nothing that says how much any particular
precision there should be for anything (there are requirements on how
equivalence must work between them, but this only affects how stored results
must align, not expression evaluation and not the values stored).

> None of the Fortran problems will be repaired by F90.

The real "problem" is: according the the F77 _standard_, there is no way to
be sure that real arithmetic is done by any particular scheme, so any program
that tries to depend on a specific result (i.e., comparing two reals for
equality) or expecting the result of any real operation to be precise within
certain limits, is nonportable.  You certainly may have your own criteria
about how things should work and pick a vendor that does things that way, but
there is no help for you in the standard.

The only reasonable "repair" I can think of is to require all arithmetic to
be done as IEEE 754, or some such, a position that clearly is not acceptable
to all vendors this year.

The new KIND mechanism in Fortran 90 will at least allow the programmer to
indicate that a variable has a certain number of digits of precision; I think
this should help some.
--
Walt Brainerd        Sun Microsystems, Inc.
wsb@eng.sun.com      MS MTV 5-40
                     Mountain View, CA 94043
                     415/336-5991
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

sjc@key.COM (Steve Correll) (10/26/90)

In article <BURLEY.90Oct24025053@world.std.com>, burley@world.std.com (James C Burley) writes:
> Also, Fortran specifically prohibits zero from being negative (or being
> significantly negative -- for example, I think, given that X is 0.0 and Y is
> -0.0, IEEE 754 specifies that (Y.LT.X) whereas Fortran specifies that
> (Y.EQ.X)).  I might be wrong about the specifics...

IEEE 754 says that -0 and +0 compare equal.

However, IEEE 754 introduces a host of problems because of NaNs. Compilers
which assume that .not.(x.gt.y) implies x.le.y, or that .not.(x.eq.y) implies
x.ne.y may founder on the IEEE requirement that NaN compares unordered (false)
with any other value. And the Fortran arithmetic IF is a particularly strange
case, since NaN is neither less than, nor equal to, nor greater than zero.
Thus IEEE technically requires three tests rather than two to implement this.
-- 
sjc@key.com or ...{sun,pyramid}!pacbell!key!sjc 		Steve Correll
[It is an interesting metaphysical problem to consider what an arithmetic IF
should do when confronted with a NaN. Fall through? Ugh. -John]
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

diamond@tkov50.enet.dec.com (diamond@tkovoa) (10/26/90)

In article <1990Oct25.010604.4796@twinsun.com> eggert@twinsun.com (Paul Eggert) writes:
>ANSI X3.9-1978 says that -0.0 and +0.0 must print the same way. IEEE 754-1985
>says that -0.0 and +0.0 are distinct, and that printing any number and reading
>it back in must yield the original number if the proper precision is used.
 
Unfortunately, it is possible.  For example, "proper precision" could be
specified as leaving room for +0.0 to print as 0.0 while -0.0 prints as
(I can hardly type it)
00.0
(oh I feel ill).
--
Norman Diamond, Nihon DEC    diamond@tkov50.enet.dec.com
                                    (tkou02 is scheduled for demolition)
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

bsy+@CS.CMU.EDU (Bennet Yee) (10/28/90)

In article <2408@charon.cwi.nl>, dik@cwi.nl (Dik T. Winter) writes:
|> ... We are working on the
|> standardization of elementary mathematical functions in Ada.  The current
|> status is that the basic package of functions like SINE, COSINE etc. is
|> very near to standardization.  One of the functions included is:
|> 	ARCTAN(Y, X)
|> which returns the arctangent of Y/X.  (Fortran users will recognize the
|> ATAN2 function.)  The specification tells us that the result is the range
|> from -PI to +PI (approximately).  The problem is, what is the result of
|> ARCTAN(Y, zero).  Does it depend on the sign of zero?  Offhand I do not
|> know what the Fortran standard tells us.

Is there something I'm missing?  As far as I knew,

	ARCTAN(Y,0) =  pi/2 if Y > 0 
	              -pi/2 if Y < 0

regardless of the sign on the zero.  If Y = 0, then it's not well defined.
(I don't recall what 754 has to say about this case.)

Perhaps the problem occurs when you compute

	ARCTAN(0,X)

when X<0?  Assuming that you want the range to be from -pi to +pi, you can
argue that the answer should depend on the sign on 0:  if the zero is
negative, then use -pi; otherwise use pi.

In some sense, this isn't really a problem with signed zeros but rather more
a problem with the discontinuity in the definition of the inverse
tangent....

Bennet S. Yee, +1 412 268-7571
School of Cucumber Science, Cranberry Melon, Pittsburgh, PA 15213-3890
Internet: bsy+@cs.cmu.edu		Uunet: ...!seismo!cs.cmu.edu!bsy+
Csnet: bsy+%cs.cmu.edu@relay.cs.net	Bitnet: bsy+%cs.cmu.edu@cmuccvma
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.

mcohen@amsaa-seer.brl.mil (Marty Cohen) (10/30/90)

Does IEEE 754 give a purpose for signed zero?
I.E., is it to result from underflow?
If so, then the -pi, +pi result would be appropriate for atan(-+0,x).
--
Marty Cohen mcohen@brl.mil  {uunet|rutgers}!brl!mcohen
Custom House Rm 800, Phila. PA 19106 (215)597-8377
[IEEE 754 doesn't give much of a rationale for anything, it's pretty short.
-John]
-- 
Send compilers articles to compilers@esegue.segue.boston.ma.us
{ima | spdcc | world}!esegue.  Meta-mail to compilers-request@esegue.