Web And Email Hosting, Jsp Web Hosting Blog

4.2.2 ACCURACY OF FLOATING POINT ARITHMETIC 227 The relations <, -, >, z introduced in this section are similar to ideas published by A. van Wijngaarden in BIT 6 (1966), 66-81. Theorems A and B above were inspired by some related work of Ole Mgller, BIT 5 (1965), 37-50, 251-255; Theorem C is due to T. J. Dekker, Numer. Math. 18 (1971), 224-242. Extensions and refinements of all three theorems have been published by S. Linnainmaa, BIT 14 (1974), 167-202. W. M. Kahan introduced Theorem D in some unpublished notes; for a complete proof and further commentary, see J. F. Reiser and D. E. Knuth, Inf. Proc. Letters 3 (1975), 84-87, 164. Unnormalized floating point arithmetic was recommended by F. L. Bauer and K. Samelson in the article cited above, and it was independently used by J. W. Carr III at the University of Michigan in 1953. Several years later, the MANIAC III computer was designed to include both kinds of arithmetic in its hardware; see R. L. Ashenhurst and N. Metropolis, JACA4 6 (1959), 415- 428, IEEE Trans. EC-12 (1963), 896-901; R. L. Ashenhurst, Proc. Spring Joint Computer Conf. 21 (1962), 195-202. See also H. L. Gray and C. Harrison, Jr., Proc. Eastern Joint Computer Conf. 16 (1959), 244-248, and W. G. Wadey, JACM 7 (1960), 129-139, for further early discussions of unnormalized arith- metic. For early developments in interval arithmetic, and some modifications, see A. Gibb, CACM 4 (1961), 319-320; B. A. Chartres, JACM 13 (1966), 386-403; and the book Interval Analysis by Ramon E. Moore (Prentice-Hall, 1966). The subsequent flourishing of this subject is described in Moore s later book, Methods and Applications of Interval Analysis (SIAM, 1979). The book Grundlagen des Numerischen Rechnens: Mathematische Begrund- ung der Rechenarithmetik by Ulrich Kulisch (Mannheim: Bibl. Inst., 1976) is entirely devoted to the study of floating point arithmetic systems; see also Kulisch s article in IEEE Trans. C-26 (1977), 610-621, and his more recent book written jointly with W. L. Miranker, entitled Computer Arithmetic in Theory and Practice (New York: Academic Press, 1980). EXERCISES Note: Normalized floating point arithmetic is assumed unless the contrary is specified. 1. [MS?] Prove that identity (7) is a consequence of (2) through (6). 2. [A&O] Use identities (2) through (8) to prove that (u @ X) @ (w @ y) 2 u @ 2, whenever x 2 0 and y 2 0. 3. [M.%?] Find eight-digit floating decimal numbers u, 21, and w such that and such that no exponent overflow or underflow occurs during the computations. 4. [IO] Is it possible to have floating point numbers u, v, and w for which exponent overflow occurs during the calculation of u @I (v @ w) but not during the calculation of (u @ ?J) @w?

226 ARITHMETIC 4.2.2 their sum and product would then turn out to be N $ h = [(24, +.60222400), (24, +.60228001)], N @ h = [(-3, +.63500305), (-3, +.63514642)]. If we try to divide by [ZIO, wi] when wo < 0 < ~1, there is a possibility of division by zero. Since the philosophy underlying interval arithmetic is to provide rigorous error estimates, a divide-by-zero error should be signalled in this case. However, overflow and underflow need not be treated as errors in interval arithmetic, if special conventions are introduced as discussed in exercise 24. Interval arithmetic takes only about twice as long as ordinary arithmetic, and it provides truly reliable error estimates. Considering the difficulty of mathematical error analyses, this is indeed a small price to pay. Since the intermediate values in a calculation often depend on each other, as explained above, the final estimates obtained with interval arithmetic will tend to be pessimistic; and iterative numerical methods often have to be redesigned if we want to deal with intervals. The prospects for effective use of interval arithmetic look very good, however, and efforts should be made to increase its availability. D. History and bibliography. Jules Tannery s classic treatise on decimal calcula- tions, LeCons d Arithm&ique (Paris: Colin, 1894), stated that positive numbers should be rounded upwards if the first discarded digit is 5 or more; since exactly half of the decimal digits are 5 or more, he felt that this rule would round up- wards exactly half of the time, on the average, so it would produce compensating errors. The idea of round to even in the ambiguous cases seems to have been mentioned first by James B. Scarborough in the first edition of his pioneering book Numerical Mathematical Analysis (Baltimore: Johns Hopkins Press, 1930), p. 2; in the second (1950) edition he amplified his earlier remarks, stating that It should be obvious to any thinking person that when a 5 is cut off, the preceding digit should be increased by 1 in only half the cases, and he recommended round-to-even in order to achieve this. The first analysis of floating point arithmetic was given by F. L. Bauer and K. Samelson, Zeitschrift fur angewandte Math. und Physik 4 (1953), 312- 316. The next publication was not until over five years later: J. W. Carr III, CACM 2, 5 (May 1959), 10-15. See also P. C. Fischer, Proc. ACM Nat. Meeting 13 (Urbana, Illinois, 1958), paper 39. The book Rounding Errors in Algebraic Processes (Englewood Cliffs: Prentice-Hall, 1963), by J. H. Wilkinson, shows how to apply error analysis of the individual arithmetic operations to the error analysis of large-scale problems; see also his treatise on The Algebraic Eigenvalue Problem (Oxford: Clarendon Press, 1965). More recent work on floating point accuracy is summarized in two important papers that can be especially recommended for further study: W. M. Kahan, Proc. IFLF Congress (1971), 2, 1214-1239; R. P. Brent, LEEE Trans. C-22 (1973), 601-607. Both papers include useful theory and demonstrate that it pays off in practice.

4.2.2 ACCURACY OF FLOATING POINT ARITHMETIC 225 find that (52) is valid when E 2 2b2-p; unnormalized addition is not as erratic as normalized addition with respekt to the associative law. It should be emphasized that unnormalized arithmetic is by no means a panacea. There are examples where it indicates greater accuracy than is present (e.g., addition of a great many small quantities of about the same magnitude, or evaluation of xn for large n); and there are many more examples when it indicates poor accuracy while normalized arithmetic actually does produce good results. There is an important reason why no straightforward one-operation-at- a-time method of error analysis can be completely satisfactory, namely the fact that operands are usually not independent of each other. This means that errors tend to cancel or reinforce each other in strange ways. For example, suppose that 2 is approximately f , and suppose that we have an approximation y = z + 6 with absolute error 6. If we now wish to compute ~(1 -x), we can form y(1 -y); if x = 4 + t we find y(1 -y) = ~(1 -2) -2t6 -b2, so the error has decreased substantially: it has been multiplied by a factor of 2~ + 6. This is just one case where multiplication of imprecise quantities can lead to a quite accurate result when the operands are not independent of each other. A more obvious example is the computation of x 8 x, which can be obtained with perfect accuracy regardless of how bad an approximation to x we begin with. The extra information that unnormalized arithmetic gives us can often be more important than the information it destroys during an extended calcula- tion, but (as usual) we mu& use it with care. Examples of the proper use of unnormalized arithmetic are discussed by R. L. Ashenhurst and N. Metropolis in Computers and Computing, AMM Slaught Memorial Papers 10 (February, 1965), 47-59; by N. Metropolis in Numer. Math. 7 (1965), 104-112; and by R. L. Ashenhurst in Error in Digital Computation 2, ed. by L. B. Rall (New York: Wiley, 1965), 3-37. Appropriate methods for computing standard mathematical functions with both input and output in unnormalized form are given by R. L. Ashenhurst in JACK 11 (1964), 168-187. An extension of unnormalized arith- metic, which remembers that certain values are known to be exact, has been discussed by N. Metropolis in IEEE Trans. C-22 (1973), 573-576. C. Interval arithmetic. Another approach to the problem of error determination is the so-called interval or range arithmetic, in which upper and lower bounds on each number are maintained during the calculations. Thus, for example, if we know that uo 2 u 5 ~1 and ~0 5 w 5 ~1, we represent this by the interval notation u. = [uo, ~11, v = [vo, WI]. The sum U@V is [ug v ~0, ~1 &, ~11, where v denotes lower floating point addition, the greatest representable number less than or equal to the true sum, and A is defined similarly (see exercise 4.2.1-13). Furthermore u 8 u = [UO 7 ul, ~1 a ~01; and if uo and ~0 are positive, we have u @ u = [ug v vo, u1 a ~~1, u @ u = [u. v wl, u1 h wo]. For example, we might represent Avogadro s number and Planck s constant as N = [(24, +.60222400), (24, +.60228000)], h = [(-26, +.10544300), (-26, +.10545700)];

224 ARITHMETIC 4.2.2 just as in Algorithm 4.2.1A, except that all scaling to the left is suppressed. Multiplication and division are performed as in Algorithm 4.2.1M, except that the answer is scaled right or left so that precisely m&x&, lV) leading zeros appear. Essentially the same rules have been used in manual calculation for many years. It follows that, for unnormalized computations, euev, euev = max(e,, ev) + (0 or 1) (48) eumv = eu + ev -4 - min(l,, lV) - (0 or 1) (49) euBv = eu -e, + q - L + L + max(L, L) + (0 or 1). (50) When the result of a calculation is zero, an unnormalized zero (often called an order of magnitude zero ) is given as the answer; this indicates that the answer may not truly be zero, we just don t know any of its significant digits. Error analysis takes a somewhat different form with unnormalized floating point arithmetic. Let us define 6, = w -q-P if u = (e,, fU). (51) This quantity depends on the representation of u, not just on the value beu-qfu. Our rounding rule tells us that These inequalities apply to normalized as well as unnormalized arithmetic; the main difference between the two types of error analysis is the definition of the exponent of the result of each operation (Eqs. (48) to (50)). We have remarked that the relations 4, -, %, and % defined earlier in this section are valid and meaningful for unnormalized numbers as well as for normalized numbers. As an example of the use of these relations, let us prove an approximate associative law for unnormalized addition, analogous to (39): for suitable E. We have A similar formula holds for [U $ (V $ w) - (U + v + w)l. Now since e(,e,)e, = max(e,,e,,e,) + (0, 1, or 2), we have 6(,~+)@,) 5 b26,a(,B,). Therefore we

4.2.2 ACCURACY OF FLOATING POINT ARITHMETIC 223 hardware, or a small percentage of the average running time of a subroutine. But the above discussion shows that such a policy is mistaken. We could save about five percent of the running time of the FADD subroutine, Program 4.2.1A, and about 25 percent of its space, if we took the liberty of rounding incorrectly in a few cases, but we are much better off leaving it as it is. The reason is not to glorify bit chasing , . a more fundamental issue is at stake here: Numerical subroutines should deliver results that satisfy simple, useful mathematical Jaws whenever possible. The crucial formula u $ v = round(u. + V) is a regularity property that makes a great deal of difference between whether mathematical analysis of computational algorithms is worth doing or worth avoiding. Without any underlying symmetry properties, the job of proving interesting results becomes extremely unpleasant. The enjoyment of one s tools is an essential ingredient of successful work. B. Unnormalized floating point arithmetic. The policy of normalizing all floating point numbers may be construed in two ways: We may look on it favorably by saying that it is an attempt to get the maximum possible accuracy obtainable with a given degree of precision, or we may consider it to be potentially dangerous since it tends to imply that. the results are more accurate than they really are. When we normalize the result of (1, +.31428571) 8 (1, +.31415927) to (-2, +.12644000), we are suppressing information about the possibly greater inaccuracy of the latter quantity. Such information would be retained if the answer were left as (1, +.00012644). The input data to a problem is frequently not known as precisely as the floating point representation allows. For example, the values of Avogadro s number and Planck s constant are not known to eight significant digits, and it might be more appropriate to denote them, respectively, by (27, +.00060225) and (-23, +.00010545) instead of by (24, +.60225200) and (-26, +.10545000). It would be nice if we could give our input data for each problem in an unnormalized form that expresses how much precison is assumed, and if the output would indicate just how much precision is known in the answer. Unfortunately, this is a terribly difficult problem, although the use of unnormalized arithmetic can help to give some indication. For example, we can say with a fair degree of certainty that the product of Avogadro s number by Planck s constant is (0, +.00063507), and that their sum is (27, +.00060225). (The purpose of this example is not to suggest that any important physical significance should be attached to the sum and product of these fundamental constants; the point is that it is possible to preserve a little of the information about precision in the result of calculations with imprecise quantities, when the original operands are independent of each other.) The rules for unnormalized arithmetic are simply this: Let 1, be the number of leading zeros in the fraction part of u = (ezL, fu), so that 1, is the largest integer 2 p with if211 < b-l-. Then addition and subtraction are performed

222 ARITHMETIC 4.2.2 we floating-subtract TJ from this result we get u = 1.0000001. Adding and subtracting w from U gives 1.0000002, and the next time we get 1.0000003, etc.; the result keeps growing although we are adding and subtracting the same value. This phenomenon, called drift, will not occur when we use a stable rounding rule based on the parity of the least significant digit. More precisely: Theorem D. (((u $ v) 8 ZJ) $ v) 8 II = (u $ v) 8 u. For example, if u = 1.2345679 and v = -0.23456785, we find u $ v = 1.0000000, (u $ V) 0 w = 1.2345678, ((u $ w) 8 w) $ w = 0.99999995, and (((u $ v) 8 V) $ w) 0 v = 1.2345678. The proof for general u and ZJ seems to require a case analysis even more detailed than that in the above theorems; see the references at the end of this section. 1 Theorem D is valid both for round to even and round to odd ; how should we choose between these possibilities? When the radix b is odd, ambiguous cases never arise except during floating point division, and the rounding in such cases is comparatively unimportant. For even radices, there is reason to prefer the following rule: Round to even when b/2 is odd, round to odd when b/2 is even. The least significant digit of a floating point fraction occurs frequently as a remainder to be rounded off in subsequent calculations, and this rule avoids generating the digit b/2 in the least significant position whenever possible; its effect is to provide some memory of an ambiguous rounding so that subsequent rounding will tend to be unambiguous. For example, if we were to round to odd in the decimal system, repeated rounding of the number 2.44445 to one less place each time leads to the sequence 2.4445, 2.445, 2.45, 2.5, 3; but if we round to even, such situations do not occur. [Roy A. Keir, hf. Proc. Letters 3 (1975), 188-189.1 On the other hand, some people prefer rounding to even in all cases, so that the remainder will tend to be 0 more often. Neither alternative conclusively dominates the other; fortunately the base is usually b = 2 or b = 10, when everyone agrees that round-to-even is best. A reader who has checked some of the details of the above proofs will realize the immense simplification that has been afforded by the simple rule u $ w = round(u + v). If our floating point addition routine would fail to give this result even in a few rare cases, the proofs would become enormously more complicated and perhaps they would even break down completely. Theorem B fails if truncation arithmetic is used in place of rounding, i.e., if we let u $ v = trunc(u + TJ) and u 8 v = trunc(u -u), where trunc(z) takes all positive real 5 into the largest floating point number 2 z. An ex- ception to Theorem B would then occur for cases such as (20, +.10000001) $ (10, -.10000001) = (20, +.lOOOOOOO), when the difference between u + v and u $ w cannot be expressed exactly as a floating point number; and also for cases such as 12345678 $ .012345678, when it can be. Many people feel that, since floating point arithmetic is inexact by nature, there is no harm in making it just a little bit less exact in certain rather rare cases, if it is convenient to do so. This policy saves a few cents in the design of computer

4.2.2 ACCURACY OF FLOATING POINT ARITHMETIC 221 Using the notation of the preceding proof, these relations are respectively equiv- alent to round(t + T) = t + r, round(t) = t, round(r) = r. (46) Exercise 12 establishes the theorem in the special case le, -e,l 2 p. Otherwise u + 2, has at most 2p significant digits and it is easy to see that round(r) = P. If now e > 0, the proof of Theorem A shows that t = –r or t = r = f a. If e 5 0 we have t + r = u and t = –2) (modulo be); this is enough to prove that t + r and r round to themselves, provided that e, 2 e and e, 2 e. But either e, < 0 or e, < 0 would contradict our hypothesis that le, -e,l < p, since e, = p. 1 Theorem B gives an explicit formula for the difference between u + w and u $21, in terms of quantities that can be calculated directly using five operations of floating point arithmetic. If the radix b is 2 or 3, we can improve on this result, obtaining the exact value of the correction term with only two floating point operations and one (fixed point) comparison of absolute values: Theorem C. If b 5 3 and IuI 2 1~1, then Proof. Following the conventions of preceding proofs again, we wish to show that w 0 w = r. It suffices to show that w = w-u, because (46) will then yield w 0 w = round(v -w ) = round(u + w - w) = round(r) = r. We shall in fact prove (47) whenever b 5 3 and eu 2 e,. If e, 2 p, then r is a tail of v modulo 1, hence w = w 8 u = w 0 r = w -r = w -u as desired. If e, < p, then we must have e, = p - 1, and w -u is a multiple of b-l; it will therefore round to itself if its magnitude is less than bp-l + b-l. Since b 5 3, wehaveindeed Iw–uI 5 Iw-u-wI+IwI 5 &+(bp- -b-l) < bPpl+bpl. This completes the proof. 1 The proofs of Theorems A, B, and C do not rely on the precise definitions of round(z) in the ambiguous cases when x is exactly midway between consecutive floating point numbers; any way of resolving the ambiguity will suffice for the validity of everything we have proved so far. No rounding rule can be best for every application. For example, we gener- ally want a special rule when computing our income tax. But for most numerical calculations the best policy appears to be the rounding scheme specified in Algorithm 4.2.1N, which insists that the least significant digit should always be made even (or always odd) when an ambiguous value is rounded. This is not a trivial technicality, of interest only to nit-pickers; it is an important practical consideration, since the ambiguous case arises surprisingly often and a biased rounding rule produces significantly poor results. For example, consider decimal arithmetic and assume that remainders of 5 are always rounded upwards. Then if u = 1.0000000 and v = 0.55555555 we have u $ w = 1.5555556; and if

220 ARITHMETIC 4.2.2 This rather cumbersome-looking identity can be rewritten in a simpler manner: Let u = (u $ ?J) 8 v , 21 = (u $ v) 8 u ; (41) U = (u $ ?J) e w , v = (u @ v) e 21 . Intuitively, u and u should be approximations to u, and v and v should be approximations to V. Theorem A tells us that u$w=u +w =u +w . (42) This is a stronger statement than the identity u $ v = u $21 = 21 $ v , (43) which follows by rounding (42). Proof. Let us say that t is a tail of z modulo b if t E z (modulo b ), (44) thus, z - round(z) is always a tail of z. The proof of Theorem A rests largely on the following simple fact proved in exercise 11: Lemma T. If t is a tail of the floating point number z, then 2 0 t = z - t. u Let w = u $ V. Theorem A holds trivially when w = 0. By multiplying all variables by a suitable power of b, we may assume without loss of generality that e, = p. Then u + v = w + T, where r is a tail of u + v modulo 1. Furthermore u = round(w -V) = round(u -r) = u -r -t, where t is a tail of u - r modulo be and e = eUl -p. If e 5 0, then t E u - T z –21 (modulo be), hence t is a tail of -V and v = round(w -u ) = round(w + t) = v + t; this proves (40). If e > 0, then Iu - rl 2 bP -4; and since Irl <_ 4, we have IuI 2 bP - 1. It follows that T is a tail of w modulo 1. If u = u, we have v = round(w -u) = round(w -r) = v - r. Otherwise the relation round(u -r) # u implies that IuJ = bP - 1, (rl = 4, (u ( = bp; we have ZI = round(w -u ) = round(w + r) = w + T, because r is also a tail of -v in this case. 1 Theorem A exhibits a regularity property of floating point addition, but it doesn t seem to be an especially useful result. The following identity is more significant: Theorem B. Under the hypotheses of Theorem A and (41), (45) Proof. In fact, we can show that u 0 u = u -u , v 8 v = v -v , and (u-u )$(v-d ) = (u-u/)+(2,-I! ), h ence (45) will follow from Theorem A.

4.2.2 ACCURACY OF FLOATING POINT ARITHMETIC 219 Moreover, we can prove without difficulty that Iu-w) 5 ~1~1 and /U-VI 5 ~1~1 implies u c=x v (6); (34) Iu-uyl 5 ElUl or Iu -VI 5 tlwl implies u - v (t); (35) and conversely, for normalized floating point numbers u and w, when E < 1, u Fx 21 (t) implies Iu-VI 2 b~luj and Iu-u[ 2 b+[; (36) u -v (6) implies Iu - 211 5 b+L or Iu -VI 5 b+l. (37) Let ~0 = bl-p be one ulp. The derivation of (17) establishes the inequality IX -round(z)1 = Ip( < &EO min(lzl, Iround(z)l), hence x zz round(s) (4~0); (38) it follows that u $ w z u + v ($E~), etc. The approximate associative law for multiplication derived above can be recast as follows: We have by (19), and the same inequality is valid with (u @ U) @ w and u @ (V @ w) interchanged. Hence by (34), whenever E 2 2~/(1 -&CO) . For example, if b = 10 and p = 8 we may take E = 0.00000021. The relations 4, -, +, and M are useful within numerical algorithms, and it is therefore a good idea to provide routines for comparing floating point numbers as well as for doing arithmetic on them. Let us now shift our attention back to the question of finding exact relations that are satisfied by the floating point operations. It is interesting to note that floating point addition and subtraction are not completely intractable from an axiomatic standpoint, since they do satisfy the nontrivial identities stated in the following theorems. Theorem A. Let u and v be normalized Aoating point numbers. Then ((u G3 VI 0 4 + ((u $ v) 8 ((u $ v) e u)) = u $21, (40) provided that no exponent overAow or underflow occurs.

218 ARITHMETIC 42.2 indicates that a relation of approximately equal would be more useful than the traditional relation of equality, when floating point computations are involved, if we could only define a suitable approximation relation. In other words, the fact that strict equality of floating point values is of little importance implies that we ought to have a new operation, Boating point comparison, which is intended to help assess the relative values of two floating point quantities. The following definitions seem to be appropriate for base b, excess 4, floating point numbers u = (eU, fU) and ZI = (e,, f,,): u < 2) (E) if and only if 2, - u > e max(,x(beu–q, Vu+); (21) u -v (6) if and only if Iv - uI 2 E max(beu-q, bevP4); (22) u t v (f) if and only if u - 2, > E max(beu-q, Vu- ); (23) u M v (E) if and only if Iz, - u[ 5 E min(beu–q, beues). (24 These definitions apply to unnormalized values as well as to normalized ones. Note that exactly one of the conditions u < v (definitely less than), u -v (approximately equal to), or u + w (definitely greater than) must always hold for any given pair of values u and w. The relation u NN 2) is somewhat stronger than u -V, and it might be read u is essentially equal to v. All of the relations are given in terms of a positive real number E that measures the degree of approximation being considered. One way to view the above definitions is to associate a neighborhood set N(u) = { 5 1 1~ - u( 5 MU--q } with each floating point number u; thus, N(u) represents a set of values near u based on the exponent of u s floating point representation. In these terms, we have u + v if and only if N(u) < v and u < N(w); u -v if and only if u E N(w) or v E N(u); u + w if and only if u > N(v) and N(u) > V; u z v if and only if u E N(v) and ZJ E N(u). (Here we are assuming that the parameter 6, which measures the degree of approximation, is a constant; a more complete notation would indicate the dependence of N(u) upon 6.) Here are some simple consequences of the above definitions: if u % 2, (E) then 2, > u (6); (25) if u sz 21 (E) then u -v (e); (26) u M u (e); (27) if u + 21 (E) then u < v; (28) if u -K v (~1) and ~1 2 ~2 then u+v (E2h (29) if u - 21 (El) and ~1 5 ~2 then u-v (c2); (30) if u 25 v (El) and 61 5 ~2 then UZV (E2h (31) if u 4 w (~1) and v -C w (~2) then u+w (cl + e2); (32) if u M w (Ed) and w z w (~2) then u-w (El + 62). (33)