Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!utgpu!water!watnot!watmath!clyde!ima!mirror!cca!g-rh
From: g-rh@cca.UUCP
Newsgroups: comp.lang.c
Subject: Re: C and Floating Point
Message-ID: <14813@cca.CCA.COM>
Date: Tue, 7-Apr-87 21:49:14 EST
Article-I.D.: cca.14813
Posted: Tue Apr  7 21:49:14 1987
Date-Received: Sat, 11-Apr-87 03:38:39 EST
References: <15958@sun.uucp> <5716@brl-smoke.ARPA> <14681@cca.CCA.COM> <1163@drivax.UUCP>
Reply-To: g-rh@CCA.UUCP (Richard Harter)
Organization: Computer Corp. of America, Cambridge, MA
Lines: 59
Keywords: C Fortran Floating Point

Ah, the ignominy of it all, being lectured on the basics by an earnest
bright eyed young lad.  Ah, well, Richard, tis you own fault for being
too subtle.

I raised the question of iteration equations and asked about thoughts
about various pseudocode alternatives for testing convergence.  Now, as
we all know, bad things often happen in iteration equations when they
are done in floating point.  When we get near a limit point we tend to
get oscillation in the last few bits with no final convergence.  I
(incorrectly) supposed that I had made it clear that these problems had
already been addressed -- that we had arranged our code so that there
would be a limit point (i.e. convergence is forced) -- so that the only
issue was whether or not we had reached the limit point.  To this end
I included the following paragraph:

"Typically we will find the the calculation of the delta in question
becomes numerically unstable when x is sufficiently close to x0.
Let us suppose that we have addressed this problem, so that we
are simply left with the convergence issue."

With this in mind, I posted some some pseudocode for dealing
with the actual detection of convergence.  William Tyler commented
on this pseudocode by noting that it would break down if the iteration
(when implemented in floating point) does not converge and said that
was unacceptable for that reason.

Well, he's right.  If we haven't minded our p's and q's, that's exactly
what does happen.  However, if we have done our job the calculation will
have the property that there will be an interval such that, for all x
in the interval, the iteration will converge.

I am assuming that we have done this.  The issue I want to deal with
is simpler.  Given that we have arranged things so that convergence
is forced, i.e. so that are guaranteed to reach a point such that,
for some n,

	x[n+1] = x[n] + delta[n] = x[n],

what is the proper way to test for this, if the calculations are
being done in floating point?  For example, is it not conceivable
that the following two pieces of pseudocode will not yield the
same results:

(1)	if ((x+delta)==x) then done = TRUE
	else                   done = FALSE

(2)	temp = x + delta
	if (temp==x) then done = TRUE
	else              done = FALSE

In (1) 'x+delta' might be stored in registers with extension bits
whereas x, also in a register, has no extension bits.  But is (2)
safe?  If not (and why not) then what is?  Again, I am assuming that
we have arranged things so that the iteration will actually converge.
This may seem picky, but it is inattention to picky little details
that gets your ass in a sling.
-- 

Richard Harter, SMDS Inc. [Disclaimers not permitted by company policy.]