Path: utzoo!attcan!utgpu!news-server.csri.toronto.edu!mailrus!uwm.edu!zaphod.mps.ohio-state.edu!samsung!munnari.oz.au!goanna!ok
From: ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe)
Newsgroups: comp.lang.fortran
Subject: Re: style--- REAL comparisons
Message-ID: <3516@goanna.cs.rmit.oz.au>
Date: 5 Aug 90 13:04:49 GMT
References: <1990Aug4.173025.15920@ux1.cso.uiuc.edu>
Organization: Comp Sci, RMIT, Melbourne, Australia
Lines: 54

In article <1990Aug4.173025.15920@ux1.cso.uiuc.edu>, mcdonald@aries.scs.uiuc.edu (Doug McDonald) writes:
> I just thought of one good - very good - use of comparisons of
> [floating-point] variables for equality:

>        real function power(a,b)
>        if (b .eq. 1.) then
>            power = a
>        else if (b .eq. 0.) then
>            power = 1.
>        else if (b .eq. 0.5) then
>            power = sqrt(a)
>        else
>            (use log and exp)
>	 end if

> There ARE good uses for such things.  People who say NEVER use
> equality comparisons on reals are just being silly.

(a) I just explained in comp.lang.c that unless the rest of your power
    function is very good indeed, this kind of non-uniformity is likely
    to make power(a,b) non-monotone in b, which can do nasty things to
    some numerical algorithms...  For example, you might get
	x .gt. 1.0 .and. power(x, 1.00000000001) .lt. power(x, 1.0)
    for some values of x.
(b) Given that Fortran already has X**N, it seems pointless to reinvent it.
(c) You are assuming that 0.5 is converted exactly.  It *may* be, but what
    if the compiler calculates it as FLOAT(5) * 10.0**(-1)?  Then there is
    going to be round-off error in the computation of 0.1, and the relative
    error in 0.5 is going to be the same (assuming FLOAT(5) is exact), and
    0.5 .eq. 1.0/2.0 will fail.  Which means that someone who knows about
    this and writes power(x, 1.0/2.0) in order to get your magic sqrt() is
    going to be disappointed.  All of this is why programs like Paranoia
    and the Elefunt tests do stuff like this:
	One = FLOAT(1)
	Two = One + One
	Half = One / Two
   This isn't guaranteed either, but it's rather safer, and they're just
   trying to get good approximations, not things they can test.

Consider also
	READ (5, *) X, Y
	WRITE (6, *) POWER(X, Y)
where the input has
	4.0 0.5
There is no guarantee that _this_ 0.5 will be converted the same way that
0.5 occurring in the program will, so despite the right magic literal
appearing, there is no reason to expect that your SQRT special case will
be selected.

People who say NEVER compare floats for equality are being realistic.
-- 
Distinguishing between a work written in Hebrew and one written in Aramaic
when we have only a Latin version made from a Greek translation is not easy.
(D.J.Harrington, discussing pseudo-Philo)