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

pmk@hall.cray.com (Peter Klausler) (03/02/88)

In article <832@unmvax.unm.edu>, mike@turing.UNM.EDU (Michael I. Bushnell) writes:
> In article <296@draken.nada.kth.se> d86-ctg@nada.kth.se (Claes Thornberg) writes:
> >... And thinking of simple things
> >like not using real variables in for-loops isn't that hard...
> 
> I see NOTHING which precludes:
> 	float x;
> 	for (x = 0; x < 20; x += .5)
> 		printf("%f ", x);
> The output would, of course, be
> 	0.0 0.5 1.0 ... 19.0 19.5
> But maybe the condescending poster didn't mean loop control variables.

No, I'm sure that Claes did, for the use of floating point in loop control
is a nasty habit that produces very unportable code. Michael's example
might be all right, as the loop increment is a power of two, but, in general,
the construct

	float start, end, incr, var;
	for (var=start; var<end; var+=incr)
		/* body */

can't be trusted to execute exactly (int)((end-start)/incr) times.


Kernighan and Plauger's "Elements of Programming Style" states
		10.0 times 0.1 is hardly ever 1.0.

I submit that
	0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1 is hardly ever
	10.0 times 0.1, and neither need be 1.0.


Let me back this up with a real example. Try this program on your
favorite C implementation:

main (argc, argv)
int argc;
char **argv;
{
        float s, e, i, v;
        int c = 0;
        sscanf (*++argv, "%f", &s);
        sscanf (*++argv, "%f", &e);
        sscanf (*++argv, "%f", &i);
        printf ("start=%f, end=%f, incr=%f\n", s, e, i);
        for (v = s; v < e; v += i)
		c++;
        printf ("%d iterations executed; expected %d; final value=%f",
                c, (int)((e-s)/i), v);
        }


I ran this on three of our development systems, and saw:

	cray-2% a.out 0 1 0.1
	start=0.000000, end=1.000000, incr=0.100000
	10 iterations executed; expected 10; final value=1.000000

	cray-ymp% a.out 0 1 0.1
	start=0.000000, end=1.000000, incr=0.100000
	iteration 11, v=1.000000
	11 iterations executed; expected 10; final value=1.100000

	sun-3% a.out 0 1 0.1
	start=0.000000, end=1.000000, incr=0.100000
	10 iterations executed; expected 9; final value=1.000000

For "a.out 0 1 0.01", the Cray-2 executed 100 iterations, while the
Y-MP and Sun-3 did 101.


This is not a bug of any sort, although compiler writers are
used to getting problem reports on it. This is a consequence of the
inexactitude of finite-precision floating-point arithmetic,
and its different hardware/software implementations.

If this scares you, good! Floating-point can be a nasty nemesis of
the numerically naive. Claes' advice is sound:

		USE INTEGERS FOR COUNTING.

Peter Klausler @ Compiler Development, Cray Research, Inc.
	...!{ihnp4!cray,sun!tundra}!hall!pmk

rob@kaa.eng.ohio-state.edu (Rob Carriere) (06/24/88)

In article <16276@brl-adm.ARPA> jrv%sdimax2@mitre-bedford.arpa writes:
> [...]
>I believe that floating point arithmetic is exact as long as all the values
>are integers with not too many bits - and they typically allow more bits than
>a long would.  If there are exceptions to this, I'd like to hear about them.

Ahem, I hope you are using tests on the loop that do not depend on
exactness.

If not, two things can go wrong:
1)  somewhere in your computation, a roundoff error creeps in, and
  your results are no longer exact (example: raise a to the power b: for
  int's this is typically done by multiplication (exact), for floats by
  logarithms or similar (not exact, even if both a and b have integer
  values)).
2)  you port to a machine with a slightly different floating point
  representation, and what used to be exact is no longer.

In fact, now that I start thinking about it, some languages don't even
guarantee that (representable number)*(representable number)=(1 of
{representable number, overflow, underflow}).  In the same vein,
nobody I know *guarantees* that integers are representable (i.e. the
closest approximation to 2 might be 1.999999)

>
>                                - Jim Van Zandt
Rob Carriere
"It's a poor computer that can think of only one way to confuse the
issue"

jrv%sdimax2@mitre-bedford.arpa (06/24/88)

> If this scares you, good! Floating-point can be a nasty nemesis of
> the numerically naive. Claes' advice is sound:
>        USE INTEGERS FOR COUNTING.
> 
> Peter Klausler @ Compiler Development, Cray Research, Inc.
>              ...!{ihnp4!cray,sun!tundra}!hall!pmk

I suggest the revision

         Use integer values for counting.
                     ^^^^^^

I believe that floating point arithmetic is exact as long as all the values
are integers with not too many bits - and they typically allow more bits than
a long would.  If there are exceptions to this, I'd like to hear about them.

                                - Jim Van Zandt

daniels@teklds.TEK.COM (Scott Daniels) (06/25/88)

In article <329@accelerator.eng.ohio-state.edu> 
	rob@kaa.eng.ohio-state.edu (Rob Carriere) writes:
>... nobody I know *guarantees* that integers are representable (i.e. 
>the closest approximation to 2 might be 1.999999)
    In fact, I have read interesting arguments that propose that many
integers be purposely mis-represented.  The argument goes as follows:

Assume (for concreteness) a binary floating point representation with 
an 4-bit mantissa ranging [0..1),and a 3-bit signed exponent field, and a 
1-bit sign.  We will not use IEEE hidden most-significant bit.  Examples:
 Notation: Exponent/ binary fraction  (only scribble positive numbers here)
	-2/.1000 = 0.125   -1/.1000 = 0.25    0/.1000 = 0.5  
	 0/.1100 = 0.75     1/.1000 = 1.0     1/.1001 = 1.125
	 1/.1010 = 1.25     1/.1100 = 1.5     2/.1000 = 2.0  
	...

Now notice that in this representation, a number E/F actually represents
the range:   E / (F - 1/32)  through E / (F + 1/32).  The size of the interval
depends only on the exponent.
	1/.1110 = 1.750    represents  1.6875 through 1.8125
	1/.1111 = 1.875    represents  1.8125 through 1.9375
	2/.1000 = 2.000    represents  1.8750 through 2.1250 
	2/.1001 = 2.250    represents  2.1250 through 2.3750 

Ok, here we go:
  To represent 2.0 exactly, we could use 2/.1000, but that represents the
interval 1.8750:2.1250.  Now, there is a tighter specification which is 
entirely within that interval: 1/.1111 (which represents 1.8125:1.9375), so
we should use that tighter interval since no poin inside it is any further
from the desired value 2.0 than the range that 2/.1000 gives.  Hence the
besty representation (the tightest) for 2.0 is an interval which does not
even include 2.0!

-Scott Daniels daniels@teklds.TEK.COM (or daniels@teklds.UUCP)

bgibbons@Apple.COM (Bill Gibbons) (06/25/88)

>>I believe that floating point arithmetic is exact as long as all the values
>>are integers with not too many bits ....
>>...

>Ahem, I hope you are using tests on the loop that do not depend on exactness.
>If not, two things can go wrong:
>1)  somewhere in your computation, a roundoff error creeps in, and
>  your results are no longer exact (example: raise a to the power b: for
>  int's this is typically done by multiplication (exact), for floats by
>  logarithms or similar (not exact, even if both a and b have integer
>  values)).
>2)  you port to a machine with a slightly different floating point
>  representation, and what used to be exact is no longer.
>  your results are no longer exact (example: raise a to the power b: for
>  int's this is typically done by multiplication (exact), for floats by
>  logarithms or similar (not exact, even if both a and b have integer
>  values)).
>...nobody I know *guarantees* that integers are representable (i.e. the
>closest approximation to 2 might be 1.999999)

If you press them for it, virtually everyone *guarantees* that integers are
representable, as long as they fit within the MANTISSA portion of a floating-
point number.  Almost everyone guarantees that add, subtract and multiply 
(of integer-valued floating-point numbers) produce exact results, as long as
the result is representable (i.e. the exact integer value still fits).  Divide
is less reliable, but everyone with hardware floating-point will give you the
right result for "trunc(x/y)".

Since add and subtract are safe, so are comparisons.

I agree that exponentiation is not safe, but raising a floating-point value to
an integer (typed as integer) power is always done with repeated multiplies
(in about log2(exponent) time), and since multiplies are accurate, so is this
type of exponentiation.

The log function is NOT accurate enough, but the usual reason for taking log
of an integer is to get the number of bits or decimal digits.  The bit length
can be computed from the floating-point exponent (using the IEEE logb()
function).  The number of decimal digits can be computed as
     (int) (log10(x)*(1.0+epsilon)) + 1
where epsilon is about 2**(-n+3), and floating-point numbers have "n" bits
in the mantissa.  This is pretty safe as long as the integer is less than (n-5)
bits long; the log() function on any given machine is almost always accurate 
enough to get this right.

As for converting between different machines, this is safe as long as both
machines have enough bits in the mantissa and a binary conversion is used.
Going through ASCII is usually safe, as long as you print enough digits and you
round the result after reading it back.

Bill Gibbons
(consulting to Apple)

cik@l.cc.purdue.edu (Herman Rubin) (06/25/88)

In article <12784@apple.Apple.COM>, bgibbons@Apple.COM (Bill Gibbons) writes:
> >>I believe that floating point arithmetic is exact as long as all the values
> >>are integers with not too many bits ....
> >>...
  
		--------------------------------------------

> I agree that exponentiation is not safe, but raising a floating-point value to
> an integer (typed as integer) power is always done with repeated multiplies
> (in about log2(exponent) time), and since multiplies are accurate, so is this
> type of exponentiation.

		---------------------------------------------

Apparently Bill missed the discussion a few months ago about exponentiation
in C.  There was considerable disagreement about what should be done about
the problem, and a rather large number of correspondents seemed unable to
consider computing powers except by the pow(x,y) function, which does not
use the reasonably quick and exact method.  

There seems to be a considerable lack of understanding of the need for good
integer arithmetic for mathematical computations.  This shows up in both
hardware and software.  If one looks at Brent's multiple precision package,
one sees that it is necessary to go to short words to handle the problem.
It is easy to simulate floating point with integer arithmetic, but clumsy
and slow the other way around.
-- 
Herman Rubin, Dept. of Statistics, Purdue Univ., West Lafayette IN47907
Phone: (317)494-6054
hrubin@l.cc.purdue.edu (ARPA or UUCP) or hrubin@purccvm.bitnet

daniels@teklds.TEK.COM (Scott Daniels) (06/28/88)

Oh do I have egg all over my face.  Ah well, I will go home and calculate 
fractions for a week.

In article <3637@teklds.TEK.COM> daniels@teklds.UUCP (Scott Daniels) writes:
--That's me--
>	... [ notation: E/F is F (a binary fraction) * pow(2.0,E) ] ...
>Now notice that in this representation, a number E/F actually represents
>the range:   E / (F - 1/32)  through E / (F + 1/32).  The size of the interval
>depends only on the exponent.
>	1/.1110 = 1.750    represents  1.6875 through 1.8125
>	1/.1111 = 1.875    represents  1.8125 through 1.9375
>	2/.1000 = 2.000    represents  1.8750 through 2.1250 
>	2/.1001 = 2.250    represents  2.1250 through 2.3750 

Assume we treat the last bit as possibly in error.  The numbers are in the 
range: E / (F - 1/16)  through E / (F + 1/16).  Again, the size of the 
interval depends only on the exponent.
	1/.1110 = 1.750    represents  1.625  through 1.875
	1/.1111 = 1.875    represents  1.750  through 2.000
	2/.1000 = 2.000    represents  1.750  through 2.250 
	2/.1001 = 2.250    represents  2.000  through 2.500

>Ok, here we go:
>  To represent 2.0 exactly, we could use 2/.1000, but that represents the
>interval 1.8750:2.1250.  Now, there is a tighter specification which is 
>entirely within that interval: 1/.1111 (which represents 1.8125:1.9375), so
>we should use that tighter interval since no poin inside it is any further
>from the desired value 2.0 than the range that 2/.1000 gives.  Hence the
>besty representation (the tightest) for 2.0 is an interval which does not
>even include 2.0!
>
	1/.1110 = 1.750    represents  1.625  through  1.875
	1/.1111 = 1.875    represents  1.750  through  2.000
	2/.1000 = 2.000    represents  1.750  through  2.250 
	2/.1001 = 2.250    represents  2.000  through  2.500
  To represent 2.0 exactly, we could use 2/.1000, but that represents the
interval 1.750:2.250.  Now, there is a tighter specification which is 
entirely within that interval: 1/.1111 (which represents 1.750:2.000), so
we should use that tighter interval since no poin inside it is any further
from the desired value 2.0 than the range that 2/.1000 gives.  Hence the
best representation (the tightest) for 2.0 is an interval which does not
even include 2.0!

>-Scott Daniels daniels@teklds.TEK.COM (or daniels@teklds.UUCP)

Basically, I remebered the argument as open intervals +/- lsbit, but I
reproduced the argument using +/- 1/2 lsbit (seemed clearer and more 
appropriate), but didn't really check that that also exhibitted the behavior.
I will now sneak off into a corner and cry myself to death.

still:
  -Scott Daniels daniels@teklds.TEK.COM (or daniels@teklds.UUCP)

michael@stb.UUCP (Michael) (07/15/88)

In article <3637@teklds.TEK.COM> daniels@teklds.UUCP (Scott Daniels) writes:
>In article <329@accelerator.eng.ohio-state.edu> 
>	rob@kaa.eng.ohio-state.edu (Rob Carriere) writes:
>>... nobody I know *guarantees* that integers are representable (i.e. 
>>the closest approximation to 2 might be 1.999999)
>    In fact, I have read interesting arguments that propose that many
>integers be purposely mis-represented.  The argument goes as follows:
>
>Assume (for concreteness) a binary floating point representation with 
>an 4-bit mantissa ranging [0..1),and a 3-bit signed exponent field, and a 

Hang on, you give a nice example for why FLOATS may not be accurate. If you
actually did that for integers, I wouldn't touch you with a ten foot pole.

Besides, in your example you give a tighter specification that does not
include 2. However, I think you'll find that a specification that is
symetrical around two (i.e., has 2 in the middle of the what it represents
range) will work better.

				Michael
: --- 
: Michael Gersten			 uunet.uu.net!denwa!stb!michael
:				sdcsvax!crash!gryphon!denwa!stb!michael
: What would have happened if we had lost World War 2. Well, the west coast
: would be owned by Japan, we would all be driving foreign cars, hmm...

michael@stb.UUCP (Michael) (07/15/88)

In article <12784@apple.Apple.COM> bgibbons@apple.apple.com.UUCP (Bill Gibbons) writes:
[ discussion of whats representable or not]
>but everyone with hardware floating-point will give you the
>right result for "trunc(x/y)".

ONLY if X and Y are the same sign.

>I agree that exponentiation is not safe, but raising a floating-point value to
>an integer (typed as integer) power is always done with repeated multiplies
>(in about log2(exponent) time), and since multiplies are accurate, so is this
>type of exponentiation.

Hmm, when I've written exponentiation routines I do it with logs and exp().
Are you sure all vendors do it as you say?
			Michael
: --- 
: Michael Gersten			 uunet.uu.net!denwa!stb!michael
:				sdcsvax!crash!gryphon!denwa!stb!michael
: What would have happened if we had lost World War 2. Well, the west coast
: would be owned by Japan, we would all be driving foreign cars, hmm...

ok@quintus.uucp (Richard A. O'Keefe) (07/19/88)

In article <10477@stb.UUCP> michael@stb.UUCP (Michael Gersten) writes:
>Hmm, when I've written exponentiation routines I do it with logs and exp().

That's the obvious way to do it, but according to Cody & Waite in their
classic "Software Manual for the Elementary Functions" it is _not_ a good
idea.  There are enough gotchas in the elementary functions that you really
need to be able to rely on them as part of the language, and even then if
you care about the results it is worth testing them.  [Which is not to say
that the elementary functions need special syntactic forms!]