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!]