ok@quintus.UUCP (Richard A. O'Keefe) (01/18/88)
Formerly, C treated parentheses according to mathematical practice:
they had only a grouping effect, and whether (a+b)-c was reordered
or not was up to the compiler.
FORTRAN and some other languages give parentheses an operational
effect as well. That is, (a+b)-c and a+b-c are different expressions
and the programmer is entitled to assume that the calculation is as
written, not (a-c)+b or any other re-ordering.
The ANSI C committee, realising that control over the evaluation of
floating-point calculations was desirable, suggested adding an
identity function to C: +(expr) would have the same value as (expr)
but would not allow re-ordering across the "+". This is less of a
kludge than the FORTRAN method, because it would let the reader see
when the ordering mattered (a unary "+" would flag this case) and
when parentheses were just for syntax. The snag with the FORTRAN
method is that you get strict ordering whether you want it or not.
Apparently, a lot of people wrote to the ANSI C committee said that
they would write numerical software in C, only if it looked like
FORTRAN. (Sorry, only if it had strict parens.)
The point of this message is that strict parentheses do not provide
sufficient control over floating-point arithmetic. Something stronger
is needed. This is just as true in Fortran as in C. (In fact I met
the problem in Fortran several years before I first heard of C.)
The problem is this:
Several machines support higher precision in their floating
point registers than in storage. The PR1ME superminis are
an example. Floating-point coprocessors with an 80-bit
format (such as the MC68881) are another.
This means that
DOUBLE PRECISION TEMP, X, Y | double temp, x, y;
TEMP = X+Y | temp = x+y;
IF (TEMP .LT. 0.0D0) CALL HOME | if (temp < 0.0) home();
and
IF (X+Y .LT. 0.0D0) CALL HOME | if (x+y < 0.0) home();
may behave differently. Putting parentheses around X+Y doesn't
cure the problem, because this tends to be compiled as
load X into FP reg [converts to high precision]
add Y to FP reg [result is in high precision]
test sign of FP reg [high precision comparison]
whereas the version which stores into TEMP would have rounded or
truncated. Basically, any floating-point comparison with an operand
which is neither a constant, nor a variable, nor a function call
beware of intrinsics!) is suspect. In fact, function calls in C
may not be safe either (they should be OK in FORTRAN, as defining
the result looks like a store into memory).
The answer is obvious: you have to say exactly what precision
you want. FORTRAN and C both permit this.
FLOAT(X) | (float)(x)
DBLE(X) | (double)(x)
specify that X is to be converted to single or double precision,
respectively. Since FLOAT(I) converts an integer UP to single
precision and FLOAT(D) converts a double DOWN to single precision,
you would think that FLOAT(X) means "convert X to *exactly*
single precision", and DBLE(X) means "convert X to *exactly*
double precision". So you'd write
IF (DBLE(X+Y) .LT. 0.0D0) CALL HOME
or if ((double)(x+y) < 0.0) home();
in the assurance that since you had now specified the precision
you wanted, that's what you'd get.
'Tisn't so.
In the compilers I've tried, if an expression is already in double
[*really* "extended"] precision, DBLE(expr) or (double)(expr)
leave the intermediate result unchanged instead of converting it as
would be done for storing in memory.
I haven't got a copy of the FORTRAN standard, so I don't know whether
this is legal in FORTRAN, and the ANSI C standard isn't finished yet.
It seems to me that FLOAT()/(float) and DBLE()/(double) should mean
"convert the operand to exactly the precision that would be given
to a REAL {DOUBLE PRECISION} variable in storage."
Note that only converting an expression which already appears to
have the desired type is changed, and then only if intermediate
results are calculated in a higher precision. Most programs should
notice no change.
Without this definition in the standard, you have to explicitly
introduce temporary variables, as shown above. As things stand, all
FLOAT(float) or DBLE(double) are good for is providing strict
parentheses...
I haven't had this problem myself, for the simple reason that a
friend of mine had so much trouble with it that it has always stuck
in my memory. But it hits, for example, the subroutine MACHAR in
the Cody & Waite ELEFUNT package.
Can anyone think of a reason why FLOAT & DBLE should *not* have the
effect I suggest?ags@j.cc.purdue.edu (Dave Seaman) (01/21/88)
In article <546@cresswell.quintus.UUCP> ok@quintus.UUCP (Richard A. O'Keefe) writes: >Formerly, C treated parentheses according to mathematical practice: >they had only a grouping effect, and whether (a+b)-c was reordered >or not was up to the compiler. >FORTRAN and some other languages give parentheses an operational >effect as well. That is, (a+b)-c and a+b-c are different expressions >and the programmer is entitled to assume that the calculation is as >written, not (a-c)+b or any other re-ordering. >... The snag with the FORTRAN >method is that you get strict ordering whether you want it or not. It is not true that FORTRAN gives you strict ordering. Taking your example, if a FORTRAN programmer writes (a+b)-c, then he is entitled to assume that it is interpreted as one of the following: (a+b)-c (b+a)-c -c+(a+b) -c+(b+a) but not as: (a-c)+b (b-c)+a (-c+a)+b (-c+b)+a b+(a-c) a+(b-c) b+(-c+a) a+(-c+b) The C programmer, on the other hand, cannot rule out any of the interpretations in the second group. The difference is crucial because floating point arithmetic is commutative but not associative. Therefore the FORTRAN programmer knows that the value of the expression cannot be changed by the compiler's evaluation order, while the C programmer has no such assurance. >It seems to me that FLOAT()/(float) and DBLE()/(double) should mean >"convert the operand to exactly the precision that would be given >to a REAL {DOUBLE PRECISION} variable in storage." >Note that only converting an expression which already appears to >have the desired type is changed, and then only if intermediate >results are calculated in a higher precision. Most programs should >notice no change. > >Without this definition in the standard, you have to explicitly >introduce temporary variables, .... This seems to be a reasonable requirement, especially in view of the fact that without it you have absolutely no guarantee, EVEN IF YOU INTRODUCE TEMPORARY VARIABLES, that the result will be as it would be stored in memory. This is because an optimizing compiler is perfectly free to assign your temporary variables to registers rather than to memory. -- Dave Seaman ags@j.cc.purdue.edu