Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!utgpu!water!watmath!clyde!rutgers!husc6!linus!philabs!pwa-b!mmintl!franka
From: franka@mmintl.UUCP
Newsgroups: comp.lang.c
Subject: Re: C and Floating Point
Message-ID: <2107@mmintl.UUCP>
Date: Tue, 14-Apr-87 20:51:05 EST
Article-I.D.: mmintl.2107
Posted: Tue Apr 14 20:51:05 1987
Date-Received: Sat, 18-Apr-87 03:30:09 EST
References: <15958@sun.uucp> <5716@brl-smoke.ARPA> <14681@cca.CCA.COM>
Reply-To: franka@mmintl.UUCP (Frank Adams)
Organization: Multimate International, E. Hartford, CT
Lines: 59
Keywords: C Fortran Floating Point

In article <14681@cca.CCA.COM> g-rh@CCA.CCA.COM.UUCP (Richard Harter) writes:
>We have a iterative process which generates a floating point number,
>x, with the following iteration equation
>
>	x <- x + delta
>
>where delta is a function of x and where it is expected that the
>iteration converges to some limit point x0.
>
>The iteration will have converged if the iteration equation does
>not change the value of x; i.e. we will have found the limit point
>of the sequence, given that we are using floating point arithmetic
>of the machine in question.  The question then is, how do we test
>for this condition?  Is the following pseudo code acceptable?
>
>	if ((x+delta)==x) then terminate_iteration
>	else continue_iteration
>
>If not, why not?  Alternatively should one use
>
>	if (delta>0.) then
>		if ((x+delta)<=x) then terminate_iteration
>		else continue_iteration
>	else if (delta<0.) then
>		if ((x+delta)>=x) then terminate_iteration
>		else continue_iteration
>	else terminate_iteration
>
>or is this also unacceptable?

As far as I can tell, based on comments made by various people, there are
machines and C compilers in existence where neither of these will work.  In
either case, x+delta may be in a register with more precision than x.

I assert that the first test should work.  When the test is for equality,
the compiler should do the extra work to ensure that there are no extra bits
floating around.  In particular, if x and y are expressions of the same
type, then x == y should be true if assignment of x and y to any variable of
that type would give the same result. There may still be cases where x<=y is
false, while x==y is true, which is ugly and potentially dangerous; but one
can't have everything.

(Actually, based on the above, the first example is still not guaranteed.
What should work is:

	y = x + delta;
	if (x == y) terminate the iteration;
	x = y;

If we compare x directly to x+delta, we will have to write the expression
(x+delta) again to update x; and the compiler may redo the computation and
get a different result.)

This sort of thing is one of the reasons for, as somebody suggested, writing
the floating point libraries in assembler.

Frank Adams                           ihnp4!philabs!pwa-b!mmintl!franka
Ashton-Tate          52 Oakland Ave North         E. Hartford, CT 06108