more comments on base conversion
Tim Peters
uunet!ksr!tim
Sat Jun 9 01:24:39 PDT 1990
> > [tim p griping about fortran defining x.eq.y as ((x)-(y)).eq.0]
> > ...
> [david hough]
> I think that this kind of problem is avoided in Fortran implementations
> since there is a catch-all clause allowing "mathematically equivalent"
> transformations, without defining what these are.
The hitch is that the various "this-or-that-ly equivalent" clauses
don't "catch" *all* (nor should they): The compiler's license is
absolutely revoked across parentheses. So, e.g., the user who writes
X = (X+Y)+Z
is *guaranteed* that the code will not evaluate it as X+(Y+Z) --
although the compiler is quite free to evaluate the "Z" part as
"LOG(EXP(Z))", or the "X+Y" part as "(X**2-Y**2)/(X-Y)" if it knows
X-Y can't be 0. Here's a practical-- and directly relevant --example:
in lots of Fortran you'll see stuff like
IF (EPS .GT. 0. .AND. 1.+EPS .EQ. 1.) THEN
C EPS is positive & smaller than machine epsilon
and this is about as portable a cheap test as there is; an experienced
Fortran programmer would have no fears about relying on it. If the
compiler were really free to do its "mathematically equivalent" tricks
here, it could act as if the user had written
IF (.FALSE.) THEN
because "1.+EPS.EQ.1" is mathematically equivalent to EPS.EQ.0., so
(EPS.GT.0..AND.1.+EPS.EQ.1.) is mathematically equivalent to
(EPS.GT.0..AND.EPS.EQ.0.) which is in turn equivalent (via trichotomy)
to plain .FALSE..
The reason a conforming compiler *can't* do this is that 1.+EPS.EQ.1.
is *defined* to mean ((1.+EPS)-1.).EQ.0., and the parens around
(1.+EPS) force the compiler to evaluate some expression equivalent to
1.+EPS before it can combine it with the "-1.". In exactly the same
way, in ((x)-(y)).eq.0 the outermost parens force the compiler to
evaluate the "(x)-(y)" part as a unit (unless, of course, the compiler
can determine that the value of the expression isn't needed -- but
that's another twist to Fortran optimization esoterica that isn't
relevant in the example).
All Fortran has to do to legitimize using a direct comparison
instruction is to drop the outermost parens in the definition; to
define "x relop y" as meaning "(x)-(y).relop.0". Then the appropriate
"equivalence" license (there's actually more than one) will allow a
compiler to evaluate that as "x.relop.y" directly. But as it stands
today, it ain't strictly conformant. The X3J3 old-timers I've talked
to about this recognize the problem, but can't remember whether they
stuck in the outermost parens to force the order of evaluation or to
make the precedence clearer, but were quite willing to change it; I
really don't understand why those who prosecuted the case in this
incarnation of X3J3 got shot down (seems like they couldn't get enough
current members to understand there *is* a problem).
> Sun's compilers have always done comparisons with comparison
> operators, and nobody complains about that.
Violating an obscure requirement of a std doesn't necessarily lead to
customer complaints (as any number of deficient 754 vendors can
testify <grin>). You're certainly within the *spirit* of Fortran to
generate the fastest comparison-like code you can, so no reasonable
Fortran programmer is going to carp. I think you're really within the
letter of the std too, but this is a little tricky & isn't the kind of
argument you're going to feel good about:
1) An interpretation to 77 ruled that a compiler can violate
parentheses if no standard-conforming program could detect the
difference. The case in question was brought up by some nit-
picker who looked at the *generated code* for something like
LOGICAL L1, L2, L3
...
L1 = L1 .AND. (L2 .AND. L3)
and whined because L1.AND.L2 was being done first (note: Fortran
is silent about whether its logical connectives are short-
circuiting or not).
Since Sun supports denorms, then, for normal x and y a std-
conforming program can't tell the difference between (x-y).eq.0 and
x.eq.y. Although on a Cray they could (pick x and y so that x-y
underflows to 0.), so a Cray *must* subtract (big deal, they don't
have a comparison instruction anyway).
2) While *a* program can be written to tell the difference between
(x-x).eq.0 and x.eq.x when, e.g., x is +inf, the very fact that
it uses +inf (or a NaN ... a non-normal) renders it non-conforming
a priori. Hence the std has nothing to say about its behavior, and
you can do whatever you like.
So, in the end, I believe the user can't write a *conforming* program
that detects Sun's violation of the parens, hence I believe you're
technically off the hook (if someone gripes that, e.g., you didn't set
the invalid operation flag after +inf.eq.+inf (hence must not have
done the subtraction), it takes some fancier sleaze to weasel out of
it, but I think it can be done ...).
Just wait 'til you try to sell "you can't change '0+x' to 'x'" to X3J3
<grin> ...
> If most of the problems relate to ambiguities in the language definitions,
> they can be fixed by clarifying the definitions.
Unclear to me (haven't worried much about it yet); am sure it won't be
easy, though (C seems to be a much easier fit to 754, partly because
it doesn't have so strong a tradition of free-wheeling "optimization").
> > What may not work perfectly is *starting* from the machine with the
> > greater precision.
> Almost all the time, you'd get a different value in the less
> precise machine.
I take it you're thinking that, e.g., if I start with w=53 bits of
precision on a Sun and move to n=48 bits on a Cray then every number I
move over with more than n significant bits must necessarily lose
information? Sure. I don't care, though -- that I can't always do it
perfectly doesn't mean I don't want (or may not *have*!) to do it. I
want to do it as well as I can, and it seems the best I can
realistically do (i.e., short of simulating a wider precision on the
Cray) is to map the w-bit x on the Sun to the closest n-bit x' on the
Cray. The question for me is just how close to that non-insane ideal
perfect conversions of 17-digit decimal strings on both ends can
guarantee to get. If x has no more than n significant bits, it seems
that the end result of the two conversions must be x' exactly; if x
has more than n, it seems that the result of the conversions must be
within a half (n-bit; Cray) ulp of x', and will in fact equal x'
exactly except for some "nearest/even" end cases. And after that, I
can move x' back and forth any number of times without losing anything
more. That's really awfully good (if true <grin>)! So I expect the
Cray people to lobby vigorously for un-exempting them from the
perfect-rounding requirement <ahem>.
> > [tgp wondering why the all-754 shops don't just transfer numbers
> > in binary or hex form]
> You can't do this in a standard portable way, unfortunately. Maybe a
> standard is needed.
If there's no std way to ship binary data on homogeneous networks, I'd
say there's a *desperate* need for a std! Amazing (not trying to be
snide here at all -- really think it is amazing -- I've avoided
networks as much as possible so far, and am starting to think that was
a wise idea <0.9 grin>).
> > configure a compiler to run with optimization off by default (no
> > kidding!).
> Doing things this way puts everybody into the optimizing compiler
> debugging business whether they've thought about it or not.
Of course. That's the true "high end" Fortran frame of mind.
> As long as Fortran-8x is a compatible upward superset of
> Fortran-77, there will be good reasons to use C instead (the
> complexity of Fortran I/O is my particular pet peeve).
I confess: I *still* feel lost whenever I have to deal with a non-
trivial FORMAT statement ...
but-i-have-to-compile-the-language-not-use-it<grin>-ly y'rs - tim
Tim Peters Kendall Square Research Corp
timaksr.com, ksr!timaharvard.harvard.edu
More information about the Numeric-interest
mailing list