Floating-point arithmetic is considered an esoteric subject by many people. This
ID: 3735189 • Letter: F
Question
Floating-point arithmetic is considered an esoteric subject by many people. This is rather surprising because floating point is ubiquitous in computer systems. Almost every language has a floating-point data type, 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. Only theoretical mathematicians care about floating point accuracy. Do you agree with this? Support your answer with examples. Search for news or online articles discuss the consequences of a lack of floating point accuracy. Also consider the practical applications of floating numbers in your discussions. Support your opinion with at least one reference.
Explanation / Answer
Dear,
Please find the Answer.
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 arithmetic (floating-point hereafter) that have a direct connection to systems building. It consists of three loosely connected parts. The first section, Rounding Error, discusses the implications of using different rounding strategies for the basic operations 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. The discussion of the standard draws on the material in the section Rounding Error. The third part discusses the connections between floating-point and the design of various aspects of computer systems. Topics include instruction set design, optimizing compilers and exception handling.
I have tried to avoid making statements about floating-point without also giving reasons why the statements are true, especially since the justifications involve nothing more complicated than elementary calculus. Those explanations that are not central to the main argument have been grouped into a section called "The Details," so that they can be skipped if desired. In particular, the proofs of many of the theorems appear in this section. The end of each proof is marked with the z symbol. When a proof is not included, the z appears immediately following the statement of the theorem.
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 fit back into its finite representation. This rounding error is the characteristic feature of floating-point computation. The section Relative Error and Ulps describes how it is measured.
Since most floating-point calculations have rounding error anyway, does it matter if the basic arithmetic operations introduce a little bit more rounding error than necessary? That question is a main theme throughout this section. The section Guard Digits 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 Exactly Rounded Operations.
Floating-point Formats
Several different representations of real numbers have been proposed, but by far the most widely used is the floating-point representation.1 Floating-point representations have a base (which is always assumed to be even) and a precision p. If = 10 and p = 3, then the number 0.1 is represented as 1.00 × 10-1. If = 2 and p = 24, then the decimal number 0.1 cannot be represented exactly, but is approximately 1.10011001100110011001101 × 2-4.
In general, a floating-point number will be represented as ± d.dd... d × e, where d.dd... d is called the significand2 and has p digits. More precisely ± d0 . d1 d2 ... dp-1 × e represents the number
(1) .
The term floating-point number will 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, emax and emin. Since there are p possible significands, and emax - emin + 1 possible exponents, a floating-point number can be encoded in
bits, 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 = 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 × or smaller than 1.0 × . Most of this paper discusses issues due to the first reason. However, numbers that are out of range will be discussed in the sections Infinity and Denormalized Numbers.
Floating-point representations are not necessarily unique. For example, both 0.01 × 101 and 1.00 × 10-1 represent 0.1. If the leading digit is nonzero (d0 0 in equation (1) above), then the representation is said to be normalized. The floating-point number 1.00 × 10-1 is normalized, while 0.01 × 101 is not. When = 2, p = 3, emin = -1 and emax = 2 there are 16 normalized floating-point numbers, as shown in FIGURE D-1. The bold hash marks correspond to numbers whose significand 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 represent 0 is with 1.0 × , since this preserves the fact that the numerical ordering of nonnegative real numbers corresponds to the lexicographic ordering of their floating-point representations.3 When the exponent is stored in a k bit field, that means that only 2k - 1 values are available for use as exponents, since one must be reserved to represent 0.
Note that the × in a floating-point number is part of the notation, and different from a floating-point multiply operation. The meaning of the × symbol should be clear from the context. For example, the expression (2.5 × 10-3) × (4.0 × 102) involves only a single floating-point multiplication.
FIGURE D-1 Normalized numbers when = 2, p = 3, emin = -1, emax = 2
Relative Error and Ulps
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 = 10 and p = 3, which will be used throughout this section. If the result of a floating-point computation is 3.12 × 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 × 10-2, then it is in error by .159 units in the last place. In general, if the floating-point number d.d...d × e is used to represent z, then it is in error by d.d...d - (z/e)p-1 units in the last place.4, 5 The term ulps will be used as shorthand for "units in the last place." If the result of a calculation is the floating-point number nearest to the correct result, it still might be in error by as much as .5 ulp. Another way to measure the difference between a floating-point number and the real number it is approximating is relative error, which is simply the difference between the two numbers divided by the real number. For example the relative error committed when approximating 3.14159 by 3.14 × 100 is .00159/3.14159 .0005.
To compute the relative error that corresponds to .5 ulp, observe that when a real number is approximated by the closest possible floating-point number d.dd...dd × e, the error can be as large as 0.00...00' × e, where ' is the digit /2, there are p units in the significand of the floating-point number, and p units of 0 in the significand of the error. This error is ((/2)-p) × e. Since numbers of the form d.dd...dd × e all have the same absolute error, but have values that range between e and × e, the relative error ranges between ((/2)-p) × e/e and ((/2)-p) × e/e+1. That is,
(2)
In particular, the relative error corresponding to .5 ulp can vary by a factor of . This factor is called the wobble. Setting = (/2)-p to the largest of the bounds in (2) above, we can say that when a real number is rounded to the closest floating-point number, the relative error is always bounded by e, which is referred to as machine epsilon.
In the example above, the relative error was .00159/3.14159 .0005. In order to avoid such small numbers, the relative error is normally written as a factor times , which in this case is = (/2)-p = 5(10)-3 = .005. Thus the relative error would be expressed as (.00159/3.14159)/.005) 0.1.
To illustrate the difference between ulps and relative error, consider the real number x = 12.35. It is approximated by = 1.24 × 101. The error is 0.5 ulps, the relative error is 0.8. Next consider the computation 8 . The exact value is 8x = 98.8, while the computed value is 8 = 9.92 × 101. The error is now 4.0 ulps, but the relative error is still 0.8. The error measured in ulps is 8 times larger, even though the relative error is the same. In general, when the base is , a fixed relative error expressed in ulps can wobble by a factor of up to . And conversely, as equation (2) above shows, a fixed error of .5 ulps results in a relative error that can wobble by .
The most natural way to measure rounding error is in ulps. For example rounding to the nearest floating-point number corresponds to an error of less than or equal to .5 ulp. However, when analyzing the rounding error caused by various formulas, relative error is a better measure. A good illustration of this is the analysis in the section Theorem 9. Since can overestimate the effect of rounding to the nearest floating-point number by the wobble factor of , error estimates of formulas will be tighter on machines with a small .
When only the order of magnitude of rounding error is of interest, ulps and 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 that the number of contaminated digits is log n. If the relative error in a computation is n, then
(3) contaminated digits log n.
Guard Digits
One method of computing the difference between two floating-point numbers is to compute the difference exactly and 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 × 1012 - 1.25 × 10-5 would be calculated as
x = 2.15 × 1012
y = .0000000000000000125 × 1012
x - y = 2.1499999999999999875 × 1012
which rounds to 2.15 × 1012. Rather than using all these digits, floating-point hardware normally operates on a fixed number of digits. Suppose that 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 × 1012 - 1.25 × 10-5 becomes
x = 2.15 × 1012
y = 0.00 × 1012
x - y = 2.15 × 1012
The answer is exactly the same as if the difference had been computed exactly and then rounded. Take another example: 10.1 - 9.93. This becomes
x = 1.01 × 101
y = 0.99 × 101
x - y = .02 × 101
The correct answer is .17, so the computed difference is off by 30 ulps and is wrong in every digit! How bad can the error be?
Theorem 1
Using a floating-point format with parameters and p, and computing differences using p digits, the relative error of the result can be as large as - 1.
Proof
A relative error of - 1 in the expression x - y occurs when x = 1.00...0 and y = ...., where = - 1. Here y has p digits (all equal to ). The exact difference is x - y = -p. However, when computing the answer using only p digits, the rightmost digit of y gets shifted off, and so the computed difference is -p+1. Thus the error is -p - -p+1 = -p ( - 1), and the relative error is -p( - 1)/-p = - 1. z
When =2, the relative error can be as large as the result, and when =10, it can be 9 times larger. Or to put it another way, when =2, equation (3) shows that the number of contaminated digits is log2(1/) = log2(2p) = p. That is, all of the p digits in the result are wrong! Suppose that one extra digit is added to guard against this situation (a guard digit). That is, the smaller number is truncated to p + 1 digits, and then the result of the subtraction is rounded to p digits. With a guard digit, the previous example becomes
x = 1.010 × 101
y = 0.993 × 101
x - y = .017 × 101
and the answer is exact. With a single guard digit, the relative error of the result may be greater than , as in 110 - 8.59.
x = 1.10 × 102
y = .085 × 102
x - y = 1.015 × 102
This rounds to 102, compared with the correct answer of 101.41, for a relative error of .006, which is greater than = .005. In general, the relative error of the result can be only slightly larger than . More precisely,
Theorem 2
If x and y are floating-point numbers in a format with parameters 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 Rounding Error. Addition is included in the above theorem since x and y can be positive or negative.
Cancellation
The last section 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 that all the digits are meaningless (Theorem 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 formula, the expression b2 - 4ac occurs. The quantities b2 and 4ac are subject to rounding errors since they are the results of floating-point multiplications. Suppose that they are rounded to the nearest floating-point number, and so are accurate to within .5 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 example, consider b = 3.34, a = 1.22, and c = 2.28. The exact value of b2 - 4ac is .0292. But b2 rounds to 11.2 and 4ac 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 to .16. 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).
A formula that exhibits catastrophic cancellation can sometimes be rearranged to eliminate the problem. Again consider the quadratic formula
(4)
When , then does not involve a cancellation and
.
But the other addition (subtraction) in one of the formulas will have a catastrophic cancellation. To avoid this, multiply the numerator and denominator of r1 by
(and similarly for r2) to obtain
(5)
If and , then computing r1 using formula (4) will involve a cancellation. Therefore, use formula (5) for computing r1 and (4) for r2. On the other hand, if b < 0, use (4) for computing r1 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).7 Unlike the quadratic 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. 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 the section Rounding Error).
In order 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 to x2 - y2, the floating-point numbers x and y might themselves be approximations to some true quantities and . For example, and 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 , and so the advantage of (x + y)(x - y) over x2 - y2 is not as dramatic. Since computing (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 eliminating 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 exact and that the results are computed as accurately as possible.
The expression x2 - y2 is more accurate when rewritten as (x - y)(x + y) because a catastrophic cancellation is 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
(6)
(Suppose the triangle is very flat; that is, a b + c. Then s a, and the term (s - a) in formula (6) subtracts two nearby numbers, one of which may have rounding error. For example, if a = 9.0, b = c = 4.53, the correct value of s is 9.03 and A is 2.342.... 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 70 ulps.
There is a way to rewrite formula (6) so that it will return accurate results even for flat triangles [Kahan 1986]. It is
(7)
If a, b, and c do not satisfy a b c, 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 (7) to compute the area of a triangle is at most 11, provided that subtraction is performed with a guard digit, e .005, and that square roots are computed to within 1/2 ulp.
The condition that e < .005 is met in virtually every actual floating-point system. For example when = 2, p 8 ensures that e < .005, and when = 10, p 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
(8) SQRT((a (b c)) (c (a b)) (c (a b)) (a (b c))) 4
Because of the cumbersome nature of (8), in the statement of theorems we will usually say the computed value of E 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 0.7, which is much less than 11. 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)n, where . This expression arises in financial calculations. Consider depositing $100 every day into a bank account that earns an annual interest rate of 6%, compounded daily. If n = 365 and i = .06, the amount of money accumulated at the end of one year is
100
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 expression 1 + i/n involves adding 1 to .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.
The troublesome expression (1 + i/n)n can be rewritten as enln(1 + i/n), where now the problem is to compute ln(1 + x) for small x. One approach is to use the approximation ln(1 + 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(1 + x) very accurately, as Theorem 4 shows [Hewlett-Packard 1982]. This formula yields $37614.07, accurate to within two cents!
Theorem 4 assumes that LN(x) approximates ln(x) to within 1/2 ulp. The problem it solves is that when x is small, LN(1 x) is not close to ln(1 + x) because 1 x has lost the information in the low order bits of x. That is, the computed value of ln(1 + x) is not close to its actual value when .
Theorem 4
If ln(1 + x) is computed using the formula
the relative error is at most 5 when 0 x < 3/4, provided subtraction is performed with a guard digit, e < 0.1, and ln is computed to within 1/2 ulp.
This formula will work for any value of x but is only interesting for , which is where catastrophic cancellation occurs in the naive formula ln(1 + x). Although the formula may seem mysterious, there is a simple explanation for why it works. Write ln(1 + x) as
.
The left hand factor can be computed exactly, but the right hand factor µ(x) = ln(1 + x)/x will suffer a large rounding error when adding 1 to x. However, µ is almost constant, since ln(1 + x) x. So changing x slightly will not introduce much error. In other words, if , computing will be a good approximation to xµ(x) = ln(1 + x). Is there a value for for which and can be computed accurately? There is; namely = (1 x) 1, because then 1 + 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 numerical 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 it merely requires making the adder one 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 ln(1 + x). Although most modern computers have a guard digit, there are a few (such as Cray systems) that do not.
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.8 The example immediately preceding Theorem 2 shows that a single guard digit will not always give exactly rounded results. The previous section 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 round to 12 or 13? One school 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 Equipment Corporation's VAX computers. 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 obtaining this 50% behavior 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 x0 = x, x1 = (x0 y) y, ..., xn = (xn-1 y) y. If and are exactly rounded using round to even, then either xn = x for all n or xn = x1 for all n 1. z
To clarify this result, consider = 10, p = 3 and let x = 1.00, y = -.555. When rounding up, the sequence becomes
x0 y = 1.56, x1 = 1.56 .555 = 1.01, x1 y = 1.01 .555 = 1.57,
and each successive value of xn increases by .01, until xn = 9.45 (n 845)9. Under round to even, xn 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 floating-point numbers using a very large significand, 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 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 = xh + xl and y = yh + yl, the exact product is
xy = xh yh + xh yl + xl yh + xl yl.
If x and y have p bit significands, the summands will also have p bit significands provided that xl, xh, yh, yl can be represented using [p/2] bits. When p is even, it is easy to find a splitting. The number x0.x1 ... xp - 1 can be written as the sum of x0.x1 ... xp/2 - 1 and 0.0 ... 0xp/2 ... xp - 1. 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 xh = .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 6
Let p be the floating-point precision, with the restriction that p is even when > 2, and assume that floating-point operations are exactly rounded. Then if k = [p/2] is half the precision (rounded up) and m = k + 1, x can be split as x = xh + xl, where
xh = (m x) (m x x), xl = x xh,
and each xi is representable using [p/2] bits of precision.
To see how this theorem works in an example, let = 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 ulps. Using Theorem 6 to write b = 3.5 - .024, a = 3.5 - .037, and c = 3.5 - .021, b2 becomes 3.52 - 2 × 3.5 × .024 + .0242. Each summand is exact, so b2 = 12.25 - .168 + .000576, where the sum is left unevaluated at this point. Similarly, ac = 3.52 - (3.5 × .037 + 3.5 × .021) + .037 × .021 = 12.25 - .2030 +.000777. Finally, subtracting these two series term by term gives an estimate for b2 - ac of 0 .0350 .000201 = .03480, which is identical to the exactly rounded result. To show that Theorem 6 really requires exact rounding, consider p = 3, = 2, and x = 7. Then m = 5, mx = 35, and m x = 32. If subtraction is performed with a single guard digit, then (m x) x = 28. Therefore, xh = 4 and xl = 3, hence xl is not representable with [p/2] = 1 bit.
As a final example of exact rounding, consider dividing m by 10. The result is a floating-point number that will in general not be equal to m/10. When = 2, multiplying m/10 by 10 will 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 The IEEE Standard.
Theorem 7
When = 2, if m and n are integers with |m| < 2p - 1 and n has the special form n = 2i + 2j, then (m n) n = m, provided floating-point operations are exactly rounded.
Proof
Scaling by a power of two is harmless, since it changes only the exponent, not the significand. If q = m/n, then scale n so that 2p - 1 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 one bit to the right of the binary point. Changing the sign of m is harmless, so assume that q > 0.
If = m n, to prove the theorem requires showing that
(9)
That is because m has at most 1 bit right of the binary point, so n will round to m. To deal with the halfway case when |n - m| = 1/4, note that since the initial unscaled m had |m| < 2p - 1, its low-order bit was 0, so the low-order bit of the scaled m is also 0. Thus, halfway cases will round to m.
Suppose that q = .q1q2 ..., and let = .q1q2 ... qp1. To estimate |n - m|, first compute
| - q| = |N/2p + 1 - m/n|,
where N is an odd integer. Since n = 2i + 2j and 2p - 1 n < 2p, it must be that n = 2p - 1 + 2k for some k p - 2, and thus
.
The numerator is an integer, and since N is odd, it is in fact an odd integer. Thus,
| - q| 1/(n2p + 1 - k).
Assume q < (the case q > is similar).10 Then n < m, and
|m-n |= m-n = n(q- ) = n(q-( -2-p-1))
=(2p-1+2k)2-p-1-2-p-1+k =
This establishes (9) and proves the theorem.11 z
The theorem holds true for any base , as long as 2i + 2j is replaced by i + j. As gets larger, however, denominators of the form i + j are farther and farther apart.
We are now in a position to answer the question, Does it matter if the basic 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. The section Cancellation 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 floating-point variable is only an approximation to some actual value.
Thank you.
Related Questions
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.