[comp.lang.fortran] style--- REAL comparisons

mcdonald@aries.scs.uiuc.edu (Doug McDonald) (08/05/90)

I just thought of one good - very good - use of compariaons of
real variables for equality:

THIS IS PSEUDOCODE, folks, - bugs guaranteed, you get the idea :

       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)

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

Doug McDonald

ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe) (08/05/90)

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)

ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe) (08/06/90)

In article <3516@goanna.cs.rmit.oz.au>, ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe) writes:
> 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:

> People who say NEVER compare floats for equality are being realistic.

I forgot to mention another problem.  Suppose you compile the subroutine
with version 18.5 of your compiler, and then compile the main program that
calls it with version 19.1.  One of the reasons for the upgrade may be that
the new version does a better version of converting floating-point literals...

A Fortran compiler which claims to conform to the IEEE 754 and/or IEEE 854
standards must get literals like 0.5, 0.625, and so on right.  But without
that claim of conformance there is no guarantee.  And in past times it just
wasn't _possible_ to send the compiler back if it didn't get 0.5 exactly
right; you often had the choice of the vendor's Fortran compiler or _none_.

I have been seriously flamed for pointing out that McDonald's example is
a good example of why comparing floats for equality is risky.  If you are
Herman Rubin, totally in control of the machine, then it is safe for you
to compare floats for equality.  (It wasn't he who flamed me,  I merely
cite him as the best known example of someone who might understand his
machine well enough to get away with this kind of thing.)  But if you
expect your code to be portable, or if you haven't got a couple of years
of numerical analysis under your hat, don't be silly.

-- 
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)