Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP Path: utzoo!utgpu!water!watnot!watmath!clyde!rutgers!mit-eddie!uw-beaver!tektronix!reed!psu-cs!omepd!jimv From: jimv@omepd.UUCP Newsgroups: comp.lang.c Subject: Re: C and Floating Point (LONG) Message-ID: <532@omepd> Date: Sat, 4-Apr-87 05:10:43 EST Article-I.D.: omepd.532 Posted: Sat Apr 4 05:10:43 1987 Date-Received: Sun, 5-Apr-87 11:37:50 EST References: <15958@sun.uucp> <5716@brl-smoke.ARPA> Reply-To: jimv@omepd.UUCP (Jim Valerio) Organization: Intel Corp., Hillsboro Lines: 263 Doug, I believe that you are mistaken in a number of points you make in your rejoinder to Dave Hough's note here. Before I dive into the details, though, let me first comment on your trailing remark: >Finally, I would urge newsgroup readers to (a) consider that others >may have something useful to say, and if you think you disagree, >try at least to understand what their reasons are for their differing >opinions, and to (b) avoid posting when you really are unlikely to >contribute significantly to the discussion. This is good advice, but it seems to me that you ignored it when you posted your article. When reading (and responding?) to this article, I hope that you keep in mind that I am posting this because I believe that I have something valuable to contribute. As you will see below, I also think that Dave's article had a lot to contribute -- in fact more than your response. Okay, on to what I want to say. In his article, Dave made a number of important points about programming in C and numerical programming. You picked up on his references to Fortran and responded rather angrily (more on that later), but didn't address or acknowledge anything about the references he suggested as important background. Those references are a veritable gold mine of important information on floating-point arithmetic. I strongly recommend that anyone and everyone who is going to do any significant amount of floating-point programming familiarize themselves with the information in these references. You say: >It is BECAUSE I care about good floating-point algorithms that I even >bothered to suggest that floating-point == was not a good idea. I'm >not going to try to get this put in the standard, but I WAS hoping that >people like Mr. Hough might learn something from the ensuing discussion. and >It is certainly true that only a few members of X3J11 feel really >comfortable concerning requirements for floating-point algorithms. >However, there ARE members who do have lots of experience in this area. I would hope that you would learn something from Dave. Dave was one of the major players in the 754 and 854 IEEE floating-point standards. He has been involved in the design and implementation of many floating-point systems over the last 10 years, including Sun's more recent offerings. I don't know exactly what Dave's PhD was in at UC Berkeley, but the gist of it was numerical programming. In short, Dave is a trained and experienced *expert* in floating-point programming. His credentials are excellent and are not to be ignored. I agree that there are many ways to misuse a floating ==, but that doesn't mean that all uses of it are misuses. I believe that there are equally as many bad uses of "goto", pointer conversions, and so on that are possible in the language. Is the possibility of misuse sufficient justification for prohibiting useful operations? I say no. Surely you're not suggesting that we need to protect the ignorant floating-point users from themselves? >I am still waiting for a really good reason for using floating-point >== in any robust algorithm. I have posted the only one that anyone >has yet come up with (== 0.0 special case switch) and have also posted >an example of the problems its general use can cause. I have many years >of experience in programming significant numerical applications in both >Fortran and C, and I am convinced that floating == is virtually NEVER >a good thing to write in one's code. It would be a Very Good Thing in >my estimation if "lint" (or your local equivalent) would warn about >any use of floating ==, so you would be likely to re-think those >sections of code. Maybe the "lint" approach is reasonable, although I am not entirely convinced. In any case, there are just too many good uses of floating == to put this warning in the compiler, or to uniformly deprecate it as an operation. Let me give you some examples of good uses of floating ==. I don't have a whole lot of floating-point code on this machine to draw on, but here are some examples taken from the code which I can easily get my hands on. (1) First, there's your example of comparing a value for == or != to 0.0. A quick textual search of Spice3 found 109 obvious instances of this usage. Usually this check is to avoid a zero divide fault, but sometimes it also is used to detect discontinuities in approximated functions (e.g. ATAN2(0,0)). Some typical examples are: if (model->MOS2substrateDoping == 0.0) xwb = 0.25e-6; if ((model->MOS2fastSurfaceStateDensity == 0.0) || (OxideCap == 0.0)) { if (dd[i] == 0.0) if (realpart(&cc[i]) == 0.0) { if (largest == 0.0) { } else if (lo == 0.0) { if (delta == 0.0) { if (mat1[i * n + i] == 0.0) else if (((num < 10.0) && (num > -10.0)) || (num == 0.0)) Sometimes the check is made for efficiency, such as in the following example: if(arg != 0) sqarg=sqrt(arg); (2) Similar to (1), there are also comparisons against other "known" values. Sometimes these values are extrema, sometimes these are initial values, and sometimes these are just critical or undefined domains in the computed function. Examples from Spice3 include: if(model->MOS1bulkJctBotGradingCoeff == .5) { } else if ((ox != HUGE) && (oy != HUGE)) { if (stack[st].e_dlist[1] != HUGE) for (d = dv; *d != HUGE; d++) Admittedly there are less of these in Spice than comparisons against zero, but they do occur. Another place I looked in the 4.3 math library. There I found more examples, of which the following are typical: if( x != -1.0) if (x == 1.0) { if ( x == zero ) if ( x == negone ) if((t=copysign(x,one))==one) return(zero/zero); else if(y==two) return(x*x); else if ( (t=drem(y,two)) == zero) return( pow_p(-x,y) ); Yet another example of this is seen in some of the following code which tests a math library for correct handling of special case operands: if(x!=pzero||s!=1.0) printf("Failed; "); else printf("O.K. ; "); if(x!=PI) printf("Failed; "); else printf("O.K. ; "); if(x!= -PI) printf("Failed; "); else printf("O.K. ; "); if(x!= PIh) {printf("Failed; "); else if(y!= PIh) {printf("Failed; "); if(x==pinf) printf("O.K. ; "); else printf("Failed; "); if(x==none) printf("O.K. ; "); else printf("Failed; "); if(x==nzero&&s== -1.0) printf("O.K. ; "); else printf("Failed; "); if(x==unft) printf("O.K. ; "); else printf("Failed; "); if(x==gunt) printf("O.K. ; "); else printf("Failed; "); if(x==pone) printf("O.K. ; "); else printf("Failed; "); if(x==pinf&&y==pinf) {printf("O.K. ; "); if(x== 5.0) printf("O.K. ; "); else printf("Failed; "); if(x== 5.0*unft) printf("O.K. ; "); else printf("Failed; "); if(x==ninf) printf("O.K. ; "); else printf("Failed; "); if(x!= 5929741.) printf("Failed; "); else printf("O.K. ; "); (3) Sometimes it is quite natural to compare two floating-point numbers for equality. Often this is done as an early guard in an expression where the difference of the two numbers must be non-zero. (This guard doesn't always guarantee that if the difference were computed that the result would be nonzero for all rounding modes, but that may be irrelevant because the function is actually computed using an algebraically equivalent formula.) Other times one is looking to see if a floating-point number has an integral value. Some examples of this from Spice3 and 4.3bsd libm are: rcheck((dd1[i] >= 0) || (floor(dd2[i]) == ceil(dd2[i])), "power"); d[i] = dd1[i] == dd2[i]; d[i] = ((realpart(&c1) == realpart(&c2)) && (imagpart(&c1) == imagpart(&c2))); if (arg == t) { if((double)k==y) { /* if y is an integer */ In yet other circumstances, a function may be defined in such a way that it's value is computed differently or in a special way when the arguments are equal. The common example of this is ATAN2(x,x). (4) In IEEE code, the construct (x != x) is often used to quickly and easily detect NaN operands. The following example from 4.3bsd libm says it all: if(x!=x) return(x); /* x is NaN */ One piece of software which I have not quoted examples from, but does make heavy use of floating ==, is Kahan's Paranoia. I argue that even though it isn't a typical floating-point program, it is performing a function that should be expressible in a language which purports to support floating-point arithmetic. On to the next point you make: >It is totally bogus to argue that mathematicians' idea of the meaning >of parentheses in expressions has anything to do with time order of >evaluation. That is a Fortran notion. C actually conforms more >closely to what a mathematician means by parentheses, which is to >group subexpressions in order to override default rules of operator >precedence. I don't believe that Dave explicitly said that the reason for a compiler to honor parentheses is a mathematical one. I think that he said that the mathematicians who did the early programming in floating-point influenced Fortran to guarantee the order of operations by obeying parentheses. All of us who have programmed floating-point know that the order that operations are done in is sometimes critical to the correct evaluation of expressions. Granted: ANSI has specified the (albeit clumsy) "+" operator to guarantee order of evaluation. Dave's point was that (a) pre-ANSI C did not have this concept, and (b) the people doing early floating-point programming (in Fortran) recognized the value of enforced order of evaluation. Sure it's possible in pre-ANSI C to guarantee order of evaluation by breaking the expressions up into separate statements. The cost there is the (un)clear code of what is being computed. Requiring order of evaluation to be enforced by separate statements is almost as clumsy requiring each integer operation to be a separate statement, disallowing side-effects, or not supporting more than 1 level of dereferencing pointers in a single expression. Sure, you can program with restrictions like these, but it sure isn't going to be your language of choice. Here are some examples of parenthesized expressions, pulled from the Elefunt tests, which want the order of evaluation guaranteed by the parentheses: x = ((one + x * a) - one) * 16.0; y = x / ((half + x * half) *((half - x) + half)); x = (x + eight) - eight; while (((b+one)-b)-one == zero); if ((a+betam1)-a != zero) if ((one - a) - one != zero) y2 = (y/two + y) - y; w = zz * xsq/(den*(den-one)); w = (x - zz) + w; w = ((half-zz) + half) * ((half + zz) + half); (I suppose some of these look pretty obscure even coded like this, eh? :-)) >Because C also has a heavily-used macro expansion >facility, we have to take that into account when deciding what the >"right" thing to do with parentheses is. I assure you that X3J11 has >heard the arguments that the Fortran-style proponents give, and we >decided differently for SOUND REASONS. We ALSO found that their >legitimate concern could be met by use of the unary + operator, >which was originally introduced for unrelated reasons. Therefore >there is no reason remaining other than "to act just like Fortran" >that can be given for what Mr. Hough apparently would prefer. I understand that there are historical reasons and compatibility issues for the ANSI spec to read the way it does about parentheses. I understand that C macro facility tends to require excess parentheses. And, I accept this decision by ANSI. But, that doesn't mean that Dave's comments are invalid. It is the "easy" coding using parentheses (as seen above) that Dave was arguing is a *good thing* for numerical programmers. What is being traded off here is clean expressibility of floating-point expressions. It's okay in some sense that this has happened, but you can't deny that it has happened, and you can't claim that somehow the resulting language is better for floating-point computations. Let me try to wrap this up with a few more miscellaneous thoughts. First, ANSI C is much better for floating-point programming than the existing C compilers. The biggest win that comes to mind is prototyped functions, which gives much-needed control over data types of passed arguments. Personally, I wish that there was a function storage class called "inline" which would allow me to get rid of a number of hairy macros. The improvements in the ANSI specification of the C math library are good if not perfect, as is the language support for extended (long double) variables. I'm even reasonably happy with the parenthesis situation, since I just overspec my compiler and require it to preserve parentheses for floating-point operations. (Sure, that's not encouraging portable code, but it's an easy dependency to document and consequently doesn't bother me much.) And who am I? Just a floating-point newcomer: only 5 years of continuous floating-point experience, and only working on my 4th floating-point processor. -- Jim Valerio {verdix,intelca!mipos3}!omepd!jimv, jimv@omepd.intel.com