What Every Computer Scientist Floating-Point Arithmetic Should Know About DAVID GOLDBERG Xerox Palo Alto Research Center, 3333 Coyote Hill Road, Palo Alto, CalLfornLa 94304 Floating-point arithmetic is considered an esotoric subject by many people. This is rather surprising, because floating-point is ubiquitous in computer systems: Almost every language has a floating-point datatype; computers from PCs to supercomputers have floating-point accelerators; most compilers will be called upon to compile floating-point algorithms from time to time; and virtually every operating system must respond to floating-point exceptions such as overflow This paper presents a tutorial on the aspects of floating-point that have a direct impact on designers of computer systems. It begins with background on floating-point representation and rounding error, continues with a discussion of the IEEE floating-point standard, and concludes with examples of how computer system builders can better support floating point, Categories and Subject Descriptors: (Primary) C.0 [Computer Systems Organization]: Languages]: General– instruction set design; D.3.4 [Programming G. 1.0 [Numerical Analysis]: General—computer Processors —compders, optirruzatzon; arithmetic, error analysis, numerzcal algorithms (Secondary) D. 2.1 [Software languages; D, 3.1 [Programming Engineering]: Requirements/Specifications– Languages]: Formal Definitions and Theory —semantZcs D ,4.1 [Operating Systems]: Process Management—synchronization General Terms: Algorithms, Design, Languages Additional Key Words and Phrases: denormalized number, exception, floating-point, floating-point standard, gradual underflow, guard digit, NaN, overflow, relative error, rounding error, rounding mode, ulp, underflow INTRODUCTION Builders of computer systems often need information about floating-point arithmetic. There are however, remarkably few sources of detailed information about it. One of the few books on the subject, Floating-Point Computation by Pat Sterbenz, is long out of print. This paper is a tutorial on those aspects of floating-point hereafter) that arithmetic ( floating-point have a direct connection to systems building. It consists of three loosely connected parts. The first (Section 1) discusses the implications of using different rounding strategies for the basic opera- tions of addition, subtraction, multiplication, and division. It also contains background information on the two methods of measuring rounding error, ulps and relative error. The second part discusses the IEEE floating-point standard, which is becoming rapidly accepted by commercial hardware manufacturers. Included in the IEEE standard is the rounding method for basic operations; therefore, the discussion of the standard draws on the material in Section 1. The third part discusses the connections between floating point and the design of various aspects of computer systems. Topics include instruction set design, Permission to copy without fee all or part of this material is granted provided that the copies are not made or distributed for direct commercial advantage, the ACM copyright notice and the title of the publication and its data appear, and notice is given that copying is by permission of the Association for Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specific permission. @ 1991 ACM 0360-0300/91/0300-0005 $01.50 ACM Computing Surveys, Vol 23, No 1, March 1991 6 1 . David Goldberg CONTENTS INTRODUCTION 1 ROUNDING ERROR 1 1 Floating-Point Formats 12 Relatlve Error and Ulps 1 3 Guard Dlglts 14 Cancellation 1 5 Exactly Rounded Operations 2 IEEE STANDARD 2 1 Formats and Operations 22 S~eclal Quantltles 23 Exceptions, Flags, and Trap Handlers SYSTEMS ASPECTS 3 1 Instruction Sets 32 Languages and Compders 33 Exception Handling 4 DETAILS 4 1 Rounding Error 3 42 Bmary-to-Decimal Conversion 4 3 Errors in Summatmn 5 SUMMARY APPENDIX ACKNOWLEDGMENTS REFERENCES optimizing compilers, and exception handling. All the statements made about floating-point are provided with justifications, but those explanations not central to the main argument are in a section called The Details and can be skipped if desired. In particular, the proofs of many of the theorems appear in this section. The end of each m-oof is marked with the H symbol; whe~ a proof is not included, the ❑ appears immediately following the statement of the theorem. 1. ROUNDING ERROR Squeezing infinitely many real numbers into a finite number of bits requires an approximate representation. Although there are infinitely many integers, in most programs the result of integer computations can be stored in 32 bits. In contrast, given any fixed number of bits, most calculations with real numbers will produce quantities that cannot be exactly represented using that many bits. Therefore, the result of a floating-point calculation must often be rounded in order to ACM Computing Surveys, Vol 23, No 1, March 1991 fit back into its finite representation. The resulting rounding error is the characteristic feature of floating-point computation. Section 1.2 describes how it is measured. Since most floating-point calculations have rounding error anyway, does it matter if the basic arithmetic operations introduce a bit more rounding error than necessary? That question is a main theme throughout Section 1. Section 1.3 discusses guard digits, a means of reducing the error when subtracting two nearby numbers. Guard digits were considered sufficiently important by IBM that in 1968 it added a guard digit to the double precision format in the System/360 architecture (single precision already had a guard digit) and retrofitted all existing machines in the field. Two examples are given to illustrate the utility of guard digits. The IEEE standard goes further than just requiring the use of a guard digit. It gives an algorithm for addition, subtraction, multiplication, division, and square root and requires that implementations produce the same result as that algorithm. Thus, when a program is moved from one machine to another, the results of the basic operations will be the same in every bit if both machines support the IEEE standard. This greatly simplifies the porting of programs. Other uses of this precise specification are given in Section 1.5. 2.1 Floating-Point Formats Several different representations of real numbers have been proposed, but by far the most widely used is the floating-point representation.’ Floating-point representations have a base O (which is always assumed to be even) and a precision p. If 6 = 10 and p = 3, the number 0.1 is represented as 1.00 x 10-1. If P = 2 and P = 24, the decimal number 0.1 cannot lExamples of other representations slas;, aud szgned logan th m [Matula 1985; Swartzlander and Alexopoulos are floatzng and Kornerup 1975] Floating-Point 100X22 [,!,1 I o Figure 1. r \ c I r 2 1 Normalized I , , i 3 4 5 6 7 numbers + ““. <~. (1) number will The term floating-point be used to mean a real number that can be exactly represented in the format under discussion. Two other parameters associated with floating-point representations are the largest and smallest allowable exponents, e~~X and e~,~. Since there are (3P possible significands and exponents, a emax — e~i. + 1 possible floating-point number can be encoded in L(1°g2 ‘ma. – ‘m,. + 1)] + [log2((3J’)] + 1 its, where the final + 1 is for the sign bit. The precise encoding is not important for now. There are two reasons why a real number might not be exactly representable as a floating-point number. The most common situation is illustrated by the decimal number 0.1. Although it has a finite decimal representation, in binary it has an infinite repeating representation. Thus, when D = 2, the number 0.1 lies strictly between two floating-point numbers and is exactly representable by neither of them. A less common situation is that a real number is out of range; that is, its absolute value is larger than f? x 2This term was introduced [196’71and has generally mantissa. by Forsythe 7 when (3 = 2, p = 3, em,n = – 1, emax = 2. +dP_l&(P-l))&, o<(il ● 11 O X221.11X22 I be represented exactly but is approximately 1.10011001100110011001101 x 2-4. In general, a floating-point number will be represented as ~ d. dd “ . . d x /3’, where d. dd . . . d is called the significand2 and has p digits. More prekdO. dld2 “.” dp_l x b’ reprecisely, sents the number + ( do + dl~-l 101X22 Arithmetic o‘m= or smaller than 1.0 x ~em~. Most of this paper discusses issues due to the first reason. Numbers that are out of range will, however, be discussed in Sections 2.2.2 and 2.2.4. Floating-point representations are not necessarily unique. For example, both represent 0.01 x 101 and 1.00 x 10-1 0.1. If the leading digit is nonzero [ do # O in eq. (1)], the representation is said to The floating-point numbe normalized. ber 1.00 x 10-1 is normalized, whereas 0.01 x 101 is not. When ~ = 2, p = 3, e~i~ = – 1, and e~~X = 2, there are 16 normalized floating-point numbers, as shown in Figure 1. The bold hash marks correspond to numbers whose significant is 1.00. Requiring that a floating-point representation be normalized makes the representation unique. Unfortunately, this restriction makes it impossible to represent zero! A natural way to repre sent O is with 1.0 x ~em~- 1, since this preserves the fact that the numerical ordering of nonnegative real numbers corresponds to the lexicographical ordering of their floating-point representations. 3 When the exponent is stored in a k bit field, that means that only 2 k – 1 values are available for use as exponents, since one must be reserved to represent O. Note that the x in a floating-point number is part of the notation and different from a floating-point multiply operation. The meaning of the x symbol should be clear from the context. For example, the expression (2.5 x 10-3, x (4.0 X 102) involves only a single floating-point multiplication. and Moler replaced the older term 3This assumes the usual arrangement where exponent is stored to the left of the significant ACM Computing the Surveys, Vol. 23, No 1, March 1991 8* 1.2 David Relative Error Goldberg and Since rounding error is inherent in floating-point computation, it is important to have a way to measure this error. Consider the floating-point format with ~ = will be used 10 and p = 3, which throughout this section. If the result of a floating-point computation is 3.12 x 10’2 and the answer when computed to infinite precision is .0314, it is clear that this is in error by 2 units in the last place. Similarly, if the real number .0314159 is represented as 3.14 x 10-2, then it is in error by .159 units in the last place. In general, if the floating-point number d. d . . . d x fle is used to represent z, it is in error by Id. d . . . d– ( z//3’) I flp - 1 units in the last place.4 The term ulps will be used as shorthand for “units in the last place. ” If the result of a calculation is the floating-point num ber nearest to the correct result, it still might be in error by as much as 1/2 ulp. Another way to measure the difference between a floating-point number and the real number it is approximating is relawhich is the difference betive error, tween the two numbers divided by the real number. For example, the relative error committed when approximating 3.14159 by 3.14 x 10° is .00159 /3.14159 = .0005. To compute the relative error that corresponds to 1/2 ulp, observe that when a real number is approximated by the closest possible floating-point number P d dd ~. dd X ~e, as large the digit /3’ Since the absolute u error can be x /3’ where & is as ‘(Y x ~/2. This error is ((~/2)&P) numb... of the form d. dd --- dd x /3e all have this same absolute error but have values that range between ~’ and O x fle, the relative error ranges bex /3’//3’ and ((~/2)&J’) tween ((&/2 )~-’) 4Un1ess the number z is larger than ~em=+ 1 or smaller than (lem~. Numbers that are out of range in this fashion will not be considered until further notice. ACM Computmg That is, x /3’/~’+1. Ulps Surveys, Vol. 23, No 1, March 1991 :(Y’ 2 s ;Ulp s ;6-’. (2) ~n particular, the relative error corre spending to 1/2 ulp can vary by a factor of O. This factor is called the wobble. Setting E = (~ /2)~-P to the largest of the bounds in (2), we can say that when a real number floating-point is rounded to the closest number, the relative error is always bounded by c, which to as machine epsilon is referred In the example above, the relative er~or was .oo159i3, ~4159 = 0005. To avoid such small numbers, the relative error is normally written as a factor times 6, which in this case is c = (~/2)P-P = error 5(10) -3 = .005. Thus, the relative would be expressed as ((.00159/ 3.14159) /.oo5)e = O.l E. To illustrate the difference between ulps and relative error, consider the real number x = 12.35. It is approximated by Z = 1.24 x 101. The error is 0.5 ulps; the relative error is 0.8 e. Next consider the computation 8x. The exact value is 8 x = 98.8, whereas, the computed value is 81 = 9.92 x 101. The error is now 4.0 ulps, but the relative error is still 0.8 e. The error measured in ulps is eight times larger, even though the relative error is the same. In general, when the base is (3, a fixed relative error expressed in ulps can wobble by a factor of up to (3. Conversely, as eq. (2) shows, a fixed error of 1/2 ulps results in a relative error that can wobble by (3. The most natural way to measure rounding error is in ulps. For example, rounding to the neared flo~ting.point number corresponds to 1/2 ulp. When analyzing the rounding error caused by various formulas, however, relative error is a better measure. A good illustration of this is the analysis immediately following the proof of Theorem 10. Since ~ can overestimate the effect of rounding to the nearest floating-point number b; the wobble factor of (3, error estimates of formulas will be tighter on machines with a small p. Floating-Point When only the order of magnitude of rounding error is of interest, ulps and e may be used interchangeably since they differ by at most a factor of ~. For example, when a floating-point number is in error by n ulps, that means the number of contaminated digits is logD n. If the relative error in a computation is ne, then contaminated 1.3 Guard digits = log,n. (3) Digits One method of computing the difference between two floating-point numbers is to compute the difference exactly, then round it to the nearest floating-point number. This is very expensive if the operands differ greatly in size. Assuming P = 3, 2,15 X 1012 – 1.25 X 10-5 would be calculated as x = 2.15 X 1012 .0000000000000000125 X 1012 X – y = 2.1499999999999999875 X 1012, y = which rounds to 2.15 x 1012. Rather than using all these digits, floating-point hardware normally operates on a fixed number of digits. Suppose the number of digits kept is p and that when the smaller operand is shifted right, digits are simply discarded (as opposed to rounding). Then, 2.15 x 1012 – 1.25 x 10-5 becomes Arithmetic wrong in every error be? digit! g 9 How bad can the Theorem 1 Using a floating-point format with parameters /3 and p and computing differences using p digits, the relative error of the result can be as large as b – 1. A relative error of 13– 1 in Proofi the expression x – y occurs when x = 1.00””” Oandy=. pp. ””p, wherep= @– 1. Here y has p digits (all equal to Q). The exact difference is x – y = P ‘p. When computing the answer using only p digits, however, the rightmost digit of y gets shifted off, so the computed difference is P–p+l. Thus, the error is p-p – @-P+l = ~-P(~ – 1), and the ror is $-P((3 – l)/O-p relative = 6 – 1. er- H When f? = 2, the absolute error can be as large as the result, and when 13= 10, it can be nine times larger. To put it another way, when (3 = 2, (3) shows that the number of contaminated digits is log2(l/~) = logJ2 J’) = p. That is, all of the p digits in the result are wrong! Suppose one extra digit is added to guard against this situation (a guard That is, the smaller number is digit). truncated to p + 1 digits, then the result of the subtraction is rounded to p digits. With a guard digit, the previous example becomes x = 1.010 x 101 x = 2.15 X 1012 ‘y = 0.00 x 1012 x–y y = 0.993 x 101 x–y= .017 x 101, =2.15x1012. The answer is exactly the same as if the difference had been computed exactly then rounded. Take another example: 10.1 – 9.93. This becomes and the answer is exact. With a single guard digit, the relative error of the re suit may be greater than ~, as in 110 – 8.59: x= x= 1.01 x 101 ‘y = 0.99 x 101 X–yz .02 x 101. The correct answer is .17, so the computed difference is off by 30 ulps and is y = z–y= 1.1OX 102 .085 X 102 1.015 x 102 This rounds to 102, compared with the correct answer of 101.41, for a relative error of .006, which is greater than ACM Computing Surveys, Vol 23, No. 1, March 1991 10 David “ Goldberg e = .005. In general, the relative error of the result can be only slightly larger than c. More precisely, we have Theorem 2. Theorem 2 If x and y are floating-point numbers in a format with 13 and p and if subtraction is done with p + 1 digits (i. e., one guard digit), then the relative rounding error in the result is less than 2 ~. This theorem will be proven in Section 4.1. Addition is included in the above theorem since x and y can be positive or negative. 1.4 –b+ Cancellation Section 1.3 can be summarized by saying that without a guard digit, the relative error committed when subtracting two nearby quantities can be very large. In other words, the evaluation of any expression containing a subtraction (or an addition of quantities with opposite signs) could result in a relative error so large (Thethat all the digits are meaningless orem 1). When subtracting nearby quantities, the most significant digits in the operands match and cancel each other. There are two kinds of cancellation: catastrophic and benign. Catastrophic cancellation occurs when the operands are subject to rounding errors. For example, in the quadratic forbz – 4 ac occurs. mula, the expression The quantities 62 and 4 ac are subject to rounding errors since they are the results of floating-point multiplications. Suppose they are rounded to the nearest floating-point number and so are accurate to within 1/2 ulp. When they are subtracted, cancellation can cause many of the accurate digits to disappear, leaving behind mainly digits contaminated by rounding error. Hence the difference might have an error of many ulps. For b = 3.34, a = 1.22, example, consider and c = 2.28. The exact value of b2 -4 ac is .0292. But b2 rounds to 11.2 and 4 ac rounds to 11.1, hence the final answer is .1, which is an error by 70 ulps even though 11.2 – 11.1 is exactly equal ACM .1. The subtraction did not introduce any error but rather exposed the error introduced in the earlier multiplications. Benign cancellation occurs when subtracting exactly known quantities. If x and y have no rounding error, then by Theorem 2 if the subtraction is done with a guard digit, the difference x – y has a very small relative error (less than 2 e). A formula that exhibits catastrophic cancellation can sometimes be rearranged to eliminate the problem. Again consider the quadratic formula to Computmg Surveys, Vol 23, No 1, March 1991 r2 = ~b2–4ac –b–~ (4) 2a When b2 P ac, then “ b2 – 4 ac does not involve a cancellation and ~ = \ b 1. But the other addition (subtraction) in one of the formulas will have a catastrophic cancellation. To avoid this, multiply the numerator and denominator of r-l by – b – ~ for r2 ) to obtain rl (and similarly 2C = –b–~’ 2C rz = (5) –b+~” If b2 % ac and b >0, then computing rl using formula (4) will involve a cancellarl tion. Therefore, use (5) for computing and (4) for rz. On the other hand, if b <0, use (4) for computing rl and (5) for r2. The expression X2 – y2 is another formula that exhibits catastrophic cancellation. It is more accurate to evaluate it as ( x – y)( x + y). 5 Unlike the quadratic 5Although the expression ( x – .Y)(x + y) does not cause a catastrophic cancellation, it IS shghtly less accurate than X2 – y2 If x > y or x < y In this case, ( x – -Y)(x + y) has three rounding errors, but X2 – y2 has only two since the rounding error committed when computing the smaller of x 2 and y 2 does not affect the final subtraction. Floating-Point formula, this improved form still has a subtraction, but it is a benign cancellation of quantities without rounding error, not a catastrophic one. By Theorem 2, the relative error in x – y is at most 2 e. The same is true of x + y. Multiplying two quantities with a small relative error results in a product with a small relative error (see Section 4.1). To avoid confusion between exact and computed values, the following notation is used. Whereas x – y denotes the exact difference of x and y, x @y denotes the computed difference (i. e., with rounding error). Similarly @, @, and @ denote computed addition, multiplication, and division, respectively. All caps indicate the computed value of a function, as in LN( x) or SQRT( x). Lowercase functions and traditional mathematical notation denote their exact values as in ln( x) and &. Although (x @y) @ (x @ y) is an excellent approximation of x 2 – y2, the floating-point numbers x and y might themselves be approximations to some true quantities 2 and j. For example, 2 and j might be exactly known decimal numbers that cannot be expressed exactly in binary. In this case, even though x ~ y is a good approximation to x – y, it can have a huge relative error compared to the true expression 2 – $, and so the advantage of ( x + y)( x – y) over X2 – y2 is not as dramatic. Since comput ing ( x + y)( x – y) is about the same amount of work as computing X2 – y2, it is clearly the preferred form in this case. In general, however, replacing a catastrophic cancellation by a benign one is not worthwhile if the expense is large because the input is often (but not always) an approximation. But eliminat ing a cancellation entirely (as in the quadratic formula) is worthwhile even if the data are not exact. Throughout this paper, it will be assumed that the floating-point inputs to an algorithm are ex are computed as .aGt and Qxat the results accurately as possible. The expression X2 – y2 is more accurate when rewritten as (x – y)( x + y) because a catastrophic cancellation is Arithmetic 8 11 replaced with a benign one. We next present more interesting examples of formulas exhibiting catastrophic cancellation that can be rewritten to exhibit only benign cancellation. The area of a triangle can be expressed directly in terms of the lengths of its sides a, b, and c as A = ~s(s - a)(s - b)(s - c) , a+b+c where s = . 2 (6) Suppose the triangle is very flat; that is, s = a, and the term a = b + c. Then (s – a) in eq. (6) subtracts two nearby numbers, one of which may have rounding error. For example, if a = 9.0, b = c = 4.53, then the correct value of s is 9.03 and A is 2.34. Even though the computed value of s (9.05) is in error by only 2 ulps, the computed value of A is 3.04, an error of 60 ulps. There is a way to rewrite formula (6) so that it will return accurate results even for flat triangles [Kahan 1986]. It is A= [(la+ X(C+ (b+c))(c (a– - b))(a+ (a-b)) (b– a? c))] ’/’/4, b?c. (7) If a, b, and c do not satisfy a > b > c, simply rename them before applying (7). It is straightforward to check that the right-hand sides of (6) and (7) are algebraically identical. Using the values of a, b, and c above gives a computed area of 2.35, which is 1 ulp in error and much more accurate than the first formula. Although formula (7) is much more accurate than (6) for this example, it would be nice to know how well (7) performs in general. Theorem 3 The rounding error incurred when using the area of a t.icqqle ie at (T) #o compuie most 11 e, provided subtraction is performed with a guard digit, e <.005, and square roots are computed to within 1/2 Ulp. ACM Computing Surveys, Vol. 23, No. 1, March 1991 12 “ David Goldberg The condition that c s .005 is met in virtually every actual floating-point system. For example, when 13= 2, p >8 ensures that e < .005, and when 6 = 10, p z 3 is enough. In statements like Theorem 3 that discuss the relative error of an expression, it is understood that the expression is computed using floating-point arithmetic. In particular, the relative error is actually of the expression (sQRT(a @(b @c))@ F3(c @(a @b))@ cents! (C @(a @b)) (a @(b @c))) (8) @4. Because of the cumbersome nature of (8), in the statement of theorems we will value of E usually say the computed rather than writing out E with circle notation. Error bounds are usually too pessimistic. In the numerical example given above, the computed value of (7) is 2.35, compared with a true value of 2.34216 for a relative error of O.7c, which is much less than 11 e. The main reason for computing error bounds is not to get precise bounds but rather to verify that the formula does not contain numerical problems. A final example of an expression that can be rewritten to use benign cancellation is (1 + x)’, where x < 1. This expression arises in financial calculations. Consider depositing $100 every day into a bank account that earns an annual interest rate of 6~o, compounded daily. If n = 365 and i = ,06, the amount of money accumulated at the end of one year is 100[(1 + i/n)” – 11/(i/n) dollars. If this is computed using ~ = 2 and P = 24, the result is $37615.45 compared to the exact answer of $37614.05, a discrepancy of $1.40. The reason for the problem is easy to see. The expresinvolves adding 1 to sion 1 + i/n .0001643836, so the low order bits of i/n are lost. This rounding error is amplified when 1 + i / n is raised to the nth power. ACM Computmg Surveys, Vol 23, No 1, March The troublesome expression (1 + i/n)’ can be rewritten as exp[ n ln(l + i / n)], where now the problem is to compute In(l + x) for small x. One approach is to use the approximation ln(l + x) = x, in which case the payment becomes $37617.26, which is off by $3.21 and even less accurate than the obvious formula. But there is a way to compute ln(l + x) accurately, as Theorem 4 shows [Hewlett-Packard 1982], This formula yields $37614.07, accurate to within 2 1991 Theorem 4 assumes that LN( x) approximate ln( x) to within 1/2 ulp. The problem it solves is that when x is small, LN(l @ x) is not close to ln(l + x) because 1 @ x has lost the information in the low order bits of x. That is, the computed value of ln(l + x) is not close to its actual value when x < 1. Theorem If ln(l mula ln(l 4 – x) is computed using the for- + x) — — Ix 1 xln(l forl~x=l + x) (1 +X)-1 forl G3x#l the relative error is at most 5 c when x < 3/4, provided subtraction is formed with a guard digit, e <0.1, in is computed to within 1/2 ulp. O< perand This formula will work for any value of x but is only interesting for x + 1, which is where catastrophic cancellation occurs in the naive formula ln(l + x) Although the formula may seem mysterious, there is a simple explanation for why it works. Write ln(l + x) as x[ln(l + x)/xl = XV(x). The left-hand factor can be computed exactly, but the right-hand factor P(x) = ln(l + x)/x will suffer a large rounding error when adding 1 to x. However, v is almost constant, since ln(l + x) = x. So changing x slightly will not introduce much error. In other words, if z= x, computing XK( 2) will be a good Floating-Point approximation to xp( x) = ln(l + x). Is there a value for 5 for which 2 and 5 + 1 can be computed accurately? There is; namely, 2 = (1 @ x) e 1, because then 1 + 2 is exactly equal to 1 @ x. The results of this section can be summarized by saying that a guard digit guarantees accuracy when nearby precisely known quantities are subtracted (benign cancellation). Sometimes a formula that gives inaccurate results can be rewritten to have much higher numeri cal accuracy by using benign cancellation; however, the procedure only works if subtraction is performed using a guard digit. The price of a guard digit is not high because is merely requires making the adder 1 bit wider. For a 54 bit double precision adder, the additional cost is less than 2%. For this price, you gain the ability to run many algorithms such as formula (6) for computing the area of a triangle and the expression in Theorem 4 for computing ln(l + ~). Although most modern computers have a guard digit, there are a few (such as Crays) that do not. 1.5 Exactly Rounded Operations When floating-point operations are done with a guard digit, they are not as accurate as if they were computed exactly then rounded to the nearest floating-point number. Operations performed in this manner will be called exactly rounded. The example immediately preceding Theorem 2 shows that a single guard digit will not always give exactly rounded results. Section 1.4 gave several examples of algorithms that require a guard digit in order to work properly. This section gives examples of algorithms that require exact rounding. So far, the definition of rounding has not been given. Rounding is straightforward, with the exception of how to round halfway cases; for example, should 12.5 mnnd to 12 OP 12? Ofie whool of thought divides the 10 digits in half, letting {0, 1,2,3,4} round down and {5,6,’7,8,9} round up; thus 12.5 would round to 13. This is how rounding works on Digital Arithmetic 8 13 Equipment Corporation’s VAXG comput ers. Another school of thought says that since numbers ending in 5 are halfway between two possible roundings, they should round down half the time and round up the other half. One way of ob taining this 50’%0behavior is to require that the rounded result have its least significant digit be even. Thus 12.5 rounds to 12 rather than 13 because 2 is even. Which of these methods is best, round up or round to even? Reiser and Knuth [1975] offer the following reason for preferring round to even. Theorem 5 Let x and y be floating-point numbers, and define X. = x, xl=(xOey)O y,...,=(x(ley)@y)If@If@ and e are exactly rounded using round to even, then either x. = x for all n or x. = xl ❑ foralln >1. To clarify this result, consider ~ = 10, and let x = 1.00, y = –.555. When rounding up, the sequence becomes X. 9 Y = 1.56, Xl = 1.56 9 .555 = 1.01, xl e y ~ LO1 Q .555 = 1.57, and each successive value of x. increases by .01. Under round to even, x. is always 1.00. This example suggests that when using the round up rule, computations can gradually drift upward, whereas when using round to even the theorem says this cannot happen. Throughout the rest of this paper, round to even will be used. One application of exact rounding occurs in multiple precision arithmetic. There are two basic approaches to higher precision. One approach represents float ing-point numbers using a very large significant, which is stored in an array of words, and codes the routines for manipulating these numbers in assembly language. The second approach represents higher precision floating-point numbers as an array of ordinary floating-point p = 3 ‘VAX is a Corporation. trademark of Digital Equipment ACM Computmg Surveys, Vol 23, No. 1, March 1991 14 “ David Goldberg numbers, where adding the elements of the array in infinite precision recovers the high precision floating-point number. It is this second approach that will be discussed here. The advantage of using an array of floating-point numbers is that it can be coded portably in a high-level language, but it requires exactly rounded arithmetic. The key to multiplication in this system is representing a product xy as a sum, where each summand has the same precision as x and y. This can be done by splitting x and y. Writing x = x~ + xl and y = y~ + yl, the exact product is xy = xhyh + xhyl + Xlyh + Xlyl. If X and y the summands have p bit significands, prowill also have p bit significands, vided XI, xh, yh? Y1 carI be represented bits. When p is even, it is using [ p/2] easy to find a splitting. The number Xo. xl ““” xp_l can be written as the sum of Xo. xl ““” xp/2–l and O.O.. .OXP,Z . . . XP ~. When p is odd, this simple splitting method will not work. An extra bit can, however, be gained by using negative numbers. For example, if ~ = 2, P = 5, and x = .10111, x can be split as x~ = .11 and xl = – .00001. There is more than one way to split a number. A splitting method that is easy to compute is due to Dekker [1971], but it requires more than a single guard digit. Theorem Let the and are half 6 p be the floating-point precision, with restriction that p is even when D >2, assume that fl;ating-point operations exactly rounded. Then if k = ~p /2~ is the precision (rounded up) and m = fik + 1, x can je split as x = Xh + xl, where xh=(m Q9x)e (m@ Xe x), xl — — x e Xh, and each x, is representable using ~p/2] bits of precision. To see how this theorem works in an example, let P = 10, p = 4, b = 3.476, a = 3.463, and c = 3.479. Then b2 – ac rounded to the nearest floating-point number is .03480, while b @ b = 12.08, a @ c = 12.05, and so the computed value of b2 – ac is .03. This is an error of 480 ACM Computmg Surveys, Vol 23, No 1, March 1991 ulps. Using Theorem 6 to write b = 3.5 – .024, a = 3.5 – .037, and c = 3.5 – .021, b2 becomes 3.52 – 2 x 3.5 x .024 + .0242. Each summand is exact, so b2 where the = 12.25 – .168 + .000576, sum is left Similarly, ac = 3.52 unevaluated – (3.5 x at this .037 + 3.5 x point. .021) + .037 x .021 = 12.25 – .2030 + .000777. Finally, subtracting these two series term by term gives an estimate for b2 – ac of O @ .0350 e .04685 = .03480, which is identical to the exactly rounded result. To show that Theorem 6 really requires exact rounding, consider p = 3, P = 2, and x = 7. Then m = 5, mx = 35, and is performed m @ x = 32. If subtraction with a single guard digit, then ( m @ x) x~ = 4 and xl = 3, e x = 28. Therefore, ~~e xl not representable with \ p/2] = As a final example of exact rounding, m by 10. The result is consider dividing a floating-point number that will in general not be equal to m /10. When P = 2, m @10 by 10 will however, multiplying miraculously restore m, provided exact rounding is being used. Actually, a more general fact (due to Kahan) is true. The proof is ingenious, but readers not interested in such details can skip ahead to Section 2. Theorem 7 When O = 2, if m and n are integers ~m ~ < 2p-1 and n has the special n=2z+2J then (m On)@n=m, provided exactly fi?~ating-point operations with form are rounded. Scaling by a power of 2 is Proof harmless, since it changes only the exponent not the significant. If q = m /n, then scale n so that 2P-1 s n < 2P and scale m so that 1/2 < q < 1. Thus, 2P–2 < m < 2P. Since m has p significant bits, it has at most 1 bit to the right of the binary point. Changing the sign of m is harmless, so assume q > 0. Floating-Point If ij = m @ n, to prove requires showing that the theorem That is because m has at most 1 bit right of the binary point, so nij will round to m. TO deal with the halfway case when I T@ – m I = 1/4, note that since the inim had I m I < 2‘- 1, its tial unscaled low-order bit was O, so the low-order bit of the scaled m is also O. Thus, halfway cases will round to m. Suppose q = .qlqz “.. , and & g = . . . qP1. To estimate I nq – m 1, compute I ~ – q I = I N/2p+1 – ifs? N is an odd integer. m/nl, where and 2P-l l/(n2P+l-k). Assume q > Q is similar). Then Im-n@l= = n(q m-nij=n(q-@) – (~ – 2-P-1)) 1 < n 2–P–1 — ( = ~p+l-km n2 p (2 P-1 +2’)2-’-’ This establishes rem. ❑ +1–k ) +2-P-’+’=:. (9) and proves the theo- The theorem holds true for any base 6, as long as 2 z + 2 J is replaced by (3L + DJ. As 6 gets larger. however, there are fewer and fewer denominators of the form ~’ + p’. We are now in a position to answer the question, Does it matter if the basic Arithmetic g 15 arithmetic operations introduce a little more rounding error than necessary? The answer is that it does matter, because accurate basic operations enable us to prove that formulas are “correct” in the sense they have a small relative error. Section 1.4 discussed several algorithms that require guard digits to produce correct results in this sense. If the input to those formulas are numbers representing imprecise measurements, however, the bounds of Theorems 3 and 4 become less interesting. The reason is that the benign cancellation x – y can become catastrophic if x and y are only approximations to some measured quantity. But accurate operations are useful even in the face of inexact data, because they enable us to establish exact relationships like those discussed in Theorems 6 and 7. These are useful even if every floatingpoint variable is only an approximation to some actual value. 2. IEEE STANDARD There are two different IEEE standards for floating-point computation. IEEE 754 is a binary standard that requires P = 2, p = 24 for single precision and p = 53 for double precision [IEEE 19871. It also specifies the precise layout of bits in a single and double precision. IEEE 854 allows either L?= 2 or P = 10 and unlike 754, does not specify how floating-point numbers are encoded into bits [Cody et al. 19841. It does not require a particular value for p, but instead it specifies constraints on the allowable values of p for single and double precision. The term IEEE Standard will be used when discussing properties common to both standards. This section provides a tour of the IEEE standard. Each subsection discusses one aspect of the standard and why it was included. It is not the purpose of this paper to argue that the IEEE standard is the best possible floating-point standard but rather to accept the standard as given and provide an introduction to its use. For full details consult the standards [Cody et al. 1984; Cody 1988; IEEE 19871. ACM Computing Surveys, Vol 23, No 1, March 1991 16 2.1 David ● Formats 2. 1.1 and Goldberg Operations Base It is clear why IEEE 854 allows ~ = 10. Base 10 is how humans exchange and think about numbers. Using (3 = 10 is especially appropriate for calculators, where the result of each operation is displayed by the calculator in decimal. There are several reasons w~y IEEE 854 requires that if the base is not 10, it must be 2. Section 1.2 mentioned one reason: The results of error analyses are much tighter when ~ is 2 because a rounding error of 1/2 ulp wobbles by a factor of fl when computed as a relative error, and error analyses are almost always simpler when based on relative error. A related reason has to do with the effective precision for large bases. Consider fi = 16, p = 1 compared to ~ = 2, p = 4. Both systems have 4 bits of significant. Consider the computation of 15/8. When ~ = 2, 15 is represented as 1.111 x 23 and 15/8 as 1.111 x 2°, So 15/8 is exact. When p = 16, however, 15 is represented as F x 160, where F is the hexadecimal digit for 15. But 15/8 is represented as 1 x 160, which has only 1 bit correct. In general, base 16 can lose up to 3 bits, so a precision of p can have an effective precision as low as 4p – 3 rather than 4p. Since large values of (3 have these problems, why did IBM choose 6 = 16 for its system/370? Only IBM knows for sure, but there are two possible reasons. The first is increased exponent range. Single precision on the system/370 has ~ = 16, p = 6. Hence the significant requires 24 bits. Since this must fit into 32 bits, this leaves 7 bits for the exponent and 1 for the sign bit. Thus, the magnitude of representable numbers ranges from about 16-2’ to about 1626 = 228. To get a similar exponent range when D = 2 would require 9 bits of exponent, leaving only 22 bits for the significant. It was just pointed out, however, that when D = 16, the effective precision can be as low as 4p – 3 = 21 bits. Even worse, when B = 2 it is possible to gain an extra bit of ACM Computing Surveys, Vol 23, No 1, March 1991 precision (as explained later in this section), so the ~ = 2 machine has 23 bits of precision to compare with a range of 21-24 bits for the ~ = 16 machine. Another possible explanation for choosing ~ = 16 bits has to do with shifting. When adding two floating-point numbers, if their exponents are different, one of the significands will have to be shifted to make the radix points line up, slowing down the operation. In the /3 = 16, p = 1 system, all the numbers between 1 and 15 have the same exponent, so no shifting is required when adding 15 = 105 possible pairs of of the () distinct numb~rs from this set. In the however, these b = 2, P = 4 system, numbers have exponents ranging from O to 3, and shifting is required for 70 of the 105 pairs. In most modern hardware, the performance gained by avoiding a shift for a subset of operands is negligible, so the small wobble of (3 = 2 makes it the preferable base. Another advantage of using ~ = 2 is that there is a way to gain an extra bit of significance .V Since floating-point numbers are always normalized, the most significant bit of the significant is always 1, and there is no reason to waste a bit of storage representing it. Formats that use this trick bit. It was are said to have a hidden pointed out in Section 1.1 that this requires a special convention for O. The method given there was that an expoof all nent of e~,~ – 1 and a significant zeros represent not 1.0 x 2 ‘mln–1 but rather O. IEEE 754 single precision is encoded in 32 bits using 1 bit for the sign, 8 bits and 23 bits for the sigfor the exponent, nificant. It uses a hidden bit, howeve~, so the significant is 24 bits (p = 24), even though it is encoded using only 23 bits. any ‘This appears to have first been published by Goldberg [1967], although Knuth [1981 page 211] attributes this Idea to Konrad Zuse Floating-Point Arithmetic ● 17 den, the calculator presents a simple model to the operator. The IEEE standard defines four different Extended precision in the IEEE standprecision: single, double, single exard serves a similar function. It enables tended, and double extended. In 754, sinlibraries to compute quantities to within gle and double precision correspond about 1/2 ulp in single (or double) preciroughly to what most floating-point sion efficiently, giving the user of those hardware provides. Single precision oclibraries a simple model, namely, that cupies a single 32 bit word, double precieach primitive operation, be it a simple sion two consecutive 32 bit words. multiply or an invocation of log, returns Extended precision is a format that offers a value accurate to within about 1/2 ulp. just a little extra precision and exponent When using extended precision, however, range (Table 1). The IEEE standard only it is important to make sure that its use specifies a lower bound on how many is transparent to the user. For example, extra bits extended precision provides. on a calculator, if the internal represenThe minimum allowable double-extended tation of a displayed value is not rounded format is sometimes referred to as 80-bit to the same precision as the display, the even though the table shows it format, result of further operations will depend using 79 bits. The reason is that hardon the hidden digits and appear unpreware implementations of extended precidictable to the user. sion normally do not use a hidden bit and To illustrate extended precision furso would use 80 rather than 79 bits.8 ther, consider the problem of converting The standard puts the most emphasis between IEEE 754 single precision and on extended precision, making no recomdecimal. Ideally, single precision nummendation concerning double precision bers will be printed with enough digits so but strongly recommending that that when the decimal number is read Implementations should support the extended back in, the single precision number can format corresponding to the widest basic format be recovered. It turns out that 9 decimal supported, digits are enough to recover a single precision binary number (see Section 4.2). One motivation for extended precision When converting a decimal number back comes from calculators, which will often to its unique binary representation, a display 10 digits but use 13 digits internally. By displaying only 10 of the 13 rounding error as small as 1 ulp is fatal because it will give the wrong answer. digits, the calculator appears to the user Here is a situation where extended preci~ } a black box that computes exponenalgorithm. tial, cosines, and so on, to 10 digits of sion is vital for an efficient When single extended is available, a accuracy. For the calculator to compute straightforward method exists for confunctions like exp, log, and cos to within verting a decimal number to a single 10 digits with reasonable efficiency, howprecision binary one. First, read in the 9 ever, it needs a few extra digits with decimal digits as an integer N, ignoring which to work. It is not hard to find a the decimal point. From Table 1, p >32, simple rational expression that approxiand since 109 < 232 = 4.3 x 109, N can mates log with an error of 500 units in be represented exactly in single exthe last place. Thus, computing with 13 tended. Next, find the appropriate power digits gives an answer correct to 10 dig10P necessary to scale N. This will be a its. By keeping these extra 3 digits hidcombination of the exponent of the decimal number, and the position of the (up until now) ignored decimal point. *According to Kahan, extended precision has 64 Compute 10 I ‘l. If \ P I s 13, this is also bits of significant because that was the widest represented exactly, because 1013 = precision across which carry propagation could be 213513 and 513<232. Finally, multiply done on the Intel 8087 without increasing the cycle time [Kahan 19881. (or divide if P < 0) N and 10’ P‘. If this 2. 1.2 Precision ACM Computmg Surveys, Vol. 23, No. 1, March 1991 18 - David Goldberg Table 1. IEEE 754 Format Parameters Format Parameter Single P emax emln Exponent width in bits Format width in bits 24 + 127 – 126 8 32 Single Extended > 32 z + 1023 < – 1022 > 11 2 43 last operation is done exactly, the closest binary number is recovered. Section 4.2 shows how to do the last multiply (or divide) exactly. Thus, for I P I s 13, the use of the single-extended format enables 9 digit decimal numbers to be converted to the closest binary number (i. e., exactly rounded). If I P I > 13, singleextended is not enough for the above algorithm to compute the exactly rounded binary equivalent always, but Coonen [1984] shows that it is enough to guarantee that the conversion of binary to decimal and back will recover the original binary number. If double precision is supported, the algorithm above would run in double precision rather than single-extended, but to convert double precision to a 17 digit decimal number and back would require the double-extended format. 2.1.3 Exponent Since the exponent can be positive or negative, some method must be chosen to represent its sign. Two common methods of representing signed numbers are sign/magnitude and two’s complement. Sign/magnitude is the system used for the sign of the significant in the IEEE formats: 1 bit is used to hold the sign; the rest of the bits represent the magnitude of the number. The two’s complement representation is often used in integer arithmetic. In this scheme, a number is represented by the smallest nonnegative number that is congruent to it modulo 2 ~. The IEEE binary standard does not use either of these methods to represent the exponent but instead uses a- biased ACM Computmg Surveys, Vol 23, No 1, March 1991 Double 53 + 1023 – 1022 11 64 Double Extended > 64 > + 16383 < – 163$32 2 15 2 79 representation. In the case of single precision, where the exponent is stored in 8 bits, the bias is 127 (for double precisiog it is 1023). What this means is that if k is the value of the exponent bits interpreted as an unsigned integer, then the exponent of the floating-point number is ~ – 127. This is often called the biased to di~tinguish from the unbiexponent ased exponent k. An advantage of’ biased representation is that nonnegative flouting-point numbers can be treated as integers for comparison purposes. Referring to Table 1, single precision has e~~, = 127 and e~,~ = – 126. The reason for having I e~l~ I < e~,X is so that the reciprocal of the smallest number (1/2 ‘mm) will not overflow. Although it is true that the reciprocal of the largest number will underflow, underflow is usually less serious than overflow. Section 2.1.1 explained that e~,~ – 1 is used for representing O, and Section 2.2 will introduce a use for e~,X + 1. In IEEE single precision, this means that the biased exponents range between e~,~ – 1 = – 127 and e~.X + 1 = 128 whereas the unbiased exponents range between O and 255, which are exactly the nonnegative numbers that can be represented using 8 bits. 2. 1.4 Operations The IEEE standard requires that the result of addition, subtraction, multiplication, and division be exactly rounded. That is, the result must be computed exactly then rounded to the nearest floating-point number (using round to even). Section 1.3 pointed out that computing the exact difference or sum of two float- Floating-Point ing-point numbers can be very expensive when their exponents are substantially different. That section introduced guard digits, which provide a practical way of computing differences while guaranteeing that the relative error is small. Computing with a single guard digit, however, will not always give the same answer as computing the exact result then rounding. By introducing a second guard digit and a third sticky bit, differences can be computed at only a little more cost than with a single guard digit, but the result is the same as if the difference were computed exactly then rounded [Goldberg 19901. Thus, the standard can be implemented efficiently. One reason for completely specifying the results of arithmetic operations is to improve the portability of software. When a .Program IS moved between two machmes and both support IEEE arithmetic, if any intermediate result differs, it must be because of software bugs not differences in arithmetic. Another advantage of precise specification is that it makes it easier to reason about floating point. Proofs about floating point are hard enough without having to deal with multiple cases arising from multiple kinds of arithmetic. Just as integer programs can be proven to be correct, so can floating-point programs, although what is proven in that case is that the rounding error of the result satisfies certain bounds. Theorem 4 is an example of such a proof. These proofs are made much easier when the operations being reasoned about are precisely specified. Once an algorithm is proven to be correct for IEEE arithmetic, it will work correctly on any machine supporting the IEEE standard. Brown [1981] has proposed axioms for floating point that include most of the existing floating-point hardware. Proofs in this system cannot, however, verify the algorithms of Sections 1.4 and 1.5, which require features not present on all hardware. Furthermore, Brown’s axioms are more complex than simply defining operations to be performed exactly then rounded. Thus, proving theorems from Brown’s axioms is usually more difficult Arithmetic “ 19 than proving them assuming operations are exactly rounded. There is not complete agreement on what operations a floating-point standard should cover. In addition to the basic operations +, –, x, and /, the IEEE standard also specifies that square root, remainder, and conversion between integer and floating point be correctly rounded. It also requires that conversion between internal formats and decimal be correctly rounded (except for very large numbers). Kulisch and Miranker [19861 have proposed adding inner product to the list of operations that are precisely specified. They note that when inner products are computed in IEEE arithmetic, the final answer can be quite wrong. For example, sums are a special case of inner products, and the sum ((2 x 10-30 + 1030) – 10--30) – 1030 is exactly equal to 10- 30 but on a machine with IEEE arithme~ic the computed result will be – 10 – 30. It is possible to compute inner products to within 1 ulp with less hardware than it takes to implement a fast multiplier [Kirchner and Kulisch 19871.9 All the operations mentioned in the standard, except conversion between decimal and binary, are required to be exactly rounded. The reason is that efficient algorithms for exactly rounding all the operations, except conversion, are known. For conversion, the best known efficient algorithms produce results that are slightly worse than exactly rounded ones [Coonen 19841. The IEEE standard does not require transcendental functions to be exactly rounded because of the table maker’s To illustrate, suppose you are dilemma. making a table of the exponential function to four places. Then exp(l.626) = 5.0835. Should this be rounded to 5.083 or 5.084? If exp(l .626) is computed more it becomes 5.08350, then carefully, ‘Some arguments against including inner product as one of the basic operations are presented by Kahan and LeBlanc [19851. ACM Computing Surveys, Vol 23, No 1, March 1991 20 “ David Goldberg then 5.0835000. Since exp is 5.083500, transcendental, this could go on arbitrarily long before distinguishing whether or exp(l.626) is 5.083500 “ “ “ O ddd 5.0834999 “ “ “ 9 ddd. Thus, it is not practical to specify that the precision of transcendental functions be the same as if the functions were computed to infinite precision then rounded. Another approach would be to specify transcendental functions algorithmically. But there does not appear to be a single algorithm that works well across all hardware architectures. Rational approximation, CORDIC,1° and large tables are three different techniques used for computing transcendental on contemporary machines. Each is appropriate for a different class of hardware, and at present no single algorithm works acceptably over the wide range of current hardware. 2.2 Special Quantities On some floating-point hardware every bit pattern represents a valid floatingpoint number. The IBM System/370 is an example of this. On the other hand, the VAX reserves some bit patterns to rerepresent special numbers called This idea goes back to served operands. the CDC 6600, which had bit patterns for the special quantities INDEFINITE and INFINITY. The IEEE standard continues in this tradition and has NaNs (Not a Number, pronounced to rhyme with plan) and infinities. Without special quantities, there is no good way to handle exceptional situations like taking the square root of a negative number other than aborting computation. Under IBM System/370 FORTRAN, the default action in response to computing the square root of a negative number like – 4 results in the printing of an error message. Since every 10CORDIC is an acronym for Coordinate Rotation Digital Computer and is a method of computing transcendental funct~ons that uses mostly shifts and adds (i. e., very few multiplications and divisions) [Walther 1971], It is the method used on both the Intel 8087 and the Motorola 68881. ACM Computmg Surveys, Vol 23. No 1, March 1991 Table 2. Exponent ~=~ nun –1 e = ‘ml. -1 e~,n 5 e 5 emax e=emay+l e=g ~ay + 1 IEEE 754 Fraction f=o f#o f:o f#o Special Values Represents *O O fx 2’mLn 1 fx2’ N%; bit pattern represents a valid number, the return value of square root must be some floating-point number. In the case of System/370 FORTRAN, ~ = 2 is returned. In IEEE arithmetic, an NaN is returned in this situation. The IEEE standard specifies the following special values (see Table 2): f O, denormalized numbers, + co and NaNs (there is more than one NaN, as explained in the next section). These are all encoded with special values exponents of either e~.X + 1 or e~,~ – 1 (it was already pointed out that O has an exponent of e~,. – 1). 2.2.1 NaNs Traditionally, the computation of 0/0 or 4 – 1 has been treated as an unrecoverable error that causes a computation to halt. There are, however, examples for which it makes sense for a computation to continue in such a situation. Consider a subroutine that finds the zeros of a function f, say zero(f). Traditionally, zero finders require the user to input an is interval [a, b] on which the function defined and over which the zero finder will search. That is, the subroutine is called as zero(f, a, b). A more useful zero finder would not require the user to input this extra information. This more general zero finder is especially appropriate for calculators, where it is natural to key in a function and awkward to then have to specify the domain. It is easy, however, to see why most zero finders require a domain. The zero finder does its work by probing the function f at various values. If it probed for a value outside the domain of f, the code for f Floating-Point Table 3. Operation + x I REM \ Operations that Produce an NaN NaN Produced by W+(–w) Oxw 0/0, cO/03 x REM O, m REM y fi(when x < O) might well compute 0/0 or ~, and the computation would halt, unnecessarily aborting the zero finding process. This problem can be avoided by introducing a special value called NaN and specifying that the computation of expressions like 0/0 and ~ produce NaN rather than halting. (A list of some of the situations that can cause a NaN is given in Table 3.) Then, when zero(f) probes outside the domain of f, the code for f will return NaN and the zero finder can continue. That is, zero(f) is not “punished” for making an incorrect guess. With this example in mind, it is easy to see what the result of combining a NaN with an ordinary floating-point number should be. Suppose the final statement off is return( – b + sqrt(d))/ (2* a). If d <0, then f should return a sqrt(d) is an NaN, NaN. Since d <0, and – b + sqrt(d) will be a NaN if the sum of an NaN and any other number is a NaN. Similarly, if one operand of a division operation is an NaN, the quotient should be a NaN. In general, whenever a NaN participates in a floating-point operation, the result is another NaN. Another approach to writing a zero solver that does not require the user to input a domain is to use signals. The zero finder could install a signal handler for floating-point exceptions. Then if f were evaluated outside its domain and raised an exception, control would be returned to the zero solver. The problem with this approach is that every language has a different method of handling signals (if it has a method at all), and so it has no hope of portability. In IEEE 754, NaNs are represented as floating-point numbers with the expo- Arithmetic 9 21 nent e~~X + 1 and nonzero significands. Implementations are free to put systemdependent information into the significant. Thus, there is not a unique NaN but rather a whole family of NaNs. When an NaN and an ordinary floating-point number are combined, the result should be the same as the NaN operand. Thus, if the result of a long computation is an NaN, the system-dependent information in the significant will be the information generated when the first NaN in the computation was generated. Actually, there is a caveat to the last statement. If both operands are NaNs, the result will be one of those NaNs but it might not be the NaN that was generated first. 2.2.2 Infinity Just as NaNs provide a way to continue a computation when expressions like 0/0 or ~ are encountered, infinities provide a way to continue when an overflow occurs. This is much safer than simply returning to the largest representable number. As an example, consider computing ~~, when b = 10, p = 3, and e~~X = 98. If x = 3 x 1070 and y = 4 X 1070, th en X2 will overflow and be replaced by 9.99 x 1098. Similarly yz and X2 + yz will each overflow in turn and be replaced by 9.99 x 1098. So the final result will be (9.99 x 1098)112 = 3.16 x 1049, which is drastically wrong. The correct answer is 5 x 1070. In IEEE arithmetic, the result of X2 is CO,as is yz, X2 + yz, and -. SO the final result is m, which is safer than returning an ordinary floating-point number that is nowhere near the correct answer.” The division of O by O results in an NaN. A nonzero number divided by O, returns infinity: 1/0 = ~, however, – 1/0 = – co. The reason for the distinction is this: If f(x) -0 and g(x) + O as llFine point: Although the default in IEEE arithmetic is to round overflowed numbers to ~, it is possible to change the default (see Section 2.3.2). ACM Computing Surveys, Vol 23, No 1, March 1991 22 “ David Goldberg x approaches some limit, then f( x)/g( x) could have any value. For example, when f’(x) = sin x and g(x) = x, then ~(x)/g(x) ~ 1 as x + O. But when ~(x) ~ O. When =l– COSX, f(x)/g(x) thinking of 0/0 as the limiting situation of a quotient of two very small numbers, 0/0 could represent anything. Thus, in the IEEE standard, 0/0 results in an NaN. But when c >0 and f(x) ~ c, g(x) ~ O, then ~(x)/g(*) ~ * m for any analytic functions f and g. If g(x) <0 for small x, then f(x)/g(x) ~ – m; otherwise the limit is + m. So the IEEE standard defines c/0 = & m as long as c # O. The sign of co depends on the signs of c and O in the usual way, so – 10/0 = – co +m. You can distinand –10/–0= guish between getting m because of overflow and getting m because of division by O by checking the status flags (which will be discussed in detail in Section 2.3.3). The overflow flag will be set in the first case, the division by O flag in the second. The rule for determining the result of an operation that has infinity as an operand is simple: Replace infinity with a finite number x and take the limit as x + m. Thus, 3/m = O, because Iim ~+~3/x = O. Similarly 4 – co = – aI and G = w. When the limit does not exist, the result is an NaN, so m/co will be an NaN (Table 3 has additional examples). This agrees with the reasoning used to conclude that 0/0 should be an NaN. When a subexpression evaluates to a NaN, the value of the entire expression is also a NaN. In the case of & w, however, the value of the expression might be an ordinary floating-point number because of rules like I/m = O. Here is a practical example that makes use of the rules for infinity arithmetic. Consider computing the function x/( X2 + 1). This is a bad formula, because not only will it than overflow when x is larger iz but infinity arithmetic will fib’”” give the &rong answer because it will yield O rather than a number near 1/x. However, x/( X2 + 1) can be rewritten as 1/( x + x- l). This improved expression will not overflow prematurely and because of infinity arithmetic will have the ACM Computmg Surveys, Vol 23, No. 1, March 1991 correct value when x = O: 1/(0 + 0-1) = 1/(0 + CO)= l/CO = O. Without infinity arithmetic, the expression 1/( x + x-1) requires a test for x = O, which not only adds extra instructions but may also disrupt a pipeline. This example illustrates that infinity a general fact; namely, arithmetic often avoids the need for spe cial case checking; however, formulas need to be carefully inspected to make sure they do not have spurious behavior at infinity [as x/(X2 + 1) did]. 2.2.3 Slgnea Zero Zero is represented by the exponent Since the emm – 1 and a zero significant. sign bit can take on two different values, there are two zeros, + O and – O. If a distinction were made when comparing -t O and – O, simple tests like if (x = O) would have unpredictable behavior, depending on the sign of x. Thus, the IEEE standard defines comparison so that +0= –O rather than –O< +0. Although it would be possible always to ignore the sign of zero, the IEEE standard does not do so. When a multiplication or division involves a signed zero, the usual sign rules apply in computing the sign of the answer. Thus, 3(+ O) = -t O and +0/– 3 = – O. If zero did not have a sign, the relation 1/(1 /x) = x would fail to hold when x = *m. The reason is that 1/– ~ and 1/+ ~ both result in O, and 1/0 results in + ~, the sign information having been lost. One way to restore the identity 1/(1 /x) = x is to have only one kind of’ infinity; however, that would result in the disastrous consequence of losing the sign of an overflowed quantity. Another example of the use of signed zero concerns underflow and functions that have a discontinuity at zero such as log. In IEEE arithmetic, it is natural to define log O = – w and log x to be an NaN whe”n x <0. Suppose” x represents a small negative number that has underflowed to zero. Thanks to signed zero, x will be negative so log can return an NaN. If there were no signed zero, however, the log function could not Floating-Point distinguish an underflowed negative number from O and would therefore have to return – m. Another example of a function with a discontinuity at zero is the signum function, which returns the sign of a number. Probably the most interesting use of signed zero occurs in complex arithmetic. As an example, consider the equation ~ = ~/&. This is certainly true when z = O. If z = —1. the obvious computation gives ~~ = ~ = i and I/i = –i. Thus, ~# I/n= 1/W ! The problem can be traced to the fact that square root is multivalued, and there is no way to select the values so they are continuous in the entire complex plane. Square root is continuous, however, if a branch cut consisting of all negative real numbers is excluded from consideration. This leaves the problem of what to do for the negative real numbers, which are of the form – x + iO, where x > 0. Signed zero provides a perfect way to resolve this problem. Numbers of the form – x + i( + O) have a square root of i&, and numbers of the form – x + i( – O) on the other side of the branch cut have a square root with the other sign In fact, the natural formulas for (– i ~). computing ~ will give these results. Let us return to ~ = l/fi. If z = –1= –l+iO, then 1/2 = 1/(-1 + iO) 1(-1 -iO) — — (-1+ iO)(-1-iO) = (-1 - iO)/(( -1)2 - 02) = –l+i(–0), so ~= – 1+ i(–0) = –i, while Thus, IEEE arithI/&= l/i = –i, metic preserves this identity for all z. Some more sophisticated examples are given by Kahan [1987]. Although distinguishing between + O and – O has advantages, it can occasionally be confusing. For example, signed zero destroys the relation x = y * I/x = l/y, which is false when x = +0 and y = –O. The Arithmetic “ 23 IEEE committee decided, however, that the advantages of using signed zero outweighed the disadvantages. 2.2.4 Denormalized Numbers Consider normalized floating-point numbers with O = 10, p = 3, and e~,. = –98. The numbers % = 6.87 x 10-97 and y = 6.81 x 10-97 appear to be perfectly ordinary floating-point numbers, which are more than a factor of 10 larger than the smallest floating-point number 1.00 x 10-98. They have a strange property, however: x 0 y = O even though x # y! The reason is that x – y = .06 x 10-97 — — 6.0 x 10- ‘g is too small to be represented as a normalized number and so must be flushed to zero. How important is it to preserve the property X=yex–y=o? (lo) It is very easy to imagine writing the z = 1/ code fragment if (x # y) then a program fail (x – y) and later having due to a spurious division by zero. Tracking down bugs like this is frustrating and time consuming. On a more philosophical level, computer science text books often point out that even though it is currently impractical to prove large programs correct, designing programs with the idea of proving them often re suits in better code. For example, introducing invariants is useful, even if they are not going to be used as part of a proof. Floating-point code is just like any other code: It helps to have provable facts on which to depend. For example, when analyzing formula (7), it will be helpful toknowthat x/2