[comp.lang.c] How not to write a loop

msb@sq.uucp (Mark Brader) (02/12/88)

> Do all bugs bite?  This short program ran under unix on a 3b1 (and Ultrix).
>	 float x;
>	 for (x=-3*PI; x < 3*PI; x += PI/8.0)  /* Change to x = -3*PI */
>		 printf ("%f %f\n", x, sin(x)/x);

The term "ran" is evidently being used rather loosely here.  Since
PI/8 is not exactly representable in any finite number of bits, the
behavior of the loop after the intended 48 iterations is in fact
undetermined.  Or as Kernighan and Plauger put it in "The Elements
of Programming Style":

	Floating point numbers should never be used for counting.
				...
		10.0 times 0.1 is hardly ever 1.0.

On the two machines where I tried this, it did in fact loop 48 times
and not 49 times, but it did the same if I changed the "<" to "<=".

The fact that it *did* loop 48 times is in fact a perfect illustration
of the problem.  If the behavior of "x" had been exactly what the author
intended, the *next* bug would have been triggered: a division by 0
on the 25th iteration.  This was avoided on the machines I tried
because "x" did not in fact ever become exactly 0.

Kernighan & Plauger again:

			Don't stop at one bug.

If you are past the introductory phase of learning to program, and
have not done so already, acquire a copy of "The Elements of Programming
Style", and *memorize* it.

One safe way to write it is:

	int i;
	for (i = -24;  i < 24;  ++i)
		if (i != 0) {
			float x = i * (PI/8);
			printf ("%f %f\n", x, sin(x)/x);
		}

Mark Brader		"TeX has found at least one bug in every Pascal
utzoo!sq!msb		 compiler it's been run on, I think, and at least
msb@sq.com		 two in every C compiler."	-- Knuth

dhesi@bsu-cs.UUCP (Rahul Dhesi) (02/15/88)

In article <1988Feb11.200149.25172@sq.uucp> msb@sq.UUCP (Mark Brader) writes:
>	Floating point numbers should never be used for counting.

Never say never.  (Well, hardly ever say never.)  At least one
widely-used implementation of Snobol4 successfully uses a floating
point register to count the number of statements executed, relying on
floating point overflow being detected by hardware so that the program
may be aborted when the user-specified statement count limit is
reached.
-- 
Rahul Dhesi         UUCP:  <backbones>!{iuvax,pur-ee,uunet}!bsu-cs!dhesi

gwyn@brl-smoke.ARPA (Doug Gwyn ) (02/17/88)

In article <2115@bsu-cs.UUCP> dhesi@bsu-cs.UUCP (Rahul Dhesi) writes:
-In article <1988Feb11.200149.25172@sq.uucp> msb@sq.UUCP (Mark Brader) writes:
->	Floating point numbers should never be used for counting.
-Never say never.  (Well, hardly ever say never.)  At least one
-widely-used implementation of Snobol4 successfully uses a floating
-point register to count the number of statements executed, relying on
-floating point overflow being detected by hardware so that the program
-may be aborted when the user-specified statement count limit is
-reached.

Mark was talking about good style, not what some crufty program
happens to do.  Note that that SNOBOL4 implementation will malfunction
on machines like the Gould PowerNodes where the environment requires
floating-point overflow trapping to be disabled (because, in the case
of the Goulds, the compiler generates code that assumes integer overflow
is benign, and the hardware design stupidly couples the two types of
overflow to the same enable bit).

d86-ctg@nada.kth.se (Claes Thornberg) (02/27/88)

Keywords:

A simple recommendation for all hackers(?) writing programs. Even
though your program is 'running' on the machine you are using now,
please think a little about compatibility. When giving away programs,
which I encourage you all to do, there may be a lot of machine
dependent bugs not doing this, and it is frustrating trying to locate
these in a program you haven't written yourself, especially when you
know the program worked when you got it. And thinking of simple things
like not using real variables in for-loops isn't that hard, I think I
was told not to do this in my very first programming class. So all you
hackers using obscure features that works on only a few machines,
please spare us the time it takes to locate bugs that shouldn't exist.

Yours in hacking,

	Claes Thornberg
	The Royal Institute of Technology, Stockholm
	Sweden

mike@turing.UNM.EDU (Michael I. Bushnell) (02/29/88)

In article <296@draken.nada.kth.se> d86-ctg@nada.kth.se (Claes Thornberg) writes:
>A simple recommendation for all hackers(?) writing programs. Even
>though your program is 'running' on the machine you are using now,
>please think a little about compatibility. When giving away programs,
>which I encourage you all to do, there may be a lot of machine
>dependent bugs not doing this, and it is frustrating trying to locate
>these in a program you haven't written yourself, especially when you
>know the program worked when you got it. And thinking of simple things
>like not using real variables in for-loops isn't that hard, I think I
>was told not to do this in my very first programming class. So all you
>hackers using obscure features that works on only a few machines,
>please spare us the time it takes to locate bugs that shouldn't exist.


From K&R I find:

	The for statement

		for (expr1 ; expr2 ; expr3)
			statement

	is equivalent to

		expr1;
		while (expr2) {
			statement
			expr3;
		}

A difference in the semantics of continue is mentioned.

I see NOTHING which precludes:

	float x;
	for (x = 0; x < 20; x += .5)
		printf("%f\n", x);

The output would, of course, be

	0.0
	0.5
	1.0
	.
	.
	.
	19.0
	19.5

In fact, the "loop" variable can be anything whatsoever, not
even a variable at all!

Try
	for (open_socket(); connection_waiting(); accept())
		process_data();

as a main loop for a server which want to die when no one wants it 
anymore.

C has no notion of "for variables" at all, really.  But maybe the
condescending poster didn't mean loop control variables.  Maybe he
meant just what he says: "using real variables in for loops".  I 
guess he just doesn't want to see:

	int x;
	float f;
	for (x = 0; x < 10; ++x) {
		scanf("%f\n", f);
		printf("%f\n", f);	
	}
as a weird way of echoing floats.

Sigh.



				Michael I. Bushnell
				Internet: mike@turing.unm.edu
				UUCP: mike@turing.unm.edu
				Bitnet: mike@turing.unm.edu
				CSnet: mike@turing.unm.edu
				YourFavoriteNet: mike@turing.unm.edu
			Golly, don't domains make everything simpler?
For peoply who run UUCP but haven't switched over to smail *yet*, you
can try {ucbvax,gatech}!unmvax!turing!mike.  

Or write:
  {Box 295, Coronado Hall} or {Computer Science, Farris Engineering Center}
  University of New Mexico
  Albuquerque, NM 87131
Or call:
  (505)277- [2992=dorm][6116=work]

I work for the CS department.  But don't blame them.

gwyn@brl-smoke.ARPA (Doug Gwyn ) (03/01/88)

In article <832@unmvax.unm.edu> mike@turing.UNM.EDU.UUCP (Michael I. Bushnell) writes:
>	float x;
>	for (x = 0; x < 20; x += .5)
>		printf("%f\n", x);
>The output would, of course, be
>	0.0
	...
>	19.5

and, perhaps,
	20.0
(ignoring the nit that this isn't quite the format specified).

That is why the original poster didn't want to see such loop control
in putative portable programs.

ok@quintus.UUCP (Richard A. O'Keefe) (03/01/88)

In article <832@unmvax.unm.edu>, mike@turing.UNM.EDU (Michael I. Bushnell) writes:
> I see NOTHING which precludes:
> 	float x;
> 	for (x = 0; x < 20; x += .5) printf("%f\n", x);
> The output would, of course, be
> 	0.0
> 	0.5
> 	...
> 	19.5

Quite right, there is nothing in K&R (or dpANS) to prohibit this.
But you have provided a good illustration of why people disparage
the use of floating-point variables in loop control.
THERE IS NO "OF COURSE" ABOUT THAT OUTPUT!

You should not be surprised to see as the last value
	19.99999
I just tried the very similar loop
	for (x = 0; x < 21; x += .3)	/* .3 divides 21 exactly */
By analogy with your "of course", the last output should obviously
be 20.7.  In fact, when I tried it just now, the last output was
	20.999994

chas@gtss.UUCP (Charles Cleveland) (03/01/88)

In article <7384@brl-smoke.ARPA> gwyn@brl.arpa (Doug Gwyn (VLD/VMB) <gwyn>) writes:
>In article <832@unmvax.unm.edu> mike@turing.UNM.EDU.UUCP (Michael I. Bushnell) writes:
>>	float x;
>>	for (x = 0; x < 20; x += .5)
>>		printf("%f\n", x);
>>The output would, of course, be
>>	0.0
>	...
>>	19.5
>
>and, perhaps,
>	20.0
>(ignoring the nit that this isn't quite the format specified).
>
>That is why the original poster didn't want to see such loop control
>in putative portable programs.

Which is why the original loop, given that floats were going to be used,
should have been written in the (portable) way that follows:
	float x;
	for (x = 0; x < 20-.25; x+= .5)
		printf("%f\n",x);
This is just a matter of subtracting half the increment from the (to be
excluded) upper limit.  Not that this is a general solution.  Perhaps
we should write something like
	#define min(x,y) = (x < y ? x : y)
	float x, start=0, end=20, inc=.5, shift;
	shift=.5*min(end-start,inc);
	for (x = start; x < end-shift; x += inc)
		printf("%f\n",x);

What a disgusting mess.  Why not use ints, and avoid all this convolution?
Because the convolution may be unavoidable.

This doesn't really have to do with loops and portability, but with the fact
that round-off errors make tests for equality between floating point
numbers unreliable.  We don't need two machines to get problems --
one machine and two sets of numbers to compare will suffice.  If the
condition you need to impose is intrinsically a floating point condition
then you may have to think hard about how to code it so that round-off can't
cause problems.  If you do this right, the code will likely be portable
across machines as well.

But if your job is to write a program to print all the integers and
half-integers between 0 and 19.5 inclusively, use ints. ;-)
-- 
-Life would be so much easier if we could just look at the source code.-

Charles Cleveland    Georgia Tech School of Physics    Atlanta, GA 30332
UUCP: ...!gatech!gtss!chas         INTERNET:  chas@ss.physics.gatech.edu

cuddy@convex.UUCP (03/02/88)

ok@quintus.Sun.COM (our feed trashes addreesses, that may not be right!)
> In article <832@unmvax.unm.edu>, mike@turing.UNM.EDU (Michael I. Bushnell) writes:
> > I see NOTHING which precludes:
> > 	float x;
> > 	for (x = 0; x < 20; x += .5) printf("%f\n", x);
> > The output would, of course, be
> > 	0.0
> > 	0.5
> > 	...
> > 	19.5
> 
> Quite right, there is nothing in K&R (or dpANS) to prohibit this.
> But you have provided a good illustration of why people disparage
> the use of floating-point variables in loop control.
> THERE IS NO "OF COURSE" ABOUT THAT OUTPUT!
> 
> You should not be surprised to see as the last value
> 	19.99999

Only if you declare things as floats!

!! (from K&R:  2.4.4, pg. 181
!!   A floating constant consists of an integer part . . .
!!   . . .  may be missing.  Every floating constant is taken to be a double.
!!			     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

every "C" compiler (at least ones that conform to K&R)  will cast floats 
to doubles and sometimes they have troubles.  this has to do with your 
compiler and machine.  the above program segment ( for(x=0; x< 20; x+=5)... )
should probably be written:
	double x;
	^^^^^^  this will help you not get round off errors
	for (x = 0.0; x < 20.0; x += 0.5) printf("%f\n",x);
    		 ^^^ 	  ^^^^       ^^^
		 these don't have to be this way, but i've allways felt it
		 good programming practice.

> THERE IS NO "OF COURSE" ABOUT THAT OUTPUT!
There's lots of "of course"s  about "C" and that output if it's done right!

> I just tried the very similar loop
> 	for (x = 0; x < 21; x += .3)	/* .3 divides 21 exactly */
> By analogy with your "of course", the last output should obviously
> be 20.7.  In fact, when I tried it just now, the last output was
> 	20.999994

more float-double round off error.  This kind of thing bites you when you
use the -lm libraries without #including <math.h>. 

------------------------------------------------------------------------------
The opinions stated above are barely my own, how can they be of my employers?
Mike Cuddy, 
CONVEX Computer Corporation
701 N. Plano Rd.
Richardson, TX  75081
uucp:{uiucdcs,ihnp4,sun,allegra}!convex!cuddy
------------------------------------------------------------------------------
  

gwyn@brl-smoke.ARPA (Doug Gwyn ) (03/04/88)

In article <64400004@convex> cuddy@convex.UUCP writes:
>Only if you declare things as floats!

This is of course wrong.  Floating-point imprecision is not due
simply to float<->double conversion, but is inherent in the
whole approach.

ok@quintus.UUCP (Richard A. O'Keefe) (03/04/88)

In article <64400004@convex>, cuddy@convex.UUCP writes:
c ok@quintus.Sun.COM (our feed trashes addreesses, that may not be right!)
c o In article <832@unmvax.unm.edu>, mike@turing.UNM.EDU (Michael I. Bushnell) writes:
c o m I see NOTHING which precludes:
c o m 	float x; for (x = 0; x < 20; x += .5) printf("%f ", x);
c o m The output would, of course, be
c o m 	0.0 0.5 ... 19.5
c o Quite right, there is nothing in K&R (or dpANS) to prohibit this.
c o But you have provided a good illustration of why people disparage
c o the use of floating-point variables in loop control.
c o THERE IS NO "OF COURSE" ABOUT THAT OUTPUT!
c o You should not be surprised to see as the last value
c o 	19.99999
c 
c Only if you declare things as floats!

If I have understood his posting, Mike Cuddy claims that the problem goes
away if you use 'double' rather than 'float'.  If you believe this,
I have a second-hand bridge you might be interested in buying...

The problem has nothing to do with casting floats to doubles.
The problem is that floating-point arithmetic is APPROXIMATE.
If you go to double precision, you get a better approximation,
perhaps even a very good one, but it is still an approximation.
If the least significant bit out of 53 is wrong, that can still
be enough to make a loop execute one time too many or one time
too few.

Note that there are perfectly sensible loops where the control variable
is float (or double), e.g.

	while (fabs(y*y - x) > eps) { compute next y }

The point at issue is that if you use floats (or doubles) for
COUNTING, you are being a fool to yourself and a burden to others,
because your counts will be approximate.

{Actually, this isn't quite true either.  Using a double-precision
variable to hold multiples of 1 is usually safe, and on far too many
machines will give you a larger range of integers than any integral
type.  64-bit integers, anyone?}

gak@mhuxm.UUCP (Vincent Hatem) (03/05/88)

In article <64400004@convex>, authorplaceholder@convex.UUCP writes:
>> In article <832@unmvax.unm.edu>, mike@turing.UNM.EDU (Michael I. Bushnell) writes:
>  [lots of junk I've deleted...]
> 
>> I just tried the very similar loop
>> 	for (x = 0; x < 21; x += .3)	/* .3 divides 21 exactly */
>> By analogy with your "of course", the last output should obviously
>> be 20.7.  In fact, when I tried it just now, the last output was
>> 	20.999994
>
> more float-double round off error.  This kind of thing bites you when you
> use the -lm libraries without #including <math.h>. 
> 
> The opinions stated above are barely my own, how can they be of my employers?
> Mike Cuddy, 
> CONVEX Computer Corporation

The errors are caused by the inherent nature of ANY IEEE floating point
number. Try it with 0.1 as an increment. It's even worse. It's all because
you simply *can't* represent the number 0.1 with IEEE floating point. You 
can come very, *very*, close - but you'll still have something like
0.099999999999999 instead of 0.10000000000... 

This occurs in *both* 32 and 64 bit floating point numbers, and is called
"rounding errors." The 64 bit numbers are more accurate, but they are still
off. When you printf() a floating point number, the function *rounds* the
float, then prints it.

Take a numerical methods course.


-- 
Vincent Hatem
AT&T International, International Systems Operations, UNIX Technical Support
                        Telefon: (201) 953-8030
		Please send all e-mail to: ihnp4!atti01!vch

850347s@aucs.UUCP (Hume Smith) (03/09/88)

ok@quintus.UUCP (Richard A. O'Keefe) wrote in <712@cresswell.quintus.UUCP>:
>In article <832@unmvax.unm.edu>, mike@turing.UNM.EDU (Michael I. Bushnell) writes:
>> I see NOTHING which precludes:
>> 	float x;
>> 	for (x = 0; x < 20; x += .5) printf("%f\n", x);
>> The output would, of course, be
>> 	0.0
>> 	0.5
>> 	...
>> 	19.5
>
>You should not be surprised to see as the last value
>	19.99999
>I just tried the very similar loop
>	for (x = 0; x < 21; x += .3)	/* .3 divides 21 exactly */
>By analogy with your "of course", the last output should obviously
>be 20.7.  In fact, when I tried it just now, the last output was
>	20.999994

I would be quite surprised to see 19.99999 come out in the former
loop.  .5 is exactly representable in floating binary, and should cause
no problem as an increment.  .3 is not representable in a finite binary
fraction (.01001 something, it repeats eventually), and cannot be summed
up to get an integral value.

If the machine did its arithmetic in decimal (Texas Intstruments?  Sharp?)
then the .3 increment would again not be a problem, since the representation
is exact again.

This, however, is the point.  The behaviour of a particular increment
is unreliable.

(Not to be insulting at all, but this may be why our Computer Science
department makes its majors take our Numerical Methods Course.  The
very first thing we do is Machine Epsilon and roundoff error.  This
comment is just to make enough text to keep inews happy.)
-- 
  Hume Smith                        UUCP:
Math Department     {uunet|watmath|utai|garfield|mnetor}!
   Acadia U                   dalcs!aucs!850347s

gae@osupyr.UUCP (Gerald Edgar) (03/10/88)

Lengthy discussion about loops such as:

>>> 	float x;
>>> 	for (x = 0; x < 20; x += .5) printf("%f\n", x);

How about this:

 	float x;
 	for (x = 0; x < 19.75; x += .5) printf("%f\n", x);

-- 
  Gerald A. Edgar                               TS1871@OHSTVMA.bitnet
  Department of Mathematics                     gae@osupyr.UUCP
  The Ohio State University  ...{akgua,gatech,ihnp4,ulysses}!cbosgd!osupyr!gae
  Columbus, OH 43210                            70715,1324  CompuServe

jep@oink.UUCP (James E. Prior) (03/11/88)

In article <402@osupyr.UUCP> gae@osupyr.mast.ohio-state.edu.UUCP (Gerald Edgar) writes:
>Lengthy discussion about loops such as:
>
>>>> 	float x;
>>>> 	for (x = 0; x < 20; x += .5) printf("%f\n", x);
                          ^
                          A decimal point would have been nice!
>
>How about this:
>
> 	float x;
> 	for (x = 0; x < 19.75; x += .5) printf("%f\n", x);

Yes, this will work reliably across many machines.  However, your
suggested solution causes another problem, obfuscation of the intended
end value of the loop.

I realize that 19.75 is half the increment less that twenty, but
someone else might not have that much insight.  You should right
your code so that it is obvious how the magic number 19.75 is
derived.

Next step in obviousness:

	for (x=0; x<20.-.5/2.; x+=.5) ...

but even now, it is not clear that the .5 in 20.-.5/2. is the same
.5 as in x+=.5.  The next step in obviousness is:

#define STEP (.5)
	for (x=0; x<20.-STEP/2.; x+=STEP) ...

This is quite handy for someone who never saw the <20 to begin with.
With the way you wrote your solution, someone needing to change the
step to .25 could inadvertently re-introduce reliability problems as
follows:

	for (x=0; x<19.75; x+=.25)

Another step to de-mystify 20. would be nice.

I stronly recommend that you read Elements of Programming Style by
P. J. Plaugher.  You have reinforced one of the reasons I quit OSU.

>  Gerald A. Edgar                               TS1871@OHSTVMA.bitnet
>  Department of Mathematics                     gae@osupyr.UUCP
>  The Ohio State University  ...{akgua,gatech,ihnp4,ulysses}!cbosgd!osupyr!gae
>  Columbus, OH 43210                            70715,1324  CompuServe
-- 
Jim Prior    {ihnp4|osu-cis}!n8emr!oink!jep    jep@oink.UUCP

ray@micomvax.UUCP (Ray Dunn) (03/19/88)

In article <7409@brl-smoke.ARPA> gwyn@brl.arpa (Doug Gwyn) writes:
>
>This is of course wrong.  Floating-point imprecision is not due
>simply to float<->double conversion, but is inherent in the
>whole approach.

Yup.  Interestingly, this is one of the very first things I can remember
being taught in "computer science" - NEVER test "reals" against absolute
values - choose a precision and test against values +/- that precision.

Mind you, that was in "autocoder", on a Ferranti Sirius, and the year was
1964.

Ah well,  yawn....

Ray Dunn.  ..{philabs, mnetor, musocs}!micomvax!ray

edw@IUS1.CS.CMU.EDU (Eddie Wyatt) (03/22/88)

> >
> >This is of course wrong.  Floating-point imprecision is not due
> >simply to float<->double conversion, but is inherent in the
> >whole approach.
> 
> Yup.  Interestingly, this is one of the very first things I can remember
> being taught in "computer science" - NEVER test "reals" against absolute
> values - choose a precision and test against values +/- that precision.

One of the "tricks" I do when optimizing code is handling typical cases
different from the general case.  In particular, I do alot with
3-d geomety and 3-d rotations.  The typical transform case has zero roll and
pitch angular rotations (x and y axis rotations).  Its to my benefit
to test for this situation and handle it differently.  The test
I use is a direct equality test to see if the roll and pitch angles
are zero and if so only perform a rotation about the z-axis.
I have no problems with this on my Sun.  Profiling has shown that this
test does decrease the time spent in the routine that performs
the rotation.  I'm hesitant to use an approximately equal test
because it is more expensive to execute. Besides if the test fails
because of numerical imprecision when it actually should have succeeded
I still get the same result, it just takes a few more instructions to
calculate it.  So the point being made, well, NEVER SAY NEVER. :-)

BTW I don't use homogeneous transforms because the transform parameters
have covariances associated with them.
-- 

Eddie Wyatt 				e-mail: edw@ius1.cs.cmu.edu