floating point comparisons illegitimate?

[Tim Peters]

>  [PLS]
>  ...
>  If a compiler sees
>    if (1.0 + EPS .EQ. 1.0)
>  is it justified in optimizing it to
>    if (EPS .EQ. 0.0)
>  If not, why not.

Since that's a Fortran fragment, you need a Fortran lawyer, and you get
three answers:

1) Technical answer #1:  If EPS is any type but REAL, no it can't.  When
   X and Y have different types, F77 defines the meaning of
	X .relational. Y
   to be the meaning of
	((X) - (Y)) .relational. 0
   So if, e.g., EPS is DOUBLE PRECISION, then the type of 1.0+EPS is
   also DOUBLE, while the type of 1.0 is REAL, so this clause applies
   and the meaning of
	(1.0+EPS .EQ. 1.)
   is the meaning of
	((1.0+EPS) -  1.)) .EQ. 0
   F77 explicitly prohibits optimizations across parentheses, so a
   compiler must not change this to
	(EPS) .EQ. 0

1a) Gloss to #1:  By most accounts, the answer to #1 was not X3J3's
    *intent* -- they were really just trying to hint at how type
    conversion gets done in comparisons without actually bothering to
    spell it out.  Nevertheless, it's what F77 says.

2) Technical answer #2:  If EPS is REAL, all the std says is that the
   compiler is free to evaluate any expression "relationally equivalent"
   to "1.+EPS.EQ.1", where "relationally equivalent" is undefined, but
   where I believe most reasonable people would take it as allowing
   evaluating "EPS.EQ.0" instead.

3) Practical answer:  Of course not -- no Fortran compiler jockey worth
   their salt would impose a marginal optimization they know will break
   countless programs.  That's a "quality of implementation" kind of
   issue that has to be decided on the basis of the Fortran culture, and
   the std is no help in such matters.  The C culture may well come to
   different answers (indeed, if that was a C fragment I believe X3J11's
   "as if" rule would explicitly prohibit the optimization -- but the
   Fortran world doesn't have-- and wouldn't abide! --an optimization- 
   inhibiting rule that severe).

>  The C standard chose to make information about precision and exponent
>  range available through the environment.  In part because tricks like
>  the above don't work on a lot of machines.

F90 too -- and I think all languages are moving in this direction.  Good
thing, since most people who think they need this info don't know how to
write code to compute it correctly anyway <grin/sigh>.

>  ...
>  A couple of the previous messages have said, in effect, that the IEEE
>  standard requires a minimum precision implementation.  Does it?

I'd say yes & no.

Yes, because 754/854 define every bit of the results required for every
combination of operation, operands, and various options for rounding
mode & "precision control", and these definitions generally require that
if X and Y are the same type, X op Y is of that type too & is fully
defined by knowing X, Y, op, and the values of the options (i.e., the
result is fully defined without reference to the context (if any) in
which "X op Y" appears).

No, because 754/854 say nothing about how language constructions map to
IEEE features.  So, e.g., if your compiler does the entire RHS of

        sometype t, x, y, z
	t = x + y + z

in an extended (wrt sometype) precision, it's not a violation of the
754/854 stds unless your compiler docs *claim* that its "+" operator
maps directly into the sometype flavor of IEEE's addition operation.
In other words, 754/854 just define an arithmetic, not a way to get at
it.

As the traffic on this list demonstrates, there's no consensus on "how
to get at it" yet; defining a language binding may be harder than
defining a good arithmetic ...

except-in-assembler<grin>-ly y'rs  - tim


Tim Peters   Kendall Square Research Corp
timaksr.com,  ksr!timaharvard.harvard.edu