contaminated digits log n.

D.3.3 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 × 10¹² - 1.25 × 10^-5 would be calculated as

x = 2.15 × 10¹²
y = .0000000000000000125 × 10¹²
x - y = 2.1499999999999999875 × 10¹²

which rounds to 2.15 × 10¹². 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 × 10¹² - 1.25 × 10^-5 becomes

x = 2.15 × 10¹²
y = 0.00 × 10¹²
x - y = 2.15 × 10¹²

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 × 10¹
y = 0.99 × 10¹
x - y = .02 × 10¹

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?

D.3.3.1 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.

D.3.3.2 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.

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 log₂(1/) = log₂(2^p) = 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 × 10¹
y = 0.993 × 10¹
x - y = .017 × 10¹

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 × 10²
y = .085 × 10²
x - y = 1.015 × 10²

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,

D.3.3.3 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.

D.3.4 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 b² - 4ac occurs. The quantities b² 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 b² - 4ac is .0292. But b² 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 .1^[6]. 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

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 r₁ by

(and similarly for r₂) to obtain

If and b> 0, then computing r₁using formula (4)will involve a cancellation. Therefore, use formula (5)for computing r₁and (4)for r₂. On the other hand, if b< 0, use (4)for computing r₁and (5)for r₂.

The expression x² - y² 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 - ydenotes the exact difference of x and y, x ydenotes 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 (xy) (xy) is an excellent approximation to x²- y², 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 yis 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 x²- y²is not as dramatic. Since computing (x+y)(x-y) is about the same amount of work as computing x²- y², 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 x² - y² 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

(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

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.

D.3.4.1 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

SQRT((a (bc)) (c(ab)) (c (a b)) (a (bc))) 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)ⁿ, 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)ⁿ can be rewritten as e^{nln(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 xhas 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 .

D.3.4.2 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)/xwill 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.

D.3.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.^[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.

D.3.5.1 Theorem 5

Let x and y be floating-point numbers, and define x₀= x, x₁= (x₀y) y, …, x_n= (x_n-1y) y. If and are exactly rounded using round to even, then either x_n= x for all n or x_n= x₁for all n 1.   

To clarify this result, consider = 10, p = 3 and let x = 1.00, y = -.555. When rounding up, the sequence becomes

x₀y = 1.56, x₁= 1.56 .555 = 1.01, x₁y = 1.01 .555 = 1.57,  

and each successive value of x_n increases by .01, until x_n = 9.45 (n 845)^[9]. Under round to even, x_n 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 = x_h + x_l and y = y_h + y_l, the exact product is

xy = x_h y_h + x_h y_l + x_l y_h + x_l y_l.

If x and y have p bit significands, the summands will also have p bit significands provided that x_l, x_h, y_h, y_l can be represented using [p/2] bits. When p is even, it is easy to find a splitting. The number x₀.x₁ … x_{p - 1} can be written as the sum of x₀.x₁ … x_{p/2 - 1} and 0.0 … 0x_p/2 … x_{p - 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 x_h = .11 and x_l = -.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.

D.3.5.2 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 = x_h + x_l, where

x_h= (m x) (m x x), x_l= x x_h,   

and each x_i 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 b²- acrounded to the nearest floating-point number is .03480, while bb= 12.08, ac= 12.05, and so the computed value of b²- acis .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, b²becomes 3.5²- 2 ×3.5 ×.024 + .024². Each summand is exact, so b²= 12.25 - .168 + .000576, where the sum is left unevaluated at this point. Similarly, ac= 3.5²- (3.5 ×.037 + 3.5 ×.021) + .037 ×.021 = 12.25 - .2030 +.000777. Finally, subtracting these two series term by term gives an estimate for b²- acof 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 mx= 32. If subtraction is performed with a single guard digit, then (mx) x= 28. Therefore, x_h= 4 and x_l= 3, hence x_lis 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.

D.3.5.3 Theorem 7

When = 2, if mand nare integers with |m| < 2^{p - 1}and nhas the special form n = 2ⁱ+ 2^j, then (m n) n = m, provided floating-point operations are exactly rounded.

D.3.5.4 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 2^{p - 1} n < 2^p and scale m so that 1/2 < q < 1. Thus, 2^{p - 2} < m < 2^p. 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  

That is because mhas at most 1 bit right of the binary point, so nwill round to m. To deal with the halfway case when |n- m| = 1/4, note that since the initial unscaled mhad |m| < 2^{p - 1}, its low-order bit was 0, so the low-order bit of the scaled mis also 0. Thus, halfway cases will round to m.  

Suppose that q= .q₁q₂…, and let = .q₁q₂…q_p1. To estimate |n- m|, first compute 

|- q| = |N/2^{p + 1}- m/n|,

where N is an odd integer. Since n = 2ⁱ + 2^j and 2^{p - 1} n < 2^p, it must be that n = 2^{p - 1} + 2^k 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/(n2^{p + 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))    

=(2^p-1+2^k)2^-p-1-2^-p-1+k=

This establishes (9) and proves the theorem.^[11]

The theorem holds true for any base , as long as 2ⁱ + 2^j is replaced by ⁱ + ^j. As gets larger, however, denominators of the form ⁱ + ^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.

---

D.4 The IEEE Standard

There are two different IEEE standards for floating-point computation. IEEE 754 is a binary standard that requires = 2, p = 24 for single precision and p = 53 for double precision [IEEE 1987]. It also specifies the precise layout of bits in a single and double precision. IEEE 854 allows either = 2 or = 10 and unlike 754, does not specify how floating-point numbers are encoded into bits [Cody et al. 1984]. 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 themselves [IEEE 1987; Cody et al. 1984].

D.4.1 Formats and Operations

D.4.1.1 Base

It is clear why IEEE 854 allows = 10. Base ten is how humans exchange and think about numbers. Using = 10 is especially appropriate for calculators, where the result of each operation is displayed by the calculator in decimal.

There are several reasons why IEEE 854 requires that if the base is not 10, it must be 2. The section Relative Error and Ulpsmentioned one reason: the results of error analyses are much tighter when is 2 because a rounding error of .5 ulp wobbles by a factor of 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 = 16, p = 1 compared to = 2, p = 4. Both systems have 4 bits of significand. Consider the computation of 15/8. When = 2, 15 is represented as 1.111 ×2³, and 15/8 as 1.111 ×2⁰. So 15/8 is exact. However, when = 16, 15 is represented as F×16⁰, where Fis the hexadecimal digit for 15. But 15/8 is represented as 1 ×16⁰, which has only one bit correct. In general, base 16 can lose up to 3 bits, so that a precision of p hexadecimal digits can have an effective precision as low as 4p- 3 rather than 4pbinary bits. Since large values of have these problems, why did IBM choose = 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 significand requires 24 bits. Since this must fit into 32 bits, this leaves 7 bits for the exponent and one for the sign bit. Thus the magnitude of representable numbers ranges from about to about = . To get a similar exponent range when = 2 would require 9 bits of exponent, leaving only 22 bits for the significand. However, it was just pointed out that when = 16, the effective precision can be as low as 4p- 3 = 21 bits. Even worse, when = 2 it is possible to gain an extra bit of 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 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 = 16, p = 1 system, all the numbers between 1 and 15 have the same exponent, and so no shifting is required when adding any of the () = 105 possible pairs of distinct numbers from this set. However, in the = 2, p = 4 system, these numbers have exponents ranging from 0 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, and so the small wobble of = 2 makes it the preferable base. Another advantage of using = 2 is that there is a way to gain an extra bit of significance.^[12]Since floating-point numbers are always normalized, the most significant bit of the significand is always 1, and there is no reason to waste a bit of storage representing it. Formats that use this trick are said to have a hiddenbit. It was already pointed out in Floating-point Formatsthat this requires a special convention for 0. The method given there was that an exponent of e_min- 1 and a significand of all zeros represents not , but rather 0.

IEEE 754 single precision is encoded in 32 bits using 1 bit for the sign, 8 bits for the exponent, and 23 bits for the significand. However, it uses a hidden bit, so the significand is 24 bits (p = 24), even though it is encoded using only 23 bits.

D.4.1.2 Precision

The IEEE standard defines four different precisions: single, double, single-extended, and double-extended. In IEEE 754, single and double precision correspond roughly to what most floating-point hardware provides. Single precision occupies a single 32 bit word, double precision two consecutive 32 bit words. Extended precision is a format that offers at least a little extra precision and exponent range (TABLE D-1).

TABLE D-1 IEEE 754 Format Parameters

Parameter	Format
	Single	Single-Extended	Double	Double-Extended
p	24	32	53	64
emax	+127	1023	+1023	> 16383
emin	-126	-1022	-1022	-16382
Exponent width in bits	8	11	11	15
Format width in bits	32	43	64	79

The IEEE standard only specifies a lower bound on how many extra bits extended precision provides. The minimum allowable double-extended format is sometimes referred to as 80-bit format, even though the table shows it using 79 bits. The reason is that hardware implementations of extended precision normally do not use a hidden bit, and so would use 80 rather than 79 bits.^[13]

The standard puts the most emphasis on extended precision, making no recommendation concerning double precision, but strongly recommending that Implementations should support the extended format corresponding to the widest basic format supported, …

One motivation for extended precision comes from calculators, which will often display 10 digits, but use 13 digits internally. By displaying only 10 of the 13 digits, the calculator appears to the user as a "black box" that computes exponentials, cosines, etc. to 10 digits of accuracy. For the calculator to compute functions like exp, log and cos to within 10 digits with reasonable efficiency, it needs a few extra digits to work with. It is not hard to find a simple rational expression that approximates log with an error of 500 units in the last place. Thus computing with 13 digits gives an answer correct to 10 digits. By keeping these extra 3 digits hidden, the calculator presents a simple model to the operator.

Extended precision in the IEEE standard serves a similar function. It enables libraries to efficiently compute quantities to within about .5 ulp in single (or double) precision, giving the user of those libraries a simple model, namely that each primitive operation, be it a simple multiply or an invocation of log, returns a value accurate to within about .5 ulp. However, when using extended precision, it is important to make sure that its use is transparent to the user. For example, on a calculator, if the internal representation of a displayed value is not rounded to the same precision as the display, then the result of further operations will depend on the hidden digits and appear unpredictable to the user.

To illustrate extended precision further, consider the problem of converting between IEEE 754 single precision and decimal. Ideally, single precision numbers will be printed with enough digits so that when the decimal number is read back in, the single precision number can be recovered. It turns out that 9 decimal digits are enough to recover a single precision binary number (see the section Binary to Decimal Conversion). When converting a decimal number back to its unique binary representation, a rounding error as small as 1 ulp is fatal, because it will give the wrong answer. Here is a situation where extended precision is vital for an efficient algorithm. When single-extended is available, a very straightforward method exists for converting a decimal number to a single precision binary one. First read in the 9 decimal digits as an integer N, ignoring the decimal point. From TABLE D-1, p 32, and since 10⁹ < 2³² 4.3 × 10⁹, N can be represented exactly in single-extended. Next find the appropriate power 10^P necessary to scale N. This will be a combination of the exponent of the decimal number, together with the position of the (up until now) ignored decimal point. Compute 10^|P|. If |P| 13, then this is also represented exactly, because 10¹³ = 2¹³5¹³, and 5¹³ < 2³². Finally multiply (or divide if p < 0) N and 10^|P|. If this last operation is done exactly, then the closest binary number is recovered. The section Binary to Decimal Conversion shows how to do the last multiply (or divide) exactly. Thus for |P| 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 |P| > 13, then single-extended is not enough for the above algorithm to always compute the exactly rounded binary equivalent, 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, then the algorithm above would be 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.

D.4.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 significand in the IEEE formats: one 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 in the range [-2^p-1, 2^p-1 - 1] is represented by the smallest nonnegative number that is congruent to it modulo 2^p.

The IEEE binary standard does not use either of these methods to represent the exponent, but instead uses a biasedrepresentation. In the case of single precision, where the exponent is stored in 8 bits, the bias is 127 (for double precision it is 1023). What this means is that if 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 unbiased exponentto distinguish from the biased exponent .   

Referring to TABLE D-1, single precision has e_max= 127 and e_min= -126. The reason for having |e_min| < e_maxis so that the reciprocal of the smallest number will not overflow. Although it is true that the reciprocal of the largest number will underflow, underflow is usually less serious than overflow. The section Baseexplained that e_min- 1 is used for representing 0, and Special Quantitieswill introduce a use for e_max+ 1. In IEEE single precision, this means that the biased exponents range between e_min- 1 = -127 and e_max+ 1 = 128, whereas the unbiased exponents range between 0 and 255, which are exactly the nonnegative numbers that can be represented using 8 bits.

D.4.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 and then rounded to the nearest floating-point number (using round to even). The section Guard Digits pointed out that computing the exact difference or sum of two floating-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. However, computing with a single guard digit will not always give the same answer as computing the exact result and 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 and then rounded [Goldberg 1990]. 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 machines and both support IEEE arithmetic, then if any intermediate result differs, it must be because of software bugs, not from 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. However, proofs in this system cannot verify the algorithms of sections Cancellation and Exactly Rounded Operations, which require features not present on all hardware. Furthermore, Brown's axioms are more complex than simply defining operations to be performed exactly and then rounded. Thus proving theorems from Brown's axioms is usually more difficult 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 +, -, × 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 [1986] 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 × 10^-30 + 10³⁰) - 10³⁰) - 10^-30 is exactly equal to 10^-30, but on a machine with IEEE arithmetic 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 Kulish 1987].^{[14] [15]}

All the operations mentioned in the standard are required to be exactly rounded except conversion between decimal and binary. The reason is that efficient algorithms for exactly rounding all the operations are known, except conversion. For conversion, the best known efficient algorithms produce results that are slightly worse than exactly rounded ones [Coonen 1984].

The IEEE standard does not require transcendental functions to be exactly rounded because of the table maker's dilemma. To illustrate, suppose you are making a table of the exponential function to 4 places. Then exp(1.626) = 5.0835. Should this be rounded to 5.083 or 5.084? If exp(1.626) is computed more carefully, it becomes 5.08350. And then 5.083500. And then 5.0835000. Since exp is transcendental, this could go on arbitrarily long before distinguishing whether exp(1.626) is 5.083500…0ddd or 5.0834999…9ddd. Thus it is not practical to specify that the precision of transcendental functions be the same as if they were computed to infinite precision and 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,^[16] and large tables are three different techniques that are used for computing transcendentals 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.

D.4.2 Special Quantities

On some floating-point hardware every bit pattern represents a valid floating-point number. The IBM System/370 is an example of this. On the other hand, the VAX reserves some bit patterns to represent special numbers called reserved operands. This idea goes back to 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) and infinities. Without any 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 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, is returned. In IEEE arithmetic, a NaN is returned in this situation.

The IEEE standard specifies the following special values (see TABLE D-2): ± 0, denormalized numbers, ± and NaNs (there is more than one NaN, as explained in the next section). These special values are all encoded with exponents of either e_max + 1 or e_min - 1 (it was already pointed out that 0 has an exponent of e_min - 1).

TABLE D-2 IEEE 754 Special Values

Exponent	Fraction	Represents
e = emin - 1	f = 0	±0
e = emin - 1	f 0
emin e emax	--	1.f × 2e
e = emax + 1	f = 0	±
e = emax + 1	f 0	NaN

D.4.3 NaNs

Traditionally, the computation of 0/0 or has been treated as an unrecoverable error which causes a computation to halt. However, th