negative zero generated under HP-UX 11.11

[Peter Cebull]

My company maintains a large legacy code, written in FORTRAN, that is used in the nuclear power industry for licensing calculations. We recently transitioned from a B132L workstation running f77 to a C240 with f90. We have a number of standard test cases we run to verify that the code is behaving correctly when we release a new version, and I ran these on the new machine. One of the test case showed significant differences from the standard data created on the old machine.

The problem was traced to the Fortran SIGN intrinsic. Specifically, SIGN(ONE,VGAS) is supposed to return a 1 when VGAS is zero. Turns out VGAS was -0, so the SIGN intrinsic returned -1, which completely altered the results of the entire transient. A simple Fortran program confirms the problem:

vgasj = 2.
void = 0.
ajunt = 4.
spvjg = -2.
wgj = vgasj * void * ajunt / spvjg
write(*,*) wgj
stop
end

When run on our new machine w/ 11.11, this returns -0.0. Our old machine w/ f77 (HP-UX 10.20) returns 0.0, as do our SUN and IBM platforms. If spvjg = 2., then 0.0 is returned. I originally thought this was a compiler problem, but an equivalent C program will show the same behavior.

A call to HP software support pointed me to the +nosignedzero flag on the f90 compiler, which fixed my original problem in the Fortran code. However, I don't buy their explanation. They say that this is normal IEEE floating point behavior, and that I am seeing the results of roundoff (a very small negative number), not a negative zero. 1) Maybe I don't fully understand floating point behavior, but I think the above expression should be exactly represented by 0.0, and not rounded off to a "small negative number". 2) If this is normal IEEE behavior, why does this occur only on HP-UX 11.11, and not HP-UX 10.20, Solaris 7, or AIX 4.3?

This rather long post boils down to my main question: What has changed which causes me to see this behavior that hasn't been there before? I wonder if there is a bug in the 11.x math libraries or some such thing. If this is correct behavior, and something was changed to correct previous bugs, what was changed and why?

Has anyone else encountered this? Thanks in advance for your input.

[A. Clay Stephenson]

Hi Peter:

By coincidence, I also worked at a nuclear power plant a long time ago in a universe far, far away. I think I will have to side with HP on this one but I don't agree that it is a small negative number and thus a rounding issue. The IEEE standard defines both a negative zero and a positive zero. In both representations both the exponent and fraction fields are all zero bits. The only difference is bit0 - 1 for negative and 0 for positive. Both are perfectly valid values and -0.0 is NOT the same as 0.0. It actually appears that the 11.11 more closely adheres to the IEEE standard. You have found the workaround to force the desired behavior with the compiler option.

Regards, Clay

[Gregory Fruth]

I'm no expert on the subtleties of IEEE math, but here's
a C program that shows what's going on. It performs
the calculations in your example step by step and shows
the floating-point results on the left and the corresponding
in-memory bit patterns on the right. On my 10.20
machines (both PA1.1 and PA2.0) the results are:

2.000000 0x4000000000000000
0.000000 0x0000000000000000
0.000000 0x0000000000000000
-0.000000 0x8000000000000000

Everything looks OK until you do a "0/-2", which results
in -0 (0x8000000000000000). Look on docs.hp.com
for the "HP-UX Floating Point Guide":

http://www.docs.hp.com/hpux/onlinedocs/B3906-90006/B3906-90006.html

They have a table which lists the IEEE-mandated results
for operations with +0 and -0, but they cop out on the
one of interest:

In multiplication and division operations with positive and negative
zero operands, the sign is determined by the usual arithmetic rules.

Perhaps they elaborate somewhere else.

Here's the program, hope it helps:

#include

union wgj_u {
double dval;
long long lval;
};
void prt(union wgj_u);

int main() {
double vgasj, voidd, ajunt, spvjg;
union wgj_u wgj;

vgasj = 2.;
voidd = 0.;
ajunt = 4.;
spvjg = -2.;

wgj.dval = vgasj;
prt(wgj);

wgj.dval = vgasj * voidd;
prt(wgj);

wgj.dval = vgasj * voidd * ajunt;
prt(wgj);

wgj.dval = vgasj * voidd * ajunt / spvjg;
prt(wgj);

return(0);
}

void prt(union wgj_u wgj) {
printf("%16lf 0x%16.16llx\n", wgj.dval, wgj.lval);
}

[Gregory Fruth]

P.S. To compile the program use "cc -Ae" or "gcc".