Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!cis.ohio-state.edu!tut.cis.ohio-state.edu!udecc.engr.udayton.edu!blackbird.afit.af.mil!dlindsle
From: dlindsle@afit.af.mil (David T. Lindsley)
Newsgroups: comp.lang.fortran
Subject: Re: What is zero, REALly??
Keywords: real numbers, old Fortran, zero comparison
Message-ID: <1991Jun12.141651.5699@afit.af.mil>
Date: 12 Jun 91 14:16:51 GMT
References: <1991Jun11.143915.11722@milton.u.washington.edu>
Distribution: usa
Organization: Air Force Institute of Technology
Lines: 59


You didn't say what kind of machine you were trying to get this to
run on, so...

Arithmetic IFs can be tricky on one's-complement machines.  (Zero can
be represented two ways, signed negative or positive.)  I was porting
a fairly large (50+K LOC) package of Fortran IV to a muliprocessing
computer.  I won't name names; but FYI it had its own proprietary
operating system, on top of which we were running Un*x.  The hardware
supposedly had zero-checking built in (so it would treat +0 and -0
identically), but we tested and verified that in certain cases, a
"negative" zero would take the first rather than the second branch.

It could have been a hardware glitch, or something in their operating
system, or a foulup in the interface between that O/S and Un*x, or
maybe something in the compiler -- maybe all of the above.  I'm inclined
to think the multiple layering fouled the compiler into generating
faulty code.

<pedagogy ON>
However, depending on what your numbers are, it _could_ be roundoff
or truncation error.  For an ordered pair of numbers (x,y) with the
mean z, then abs(x-z) == abs(y-z), identically.  But the computer
represents numbers as a finite sequence of bits, thus what is stored
as a binary approximation to a number (call this B(a) for any real
number a).  Since x & y are not necessarily equal to B(x) and B(y),
respectively, abs(B(x)-B(z)) is not necessarily equal to abs(B(y)-B(z)).

Now if I remember my assembler correctly (it's been a while), 0.5
(one-half) will generate in binary form an infinite repeating decimal
fraction which will, of course, be truncated.
<pedagogy OFF>

I'm sorry if I patronize by addressing such an elementary point.
But from your article, it seems as though the problem may actually
be this simple.

How about this: (assuming 64-bit floating point storage, 80x87/68881/
G_FLOATING compatible):

	common /toler/ ztol
	data ztol /1.d-300/
...
	if (abs(x).le.ztol) x=0.
	if (x) ...

An ugly kludge, admittedly, but it works.  It also has the great advantage
of requiring only minimal editing of the code.

If your floating-point is IEEE-compatible, you could also check for
underflow, etc.

Hope this helps.


-- 
Dave Lindsley	#24601#					   OPINIONS. MINE. 
dlindsle@blackbird.afit.af.mil		(The words don't come no smaller.)
			?? lamroN eb yhW ??