Serious bug

[Craig Burley]

>I don't want to replace (x == 0.) by (x < 0. ? (x > -DBL_MIN) : (x <
>DBL_MIN)) for the simple cases (e.g., I assigned 0. to x earlier, did
>some other stuff which may have altered x, and now I want to see if x
>was altered).  I want the compiler to take care of that for me.  My
>experience says that most compilers do.  Don't bother telling me my
>expectations are too high.

I believe your experience is wrong -- I don't recall *ever* hearing of
a compiler for Fortran, C, or C++ that does the kind of substitution
(of range-checking for direct comparison) that you're talking about.

But, your expectations are *not* too high -- if you want Java, or
some other such language that provides these guarantees, go find
and use them, as there's quite a long list of computer languages that
have been developed over the years.  Also, there's always assembly
language, if you want incredibly fine control over the precision
of intermediate operations.  (Watch out for "status flags" that
affect how a given operation is implemented -- the x86 has those
in its floating-point unit.)

However, if you also want the compiler to implement your code at
speeds comparable to what *properly written* numerical code obtains
on published benchmarks, then, yes, your expectations are probably
too high.

>> This means that, in languages like Fortran, C, and, I would guess,
>> C++, it is permitted for evaluation of
>> 
>>   3.1 - 3.1
>> 
>> to return a non-zero number, because the approximations of the two
>> `3.1' constants need not be identical.
>
>Well, that's just silly.  If everyone agrees with that, I'm glad I
>didn't bother to go back to school for any "Computer Science" courses.

I agree it's "silly", but it's those doggone CompSci courses (of
which I attended precisely one, over 20 years ago, before dropping
out ;-) that might also make it clear how difficult it is to
specify much of anything else.  Not that I don't think it'd be
wise for future incarnations of standards for Fortran, C, C++, and
so on, to *try*, at least....

>>         DOUBLE PRECISION R
>>         DOUBLE PRECISION D
>> 
>>         R = 3.1
>>         D = 3.1
>>         IF (R .EQ. D) PRINT *, 'Equal!'
>>         END
>
>Get it?  I'm sorry if my statements didn't make it clear that I am only
>talking about numbers of exactly the same type (including precision).

Your statements included *some* info on that, but not nearly enough.

For example, I still can't tell from the above what you mean by
"numbers of the same type".  My *original* example used numbers of
the same type -- in Fortran, "3.1" is *always* a single-precision
floating-point number.  Therefore, the relevance of the type of
the variable into which that constant is stored *should* be precisely
zero, as long as the variable holds at least as much precision as
the constant (which, in my example, was indeed the case).

Nevertheless, while g77 obeys the implications of the above, many
Fortran compilers don't, because *users* are accustomed to compilers
performing what amount to algebraic, or mathematically correct,
transformations being performed -- transformations that are *not*
computationally consistent.

In other words, "D = 3.1" becomes, in some dialects of Fortran, equivalent
to "D = 3.1D0", so the numbers being compared are *not* of the same
precision, so "Equal!" is not printed.  (This violates the spirit of
the Fortran standards, IMO, but not the letter -- and many people are
"mad" at g77 for not [yet] supporting this sort of violation!)

As a further example, do you believe that "3.1 * 1." and "31. * .1"
are "numbers of the same type [and precision]"?  If so, then you
really are going to continue to run into trouble with all sorts of
compilers, though more immediately with some than with others,
because those two expressions do *not* necessarily produce computationally
identical results, even given "helpful" compilers.

Further yet, is "3.1 * 1." the exact same number as "3.1 * R", where
R is the same type/precision as 3.1 and contains "1."?  You might say
yes, but the answer is "not necessarily" if the means by which R
is said to contain "1." results from a series of computations that
produces an *approximation* of 1..  Even if printing out R just before
the comparison printed "1.", the comparison might still fail, because
a slightly different approximation method might apply to the version
of R used for the comparison than for the printing in decimal!
(E.g. "R = 3. * (1. / 3.)" might produce this effect.)

No widely used computer language I'm aware of provides precise control
over the intermediate results of subexpressions, perhaps even expressions,
perhaps even variables, and so on.

And, numeric constants are generally best viewed as highly compact
expressions of operations such as multiplication, addition, and
division, in this context, though some languages might so tightly
nail down exactly how constants are evaluated that this view is
not necessary.

So, even "3.1" is most "safely" viewed as an in-line expression
such as "(3. + (1. * .1))", which the compiler is (usually) free
to optimize however else it'd be free optimize such an expression.

Further, most computer languages either *allow* algebraic substitution
by the compiler, or, if popular (like Fortran), are implemented *as
if they did* (e.g. one or more Sun Fortran's do this, from what I
have heard).

Thus, the statements

	D1 = 3.1
	D2 = 3.1
	IF (D1 .EQ. D2) ...

can be substituted by such implementations *algebraically* as:

	D1 = 3.1
	IF (D1 .EQ. 3.1) ...

Further, the 3.1 in the second statement can be implemented using
more precision than is contained by D1.  That would make the
comparison fail every time.  (It seems this would violate the
spirit of the Fortran standard as well, but I'm not so sure.)

That some compilers might happen to employ some ad-hoc rules to
avoid "obvious" problems like the above does not necessarily mean
they'll consistently do what you want all the time.  Which means
that, the more you think you're using a compiler that helps you
avoid such problems without your having to "think" about how
to really use floating-point, the more likely you'll write code
that exhibits subtle numerical bugs when it runs, especially when
compiled by that compiler.  (Other compilers might tend to expose
the bugs in your code in ways that aren't so subtle.)

It is *possible* to architect a language standard, or a particular
implementation (compiler &c), to avoid *most* of these problems,
but the result will be a system that runs most numerical codes much
slower than any other system, especially on hardware like Intel x86.

Note that Fortran compilers supporting interval arithmetic might
become available soon.  Maybe C/C++ compilers will, someday, as well.
You might want to wait for those, if you want the compiler to do your
range-checking for you.  E.g., in Interval Fortran, [.1] is an
interval around the real constant .1, in that the lower and upper
bounds, viewed as precise floating-point values, contain the real
number .1 in the range they form.  (But note that this means
[3.1415926535] does not necessarily contain pi -- to enclose an
irrational using this form, one must hand-twiddle the final digits,
as in [3.1415926535,3.1415926536].)

There's all sorts of literature available on the dangers of using
floating-point, the vagaries of various language and compiler
implementations, and so on.  "Kahan" is a name that comes to mind.

As far as egcs, gcc, and so on -- there *is* a tension between providing
seamless, consistent behavior across platforms (which GNU, as a
project and UNIX offshoot, generally tries hard to provide) and
providing performance and behavior consistent with the expectations
of users of a particular platform.

And, for floating-point arithmetic, the "community" as a whole tends
to want this tension resolved in favor of performance over consistency.
That's consistent with how things like integer exceptions and pointers
(widths, etc.) are handled, but inconsistent with most everything else
(like what extensions are supported, which, prior to GNU, was usually
based on the particular machine/OS combination you were using).

Personally, I'm not entirely happy with that.  I'd like g77 to support
full IEEE 754 behavior on all machines regardless of underlying iron,
if I had my druthers, but that'd make g77 really slow on many machines,
and still wouldn't (necessarily) make it meet *your* requirements.
Still, I prefer consistency and correctness over speed, especially as
someone maintaining the product (so the "experts" who know they want
performance end up being self-selecting, in that they'd know to use
the options that meant "I know what I'm doing, go for speed instead
of the canonical GNU floating-point model").

Until this changes, however, you do have to be aware of the details
of floating-point behavior of the implementation of the language
you're using at any given time, and this includes what options to
use (e.g. on the Alpha, you might want `-mieee') when compiling,
linking, and so on, plus what the underlying hardware does (to the
extent the software exposes the details of that implementation).

And, I'm not sure offhand that any prominent vendor of C/C++ compilers
provides more consistency across a range of hardware architectures
than do egcs' gcc/g++.  (I think one or two Fortran compilers might
do better than g77 in that respect, but I don't think Digital Fortran
is an example of this -- I get the impression Digital Fortran defaults
to giving up full IEEE 754 on Alphas, moreso than it does on x86,
but probably still silently uses excess precision on x86 where it
won't on Alphas, all to get good performance by default.)

In the meantime, storing every single intermediate calculation through
a variable, and compiling with `-ffloat-store' on systems like x86
and m68k, might *effectively*, if not definitively, give you some of the
behaviors you seem to want with the egcs/gcc product line.

        tq vm, (burley)