Trăn trừ nhị phân

Số học dấu phẩy động được coi là một chủ đề bí truyền của nhiều người. Điều này khá ngạc nhiên vì dấu chấm động phổ biến trong các hệ thống máy tính. Hầu hết mọi ngôn ngữ đều có kiểu dữ liệu dấu chấm động; . Bài báo này trình bày một hướng dẫn về các khía cạnh của dấu phẩy động có tác động trực tiếp đến các nhà thiết kế hệ thống máy tính. Nó bắt đầu với thông tin cơ bản về biểu diễn dấu phẩy động và lỗi làm tròn, tiếp tục với phần thảo luận về tiêu chuẩn dấu phẩy động IEEE và kết thúc bằng nhiều ví dụ về cách các nhà chế tạo máy tính có thể hỗ trợ dấu phẩy động tốt hơn.

Danh mục và mô tả chủ đề. (Chính) C. 0 [Computer Systems Organization]. General -- instruction set design; D. 3. 4 [Programming Languages]. Processors -- compilers, optimization; G. 1. 0 [Numerical Analysis]. General -- computer arithmetic, error analysis, numerical algorithms (Secondary)

D. 2. 1 [Software Engineering]. Requirements/Specifications -- languages; D. 3. 4 Programming Languages]. Formal Definitions and Theory -- semantics; D. 4. 1 Operating Systems]. Quản lý quy trình -- đồng bộ hóa

General Terms. Algorithms, Design, Languages

Additional Key Words and Phrases. Denormalized number, exception, floating-point, floating-point standard, gradual underflow, guard digit, NaN, overflow, relative error, rounding error, rounding mode, ulp, underflow

Introduction

Builders of computer systems often need information about floating-point arithmetic. Tuy nhiên, có rất ít nguồn thông tin chi tiết về nó. Một trong số ít cuốn sách về chủ đề này, Tính toán dấu phẩy động của Pat Sterbenz, đã hết bản in từ lâu. Bài viết này là một hướng dẫn về các khía cạnh của số học dấu phẩy động (dấu phẩy động sau đây) có kết nối trực tiếp với việc xây dựng hệ thống. Nó bao gồm ba phần kết nối lỏng lẻo. Phần đầu tiên, Lỗi làm tròn, thảo luận về ý nghĩa của việc sử dụng các chiến lược làm tròn khác nhau cho các phép tính cơ bản cộng, trừ, nhân và chia. It also contains background information on the two methods of measuring rounding error, ulps and

     n = n/2
2
     n = n/2
3. 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

Trăn trừ nhị phân
(which is always assumed to be even) and a precision p. If
Trăn trừ nhị phân
= 10 and p = 3, then the number 0. 1 is represented as 1. 00 × 10-1. If
Trăn trừ nhị phân
= 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 ×

Trăn trừ nhị phân
e, where d. dd. d is called the significand2 and has p digits. More precisely ± d0 . d1 d2 . dp-1 ×
Trăn trừ nhị phân
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

Trăn trừ nhị phân
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

Trăn trừ nhị phân
= 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
Trăn trừ nhị phân
×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

Trăn trừ nhị phân
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
Trăn trừ nhị phân
 = 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

Trăn trừ nhị phân

FIGURE D-1 Normalized numbers when
Trăn trừ nhị phân
= 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

Trăn trừ nhị phân
 = 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 ×
Trăn trừ nhị phân
e is used to represent z, then it is in error by
Trăn trừ nhị phân
d. d. d - (z/
Trăn trừ nhị phân
e)
Trăn trừ nhị phân
Trăn trừ nhị phân
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 là. 00159/3. 14159 
Trăn trừ nhị phân
. 0005.

Để tính toán lỗi tương đối tương ứng với. 5 ulp, observe that when a real number is approximated by the closest possible floating-point number d. dd. dd ×

Trăn trừ nhị phân
e, the error can be as large as 0. 00. 00
Trăn trừ nhị phân
' ×
Trăn trừ nhị phân
e, where
Trăn trừ nhị phân
' is the digit
Trăn trừ nhị phân
/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 ((
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
-p) ×
Trăn trừ nhị phân
e. Since numbers of the form d. dd. dd ×
Trăn trừ nhị phân
e all have the same absolute error, but have values that range between
Trăn trừ nhị phân
e and
Trăn trừ nhị phân
×
Trăn trừ nhị phân
e, the relative error ranges between ((
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
-p) ×
Trăn trừ nhị phân
e/
Trăn trừ nhị phân
e and ((
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
-p) ×
Trăn trừ nhị phân
e/
Trăn trừ nhị phân
e+1. That is,

(2)

Đặc biệt, lỗi tương đối tương ứng với. 5 ulp có thể thay đổi theo hệ số

Trăn trừ nhị phân
. Yếu tố này được gọi là dao động. Đặt
Trăn trừ nhị phân
= (
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
-p thành giới hạn lớn nhất trong (2) ở trên, chúng ta có thể nói .

Trong ví dụ trên, lỗi tương đối là. 00159/3. 14159

Trăn trừ nhị phân
. 0005. Để tránh những con số nhỏ như vậy, sai số tương đối thường được viết dưới dạng thừa số lần
Trăn trừ nhị phân
, trong trường hợp này là
Trăn trừ nhị phân
= (
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
-p = 5(10)-3 = .005. Thus the relative error would be expressed as (.00159/3.14159)/.005)
Trăn trừ nhị phân
Trăn trừ nhị phân
0. 1
Trăn trừ nhị phân
.

Để minh họa sự khác biệt giữa ulps và lỗi tương đối, hãy xem xét số thực x = 12. 35. Nó được xấp xỉ bởi = 1. 24 × 101. lỗi là 0. 5 ulps, lỗi tương đối là 0. 8

Trăn trừ nhị phân
. Tiếp theo xét phép tính 8. Giá trị chính xác là 8x = 98. 8, trong khi giá trị tính toán là 8= 9. 92 × 101. The error is now 4. 0 ulps, but the relative error is still 0. 8
Trăn trừ nhị phân
. The error measured in ulps is 8 times larger, even though the relative error is the same. In general, when the base is
Trăn trừ nhị phân
, a fixed relative error expressed in ulps can wobble by a factor of up to
Trăn trừ nhị phân
. And conversely, as equation (2) above shows, a fixed error of . 5 ulps results in a relative error that can wobble by
Trăn trừ nhị phân
.

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. Vì

Trăn trừ nhị phân
có thể đánh giá quá cao ảnh hưởng của việc làm tròn đến số dấu phẩy động gần nhất theo hệ số dao động của
Trăn trừ nhị phân
, ước tính lỗi của công thức sẽ chặt chẽ hơn trên các máy có kích thước nhỏ .
Trăn trừ nhị phân
.

When only the order of magnitude of rounding error is of interest, ulps and

Trăn trừ nhị phân
may be used interchangeably, since they differ by at most a factor of
Trăn trừ nhị phân
. For example, when a floating-point number is in error by n ulps, that means that the number of contaminated digits is log
Trăn trừ nhị phân
n. If the relative error in a computation is n
Trăn trừ nhị phân
, then

(3) chữ số bị ô nhiễm
Trăn trừ nhị phân
nhật ký
Trăn trừ nhị phân
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
Trăn trừ nhị phân
and p, and computing differences using p digits, the relative error of the result can be as large as
Trăn trừ nhị phân
- 1.

Proof

A relative error of
Trăn trừ nhị phân
- 1 in the expression x - y occurs when x = 1. 00. 0 and y = .
Trăn trừ nhị phân
Trăn trừ nhị phân
.
Trăn trừ nhị phân
, where
Trăn trừ nhị phân
=
Trăn trừ nhị phân
- 1. Here y has p digits (all equal to
Trăn trừ nhị phân
). The exact difference is x - y =
Trăn trừ nhị phân
-p. However, when computing the answer using only p digits, the rightmost digit of y gets shifted off, and so the computed difference is
Trăn trừ nhị phân
-p+1. Do đó, lỗi là
Trăn trừ nhị phân
-p -
Trăn trừ nhị phân
-p+1 =
Trăn trừ nhị phân
-p (
Trăn trừ nhị phân
- 1), and the relative error is
Trăn trừ nhị phân
-p(
Trăn trừ nhị phân
- 1)/
Trăn trừ nhị phân
-p =
Trăn trừ nhị phân
- 1. z

When

Trăn trừ nhị phân
=2, the relative error can be as large as the result, and when
Trăn trừ nhị phân
=10, it can be 9 times larger. Or to put it another way, when
Trăn trừ nhị phân
=2, equation (3) shows that the number of contaminated digits is log2(1/
Trăn trừ nhị phân
) = 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

và câu trả lời là chính xác. Với một chữ số bảo vệ, sai số tương đối của kết quả có thể lớn hơn

Trăn trừ nhị phân
, như trong 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

Trăn trừ nhị phân
= . 005. In general, the relative error of the result can be only slightly larger than
Trăn trừ nhị phân
. More precisely,

Theorem 2

If x and y are floating-point numbers in a format with parameters
Trăn trừ nhị phân
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
Trăn trừ nhị phân
.

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. Giá trị chính xác của b2 - 4ac là. 0292. But b2 rounds to 11. 2 và 4ac làm tròn thành 11. 1, do đó câu trả lời cuối cùng là. 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

Trăn trừ nhị phân
).

A formula that exhibits catastrophic cancellation can sometimes be rearranged to eliminate the problem. Again consider the quadratic formula

(4)
Trăn trừ nhị phân

When, thendoes 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)
Trăn trừ nhị phân

Ifand, 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

Trăn trừ nhị phân
. 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, xy denotes the computed difference (i. e. , with rounding error). Similarly

Trăn trừ nhị phân
,
Trăn trừ nhị phân
, anddenote computed addition, multiplication, and division, respectively. All caps indicate the computed value of a function, as in
     n = n/2
4 or
     n = n/2
5. Lowercase functions and traditional mathematical notation denote their exact values as in ln(x) and.

Although (xy)

Trăn trừ nhị phân
(x
Trăn trừ nhị phân
y) is an excellent approximation to x2 - y2, the floating-point numbers x and y might themselves be approximations to some true quantitiesand. For example,andmight 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

Trăn trừ nhị phân
b + c. Then s
Trăn trừ nhị phân
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)
Trăn trừ nhị phân

If a, b, and c do not satisfy a

Trăn trừ nhị phân
b
Trăn trừ nhị phân
c, rename them before applying (7). It is straightforward to check that the right-hand sides of (6) and (7) are algebraically identical. Sử dụng các giá trị của a, b và c ở trên sẽ cho diện tích tính toán là 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
Trăn trừ nhị phân
, provided that subtraction is performed with a guard digit, e 
Trăn trừ nhị phân
 . 005, and that square roots are computed to within 1/2 ulp.

The condition that e < . 005 được đáp ứng trong hầu hết mọi hệ thống dấu phẩy động thực tế. Ví dụ: khi

Trăn trừ nhị phân
= 2, p
Trăn trừ nhị phân
8 đảm bảo rằng e
Trăn trừ nhị phân
= 10, p 
Trăn trừ nhị phân
 3 là đủ.

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)
     n = n/2
6((a
Trăn trừ nhị phân
(b
Trăn trừ nhị phân
c))
Trăn trừ nhị phân
(c(ab))
Trăn trừ nhị phân
(c
Trăn trừ nhị phân
(ab))
Trăn trừ nhị phân
(a
Trăn trừ nhị phân
(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

Trăn trừ nhị phân
, which is much less than 11
Trăn trừ nhị phân
. 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

Trăn trừ nhị phân
= 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)

Trăn trừ nhị phân
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

     n = n/2
4 approximates ln(x) to within 1/2 ulp. The problem it solves is that when x is small,
     n = n/2
8(1
Trăn trừ nhị phân
x) is not close to ln(1 + x) because 1
Trăn trừ nhị phân
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
Trăn trừ nhị phân

the relative error is at most 5
Trăn trừ nhị phân
when 0
Trăn trừ nhị phân
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

Trăn trừ nhị phân
.

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)

Trăn trừ nhị phân
x. So changing x slightly will not introduce much error. In other words, if, computingwill be a good approximation to xµ(x) = ln(1 + x). Is there a value forfor whichandcan be computed accurately? There is; namely= (1
Trăn trừ nhị phân
x)1, because then 1 +is exactly equal to 1
Trăn trừ nhị phân
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 = (x0y)
Trăn trừ nhị phân
y, . , xn = (xn-1  y)
Trăn trừ nhị phân
y. If
Trăn trừ nhị phân
andare exactly rounded using round to even, then either xn = x for all n or xn = x1 for all n
Trăn trừ nhị phân
1. z

Để làm rõ kết quả này, hãy xem xét

Trăn trừ nhị phân
= 10, p = 3 và đặt x = 1. 00, y = -. 555. When rounding up, the sequence becomes

x0y = 1. 56, x1 = 1. 56. 555 = 1. 01, x1  y = 1. 01
Trăn trừ nhị phân
. 555 = 1. 57,

and each successive value of xn increases by . 01, until xn = 9. 45 (n

Trăn trừ nhị phân
845)9. Dưới tròn đến chẵn, xn luôn là 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

Chìa khóa của phép nhân trong hệ thống này là biểu diễn tích xy dưới dạng tổng, trong đó mỗi phép cộng có cùng độ chính xác với x và y. Điều này có thể được thực hiện bằng cách tách x và y. Viết x = xh + xl và y = yh + yl, tích chính xác là

xy = xh yh + xh yl + xl yh + xl yl

Nếu x và y có p bit có nghĩa, thì các phép cộng cũng sẽ có p bit có nghĩa với điều kiện là xl, xh, yh, yl có thể được biểu diễn bằng cách sử dụng [p/2] bit. Khi p chẵn, dễ dàng tìm được một phép chia. số x0. x1. xp - 1 có thể được viết dưới dạng tổng của x0. x1. xp/2 - 1 và 0. 0. 0xp/2. xp - 1. Khi p là số lẻ, phương pháp tách đơn giản này sẽ không hoạt động. Tuy nhiên, có thể thu được thêm một bit bằng cách sử dụng các số âm. Ví dụ: nếu

Trăn trừ nhị phân
= 2, p = 5 và x =. 10111, x có thể được chia thành xh =. 11 và xl = -. 00001. Có nhiều hơn một cách để tách một số. Dekker [1971] đã đưa ra một phương pháp phân tách dễ tính toán, nhưng nó yêu cầu nhiều hơn một chữ số bảo vệ.

Định lý 6

Cho p là độ chính xác của dấu phẩy động, với hạn chế là p chẵn khi
Trăn trừ nhị phân
 > 2 và giả sử rằng các phép tính của dấu phẩy động được làm tròn chính xác. Sau đó, nếu k = [p/2] bằng một nửa độ chính xác (làm tròn lên) và m =
Trăn trừ nhị phân
k + 1, thì x có thể được chia thành x = xh + xl, trong đó xh = (m Trăn trừ nhị phân x)  (m
Trăn trừ nhị phân
xx), xl = xxh,

và mỗi xi có thể biểu diễn bằng cách sử dụng [p/2] bit chính xác

Để xem định lý này hoạt động như thế nào trong một ví dụ, hãy

Trăn trừ nhị phân
= 10, p = 4, b = 3. 476, một = 3. 463 và c = 3. 479. Sau đó, b2 - ac được làm tròn thành số dấu phẩy động gần nhất là. 03480, trong khi b
Trăn trừ nhị phân
b = 12. 08, a
Trăn trừ nhị phân
c = 12. 05, và do đó, giá trị tính toán của b2 - ac là. 03. Đây là lỗi 480 ulps. Sử dụng Định lý 6 để viết b = 3. 5 -. 024, a = 3. 5 -. 037 và c = 3. 5 -. 021, b2 trở thành 3. 52 - 2 × 3. 5 ×. 024 +. 0242. Mỗi mệnh đề đều chính xác, vì vậy b2 = 12. 25 -. 168 +. 000576, tại thời điểm này tổng vẫn chưa được đánh giá. Tương tự, ac = 3. 52 - (3. 5 ×. 037 + 3. 5 ×. 021) +. 037 ×. 021 = 12. 25 -. 2030 +. 000777. Cuối cùng, trừ hai chuỗi này theo số hạng sẽ đưa ra ước tính cho b2 - ac là 0 
Trăn trừ nhị phân
. 0350. 000201 =. 03480, giống với kết quả được làm tròn chính xác. Để chứng minh rằng Định lý 6 thực sự yêu cầu làm tròn chính xác, hãy xem xét p = 3,
Trăn trừ nhị phân
= 2 và x = 7. Khi đó m = 5, mx = 35 và m 
Trăn trừ nhị phân
 x = 32. Nếu phép trừ được thực hiện với một chữ số bảo vệ, thì (m 
Trăn trừ nhị phân
 x) x = 28. Do đó, xh = 4 và xl = 3, do đó xl không thể biểu diễn với [p/2] = 1 bit.

Ví dụ cuối cùng về làm tròn chính xác, xem xét chia m cho 10. Kết quả là một số dấu phẩy động nói chung sẽ không bằng m/10. Khi

Trăn trừ nhị phân
= 2, nhân m/10 với 10 sẽ khôi phục lại m, với điều kiện là sử dụng phương pháp làm tròn chính xác. Trên thực tế, một sự thật tổng quát hơn (do Kahan) là đúng. Bằng chứng là khéo léo, nhưng độc giả không quan tâm đến những chi tiết như vậy có thể bỏ qua phần Tiêu chuẩn IEEE.

Định lý 7

Khi
Trăn trừ nhị phân
= 2, nếu m và n là các số nguyên với. m. < 2p - 1 và n có dạng đặc biệt n = 2i + 2j, khi đó (mn)
Trăn trừ nhị phân
n = m, miễn là các phép toán dấu phẩy động được làm tròn chính xác.

Proof

Chia tỷ lệ theo lũy thừa hai là vô hại, vì nó chỉ thay đổi số mũ, không thay đổi ý nghĩa. Nếu q = m/n, thì chia tỷ lệ n sao cho 2p - 1
Trăn trừ nhị phân
n < 2p và chia tỷ lệ m sao cho 1/2 < q < 1. Do đó, 2p - 2 < m < 2p. Vì m có p bit có nghĩa nên nó có nhiều nhất một bit ở bên phải điểm nhị phân. Đổi dấu của m là vô hại, vì vậy giả sử rằng q > 0. Nếu = mn, để chứng minh định lý cần chứng minh rằng(9)
Đó là do m có nhiều nhất 1 bit bên phải điểm nhị phân nên n sẽ làm tròn thành m. Để đối phó với trường hợp giữa chừng khi. n-m. = 1/4, lưu ý rằng vì m chưa chia tỷ lệ ban đầu có. m. < 2p - 1, bit bậc thấp của nó là 0, vì vậy bit bậc thấp của m được chia tỷ lệ cũng là 0. Do đó, các trường hợp nửa chừng sẽ làm tròn thành m. Giả sử rằng q =. q1q2. , và để =. q1q2. qp1. Ước tính, ước lượng. n-m. , tính toán đầu tiên. - q. =. N/2p + 1 - m/n. ,
trong đó N là số nguyên lẻ. Vì n = 2i + 2j và 2p - 1
Trăn trừ nhị phân
n < 2p nên n = 2p - 1 + 2k đối với một số k
Trăn trừ nhị phân
p - 2, .
Trăn trừ nhị phân
.
Tử số là một số nguyên, và vì N là số lẻ nên thực tế nó là một số nguyên lẻ. Như vậy,. - q.
Trăn trừ nhị phân
1/(n2p + 1 - k).
Giả sử q
Trăn trừ nhị phân
Trăn trừ nhị phân

=(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

Trăn trừ nhị phân
, as long as 2i + 2j is replaced by
Trăn trừ nhị phân
i +
Trăn trừ nhị phân
j. As
Trăn trừ nhị phân
gets larger, however, denominators of the form
Trăn trừ nhị phân
i +
Trăn trừ nhị phân
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

The IEEE Standard

There are two different IEEE standards for floating-point computation. IEEE 754 is a binary standard that requires

Trăn trừ nhị phân
= 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
Trăn trừ nhị phân
= 2 or
Trăn trừ nhị phân
= 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. Mục đích của bài báo này không phải là để tranh luận rằng tiêu chuẩn IEEE là tiêu chuẩn dấu phẩy động tốt nhất có thể mà là để chấp nhận tiêu chuẩn như đã cho và giới thiệu về việc sử dụng nó. For full details consult the standards themselves [IEEE 1987; Cody et al. 1984]

Formats and Operations

Base

It is clear why IEEE 854 allows

Trăn trừ nhị phân
= 10. Base ten is how humans exchange and think about numbers. Using
Trăn trừ nhị phân
= 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 Ulps mentioned one reason. the results of error analyses are much tighter when

Trăn trừ nhị phân
is 2 because a rounding error of . 5 ulp wobbles by a factor of
Trăn trừ nhị phân
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
Trăn trừ nhị phân
= 16, p = 1 compared to
Trăn trừ nhị phân
= 2, p = 4. Both systems have 4 bits of significand. Consider the computation of 15/8. When
Trăn trừ nhị phân
= 2, 15 is represented as 1. 111 × 23, and 15/8 as 1. 111 × 20. So 15/8 is exact. However, when
Trăn trừ nhị phân
= 16, 15 is represented as F × 160, where F is the hexadecimal digit for 15. But 15/8 is represented as 1 × 160, 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 4p binary bits. Since large values of
Trăn trừ nhị phân
have these problems, why did IBM choose
Trăn trừ nhị phân
= 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
Trăn trừ nhị phân
= 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 aboutto about=. To get a similar exponent range when
Trăn trừ nhị phân
= 2 would require 9 bits of exponent, leaving only 22 bits for the significand. However, it was just pointed out that when
Trăn trừ nhị phân
= 16, the effective precision can be as low as 4p - 3 = 21 bits. Even worse, when
Trăn trừ nhị phân
= 2 it is possible to gain an extra bit of precision (as explained later in this section), so the
Trăn trừ nhị phân
= 2 machine has 23 bits of precision to compare with a range of 21 - 24 bits for the
Trăn trừ nhị phân
= 16 machine.

Another possible explanation for choosing

Trăn trừ nhị phân
= 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
Trăn trừ nhị phân
= 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
Trăn trừ nhị phân
= 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

Trăn trừ nhị phân
= 2 makes it the preferable base. Another advantage of using
Trăn trừ nhị phân
= 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 hidden bit. It was already pointed out in Floating-point Formats that this requires a special convention for 0. The method given there was that an exponent of emin - 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. Tuy nhiên, nó sử dụng một bit ẩn, do đó, ý nghĩa là 24 bit (p = 24), mặc dù nó được mã hóa chỉ bằng 23 bit

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 ParametersSingleSingle-ExtendedDoubleDouble-Extendedp24
Trăn trừ nhị phân
3253
Trăn trừ nhị phân
64emax+127
Trăn trừ nhị phân
1023+1023> 16383emin-126
Trăn trừ nhị phân
-1022-1022
Trăn trừ nhị phân
-16382Exponent width in bits8
Trăn trừ nhị phân
1111
Trăn trừ nhị phân
15Format width in bits32
Trăn trừ nhị phân
4364
Trăn trừ nhị phân
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. Để máy tính tính toán các hàm như exp, log và cos trong phạm vi 10 chữ số với hiệu quả hợp lý, nó cần thêm một vài chữ số để hoạt động với. 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

Độ chính xác mở rộng trong tiêu chuẩn IEEE phục vụ chức năng tương tự. Nó cho phép các thư viện tính toán số lượng một cách hiệu quả trong khoảng. 5 ulp với độ chính xác đơn (hoặc kép), cung cấp cho người dùng các thư viện đó một mô hình đơn giản, cụ thể là mỗi thao tác nguyên thủy, có thể là một phép nhân đơn giản hoặc một lệnh gọi nhật ký, trả về một giá trị chính xác trong khoảng. 5 lần. Tuy nhiên, khi sử dụng độ chính xác mở rộng, điều quan trọng là đảm bảo rằng việc sử dụng nó minh bạch cho người dùng. 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. Lý tưởng nhất là các số có độ chính xác đơn sẽ được in với đủ chữ số để khi đọc lại số thập phân, số có độ chính xác đơn có thể được phục hồi. Hóa ra 9 chữ số thập phân là đủ để khôi phục một số nhị phân chính xác duy nhất (xem phần Chuyển đổi nhị phân sang thập phân). Khi chuyển đổi một số thập phân trở lại biểu diễn nhị phân duy nhất của nó, một lỗi làm tròn nhỏ tới 1 ulp là nghiêm trọng, vì nó sẽ đưa ra câu trả lời sai. 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 

Trăn trừ nhị phân
 32, and since 109 
Trăn trừ nhị phân
4. 3 × 109, N can be represented exactly in single-extended. Next find the appropriate power 10P 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.  
Trăn trừ nhị phân
 13, then this is also represented exactly, because 1013 = 213513, and 513 
Trăn trừ nhị phân
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

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. Dấu hiệu/độ lớn là hệ thống được sử dụng cho dấu hiệu của ý nghĩa và trong các định dạng IEEE. 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 [-2p-1, 2p-1 - 1] is represented by the smallest nonnegative number that is congruent to it modulo 2p

The IEEE binary standard does not use either of these methods to represent the exponent, but instead uses a biased representation. 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 ifis 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 exponent to distinguish from the biased exponent

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

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 + 1030) - 1030) - 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. Điều này có nên được làm tròn thành 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

Special Quantities

Trên một số phần cứng dấu phẩy động, mọi mẫu bit biểu thị một số dấu phẩy động hợp lệ. The IBM System/370 is an example of this. On the other hand, the VAXTM 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

     n = n/2
9 and
     if (n==0) return u
0

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, số không chuẩn hóa, ±

Trăn trừ nhị phân
và NaN (có nhiều hơn một NaN, như được giải thích trong phần tiếp theo). Tất cả các giá trị đặc biệt này đều được mã hóa bằng số mũ của emax + 1 hoặc emin - 1 (đã chỉ ra rằng 0 có số mũ là emin - 1).

TABLE D-2   IEEE 754 Special ValuesExponentFractionRepresentse = emin - 1f = 0±0e = emin - 1f
Trăn trừ nhị phân
0emin
Trăn trừ nhị phân
e
Trăn trừ nhị phân
emax--1. f × 2ee = emax + 1f = 0±
Trăn trừ nhị phân
e = emax + 1f
Trăn trừ nhị phân
0NaN

NaNs

Traditionally, the computation of 0/0 orhas been treated as an unrecoverable error which causes a computation to halt. However, there are examples where it makes sense for a computation to continue in such a situation. Consider a subroutine that finds the zeros of a function f, say

     if (n==0) return u
1. Traditionally, zero finders require the user to input an interval [a, b] on which the function is defined and over which the zero finder will search. That is, the subroutine is called as
     if (n==0) return u
2,
     if (n==0) return u
3,
     if (n==0) return u
4. A more useful zero finder would not require the user to input this extra information. This more general zero finder is especially appropriate for calculators, where it is natural to simply key in a function, and awkward to then have to specify the domain. However, it is easy to see why most zero finders require a domain. The zero finder does its work by probing the function
     if (n==0) return u
5 at various values. If it probed for a value outside the domain of
     if (n==0) return u
5, the code for
     if (n==0) return u
5 might well compute 0/0 or, and the computation would halt, unnecessarily aborting the zero finding process

This problem can be avoided by introducing a special value called NaN, and specifying that the computation of expressions like 0/0 andproduce NaN, rather than halting. A list of some of the situations that can cause a NaN are given in TABLE D-3. Then when

     if (n==0) return u
1 probes outside the domain of
     if (n==0) return u
5, the code for
     if (n==0) return u
5 will return NaN, and the zero finder can continue. Tức là,
     if (n==0) return u
1 không bị "trừng phạt" vì đoán sai. With this example in mind, it is easy to see what the result of combining a NaN with an ordinary floating-point number should be. Suppose that the final statement of
     if (n==0) return u
5 is
     n = n/2
13 
     n = n/2
14. If d < 0, then
     if (n==0) return u
5 should return a NaN. Since d < 0,
     n = n/2
16 is a NaN, and
     n = n/2
17 will be a NaN, if the sum of a NaN and any other number is a NaN. Similarly if one operand of a division operation is a NaN, the quotient should be a NaN. In general, whenever a NaN participates in a floating-point operation, the result is another NaN

TABLE D-3   Operations That Produce a NaN+
Trăn trừ nhị phân
+ (-
Trăn trừ nhị phân
)×0 ×
Trăn trừ nhị phân
/0/0,
Trăn trừ nhị phân
/
Trăn trừ nhị phân
     n = n/2
18x
     n = n/2
18 0,
Trăn trừ nhị phân
     n = n/2
18 y(when x < 0)

Another approach to writing a zero solver that doesn't require the user to input a domain is to use signals. The zero-finder could install a signal handler for floating-point exceptions. Then if

     if (n==0) return u
5 was evaluated outside its domain and raised an exception, control would be returned to the zero solver. The problem with this approach is that every language has a different method of handling signals (if it has a method at all), and so it has no hope of portability

In IEEE 754, NaNs are often represented as floating-point numbers with the exponent emax + 1 and nonzero significands. Implementations are free to put system-dependent information into the significand. Thus there is not a unique NaN, but rather a whole family of NaNs. When a NaN and an ordinary floating-point number are combined, the result should be the same as the NaN operand. Thus if the result of a long computation is a NaN, the system-dependent information in the significand will be the information that was generated when the first NaN in the computation was generated. Actually, there is a caveat to the last statement. If both operands are NaNs, then the result will be one of those NaNs, but it might not be the NaN that was generated first

Infinity

Just as NaNs provide a way to continue a computation when expressions like 0/0 orare encountered, infinities provide a way to continue when an overflow occurs. This is much safer than simply returning the largest representable number. As an example, consider computing, when

Trăn trừ nhị phân
 = 10, p = 3, and emax = 98. If x = 3 × 1070 and y = 4 × 1070, then x2 will overflow, and be replaced by 9. 99 × 1098. Similarly y2, and x2 + y2 will each overflow in turn, and be replaced by 9. 99 × 1098. So the final result will be, which is drastically wrong. the correct answer is 5 × 1070. In IEEE arithmetic, the result of x2 is
Trăn trừ nhị phân
, as is y2, x2 + y2 and. So the final result is
Trăn trừ nhị phân
, which is safer than returning an ordinary floating-point number that is nowhere near the correct answer. 17

The division of 0 by 0 results in a NaN. A nonzero number divided by 0, however, returns infinity. 1/0 =

Trăn trừ nhị phân
, -1/0 = -
Trăn trừ nhị phân
. The reason for the distinction is this. if f(x)
Trăn trừ nhị phân
0 and g(x)
Trăn trừ nhị phân
0 as x approaches some limit, then f(x)/g(x) could have any value. For example, when f(x) = sin x and g(x) = x, then f(x)/g(x)
Trăn trừ nhị phân
1 as x
Trăn trừ nhị phân
0. But when f(x) = 1 - cos x, f(x)/g(x)
Trăn trừ nhị phân
0. When thinking of 0/0 as the limiting situation of a quotient of two very small numbers, 0/0 could represent anything. Thus in the IEEE standard, 0/0 results in a NaN. But when c > 0, f(x)
Trăn trừ nhị phân
c, and g(x)
Trăn trừ nhị phân
0, then f(x)/g(x) 
Trăn trừ nhị phân
 ±
Trăn trừ nhị phân
, for any analytic functions f and g. If g(x) < 0 for small x, then f(x)/g(x)
Trăn trừ nhị phân
-
Trăn trừ nhị phân
, otherwise the limit is +
Trăn trừ nhị phân
. So the IEEE standard defines c/0 = ±
Trăn trừ nhị phân
, as long as c
Trăn trừ nhị phân
0. The sign of
Trăn trừ nhị phân
depends on the signs of c and 0 in the usual way, so that -10/0 = -
Trăn trừ nhị phân
, and -10/-0 = +
Trăn trừ nhị phân
. You can distinguish between getting
Trăn trừ nhị phân
because of overflow and getting
Trăn trừ nhị phân
because of division by zero by checking the status flags (which will be discussed in detail in section Flags). The overflow flag will be set in the first case, the division by zero flag in the second.

The rule for determining the result of an operation that has infinity as an operand is simple. replace infinity with a finite number x and take the limit as x

Trăn trừ nhị phân
Trăn trừ nhị phân
. Thus 3/
Trăn trừ nhị phân
 = 0, because

.

Similarly, 4 -

Trăn trừ nhị phân
= -
Trăn trừ nhị phân
, and = 
Trăn trừ nhị phân
. When the limit doesn't exist, the result is a NaN, so
Trăn trừ nhị phân
/
Trăn trừ nhị phân
will be a NaN (TABLE D-3 has additional examples). This agrees with the reasoning used to conclude that 0/0 should be a NaN.

Khi một biểu thức con ước tính thành NaN, giá trị của toàn bộ biểu thức cũng là một NaN. In the case of ±

Trăn trừ nhị phân
however, the value of the expression might be an ordinary floating-point number because of rules like 1/
Trăn trừ nhị phân
= 0. Here is a practical example that makes use of the rules for infinity arithmetic. Consider computing the function x/(x2 + 1). This is a bad formula, because not only will it overflow when x is larger than, but infinity arithmetic will give the wrong answer because it will yield 0, rather than a number near 1/x. However, x/(x2 + 1) can be rewritten as 1/(x + x-1). This improved expression will not overflow prematurely and because of infinity arithmetic will have the correct value when x = 0. 1/(0 + 0-1) = 1/(0 +
Trăn trừ nhị phân
) = 1/
Trăn trừ nhị phân
= 0. Without infinity arithmetic, the expression 1/(x + x-1) requires a test for x = 0, which not only adds extra instructions, but may also disrupt a pipeline. This example illustrates a general fact, namely that infinity arithmetic often avoids the need for special case checking; however, formulas need to be carefully inspected to make sure they do not have spurious behavior at infinity (as x/(x2 + 1) did).

Signed Zero

Zero is represented by the exponent emin - 1 and a zero significand. Since the sign bit can take on two different values, there are two zeros, +0 and -0. If a distinction were made when comparing +0 and -0, simple tests like

     n = n/2
002 
     n = n/2
003 
     n = n/2
004 
     n = n/2
005 would have very unpredictable behavior, depending on the sign of
     n = n/2
006. Thus the IEEE standard defines comparison so that +0 = -0, rather than -0 < +0. Although it would be possible always to ignore the sign of zero, the IEEE standard does not do so. When a multiplication or division involves a signed zero, the usual sign rules apply in computing the sign of the answer. Thus 3·(+0) = +0, and +0/-3 = -0. If zero did not have a sign, then the relation 1/(1/x) = x would fail to hold when x = ±
Trăn trừ nhị phân
. The reason is that 1/-
Trăn trừ nhị phân
and 1/+
Trăn trừ nhị phân
both result in 0, and 1/0 results in +
Trăn trừ nhị phân
, the sign information having been lost. One way to restore the identity 1/(1/x) = x is to only have one kind of infinity, however that would result in the disastrous consequence of losing the sign of an overflowed quantity.

Một ví dụ khác về việc sử dụng luồng dưới mối quan tâm bằng 0 có dấu và các hàm có gián đoạn ở 0, chẳng hạn như nhật ký. In IEEE arithmetic, it is natural to define log 0 = -

Trăn trừ nhị phân
and log x to be a NaN when x < 0. Suppose that x represents a small negative number that has underflowed to zero. Thanks to signed zero, x will be negative, so log can return a NaN. However, if there were no signed zero, the log function could not distinguish an underflowed negative number from 0, and would therefore have to return -
Trăn trừ nhị phân
. Another example of a function with a discontinuity at zero is the signum function, which returns the sign of a number.

Probably the most interesting use of signed zero occurs in complex arithmetic. To take a simple example, consider the equation. Điều này chắc chắn đúng khi z

Trăn trừ nhị phân
0. Nếu z = -1, phép tính rõ ràng cho và. Thus,. The problem can be traced to the fact that square root is multi-valued, and there is no way to select the values so that it is continuous in the entire complex plane. However, square root is continuous if a branch cut consisting of all negative real numbers is excluded from consideration. This leaves the problem of what to do for the negative real numbers, which are of the form -x + i0, where x > 0. Signed zero provides a perfect way to resolve this problem. Numbers of the form x + i(+0) have one signand numbers of the form x + i(-0) on the other side of the branch cut have the other sign. In fact, the natural formulas for computingwill give these results.

Back to. If z =1 = -1 + i0, then

1/z = 1/(-1 + i0) = [(-1- i0)]/[(-1 + i0)(-1 - i0)] = (-1 -- i0)/((-1)2 - 02) = -1 + i(-0),

and so, while. Thus IEEE arithmetic preserves this identity for all z. Some more sophisticated examples are given by Kahan [1987]. Although distinguishing between +0 and -0 has advantages, it can occasionally be confusing. For example, signed zero destroys the relation x = y 

Trăn trừ nhị phân
 1/x = 1/y, which is false when x = +0 and y = -0. However, the IEEE committee decided that the advantages of utilizing the sign of zero outweighed the disadvantages.

số không chuẩn hóa

Consider normalized floating-point numbers with

Trăn trừ nhị phân
= 10, p = 3, and emin = -98. The numbers x = 6. 87 × 10-97 and y = 6. 81 × 10-97 appear to be perfectly ordinary floating-point numbers, which are more than a factor of 10 larger than the smallest floating-point number 1. 00 × 10-98. They have a strange property, however. xy = 0 even though x
Trăn trừ nhị phân
y. The reason is that x - y = . 06 × 10 -97  = 6. 0 × 10-99 is too small to be represented as a normalized number, and so must be flushed to zero. How important is it to preserve the property

(10) x = y
Trăn trừ nhị phân
x - y = 0 ?

It's very easy to imagine writing the code fragment,

     n = n/2
002 
     n = n/2
003 
Trăn trừ nhị phân
 
     n = n/2
009 
     n = n/2
010 
     n = n/2
011 
     n = n/2
004 
     n = n/2
013, and much later having a program fail due to a spurious division by zero. Tracking down bugs like this is frustrating and time consuming. On a more philosophical level, computer science textbooks often point out that even though it is currently impractical to prove large programs correct, designing programs with the idea of proving them often results in better code. For example, introducing invariants is quite useful, even if they aren't going to be used as part of a proof. Floating-point code is just like any other code. it helps to have provable facts on which to depend. For example, when analyzing formula (6), it was very helpful to know that x/2 
Trăn trừ nhị phân
 x y = x - y. Similarly, knowing that (10) is true makes writing reliable floating-point code easier. If it is only true for most numbers, it cannot be used to prove anything.

The IEEE standard uses denormalized18 numbers, which guarantee (10), as well as other useful relations. They are the most controversial part of the standard and probably accounted for the long delay in getting 754 approved. Most high performance hardware that claims to be IEEE compatible does not support denormalized numbers directly, but rather traps when consuming or producing denormals, and leaves it to software to simulate the IEEE standard. 19 The idea behind denormalized numbers goes back to Goldberg [1967] and is very simple. When the exponent is emin, the significand does not have to be normalized, so that when

Trăn trừ nhị phân
= 10, p = 3 and emin = -98, 1. 00 × 10-98 is no longer the smallest floating-point number, because 0. 98 × 10-98 is also a floating-point number.

There is a small snag when

Trăn trừ nhị phân
= 2 and a hidden bit is being used, since a number with an exponent of emin will always have a significand greater than or equal to 1. 0 because of the implicit leading bit. The solution is similar to that used to represent 0, and is summarized in TABLE D-2. The exponent emin is used to represent denormals. More formally, if the bits in the significand field are b1, b2, . , bp -1, and the value of the exponent is e, then when e > emin - 1, the number being represented is 1. b1b2. bp - 1 × 2e whereas when e = emin - 1, the number being represented is 0. b1b2. bp - 1 × 2e + 1. The +1 in the exponent is needed because denormals have an exponent of emin, not emin - 1.

Nhắc lại ví dụ về

Trăn trừ nhị phân
= 10, p = 3, emin = -98, x = 6. 87 × 10-97 and y = 6. 81 × 10-97 presented at the beginning of this section. With denormals, x - y does not flush to zero but is instead represented by the denormalized number . 6 × 10-98. This behavior is called gradual underflow. It is easy to verify that (10) always holds when using gradual underflow.

Trăn trừ nhị phân

FIGURE D-2 Flush To Zero Compared With Gradual Underflow

FIGURE D-2 illustrates denormalized numbers. The top number line in the figure shows normalized floating-point numbers. Notice the gap between 0 and the smallest normalized number. If the result of a floating-point calculation falls into this gulf, it is flushed to zero. The bottom number line shows what happens when denormals are added to the set of floating-point numbers. The "gulf" is filled in, and when the result of a calculation is less than, it is represented by the nearest denormal. When denormalized numbers are added to the number line, the spacing between adjacent floating-point numbers varies in a regular way. adjacent spacings are either the same length or differ by a factor of

Trăn trừ nhị phân
. Without denormals, the
spacing abruptly changes fromto, which is a factor of, rather than the orderly change by a factor of
Trăn trừ nhị phân
. Because of this, many algorithms that can have large relative error for normalized numbers close to the underflow threshold are well-behaved in this range when gradual underflow is used.

Without gradual underflow, the simple expression x - y can have a very large relative error for normalized inputs, as was seen above for x = 6. 87 × 10-97 and y = 6. 81 × 10-97. Large relative errors can happen even without cancellation, as the following example shows [Demmel 1984]. Consider dividing two complex numbers, a + ib and c + id. The obvious formula

Trăn trừ nhị phân
· i

suffers from the problem that if either component of the denominator c + id is larger than, the formula will overflow, even though the final result may be well within range. A better method of computing the quotients is to use Smith's formula

(11)
Trăn trừ nhị phân

Applying Smith's formula to (2 · 10-98 + i10-98)/(4 · 10-98 + i(2 · 10-98)) gives the correct answer of 0. 5 with gradual underflow. It yields 0. 4 với flush to zero, lỗi 100 ulps. It is typical for denormalized numbers to guarantee error bounds for arguments all the way down to 1. 0 x.

Exceptions, Flags and Trap Handlers

Khi một điều kiện đặc biệt như chia cho 0 hoặc tràn xảy ra trong số học IEEE, mặc định là cung cấp kết quả và tiếp tục. Typical of the default results are NaN for 0/0 and, and

Trăn trừ nhị phân
for 1/0 and overflow. The preceding sections gave examples where proceeding from an exception with these default values was the reasonable thing to do. When any exception occurs, a status flag is also set. Implementations of the IEEE standard are required to provide users with a way to read and write the status flags. The flags are "sticky" in that once set, they remain set until explicitly cleared. Testing the flags is the only way to distinguish 1/0, which is a genuine infinity from an overflow.

Sometimes continuing execution in the face of exception conditions is not appropriate. The section Infinity gave the example of x/(x2 + 1). When x >, the denominator is infinite, resulting in a final answer of 0, which is totally wrong. Although for this formula the problem can be solved by rewriting it as 1/(x + x-1), rewriting may not always solve the problem. The IEEE standard strongly recommends that implementations allow trap handlers to be installed. Then when an exception occurs, the trap handler is called instead of setting the flag. The value returned by the trap handler will be used as the result of the operation. Trình xử lý bẫy có trách nhiệm xóa hoặc đặt cờ trạng thái;

The IEEE standard divides exceptions into 5 classes. overflow, underflow, division by zero, invalid operation and inexact. There is a separate status flag for each class of exception. The meaning of the first three exceptions is self-evident. Invalid operation covers the situations listed in TABLE D-3, and any comparison that involves a NaN. The default result of an operation that causes an invalid exception is to return a NaN, but the converse is not true. When one of the operands to an operation is a NaN, the result is a NaN but no invalid exception is raised unless the operation also satisfies one of the conditions in TABLE D-3. 20

TABLE D-4   Exceptions in IEEE 754*ExceptionResult when traps disabledArgument to trap handleroverflow±
Trăn trừ nhị phân
or ±xmaxround(x2-
Trăn trừ nhị phân
)underflow0,or denormalround(x2
Trăn trừ nhị phân
)divide by zero±
Trăn trừ nhị phân
operandsinvalidNaNoperandsinexactround(x)round(x)

*x is the exact result of the operation,

Trăn trừ nhị phân
= 192 for single precision, 1536 for double, and xmax = 1. 11 . 11 ×.

The inexact exception is raised when the result of a floating-point operation is not exact. In the

Trăn trừ nhị phân
= 10, p = 3 system, 3. 5
Trăn trừ nhị phân
4. 2 = 14. 7 is exact, but 3. 5 
Trăn trừ nhị phân
 4. 3 = 15. 0 is not exact (since 3. 5 · 4. 3 = 15. 05), and raises an inexact exception. Binary to Decimal Conversion discusses an algorithm that uses the inexact exception. A summary of the behavior of all five exceptions is given in TABLE D-4.

There is an implementation issue connected with the fact that the inexact exception is raised so often. If floating-point hardware does not have flags of its own, but instead interrupts the operating system to signal a floating-point exception, the cost of inexact exceptions could be prohibitive. This cost can be avoided by having the status flags maintained by software. The first time an exception is raised, set the software flag for the appropriate class, and tell the floating-point hardware to mask off that class of exceptions. Then all further exceptions will run without interrupting the operating system. When a user resets that status flag, the hardware mask is re-enabled

Trap Handlers

One obvious use for trap handlers is for backward compatibility. Old codes that expect to be aborted when exceptions occur can install a trap handler that aborts the process. This is especially useful for codes with a loop like

     n = n/2
014 
     n = n/2
015 
     n = n/2
016 
     n = n/2
003
     n = n/2
018
     n = n/2
019. Since comparing a NaN to a number with <,
, >, , or = (but not ) always returns false, this code will go into an infinite loop if
     n = n/2
006 ever becomes a NaN.
Trăn trừ nhị phân
, >,
Trăn trừ nhị phân
, or = (but not
Trăn trừ nhị phân
) always returns false, this code will go into an infinite loop if
     n = n/2
006 ever becomes a NaN.

There is a more interesting use for trap handlers that comes up when computing products such asthat could potentially overflow. One solution is to use logarithms, and compute expinstead. Vấn đề với cách tiếp cận này là nó kém chính xác hơn và tốn kém hơn biểu thức đơn giản, ngay cả khi không có tràn. There is another solution using trap handlers called over/underflow counting that avoids both of these problems [Sterbenz 1974]

The idea is as follows. There is a global counter initialized to zero. Whenever the partial productoverflows for some k, the trap handler increments the counter by one and returns the overflowed quantity with the exponent wrapped around. In IEEE 754 single precision, emax = 127, so if pk = 1. 45 × 2130, it will overflow and cause the trap handler to be called, which will wrap the exponent back into range, changing pk to 1. 45 × 2-62 (see below). Similarly, if pk underflows, the counter would be decremented, and negative exponent would get wrapped around into a positive one. When all the multiplications are done, if the counter is zero then the final product is pn. If the counter is positive, the product overflowed, if the counter is negative, it underflowed. If none of the partial products are out of range, the trap handler is never called and the computation incurs no extra cost. Even if there are over/underflows, the calculation is more accurate than if it had been computed with logarithms, because each pk was computed from pk - 1 using a full precision multiply. Barnett [1987] discusses a formula where the full accuracy of over/underflow counting turned up an error in earlier tables of that formula

IEEE 754 specifies that when an overflow or underflow trap handler is called, it is passed the wrapped-around result as an argument. The definition of wrapped-around for overflow is that the result is computed as if to infinite precision, then divided by 2

Trăn trừ nhị phân
, and then rounded to the relevant precision. For underflow, the result is multiplied by 2
Trăn trừ nhị phân
. The exponent
Trăn trừ nhị phân
is 192 for single precision and 1536 for double precision. This is why 1. 45 x 2130 was transformed into 1. 45 × 2-62 in the example above.

Rounding Modes

In the IEEE standard, rounding occurs whenever an operation has a result that is not exact, since (with the exception of binary decimal conversion) each operation is computed exactly and then rounded. By default, rounding means round toward nearest. Tiêu chuẩn yêu cầu cung cấp ba chế độ làm tròn khác, đó là làm tròn về 0, làm tròn về phía +

Trăn trừ nhị phân
và làm tròn về phía -
Trăn trừ nhị phân
. When used with the convert to integer operation, round toward -
Trăn trừ nhị phân
causes the convert to become the floor function, while round toward +
Trăn trừ nhị phân
is ceiling. The rounding mode affects overflow, because when round toward 0 or round toward -
Trăn trừ nhị phân
is in effect, an overflow of positive magnitude causes the default result to be the largest representable number, not +
Trăn trừ nhị phân
. Similarly, overflows of negative magnitude will produce the largest negative number when round toward +
Trăn trừ nhị phân
or round toward 0 is in effect.

One application of rounding modes occurs in interval arithmetic (another is mentioned in Binary to Decimal Conversion). When using interval arithmetic, the sum of two numbers x and y is an interval, whereis x

Trăn trừ nhị phân
y rounded toward -
Trăn trừ nhị phân
, andis x
Trăn trừ nhị phân
y rounded toward +
Trăn trừ nhị phân
. The exact result of the addition is contained within the interval. Without rounding modes, interval arithmetic is usually implemented by computingand, whereis machine epsilon. 21 Điều này dẫn đến đánh giá quá cao kích thước của các khoảng. Since the result of an operation in interval arithmetic is an interval, in general the input to an operation will also be an interval. If two intervals, and, are added, the result is, whereiswith the rounding mode set to round toward -
Trăn trừ nhị phân
, andiswith the rounding mode set to round toward +
Trăn trừ nhị phân
.

When a floating-point calculation is performed using interval arithmetic, the final answer is an interval that contains the exact result of the calculation. This is not very helpful if the interval turns out to be large (as it often does), since the correct answer could be anywhere in that interval. Interval arithmetic makes more sense when used in conjunction with a multiple precision floating-point package. The calculation is first performed with some precision p. If interval arithmetic suggests that the final answer may be inaccurate, the computation is redone with higher and higher precisions until the final interval is a reasonable size

Flags

The IEEE standard has a number of flags and modes. As discussed above, there is one status flag for each of the five exceptions. underflow, overflow, division by zero, invalid operation and inexact. There are four rounding modes. round toward nearest, round toward +

Trăn trừ nhị phân
, round toward 0, and round toward -
Trăn trừ nhị phân
. It is strongly recommended that there be an enable mode bit for each of the five exceptions. This section gives some simple examples of how these modes and flags can be put to good use. A more sophisticated example is discussed in the section Binary to Decimal Conversion.

Consider writing a subroutine to compute xn, where n is an integer. When n > 0, a simple routine like

     n = n/2
9
     n = n/2
0
     n = n/2
1
     n = n/2
     n = n/2
3
     n = n/2
4
     n = n/2
5
     n = n/2
     if (n==0) return u
     n = n/2
1
     n = n/2
00
     n = n/2
3

If n < 0, then a more accurate way to compute xn is not to call

     n = n/2
021
     n = n/2
022 but rather
     n = n/2
023
     n = n/2
022, because the first expression multiplies n quantities each of which have a rounding error from the division (i. e. , 1/x). In the second expression these are exact (i. e. , x), and the final division commits just one additional rounding error. Unfortunately, these is a slight snag in this strategy. If
     n = n/2
025
     n = n/2
022 underflows, then either the underflow trap handler will be called, or else the underflow status flag will be set. This is incorrect, because if x-n underflows, then xn will either overflow or be in range. 22 But since the IEEE standard gives the user access to all the flags, the subroutine can easily correct for this. It simply turns off the overflow and underflow trap enable bits and saves the overflow and underflow status bits. It then computes
     n = n/2
023
     n = n/2
022. If neither the overflow nor underflow status bit is set, it restores them together with the trap enable bits. If one of the status bits is set, it restores the flags and redoes the calculation using
     n = n/2
021
     n = n/2
022, which causes the correct exceptions to occur

Another example of the use of flags occurs when computing arccos via the formula

arccos x = 2 arctan
Trăn trừ nhị phân
.

If arctan(

Trăn trừ nhị phân
) evaluates to
Trăn trừ nhị phân
/2, then arccos(-1) will correctly evaluate to 2·arctan(
Trăn trừ nhị phân
) =
Trăn trừ nhị phân
, because of infinity arithmetic. However, there is a small snag, because the computation of (1 - x)/(1 + x) will cause the divide by zero exception flag to be set, even though arccos(-1) is not exceptional. The solution to this problem is straightforward. Simply save the value of the divide by zero flag before computing arccos, and then restore its old value after the computation.

Systems Aspects

The design of almost every aspect of a computer system requires knowledge about floating-point. Computer architectures usually have floating-point instructions, compilers must generate those floating-point instructions, and the operating system must decide what to do when exception conditions are raised for those floating-point instructions. Computer system designers rarely get guidance from numerical analysis texts, which are typically aimed at users and writers of software, not at computer designers. As an example of how plausible design decisions can lead to unexpected behavior, consider the following BASIC program

     n = n/2
02
     n = n/2
03
     n = n/2
04

When compiled and run using Borland's Turbo Basic on an IBM PC, the program prints

     n = n/2
031
     n = n/2
032. This example will be analyzed in the next section

Incidentally, some people think that the solution to such anomalies is never to compare floating-point numbers for equality, but instead to consider them equal if they are within some error bound E. This is hardly a cure-all because it raises as many questions as it answers. What should the value of E be? If x < 0 and y > 0 are within E, should they really be considered to be equal, even though they have different signs? Furthermore, the relation defined by this rule, a ~ b 

Trăn trừ nhị phân
. a - b. < E, is not an equivalence relation because a ~ b and b ~ c does not imply that a ~ c.

Instruction Sets

It is quite common for an algorithm to require a short burst of higher precision in order to produce accurate results. One example occurs in the quadratic formula ()/2a. As discussed in the section Proof of Theorem 4, when b2

Trăn trừ nhị phân
4ac, rounding error can contaminate up to half the digits in the roots computed with the quadratic formula. By performing the subcalculation of b2 - 4ac in double precision, half the double precision bits of the root are lost, which means that all the single precision bits are preserved.

The computation of b2 - 4ac in double precision when each of the quantities a, b, and c are in single precision is easy if there is a multiplication instruction that takes two single precision numbers and produces a double precision result. In order to produce the exactly rounded product of two p-digit numbers, a multiplier needs to generate the entire 2p bits of product, although it may throw bits away as it proceeds. Thus, hardware to compute a double precision product from single precision operands will normally be only a little more expensive than a single precision multiplier, and much cheaper than a double precision multiplier. Despite this, modern instruction sets tend to provide only instructions that produce a result of the same precision as the operands. 23

If an instruction that combines two single precision operands to produce a double precision product was only useful for the quadratic formula, it wouldn't be worth adding to an instruction set. However, this instruction has many other uses. Consider the problem of solving a system of linear equations,

a11x1 + a12x2 + · · · + a1nxn= b1
a21x1 + a22x2 + · · · + a2nxn= b2
· · ·
an1x1 + an2x2 + · · ·+ annxn= bn

which can be written in matrix form as Ax = b, where

Trăn trừ nhị phân

Suppose that a solution x(1) is computed by some method, perhaps Gaussian elimination. There is a simple way to improve the accuracy of the result called iterative improvement. First compute

(12)
Trăn trừ nhị phân
= Ax(1) - b

and then solve the system

(13) Ay =
Trăn trừ nhị phân

Note that if x(1) is an exact solution, then

Trăn trừ nhị phân
is the zero vector, as is y. In general, the computation of
Trăn trừ nhị phân
and y will incur rounding error, so Ay 
Trăn trừ nhị phân
 
Trăn trừ nhị phân
 
Trăn trừ nhị phân
 Ax(1) - b = A(x(1) - x), where x is the (unknown) true solution. Then y 
Trăn trừ nhị phân
 x(1) - x, so an improved estimate for the solution is

(14) x(2) = x(1) - y

The three steps (12), (13), and (14) can be repeated, replacing x(1) with x(2), and x(2) with x(3). This argument that x(i + 1) is more accurate than x(i) is only informal. For more information, see [Golub and Van Loan 1989]

When performing iterative improvement,

Trăn trừ nhị phân
is a vector whose elements are the difference of nearby inexact floating-point numbers, and so can suffer from catastrophic cancellation. Thus iterative improvement is not very useful unless
Trăn trừ nhị phân
= Ax(1) - b is computed in double precision. Once again, this is a case of computing the product of two single precision numbers (A and x(1)), where the full double precision result is needed.

To summarize, instructions that multiply two floating-point numbers and return a product with twice the precision of the operands make a useful addition to a floating-point instruction set. Some of the implications of this for compilers are discussed in the next section

Languages and Compilers

The interaction of compilers and floating-point is discussed in Farnum [1988], and much of the discussion in this section is taken from that paper

Ambiguity

Ideally, a language definition should define the semantics of the language precisely enough to prove statements about programs. While this is usually true for the integer part of a language, language definitions often have a large grey area when it comes to floating-point. Perhaps this is due to the fact that many language designers believe that nothing can be proven about floating-point, since it entails rounding error. Nếu vậy, các phần trước đã chứng minh sự ngụy biện trong lập luận này. This section discusses some common grey areas in language definitions, including suggestions about how to deal with them

Remarkably enough, some languages don't clearly specify that if

     n = n/2
006 is a floating-point variable (with say a value of
     n = n/2
034), then every occurrence of (say)
     n = n/2
035 must have the same value. Ví dụ Ada, dựa trên mô hình của Brown, dường như ngụ ý rằng số học dấu phẩy động chỉ phải thỏa mãn các tiên đề của Brown, và do đó các biểu thức có thể có một trong nhiều giá trị có thể. Suy nghĩ về dấu phẩy động theo cách mờ nhạt này trái ngược hoàn toàn với mô hình IEEE, trong đó kết quả của mỗi phép toán dấu phẩy động được xác định chính xác. Trong mô hình IEEE, chúng ta có thể chứng minh rằng
     n = n/2
036 ước lượng thành
     n = n/2
037 (Định lý 7). Trong mô hình của Brown, chúng ta không thể

Một sự mơ hồ khác trong hầu hết các định nghĩa ngôn ngữ liên quan đến những gì xảy ra khi tràn, tràn và các ngoại lệ khác. The IEEE standard precisely specifies the behavior of exceptions, and so languages that use the standard as a model can avoid any ambiguity on this point

Another grey area concerns the interpretation of parentheses. Due to roundoff errors, the associative laws of algebra do not necessarily hold for floating-point numbers. For example, the expression

     n = n/2
038 has a totally different answer than
     n = n/2
039 when x = 1030, y = -1030 and z = 1 (it is 1 in the former case, 0 in the latter). The importance of preserving parentheses cannot be overemphasized. The algorithms presented in theorems 3, 4 and 6 all depend on it. For example, in Theorem 6, the formula xh = mx - (mx - x) would reduce to xh = x if it weren't for parentheses, thereby destroying the entire algorithm. A language definition that does not require parentheses to be honored is useless for floating-point calculations

Subexpression evaluation is imprecisely defined in many languages. Suppose that

     n = n/2
040 is double precision, but
     n = n/2
006 and
     n = n/2
042 are single precision. Then in the expression
     n = n/2
040 
     n = n/2
044 
     n = n/2
045 is the product performed in single or double precision? Another example. in
     n = n/2
006
     n = n/2
044
     n = n/2
048 where
     n = n/2
049 and
     n = n/2
050 are integers, is the division an integer operation or a floating-point one? There are two ways to deal with this problem, neither of which is completely satisfactory. The first is to require that all variables in an expression have the same type. This is the simplest solution, but has some drawbacks. First of all, languages like Pascal that have subrange types allow mixing subrange variables with integer variables, so it is somewhat bizarre to prohibit mixing single and double precision variables. Another problem concerns constants. In the expression
     n = n/2
051, most languages interpret 0. 1 to be a single precision constant. Now suppose the programmer decides to change the declaration of all the floating-point variables from single to double precision. If 0. 1 is still treated as a single precision constant, then there will be a compile time error. The programmer will have to hunt down and change every floating-point constant

The second approach is to allow mixed expressions, in which case rules for subexpression evaluation must be provided. There are a number of guiding examples. The original definition of C required that every floating-point expression be computed in double precision [Kernighan and Ritchie 1978]. This leads to anomalies like the example at the beginning of this section. The expression

     n = n/2
052 is computed in double precision, but if
     n = n/2
053 is a single-precision variable, the quotient is rounded to single precision for storage. Since 3/7 is a repeating binary fraction, its computed value in double precision is different from its stored value in single precision. Thus the comparison q = 3/7 fails. This suggests that computing every expression in the highest precision available is not a good rule

Another guiding example is inner products. If the inner product has thousands of terms, the rounding error in the sum can become substantial. One way to reduce this rounding error is to accumulate the sums in double precision (this will be discussed in more detail in the section Optimizers). If

     n = n/2
054 is a double precision variable, and
     n = n/2
055 and
     n = n/2
056 are single precision arrays, then the inner product loop will look like
     n = n/2
054 
     n = n/2
004
     n = n/2
054
     n = n/2
044
     n = n/2
061. If the multiplication is done in single precision, than much of the advantage of double precision accumulation is lost, because the product is truncated to single precision just before being added to a double precision variable

A rule that covers both of the previous two examples is to compute an expression in the highest precision of any variable that occurs in that expression. Then

     n = n/2
053 
     n = n/2
004 
     n = n/2
052 will be computed entirely in single precision24 and will have the boolean value true, whereas
     n = n/2
054
     n = n/2
004
     n = n/2
054
     n = n/2
044
     n = n/2
061 will be computed in double precision, gaining the full advantage of double precision accumulation. However, this rule is too simplistic to cover all cases cleanly. If
     n = n/2
070 and
     n = n/2
071 are double precision variables, the expression
     n = n/2
042 
     n = n/2
004 
     n = n/2
006 
     n = n/2
044 
     n = n/2
076 contains a double precision variable, but performing the sum in double precision would be pointless, because both operands are single precision, as is the result

A more sophisticated subexpression evaluation rule is as follows. First assign each operation a tentative precision, which is the maximum of the precisions of its operands. This assignment has to be carried out from the leaves to the root of the expression tree. Then perform a second pass from the root to the leaves. In this pass, assign to each operation the maximum of the tentative precision and the precision expected by the parent. In the case of

     n = n/2
053 
     n = n/2
004 
     n = n/2
052, every leaf is single precision, so all the operations are done in single precision. In the case of
     n = n/2
054 
     n = n/2
004 
     n = n/2
054 
     n = n/2
044 
     n = n/2
061, the tentative precision of the multiply operation is single precision, but in the second pass it gets promoted to double precision, because its parent operation expects a double precision operand. And in
     n = n/2
042 
     n = n/2
004 
     n = n/2
006 
     n = n/2
044 
     n = n/2
076, the addition is done in single precision. Farnum [1988] presents evidence that this algorithm in not difficult to implement

The disadvantage of this rule is that the evaluation of a subexpression depends on the expression in which it is embedded. This can have some annoying consequences. For example, suppose you are debugging a program and want to know the value of a subexpression. You cannot simply type the subexpression to the debugger and ask it to be evaluated, because the value of the subexpression in the program depends on the expression it is embedded in. A final comment on subexpressions. since converting decimal constants to binary is an operation, the evaluation rule also affects the interpretation of decimal constants. This is especially important for constants like

     n = n/2
090 which are not exactly representable in binary

Another potential grey area occurs when a language includes exponentiation as one of its built-in operations. Unlike the basic arithmetic operations, the value of exponentiation is not always obvious [Kahan and Coonen 1982]. If

     n = n/2
091 is the exponentiation operator, then
     n = n/2
092 certainly has the value -27. However,
     n = n/2
093 is problematical. If the
     n = n/2
091 operator checks for integer powers, it would compute
     n = n/2
093 as -3. 03 = -27. On the other hand, if the formula xy = eylogx is used to define
     n = n/2
091 for real arguments, then depending on the log function, the result could be a NaN (using the natural definition of log(x) =
     n = n/2
097 when x < 0). If the FORTRAN
     n = n/2
098 function is used however, then the answer will be -27, because the ANSI FORTRAN standard defines
     n = n/2
099 to be i
Trăn trừ nhị phân
+ log 3 [ANSI 1978]. The programming language Ada avoids this problem by only defining exponentiation for integer powers, while ANSI FORTRAN prohibits raising a negative number to a real power.

Trên thực tế, tiêu chuẩn FORTRAN nói rằng

Any arithmetic operation whose result is not mathematically defined is prohibited

Unfortunately, with the introduction of ±

Trăn trừ nhị phân
by the IEEE standard, the meaning of not mathematically defined is no longer totally clear cut. One definition might be to use the method shown in section Infinity. For example, to determine the value of ab, consider non-constant analytic functions f and g with the property that f(x)
Trăn trừ nhị phân
a and g(x)
Trăn trừ nhị phân
b as x
Trăn trừ nhị phân
0. If f(x)g(x) always approaches the same limit, then this should be the value of ab. This definition would set 2
Trăn trừ nhị phân
 = 
Trăn trừ nhị phân
which seems quite reasonable. In the case of 1. 0
Trăn trừ nhị phân
, when f(x) = 1 and g(x) = 1/x the limit approaches 1, but when f(x) = 1 - x and g(x) = 1/x the limit is e-1. So 1. 0
Trăn trừ nhị phân
, should be a NaN. In the case of 00, f(x)g(x) = eg(x)log f(x). Since f and g are analytic and take on the value 0 at 0, f(x) = a1x1 + a2x2 + . and g(x) = b1x1 + b2x2 + . Thus limx
Trăn trừ nhị phân
0g(x) log f(x) = limx 
Trăn trừ nhị phân
 0x log(x(a1 + a2x + . )) = limx
Trăn trừ nhị phân
0x log(a1x) = 0. So f(x)g(x)
Trăn trừ nhị phân
e0 = 1 for all f and g, which means that 00 = 1. 25 26 Using this definition would unambiguously define the exponential function for all arguments, and in particular would define
     n = n/2
093 to be -27.

The IEEE Standard

The section The IEEE Standard," discussed many of the features of the IEEE standard. However, the IEEE standard says nothing about how these features are to be accessed from a programming language. Thus, there is usually a mismatch between floating-point hardware that supports the standard and programming languages like C, Pascal or FORTRAN. Some of the IEEE capabilities can be accessed through a library of subroutine calls. For example the IEEE standard requires that square root be exactly rounded, and the square root function is often implemented directly in hardware. This functionality is easily accessed via a library square root routine. However, other aspects of the standard are not so easily implemented as subroutines. For example, most computer languages specify at most two floating-point types, while the IEEE standard has four different precisions (although the recommended configurations are single plus single-extended or single, double, and double-extended). Infinity provides another example. Các hằng để biểu diễn ±

Trăn trừ nhị phân
có thể được cung cấp bởi một chương trình con. But that might make them unusable in places that require constant expressions, such as the initializer of a constant variable.

Một tình huống tinh vi hơn là thao túng trạng thái liên quan đến tính toán, trong đó trạng thái bao gồm các chế độ làm tròn, bit kích hoạt bẫy, trình xử lý bẫy và cờ ngoại lệ. Một cách tiếp cận là cung cấp các chương trình con để đọc và ghi trạng thái. Ngoài ra, một cuộc gọi duy nhất có thể đặt nguyên tử một giá trị mới và trả về giá trị cũ thường hữu ích. Như các ví dụ trong phần Các cờ cho thấy, một kiểu sửa đổi trạng thái IEEE rất phổ biến là chỉ thay đổi nó trong phạm vi của một khối hoặc chương trình con. Do đó, lập trình viên phải tìm từng lối ra khỏi khối và đảm bảo trạng thái được khôi phục. Hỗ trợ ngôn ngữ để đặt trạng thái chính xác trong phạm vi của một khối sẽ rất hữu ích ở đây. Modula-3 là một ngôn ngữ triển khai ý tưởng này cho trình xử lý bẫy [Nelson 1991]

There are a number of minor points that need to be considered when implementing the IEEE standard in a language. Vì x - x = +0 với mọi x,27 (+0) - (+0) = +0. Tuy nhiên, -(+0) = -0, do đó -x không được định nghĩa là 0 - x. Việc giới thiệu NaN có thể gây nhầm lẫn, bởi vì một NaN không bao giờ bằng bất kỳ số nào khác (bao gồm cả NaN khác), vì vậy x = x không còn đúng nữa. Trên thực tế, biểu thức x

Trăn trừ nhị phân
x là cách đơn giản nhất để kiểm tra NaN nếu chức năng đề xuất của IEEE
     n = n/2
101 không được cung cấp. Hơn nữa, NaN không có thứ tự đối với tất cả các số khác, vì vậy x
Trăn trừ nhị phân
y không thể được định nghĩa là không phải x > y. Do việc sử dụng NaN làm cho các số dấu phẩy động trở nên có thứ tự một phần, nên hàm
     n = n/2
102 trả về một trong các

Mặc dù tiêu chuẩn IEEE xác định các phép toán dấu phẩy động cơ bản để trả về NaN nếu bất kỳ toán hạng nào là NaN, đây có thể không phải lúc nào cũng là định nghĩa tốt nhất cho các phép toán phức hợp. Ví dụ: khi tính toán hệ số tỷ lệ thích hợp để sử dụng trong việc vẽ đồ thị, giá trị lớn nhất của một tập hợp giá trị phải được tính toán. Trong trường hợp này, điều hợp lý là thao tác tối đa chỉ cần bỏ qua NaN

Cuối cùng, làm tròn có thể là một vấn đề. Tiêu chuẩn IEEE xác định làm tròn rất chính xác và nó phụ thuộc vào giá trị hiện tại của các chế độ làm tròn. Điều này đôi khi xung đột với định nghĩa làm tròn ẩn trong chuyển đổi loại hoặc hàm

     n = n/2
103 rõ ràng trong các ngôn ngữ. Điều này có nghĩa là các chương trình muốn sử dụng phương pháp làm tròn IEEE không thể sử dụng các ngôn ngữ gốc của ngôn ngữ tự nhiên và ngược lại, các ngôn ngữ gốc sẽ không hiệu quả để triển khai trên số lượng máy IEEE ngày càng tăng

Trình tối ưu hóa

Compiler texts tend to ignore the subject of floating-point. For example Aho et al. [1986] mentions replacing

     n = n/2
104 with
     n = n/2
105, leading the reader to assume that
     n = n/2
106 should be replaced by
     n = n/2
051. However, these two expressions do not have the same semantics on a binary machine, because 0. 1 cannot be represented exactly in binary. This textbook also suggests replacing
     n = n/2
108 by
     n = n/2
109, even though we have seen that these two expressions can have quite different values when y
Trăn trừ nhị phân
z. Although it does qualify the statement that any algebraic identity can be used when optimizing code by noting that optimizers should not violate the language definition, it leaves the impression that floating-point semantics are not very important. Whether or not the language standard specifies that parenthesis must be honored,
     n = n/2
038 can have a totally different answer than
     n = n/2
039, as discussed above. There is a problem closely related to preserving parentheses that is illustrated by the following code

     n = n/2
05
     n = n/2
06


:

This is designed to give an estimate for machine epsilon. If an optimizing compiler notices that eps + 1 > 1

Trăn trừ nhị phân
eps > 0, the program will be changed completely. Instead of computing the smallest number x such that 1
Trăn trừ nhị phân
x is still greater than x (x
Trăn trừ nhị phân
e
Trăn trừ nhị phân
), it will compute the largest number x for which x/2 is rounded to 0 (x
Trăn trừ nhị phân
). Avoiding this kind of "optimization" is so important that it is worth presenting one more very useful algorithm that is totally ruined by it.

Many problems, such as numerical integration and the numerical solution of differential equations involve computing sums with many terms. Because each addition can potentially introduce an error as large as . 5 ulp, a sum involving thousands of terms can have quite a bit of rounding error. A simple way to correct for this is to store the partial summand in a double precision variable and to perform each addition using double precision. If the calculation is being done in single precision, performing the sum in double precision is easy on most computer systems. However, if the calculation is already being done in double precision, doubling the precision is not so simple. One method that is sometimes advocated is to sort the numbers and add them from smallest to largest. However, there is a much more efficient method which dramatically improves the accuracy of sums, namely

Theorem 8 (Kahan Summation Formula)

Suppose thatis computed using the following algorithm
     n = n/2
07
     n = n/2
08
     n = n/2
09
     n = n/2
10
     n = n/2
11
     n = n/2
12
     n = n/2
13
     n = n/2
14
Then the computed sum S is equal towhere.

Using the naive formula, the computed sum is equal towhere .

Trăn trừ nhị phân
j.  

An optimizer that believed floating-point arithmetic obeyed the laws of algebra would conclude that C = [T-S] - Y = [(S+Y)-S] - Y = 0, rendering the algorithm completely useless. These examples can be summarized by saying that optimizers should be extremely cautious when applying algebraic identities that hold for the mathematical real numbers to expressions involving floating-point variables

Another way that optimizers can change the semantics of floating-point code involves constants. In the expression

     n = n/2
112, there is an implicit decimal to binary conversion operation that converts the decimal number to a binary constant. Because this constant cannot be represented exactly in binary, the inexact exception should be raised. In addition, the underflow flag should to be set if the expression is evaluated in single precision. Since the constant is inexact, its exact conversion to binary depends on the current value of the IEEE rounding modes. Thus an optimizer that converts
     n = n/2
113 to binary at compile time would be changing the semantics of the program. However, constants like 27. 5 which are exactly representable in the smallest available precision can be safely converted at compile time, since they are always exact, cannot raise any exception, and are unaffected by the rounding modes. Constants that are intended to be converted at compile time should be done with a constant declaration, such as
     n = n/2
114
     n = n/2
115
     n = n/2
004
     n = n/2
117

Common subexpression elimination is another example of an optimization that can change floating-point semantics, as illustrated by the following code

     n = n/2
15
     n = n/2
16
     n = n/2
17

Although

     n = n/2
118 can appear to be a common subexpression, it is not because the rounding mode is different at the two evaluation sites. Three final examples. x = x cannot be replaced by the boolean constant
     n = n/2
119, because it fails when x is a NaN; -x = 0 - x fails for x = +0; and x < y is not the opposite of x
Trăn trừ nhị phân
y, because NaNs are neither greater than nor less than ordinary floating-point numbers.

Despite these examples, there are useful optimizations that can be done on floating-point code. First of all, there are algebraic identities that are valid for floating-point numbers. Some examples in IEEE arithmetic are x + y = y + x, 2 ×  x = x + x, 1 × x = x, and 0. 5× x = x/2. However, even these simple identities can fail on a few machines such as CDC and Cray supercomputers. Instruction scheduling and in-line procedure substitution are two other potentially useful optimizations. 28

As a final example, consider the expression

     n = n/2
070 
     n = n/2
004 
     n = n/2
045, where
     n = n/2
006 and
     n = n/2
042 are single precision variables, and
     n = n/2
070 is double precision. On machines that have an instruction that multiplies two single precision numbers to produce a double precision number,
     n = n/2
070 
     n = n/2
004 
     n = n/2
045 can get mapped to that instruction, rather than compiled to a series of instructions that convert the operands to double and then perform a double to double precision multiply

Some compiler writers view restrictions which prohibit converting (x + y) + z to x + (y + z) as irrelevant, of interest only to programmers who use unportable tricks. Perhaps they have in mind that floating-point numbers model real numbers and should obey the same laws that real numbers do. The problem with real number semantics is that they are extremely expensive to implement. Every time two n bit numbers are multiplied, the product will have 2n bits. Every time two n bit numbers with widely spaced exponents are added, the number of bits in the sum is n + the space between the exponents. The sum could have up to (emax - emin) + n bits, or roughly 2·emax + n bits. An algorithm that involves thousands of operations (such as solving a linear system) will soon be operating on numbers with many significant bits, and be hopelessly slow. The implementation of library functions such as sin and cos is even more difficult, because the value of these transcendental functions aren't rational numbers. Exact integer arithmetic is often provided by lisp systems and is handy for some problems. However, exact floating-point arithmetic is rarely useful

The fact is that there are useful algorithms (like the Kahan summation formula) that exploit the fact that (x + y) + z

Trăn trừ nhị phân
x + (y + z), and work whenever the bound

a
Trăn trừ nhị phân
b = (a + b)(1 +
Trăn trừ nhị phân
)

holds (as well as similar bounds for -, × and /). Since these bounds hold for almost all commercial hardware, it would be foolish for numerical programmers to ignore such algorithms, and it would be irresponsible for compiler writers to destroy these algorithms by pretending that floating-point variables have real number semantics

Exception Handling

The topics discussed up to now have primarily concerned systems implications of accuracy and precision. Trap handlers also raise some interesting systems issues. The IEEE standard strongly recommends that users be able to specify a trap handler for each of the five classes of exceptions, and the section Trap Handlers, gave some applications of user defined trap handlers. In the case of invalid operation and division by zero exceptions, the handler should be provided with the operands, otherwise, with the exactly rounded result. Depending on the programming language being used, the trap handler might be able to access other variables in the program as well. For all exceptions, the trap handler must be able to identify what operation was being performed and the precision of its destination

The IEEE standard assumes that operations are conceptually serial and that when an interrupt occurs, it is possible to identify the operation and its operands. On machines which have pipelining or multiple arithmetic units, when an exception occurs, it may not be enough to simply have the trap handler examine the program counter. Hardware support for identifying exactly which operation trapped may be necessary

Another problem is illustrated by the following program fragment

     n = n/2
18
     n = n/2
19
     n = n/2
0
     n = n/2
1

Suppose the second multiply raises an exception, and the trap handler wants to use the value of

     if (n==0) return u
3. On hardware that can do an add and multiply in parallel, an optimizer would probably move the addition operation ahead of the second multiply, so that the add can proceed in parallel with the first multiply. Thus when the second multiply traps,
     if (n==0) return u
     n = n/2
004 
     n = n/2
132 
     n = n/2
044 
     n = n/2
134 has already been executed, potentially changing the result of
     if (n==0) return u
3. It would not be reasonable for a compiler to avoid this kind of optimization, because every floating-point operation can potentially trap, and thus virtually all instruction scheduling optimizations would be eliminated. This problem can be avoided by prohibiting trap handlers from accessing any variables of the program directly. Thay vào đó, trình xử lý có thể được cung cấp toán hạng hoặc kết quả dưới dạng đối số

But there are still problems. In the fragment

the two instructions might well be executed in parallel. If the multiply traps, its argument

     n = n/2
011 could already have been overwritten by the addition, especially since addition is usually faster than multiply. Computer systems that support the IEEE standard must provide some way to save the value of
     n = n/2
011, either in hardware or by having the compiler avoid such a situation in the first place

W. Kahan has proposed using presubstitution instead of trap handlers to avoid these problems. In this method, the user specifies an exception and the value he wants to be used as the result when the exception occurs. As an example, suppose that in code for computing (sin x)/x, the user decides that x = 0 is so rare that it would improve performance to avoid a test for x = 0, and instead handle this case when a 0/0 trap occurs. Using IEEE trap handlers, the user would write a handler that returns a value of 1 and install it before computing sin x/x. Using presubstitution, the user would specify that when an invalid operation occurs, the value 1 should be used. Kahan gọi đây là sự thay thế trước, bởi vì giá trị được sử dụng phải được chỉ định trước khi xảy ra ngoại lệ. Khi sử dụng trình xử lý bẫy, giá trị được trả về có thể được tính khi bẫy xảy ra

The advantage of presubstitution is that it has a straightforward hardware implementation. 29 As soon as the type of exception has been determined, it can be used to index a table which contains the desired result of the operation. Although presubstitution has some attractive attributes, the widespread acceptance of the IEEE standard makes it unlikely to be widely implemented by hardware manufacturers

The Details

A number of claims have been made in this paper concerning properties of floating-point arithmetic. We now proceed to show that floating-point is not black magic, but rather is a straightforward subject whose claims can be verified mathematically. This section is divided into three parts. The first part presents an introduction to error analysis, and provides the details for the section Rounding Error. The second part explores binary to decimal conversion, filling in some gaps from the section The IEEE Standard. The third part discusses the Kahan summation formula, which was used as an example in the section Systems Aspects

Rounding Error

In the discussion of rounding error, it was stated that a single guard digit is enough to guarantee that addition and subtraction will always be accurate (Theorem 2). We now proceed to verify this fact. Theorem 2 has two parts, one for subtraction and one for addition. The part for subtraction is

Theorem 9

If x and y are positive floating-point numbers in a format with parameters

Trăn trừ nhị phân
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

Trăn trừ nhị phân
e
Trăn trừ nhị phân
2e.

Proof

Interchange x and y if necessary so that x > y. It is also harmless to scale x and y so that x is represented by x0. x1 . xp - 1 ×
Trăn trừ nhị phân
0. If y is represented as y0. y1 . yp-1, then the difference is exact. If y is represented as 0. y1 . yp, then the guard digit ensures that the computed difference will be the exact difference rounded to a floating-point number, so the rounding error is at most e. In general, let y = 0. 0 . 0yk + 1 . yk + p andbe y truncated to p + 1 digits. Then(15) y -< (
Trăn trừ nhị phân
- 1)(
Trăn trừ nhị phân
-p - 1 +
Trăn trừ nhị phân
-p - 2 + . +
Trăn trừ nhị phân
-p - k).
From the definition of guard digit, the computed value of x - y is x -rounded to be a floating-point number, that is, (x -) +
Trăn trừ nhị phân
, where the rounding error
Trăn trừ nhị phân
satisfies(16) .
Trăn trừ nhị phân
.
Trăn trừ nhị phân
(
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
-p.
The exact difference is x - y, so the error is (x - y) - (x -+
Trăn trừ nhị phân
) =- y +
Trăn trừ nhị phân
. There are three cases. If x - y
Trăn trừ nhị phân
1 then the relative error is bounded by(17)
Trăn trừ nhị phân
Trăn trừ nhị phân
-p [(
Trăn trừ nhị phân
- 1)(
Trăn trừ nhị phân
-1 + . +
Trăn trừ nhị phân
-k) +
Trăn trừ nhị phân
/2] <
Trăn trừ nhị phân
-p(1 +
Trăn trừ nhị phân
/2) .
Secondly, if x -< 1, then
Trăn trừ nhị phân
= 0. Since the smallest that x - y can be is
Trăn trừ nhị phân
> (
Trăn trừ nhị phân
- 1)(
Trăn trừ nhị phân
-1 + .  +
Trăn trừ nhị phân
-k), where
Trăn trừ nhị phân
=
Trăn trừ nhị phân
- 1,
in this case the relative error is bounded by(18)
Trăn trừ nhị phân
.
The final case is when x - y < 1 but x -
Trăn trừ nhị phân
1. The only way this could happen is if x - = 1, in which case
Trăn trừ nhị phân
= 0. But if
Trăn trừ nhị phân
= 0, then (18) applies, so that again the relative error is bounded by
Trăn trừ nhị phân
-p <
Trăn trừ nhị phân
-p(1 +
Trăn trừ nhị phân
/2). z

When

Trăn trừ nhị phân
= 2, the bound is exactly 2e, and this bound is achieved for x= 1 + 22 - p and y = 21 - p - 21 - 2p in the limit as p
Trăn trừ nhị phân
Trăn trừ nhị phân
. When adding numbers of the same sign, a guard digit is not necessary to achieve good accuracy, as the following result shows.

Định lý 10

If x

Trăn trừ nhị phân
0 and y
Trăn trừ nhị phân
0, then the relative error in computing x + y is at most 2
Trăn trừ nhị phân
, even if no guard digits are used.

Proof

Thuật toán cộng với k chữ số bảo vệ tương tự như thuật toán trừ. Nếu x 
Trăn trừ nhị phân
 y, dịch chuyển y sang phải cho đến khi các điểm cơ số của x và y thẳng hàng. Loại bỏ mọi chữ số đã dịch chuyển qua vị trí p + k. Tính tổng của hai số p + k chữ số này chính xác. Then round to p digits. Chúng tôi sẽ xác minh định lý khi không có số bảo vệ nào được sử dụng; . Không mất tính tổng quát khi giả sử rằng x
Trăn trừ nhị phân
y
Trăn trừ nhị phân
0 và x có dạng d. đ. d ×
Trăn trừ nhị phân
0. Đầu tiên, giả sử không có thực hiện. Then the digits shifted off the end of y have a value less than
Trăn trừ nhị phân
-p + 1, and the sum is at least 1, so the relative error is less than
Trăn trừ nhị phân
-p+1/1 = 2e. If there is a carry out, then the error from shifting must be added to the rounding error of.

The sum is at least

Trăn trừ nhị phân
, so the relative error is less than

Trăn trừ nhị phân
Trăn trừ nhị phân
2
Trăn trừ nhị phân
. z

It is obvious that combining these two theorems gives Theorem 2. Theorem 2 gives the relative error for performing one operation. Comparing the rounding error of x2 - y2 and (x + y) (x - y) requires knowing the relative error of multiple operations. The relative error of xy is

Trăn trừ nhị phân
1 = [(xy) - (x - y)] / (x - y), which satisfies .
Trăn trừ nhị phân
1.  
Trăn trừ nhị phân
 2e. Or to write it another way

(19) xy = (x - y) (1 +
Trăn trừ nhị phân
1), .
Trăn trừ nhị phân
1.
Trăn trừ nhị phân
2e

Similarly

(20) x
Trăn trừ nhị phân
y = (x + y) (1 +
Trăn trừ nhị phân
2), .
Trăn trừ nhị phân
2.
Trăn trừ nhị phân
2e

Assuming that multiplication is performed by computing the exact product and then rounding, the relative error is at most . 5 ulp, so

(21) u
Trăn trừ nhị phân
v = uv (1 +
Trăn trừ nhị phân
3), .
Trăn trừ nhị phân
3.
Trăn trừ nhị phân
e

for any floating-point numbers u and v. Putting these three equations together (letting u = xy and v = x

Trăn trừ nhị phân
y) gives

(22) (xy)
Trăn trừ nhị phân
(x
Trăn trừ nhị phân
y) = (x - y) (1 +
Trăn trừ nhị phân
1) (x + y) (1 +
Trăn trừ nhị phân
2) (1 +
Trăn trừ nhị phân
3)

So the relative error incurred when computing (x - y) (x + y) is

(23)
Trăn trừ nhị phân

This relative error is equal to

Trăn trừ nhị phân
1 +
Trăn trừ nhị phân
2 +
Trăn trừ nhị phân
3 +
Trăn trừ nhị phân
1
Trăn trừ nhị phân
2 +
Trăn trừ nhị phân
1
Trăn trừ nhị phân
3 +
Trăn trừ nhị phân
2
Trăn trừ nhị phân
3 +
Trăn trừ nhị phân
1
Trăn trừ nhị phân
2
Trăn trừ nhị phân
3, which is bounded by 5
Trăn trừ nhị phân
+ 8
Trăn trừ nhị phân
2. In other words, the maximum relative error is about 5 rounding errors (since e is a small number, e2 is almost negligible).

A similar analysis of (x

Trăn trừ nhị phân
x)(y
Trăn trừ nhị phân
y) cannot result in a small value for the relative error, because when two nearby values of x and y are plugged into x2 - y2, the relative error will usually be quite large. Another way to see this is to try and duplicate the analysis that worked on (xy)
Trăn trừ nhị phân
(x
Trăn trừ nhị phân
y), yielding

(x
Trăn trừ nhị phân
x)(y
Trăn trừ nhị phân
y) = [x2(1 +
Trăn trừ nhị phân
1) - y2(1 +
Trăn trừ nhị phân
2)] (1 +
Trăn trừ nhị phân
3)
= ((x2 - y2) (1 +
Trăn trừ nhị phân
1) + (
Trăn trừ nhị phân
1 -
Trăn trừ nhị phân
2)y2) (1 +
Trăn trừ nhị phân
3)

When x and y are nearby, the error term (

Trăn trừ nhị phân
1 -
Trăn trừ nhị phân
2)y2 can be as large as the result x2 - y2. These computations formally justify our claim that (x - y) (x + y) is more accurate than x2 - y2.

We next turn to an analysis of the formula for the area of a triangle. In order to estimate the maximum error that can occur when computing with (7), the following fact will be needed

Theorem 11

If subtraction is performed with a guard digit, and y/2
Trăn trừ nhị phân
x
Trăn trừ nhị phân
2y, then x - y is computed exactly.

Proof

Note that if x and y have the same exponent, then certainly xy is exact. Otherwise, from the condition of the theorem, the exponents can differ by at most 1. Scale and interchange x and y if necessary so that 0
Trăn trừ nhị phân
y
Trăn trừ nhị phân
x, and x is represented as x0. x1 . xp - 1 and y as 0. y1 . yp. Then the algorithm for computing xy will compute x - y exactly and round to a floating-point number. If the difference is of the form 0. d1 . dp, the difference will already be p digits long, and no rounding is necessary. Since x
Trăn trừ nhị phân
2y, x - y
Trăn trừ nhị phân
 y, and since y is of the form 0. d1 . dp, so is x - y. z

When

Trăn trừ nhị phân
> 2, the hypothesis of Theorem 11 cannot be replaced by y/
Trăn trừ nhị phân
Trăn trừ nhị phân
Trăn trừ nhị phân
Trăn trừ nhị phân
y; the stronger condition y/2
Trăn trừ nhị phân
x
Trăn trừ nhị phân
2y is still necessary. The analysis of the error in (x - y) (x + y), immediately following the proof of Theorem 10, used the fact that the relative error in the basic operations of addition and subtraction is small (namely equations (19) and (20)). This is the most common kind of error analysis. However, analyzing formula (7) requires something more, namely Theorem 11, as the following proof will show.

Theorem 12

If subtraction uses a guard digit, and if a,b and c are the sides of a triangle (a 
Trăn trừ nhị phân
 b
Trăn trừ nhị phân
 c), then the relative error in computing (a + (b + c))(c - (a - b))(c + (a - b))(a +(b - c)) is at most 16
Trăn trừ nhị phân
, provided e < . 005.

Proof

Let's examine the factors one by one. From Theorem 10, b 
Trăn trừ nhị phân
 c = (b + c) (1 + 
Trăn trừ nhị phân
1), where
Trăn trừ nhị phân
1 is the relative error, and .
Trăn trừ nhị phân
1.
Trăn trừ nhị phân
2
Trăn trừ nhị phân
. Then the value of the first factor is(a
Trăn trừ nhị phân
(b
Trăn trừ nhị phân
c)) = (a + (b
Trăn trừ nhị phân
c)) (1 +
Trăn trừ nhị phân
2) = (a + (b + c) (1 +
Trăn trừ nhị phân
1))(1 +
Trăn trừ nhị phân
2),
and thus(a + b + c) (1 - 2
Trăn trừ nhị phân
)2
Trăn trừ nhị phân
[a + (b + c) (1 - 2
Trăn trừ nhị phân
)] · (1-2
Trăn trừ nhị phân
)
Trăn trừ nhị phân
a
Trăn trừ nhị phân
(b
Trăn trừ nhị phân
c)
Trăn trừ nhị phân
[a + (b + c) (1 + 2
Trăn trừ nhị phân
)] (1 + 2
Trăn trừ nhị phân
)
Trăn trừ nhị phân
(a + b + c) (1 + 2
Trăn trừ nhị phân
)2
This means that there is an
Trăn trừ nhị phân
1 so that(24) (a
Trăn trừ nhị phân
(b
Trăn trừ nhị phân
c)) = (a + b + c) (1 +
Trăn trừ nhị phân
1)2, .
Trăn trừ nhị phân
1.
Trăn trừ nhị phân
2
Trăn trừ nhị phân
.
The next term involves the potentially catastrophic subtraction of c and a  
     n = n/2
132, because ab may have rounding error. Because a, b and c are the sides of a triangle, a
Trăn trừ nhị phân
b+ c, and combining this with the ordering c
Trăn trừ nhị phân
b
Trăn trừ nhị phân
a gives a
Trăn trừ nhị phân
b + c
Trăn trừ nhị phân
2b
Trăn trừ nhị phân
2a. So a - b satisfies the conditions of Theorem 11. This means that a - b = ab is exact, hence c(a - b) is a harmless subtraction which can be estimated from Theorem 9 to be(25) (c(ab)) = (c - (a - b)) (1 +
Trăn trừ nhị phân
2), .
Trăn trừ nhị phân
2.
Trăn trừ nhị phân
2
Trăn trừ nhị phân

The third term is the sum of two exact positive quantities, so(26) (c
Trăn trừ nhị phân
(ab)) = (c + (a - b)) (1 +
Trăn trừ nhị phân
3), .
Trăn trừ nhị phân
3.
Trăn trừ nhị phân
2
Trăn trừ nhị phân

Finally, the last term is(27) (a
Trăn trừ nhị phân
(bc)) = (a + (b - c)) (1 +
Trăn trừ nhị phân
4)2, .
Trăn trừ nhị phân
4.
Trăn trừ nhị phân
2
Trăn trừ nhị phân
,
using both Theorem 9 and Theorem 10. If multiplication is assumed to be exactly rounded, so that x
Trăn trừ nhị phân
y = xy(1 +
Trăn trừ nhị phân
) with .
Trăn trừ nhị phân
.
Trăn trừ nhị phân
Trăn trừ nhị phân
, then combining (24), (25), (26) and (27) gives(a
Trăn trừ nhị phân
(b
Trăn trừ nhị phân
c)) (c(ab)) (c
Trăn trừ nhị phân
(ab)) (a
Trăn trừ nhị phân
(bc))
Trăn trừ nhị phân
(a + (b + c)) (c - (a - b)) (c + (a - b)) (a + (b - c)) E
whereE = (1 +
Trăn trừ nhị phân
1)2 (1 +
Trăn trừ nhị phân
2) (1 +
Trăn trừ nhị phân
3) (1 +
Trăn trừ nhị phân
4)2 (1 +
Trăn trừ nhị phân
1)(1 +
Trăn trừ nhị phân
2) (1 +
Trăn trừ nhị phân
3)
An upper bound for E is (1 + 2
Trăn trừ nhị phân
)6(1 +
Trăn trừ nhị phân
)3, which expands out to 1 + 15
Trăn trừ nhị phân
 + O(
Trăn trừ nhị phân
2). Some writers simply ignore the O(e2) term, but it is easy to account for it. Writing (1 + 2
Trăn trừ nhị phân
)6(1 +
Trăn trừ nhị phân
)3 = 1 + 15
Trăn trừ nhị phân
+
Trăn trừ nhị phân
R(
Trăn trừ nhị phân
), R(
Trăn trừ nhị phân
) is a polynomial in e with positive coefficients, so it is an increasing function of
Trăn trừ nhị phân
. Since R(. 005) = . 505, R(
Trăn trừ nhị phân
) < 1 for all
Trăn trừ nhị phân
< . 005, and hence E 
Trăn trừ nhị phân
 (1 + 2
Trăn trừ nhị phân
)6(1 + 
Trăn trừ nhị phân
)3 
Trăn trừ nhị phân
. To get a lower bound on E, note that 1 - 15
Trăn trừ nhị phân
 - 
Trăn trừ nhị phân
R(
Trăn trừ nhị phân
) < E, and so when
Trăn trừ nhị phân
< . 005, 1 - 16
Trăn trừ nhị phân
< (1 - 2
Trăn trừ nhị phân
)6(1 -
Trăn trừ nhị phân
)3. Combining these two bounds yields 1 - 16
Trăn trừ nhị phân
< E
Trăn trừ nhị phân
. Thus the relative error is at most 16
Trăn trừ nhị phân
. z

Theorem 12 certainly shows that there is no catastrophic cancellation in formula (7). So although it is not necessary to show formula (7) is numerically stable, it is satisfying to have a bound for the entire formula, which is what Theorem 3 of Cancellation gives

Proof of Theorem 3

Letq = (a + (b + c)) (c - (a - b)) (c + (a - b)) (a + (b - c))
andQ = (a
Trăn trừ nhị phân
(b
Trăn trừ nhị phân
c))
Trăn trừ nhị phân
(c(ab))
Trăn trừ nhị phân
(c
Trăn trừ nhị phân
(ab))
Trăn trừ nhị phân
(a
Trăn trừ nhị phân
(bc)).
Then, Theorem 12 shows that Q = q(1 +
Trăn trừ nhị phân
), with
Trăn trừ nhị phân
Trăn trừ nhị phân
16
Trăn trừ nhị phân
. Dễ dàng kiểm tra rằng(28)
được cung cấp
Trăn trừ nhị phân
Trăn trừ nhị phân
. 04/(. 52)2
Trăn trừ nhị phân
. 15, and since .
Trăn trừ nhị phân
.
Trăn trừ nhị phân
16
Trăn trừ nhị phân
Trăn trừ nhị phân
16(. 005) = . 08,
Trăn trừ nhị phân
does satisfy the condition. Thus,
with .
Trăn trừ nhị phân
1.
Trăn trừ nhị phân
. 52.
Trăn trừ nhị phân
.
Trăn trừ nhị phân
8. 5
Trăn trừ nhị phân
. If square roots are computed to within . 5 ulp, then the error when computingis (1 + 
Trăn trừ nhị phân
1)(1 +
Trăn trừ nhị phân
2), with .
Trăn trừ nhị phân
2.
Trăn trừ nhị phân
Trăn trừ nhị phân
. If
Trăn trừ nhị phân
= 2, then there is no further error committed when dividing by 4. Otherwise, one more factor 1 +
Trăn trừ nhị phân
3 with .
Trăn trừ nhị phân
3.  
Trăn trừ nhị phân
 
Trăn trừ nhị phân
is necessary for the division, and using the method in the proof of Theorem 12, the final error bound of (1 +
Trăn trừ nhị phân
1) (1 +
Trăn trừ nhị phân
2) (1 +
Trăn trừ nhị phân
3) is dominated by 1 +
Trăn trừ nhị phân
4, with .
Trăn trừ nhị phân
4.
Trăn trừ nhị phân
11
Trăn trừ nhị phân
. z

To make the heuristic explanation immediately following the statement of Theorem 4 precise, the next theorem describes just how closely µ(x) approximates a constant

Theorem 13

If µ(x) = ln(1 + x)/x, then for 0
Trăn trừ nhị phân
x
Trăn trừ nhị phân
,
Trăn trừ nhị phân
µ(x)
Trăn trừ nhị phân
1 and the derivative satisfies . µ'(x).  
Trăn trừ nhị phân
 .

Proof

Note that µ(x) = 1 - x/2 + x2/3 - . is an alternating series with decreasing terms, so for x
Trăn trừ nhị phân
1, µ(x)
Trăn trừ nhị phân
1 - x/2
Trăn trừ nhị phân
1/2. It is even easier to see that because the series for µ is alternating, µ(x)
Trăn trừ nhị phân
1. The Taylor series of µ'(x) is also alternating, and if x
Trăn trừ nhị phân
has decreasing terms, so -
Trăn trừ nhị phân
µ'(x)
Trăn trừ nhị phân
-+ 2x/3, or -
Trăn trừ nhị phân
µ'(x)
Trăn trừ nhị phân
0, thus . µ'(x).
Trăn trừ nhị phân
. z

Proof of Theorem 4

Since the Taylor series for ln
Trăn trừ nhị phân

is an alternating series, 0 < x - ln(1 + x) < x2/2, the relative error incurred when approximating ln(1 + x) by x is bounded by x/2. If 1
Trăn trừ nhị phân
x = 1, then . x.  
Trăn trừ nhị phân
, so the relative error is bounded by
Trăn trừ nhị phân
/2. When 1
Trăn trừ nhị phân
x
Trăn trừ nhị phân
1, definevia 1
Trăn trừ nhị phân
x = 1 +. Then since 0
Trăn trừ nhị phân
x < 1, (1
Trăn trừ nhị phân
x)1 =. If division and logarithms are computed to within ulp, then the computed value of the expression ln(1 + x)/((1 + x) - 1) is(29)(1 +
Trăn trừ nhị phân
1) (1 +
Trăn trừ nhị phân
2) =(1 +
Trăn trừ nhị phân
1) (1 +
Trăn trừ nhị phân
2) = µ() (1 +
Trăn trừ nhị phân
1) (1 +
Trăn trừ nhị phân
2)

where .

Trăn trừ nhị phân
1.
Trăn trừ nhị phân
Trăn trừ nhị phân
and .
Trăn trừ nhị phân
2.
Trăn trừ nhị phân
Trăn trừ nhị phân
. To estimate µ(), use the mean value theorem, which says that

(30) µ() - µ(x) = (- x)µ'(
Trăn trừ nhị phân
)
for some
Trăn trừ nhị phân
between x and. From the definition of, it follows that . - x.
Trăn trừ nhị phân
Trăn trừ nhị phân
, and combining this with Theorem 13 gives . µ() - µ(x).
Trăn trừ nhị phân
Trăn trừ nhị phân
/2, or . µ()/µ(x) - 1.
Trăn trừ nhị phân
Trăn trừ nhị phân
/(2. µ(x). )
Trăn trừ nhị phân
Trăn trừ nhị phân
which means that µ() = µ(x) (1 +
Trăn trừ nhị phân
3), with .
Trăn trừ nhị phân
3.
Trăn trừ nhị phân
Trăn trừ nhị phân
. Finally, multiplying by x introduces a final
Trăn trừ nhị phân
4, so the computed value ofx·ln(1 
Trăn trừ nhị phân
x)/((1
Trăn trừ nhị phân
x)1)

is


It is easy to check that if
Trăn trừ nhị phân
< 0. 1, then(1 +
Trăn trừ nhị phân
1) (1 +
Trăn trừ nhị phân
2) (1 +
Trăn trừ nhị phân
3) (1 +
Trăn trừ nhị phân
4) = 1 + 
Trăn trừ nhị phân
,

with .

Trăn trừ nhị phân
.  
Trăn trừ nhị phân
5
Trăn trừ nhị phân
. z

An interesting example of error analysis using formulas (19), (20), and (21) occurs in the quadratic formula. The section Cancellation, explained how rewriting the equation will eliminate the potential cancellation caused by the ± operation. But there is another potential cancellation that can occur when computing d = b2 - 4ac. This one cannot be eliminated by a simple rearrangement of the formula. Roughly speaking, when b2

Trăn trừ nhị phân
4ac, rounding error can contaminate up to half the digits in the roots computed with the quadratic formula. Here is an informal proof (another approach to estimating the error in the quadratic formula appears in Kahan [1972]).

If b2

Trăn trừ nhị phân
4ac, rounding error can contaminate up to half the digits in the roots computed with the quadratic formula.

Proof. Write (b

Trăn trừ nhị phân
b)(4a
Trăn trừ nhị phân
c) = (b2(1 +
Trăn trừ nhị phân
1) - 4ac(1 +
Trăn trừ nhị phân
2)) (1 +
Trăn trừ nhị phân
3), where .
Trăn trừ nhị phân
i.
Trăn trừ nhị phân
 
Trăn trừ nhị phân
. 30 Using d = b2 - 4ac, this can be rewritten as (d(1 +
Trăn trừ nhị phân
1) - 4ac(
Trăn trừ nhị phân
2 -
Trăn trừ nhị phân
1)) (1 +
Trăn trừ nhị phân
3). To get an estimate for the size of this error, ignore second order terms in
Trăn trừ nhị phân
i, in which case the absolute error is d(
Trăn trừ nhị phân
1 +
Trăn trừ nhị phân
3) - 4ac
Trăn trừ nhị phân
4, where .
Trăn trừ nhị phân
4. = .
Trăn trừ nhị phân
1 - 
Trăn trừ nhị phân
2.
Trăn trừ nhị phân
2
Trăn trừ nhị phân
. Since, the first term d(
Trăn trừ nhị phân
1 +
Trăn trừ nhị phân
3) can be ignored. To estimate the second term, use the fact that ax2 + bx + c = a(x - r1) (x - r2), so ar1r2 = c. Since b2
Trăn trừ nhị phân
4ac, then r1
Trăn trừ nhị phân
r2, so the second error term is. Thus the computed value ofis

.

The inequality

shows that

,

where

,

so the absolute error ina is about. Since

Trăn trừ nhị phân
4
Trăn trừ nhị phân
Trăn trừ nhị phân
-p,, and thus the absolute error ofdestroys the bottom half of the bits of the roots r1
Trăn trừ nhị phân
r2. In other words, since the calculation of the roots involves computing with, and this expression does not have meaningful bits in the position corresponding to the lower order half of ri, then the lower order bits of ri cannot be meaningful. z

Finally, we turn to the proof of Theorem 6. It is based on the following fact, which is proven in the section Theorem 14 and Theorem 8

Theorem 14

Let 0 < k < p, and set m =
Trăn trừ nhị phân
k + 1, and assume that floating-point operations are exactly rounded. Then (m
Trăn trừ nhị phân
x)(m
Trăn trừ nhị phân
xx) is exactly equal to x rounded to p - k significant digits. More precisely, x is rounded by taking the significand of x, imagining a radix point just left of the k least significant digits and rounding to an integer.

Proof of Theorem 6

By Theorem 14, xh is x rounded to p - k =places. If there is no carry out, then certainly xh can be represented withsignificant digits. Suppose there is a carry-out. If x = x0. x1 . xp - 1 ×
Trăn trừ nhị phân
e, then rounding adds 1 to xp - k - 1, and the only way there can be a carry-out is if xp - k - 1 =
Trăn trừ nhị phân
- 1, but then the low order digit of xh is 1 + xp - k- 1 = 0, and so again xh is representable indigits. Để xử lý xl, chia tỷ lệ x thành một số nguyên thỏa mãn
Trăn trừ nhị phân
p - 1
Trăn trừ nhị phân
x
Trăn trừ nhị phân
Trăn trừ nhị phân
p - 1. Letwhereis the p - k high order digits of x, andis the k low order digits. Có ba trường hợp cần xem xét. Nếu, thì làm tròn x đến p - k vị trí giống như cắt và, và. Sincehas at most k digits, if p is even, thenhas at most k ==digits. Otherwise,
Trăn trừ nhị phân
= 2 andis representable with k - 1
Trăn trừ nhị phân
significant bits. The second case is when, and then computing xh involves rounding up, so xh =+
Trăn trừ nhị phân
k, and xl = x - xh = x - -
Trăn trừ nhị phân
k =-
Trăn trừ nhị phân
k. Once again,has at most k digits, so is representable with
Trăn trừ nhị phân
p/2
Trăn trừ nhị phân
digits. Finally, if= (
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
k - 1, then xh =or + 
Trăn trừ nhị phân
k depending on whether there is a round up. So xl is either (
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
k - 1 or (
Trăn trừ nhị phân
/2)
Trăn trừ nhị phân
k - 1 - 
Trăn trừ nhị phân
k = -
Trăn trừ nhị phân
k/2, both of which are represented with 1 digit. z

Theorem 6 gives a way to express the product of two working precision numbers exactly as a sum. There is a companion formula for expressing a sum exactly. If . x.

Trăn trừ nhị phân
. y. then x + y = (x
Trăn trừ nhị phân
y) + (x(x
Trăn trừ nhị phân
y))
Trăn trừ nhị phân
y [Dekker 1971; Knuth 1981, Theorem C in section 4. 2. 2]. However, when using exactly rounded operations, this formula is only true for
Trăn trừ nhị phân
= 2, and not for
Trăn trừ nhị phân
= 10 as the example x = . 99998, y = . 99997 shows.

Binary to Decimal Conversion

Since single precision has p = 24, and 224 < 108, you might expect that converting a binary number to 8 decimal digits would be sufficient to recover the original binary number. However, this is not the case

Theorem 15

When a binary IEEE single precision number is converted to the closest eight digit decimal number, it is not always possible to uniquely recover the binary number from the decimal one. However, if nine decimal digits are used, then converting the decimal number to the closest binary number will recover the original floating-point number

Proof

Binary single precision numbers lying in the half open interval [103, 210) = [1000, 1024) have 10 bits to the left of the binary point, and 14 bits to the right of the binary point. Thus there are (210 - 103)214 = 393,216 different binary numbers in that interval. If decimal numbers are represented with 8 digits, then there are (210 - 103)104 = 240,000 decimal numbers in the same interval. There is no way that 240,000 decimal numbers could represent 393,216 different binary numbers. So 8 decimal digits are not enough to uniquely represent each single precision binary number. To show that 9 digits are sufficient, it is enough to show that the spacing between binary numbers is always greater than the spacing between decimal numbers. This will ensure that for each decimal number N, the interval[N -ulp, N +ulp]
contains at most one binary number. Thus each binary number rounds to a unique decimal number which in turn rounds to a unique binary number. To show that the spacing between binary numbers is always greater than the spacing between decimal numbers, consider an interval [10n, 10n + 1]. On this interval, the spacing between consecutive decimal numbers is 10(n + 1) - 9. On [10n, 2m], where m is the smallest integer so that 10n < 2m, the spacing of binary numbers is 2m - 24, and the spacing gets larger further on in the interval. Thus it is enough to check that 10(n + 1) - 9 < 2m - 24. But in fact, since 10n < 2m, then 10(n + 1) - 9 = 10n10-8 < 2m10-8 < 2m2-24. z

The same argument applied to double precision shows that 17 decimal digits are required to recover a double precision number

Binary-decimal conversion also provides another example of the use of flags. Recall from the section Precision, that to recover a binary number from its decimal expansion, the decimal to binary conversion must be computed exactly. That conversion is performed by multiplying the quantities N and 10. P. (which are both exact if p < 13) in single-extended precision and then rounding this to single precision (or dividing if p < 0; both cases are similar). Of course the computation of N · 10. P. cannot be exact; it is the combined operation round(N · 10. P. ) that must be exact, where the rounding is from single-extended to single precision. To see why it might fail to be exact, take the simple case of

Trăn trừ nhị phân
= 10, p = 2 for single, and p = 3 for single-extended. If the product is to be 12. 51, then this would be rounded to 12. 5 as part of the single-extended multiply operation. Rounding to single precision would give 12. But that answer is not correct, because rounding the product to single precision should give 13. The error is due to double rounding.

By using the IEEE flags, double rounding can be avoided as follows. Save the current value of the inexact flag, and then reset it. Set the rounding mode to round-to-zero. Then perform the multiplication N · 10. P. Lưu trữ giá trị mới của cờ không chính xác trong

     n = n/2
139 và khôi phục chế độ làm tròn và cờ không chính xác. If
     n = n/2
139 is 0, then N · 10. P. is exact, so round(N · 10. P. ) will be correct down to the last bit. If
     n = n/2
139 is 1, then some digits were truncated, since round-to-zero always truncates. The significand of the product will look like 1. b1. b22b23. b31. A double rounding error may occur if b23 . b31 = 10. 0. A simple way to account for both cases is to perform a logical
     n = n/2
142 of
     n = n/2
139 with b31. Then round(N · 10. P. ) will be computed correctly in all cases

Errors In Summation

The section Optimizers, mentioned the problem of accurately computing very long sums. The simplest approach to improving accuracy is to double the precision. To get a rough estimate of how much doubling the precision improves the accuracy of a sum, let s1 = x1, s2 = s1 

Trăn trừ nhị phân
 x2. , si = si - 1
Trăn trừ nhị phân
xi. Then si = (1 + 
Trăn trừ nhị phân
i) (si - 1 + xi), where
Trăn trừ nhị phân
Trăn trừ nhị phân
i
Trăn trừ nhị phân
Trăn trừ nhị phân
Trăn trừ nhị phân
, and ignoring second order terms in
Trăn trừ nhị phân
i gives

(31)
Trăn trừ nhị phân

The first equality of (31) shows that the computed value ofis the same as if an exact summation was performed on perturbed values of xj. The first term x1 is perturbed by n

Trăn trừ nhị phân
, the last term xn by only
Trăn trừ nhị phân
. The second equality in (31) shows that error term is bounded by. Doubling the precision has the effect of squaring
Trăn trừ nhị phân
. If the sum is being done in an IEEE double precision format, 1/
Trăn trừ nhị phân
Trăn trừ nhị phân
1016, so thatfor any reasonable value of n. Thus, doubling the precision takes the maximum perturbation of n
Trăn trừ nhị phân
and changes it to. Thus the 2
Trăn trừ nhị phân
error bound for the Kahan summation formula (Theorem 8) is not as good as using double precision, even though it is much better than single precision.

For an intuitive explanation of why the Kahan summation formula works, consider the following diagram of the procedure

Trăn trừ nhị phân

Each time a summand is added, there is a correction factor C which will be applied on the next loop. So first subtract the correction C computed in the previous loop from Xj, giving the corrected summand Y. Then add this summand to the running sum S. The low order bits of Y (namely Yl) are lost in the sum. Next compute the high order bits of Y by computing T - S. When Y is subtracted from this, the low order bits of Y will be recovered. These are the bits that were lost in the first sum in the diagram. They become the correction factor for the next loop. A formal proof of Theorem 8, taken from Knuth [1981] page 572, appears in the section Theorem 14 and Theorem 8. "

Summary

It is not uncommon for computer system designers to neglect the parts of a system related to floating-point. This is probably due to the fact that floating-point is given very little (if any) attention in the computer science curriculum. This in turn has caused the apparently widespread belief that floating-point is not a quantifiable subject, and so there is little point in fussing over the details of hardware and software that deal with it

This paper has demonstrated that it is possible to reason rigorously about floating-point. For example, floating-point algorithms involving cancellation can be proven to have small relative errors if the underlying hardware has a guard digit, and there is an efficient algorithm for binary-decimal conversion that can be proven to be invertible, provided that extended precision is supported. The task of constructing reliable floating-point software is made much easier when the underlying computer system is supportive of floating-point. In addition to the two examples just mentioned (guard digits and extended precision), the section Systems Aspects of this paper has examples ranging from instruction set design to compiler optimization illustrating how to better support floating-point

The increasing acceptance of the IEEE floating-point standard means that codes that utilize features of the standard are becoming ever more portable. The section The IEEE Standard, gave numerous examples illustrating how the features of the IEEE standard can be used in writing practical floating-point codes

Acknowledgments

This article was inspired by a course given by W. Kahan at Sun Microsystems from May through July of 1988, which was very ably organized by David Hough of Sun. My hope is to enable others to learn about the interaction of floating-point and computer systems without having to get up in time to attend 8. 00 a. m. lectures. Thanks are due to Kahan and many of my colleagues at Xerox PARC (especially John Gilbert) for reading drafts of this paper and providing many useful comments. Nhận xét từ Paul Hilfinger và một trọng tài ẩn danh cũng giúp cải thiện phần trình bày

References

Aho, Alfred V. , Sethi, R. , và Ullman J. Đ. 1986. Trình biên dịch. Principles, Techniques and Tools, Addison-Wesley, Reading, MA

ANSI 1978. American National Standard Programming Language FORTRAN, ANSI Standard X3. 9-1978, American National Standards Institute, New York, NY

Barnett, David 1987. A Portable Floating-Point Environment, unpublished manuscript

Brown, W. S. 1981. A Simple but Realistic Model of Floating-Point Computation, ACM Trans. on Math. Software 7(4), pp. 445-480

Cody, W. J et. al. 1984. A Proposed Radix- and Word-length-independent Standard for Floating-point Arithmetic, IEEE Micro 4(4), pp. 86-100

Cody, W. J. 1988. Floating-Point Standards -- Theory and Practice, in "Reliability in Computing. the role of interval methods in scientific computing", ed. by Ramon E. Moore, pp. 99-107, Academic Press, Boston, MA

Coonen, Jerome 1984. Contributions to a Proposed Standard for Binary Floating-Point Arithmetic, PhD Thesis, Univ. of California, Berkeley

Dekker, T. J. 1971. A Floating-Point Technique for Extending the Available Precision, Numer. Math. 18(3), pp. 224-242

Demmel, James 1984. Underflow and the Reliability of Numerical Software, SIAM J. Sci. Stat. Comput. 5(4), pp. 887-919

Farnum, Charles 1988. Compiler Support for Floating-point Computation, Software-Practice and Experience, 18(7), pp. 701-709

Forsythe, G. E. and Moler, C. B. 1967. Computer Solution of Linear Algebraic Systems, Prentice-Hall, Englewood Cliffs, NJ

Goldberg, I. Bennett 1967. 27 Bits Are Not Enough for 8-Digit Accuracy, Comm. of the ACM. 10(2), pp 105-106

Goldberg, David 1990. Computer Arithmetic, in "Computer Architecture. A Quantitative Approach", by David Patterson and John L. Hennessy, Appendix A, Morgan Kaufmann, Los Altos, CA

Golub, Gene H. and Van Loan, Charles F. 1989. Matrix Computations, 2nd edition,The Johns Hopkins University Press, Baltimore Maryland

Graham, Ronald L. , Knuth, Donald E. and Patashnik, Oren. 1989. Concrete Mathematics, Addison-Wesley, Reading, MA, p. 162

Hewlett Packard 1982. HP-15C Advanced Functions Handbook

IEEE 1987. IEEE Standard 754-1985 for Binary Floating-point Arithmetic, IEEE, (1985). Reprinted in SIGPLAN 22(2) pp. 9-25

Kahan, W. 1972. A Survey Of Error Analysis, in Information Processing 71, Vol 2, pp. 1214 - 1239 (Ljubljana, Yugoslavia), North Holland, Amsterdam

Kahan, W. 1986. Calculating Area and Angle of a Needle-like Triangle, unpublished manuscript

Kahan, W. 1987. Branch Cuts for Complex Elementary Functions, in "The State of the Art in Numerical Analysis", ed. by M. J. D. Powell and A. Iserles (Univ of Birmingham, England), Chapter 7, Oxford University Press, New York

Kahan, W. 1988. Unpublished lectures given at Sun Microsystems, Mountain View, CA

Kahan, W. and Coonen, Jerome T. 1982. The Near Orthogonality of Syntax, Semantics, and Diagnostics in Numerical Programming Environments, in "The Relationship Between Numerical Computation And Programming Languages", ed. by J. K. Reid, pp. 103-115, North-Holland, Amsterdam

Kahan, W. and LeBlanc, E. 1985. Anomalies in the IBM Acrith Package, Proc. 7th IEEE Symposium on Computer Arithmetic (Urbana, Illinois), pp. 322-331

Kernighan, Brian W. and Ritchie, Dennis M. 1978. The C Programming Language, Prentice-Hall, Englewood Cliffs, NJ

Kirchner, R. và Kulisch, U. 1987. Arithmetic for Vector Processors, Proc. Hội nghị chuyên đề IEEE lần thứ 8 về số học máy tính (Como, Ý), trang. 256-269

Knuth, Donald E. , 1981. Nghệ thuật lập trình máy tính, Tập II, Phiên bản thứ hai, Addison-Wesley, Reading, MA

Kulisch, U. W. , và Miranker, W. L. 1986. Số học của máy tính kỹ thuật số. Một cách tiếp cận mới, Đánh giá SIAM 28(1), trang 1-36

Matula, D. W. và Kornerup, P. 1985. Số học hợp lý chính xác hữu hạn. Hệ thống số gạch chéo, IEEE Trans. trên may tinh. C-34(1), trang 3-18

Nelson, G. 1991. Systems Programming With Modula-3, Prentice-Hall, Englewood Cliffs, NJ

Reiser, John F. and Knuth, Donald E. 1975. Evading the Drift in Floating-point Addition, Information Processing Letters 3(3), pp 84-87

Sterbenz, Pat H. 1974. Floating-Point Computation, Prentice-Hall, Englewood Cliffs, NJ

Swartzlander, Earl E. and Alexopoulos, Aristides G. 1975. The Sign/Logarithm Number System, IEEE Trans. Comput. C-24(12), pp. 1238-1242

Walther, J. S. , 1971. A unified algorithm for elementary functions, Proceedings of the AFIP Spring Joint Computer Conf. 38, trang. 379-385

Theorem 14 and Theorem 8

This section contains two of the more technical proofs that were omitted from the text

Theorem 14

Let 0 < k < p, and set m =
Trăn trừ nhị phân
k + 1, and assume that floating-point operations are exactly rounded. Then (m
Trăn trừ nhị phân
x)(m
Trăn trừ nhị phân
xx) is exactly equal to x rounded to p - k significant digits. More precisely, x is rounded by taking the significand of x, imagining a radix point just left of the k least significant digits, and rounding to an integer.

Proof

The proof breaks up into two cases, depending on whether or not the computation of mx =
Trăn trừ nhị phân
kx + x has a carry-out or not. Assume there is no carry out. It is harmless to scale x so that it is an integer. Then the computation of mx = x +
Trăn trừ nhị phân
kx looks like this.
     n = n/2
144aa. aabb. bb
     n = n/2
145where x has been partitioned into two parts. The low order k digits are marked
     n = n/2
132 and the high order p - k digits are marked
     if (n==0) return u
3. To compute m
Trăn trừ nhị phân
x from mx involves rounding off the low order k digits (the ones marked with
     n = n/2
132) so(32) m
Trăn trừ nhị phân
x = mx - x mod(
Trăn trừ nhị phân
k) + r
Trăn trừ nhị phân
k
The value of r is 1 if
     n = n/2
149 is greater thanand 0 otherwise. More precisely(33) r = 1 if
     n = n/2
150 rounds to a + 1, r = 0 otherwise.
Next compute m
Trăn trừ nhị phân
x - x = mx - x mod(
Trăn trừ nhị phân
k) + r
Trăn trừ nhị phân
k - x =
Trăn trừ nhị phân
k(x + r) - x mod(
Trăn trừ nhị phân
k). The picture below shows the computation of m
Trăn trừ nhị phân
x - x rounded, that is, (m 
Trăn trừ nhị phân
x)x. The top line is
Trăn trừ nhị phân
k(x + r), where
     n = n/2
151 is the digit that results from adding
     n = n/2
152 to the lowest order digit
     n = n/2
132.
     n = n/2
154bb. bb
     n = n/2
155If
     n = n/2
149
Trăn trừ nhị phân
kx. Nếu
     n = n/2
149 > thì r = 1 và 1 bị trừ khỏi
     n = n/2
151 vì khoản vay, nên kết quả là
Trăn trừ nhị phân
kx. Finally consider the case
     n = n/2
149 =. If r = 0 then
     n = n/2
151 is even,
     n = n/2
162 is odd, and the difference is rounded up, giving
Trăn trừ nhị phân
kx. Similarly when r = 1,
     n = n/2
151 is odd,
     n = n/2
162 is even, the difference is rounded down, so again the difference is
Trăn trừ nhị phân
kx. To summarize(34) (m
Trăn trừ nhị phân
x)x =
Trăn trừ nhị phân
kx
Combining equations (32) and (34) gives (m
Trăn trừ nhị phân
x) - (m
Trăn trừ nhị phân
xx) = x - x mod(
Trăn trừ nhị phân
k) +
Trăn trừ nhị phân
·
Trăn trừ nhị phân
k. The result of performing this computation is
     n = n/2
165 bb. bb
     n = n/2
166The rule for computing r, equation (33), is the same as the rule for rounding
     n = n/2
167
     n = n/2
168 to p - k places. Thus computing mx - (mx - x) in floating-point arithmetic precision is exactly equal to rounding x to p - k places, in the case when x +
Trăn trừ nhị phân
kx does not carry out. When x +
Trăn trừ nhị phân
kx does carry out, then mx =
Trăn trừ nhị phân
kx + x looks like this.
     n = n/2
144aa. aabb. bb
     n = n/2
170Thus, m
Trăn trừ nhị phân
x = mx - x mod(
Trăn trừ nhị phân
k) + w
Trăn trừ nhị phân
k, where w = -Z if Z <
Trăn trừ nhị phân
/2, but the exact value of w is unimportant. Next, m
Trăn trừ nhị phân
x - x =
Trăn trừ nhị phân
kx - x mod(
Trăn trừ nhị phân
k) + w
Trăn trừ nhị phân
k. In a picture
     n = n/2
171 w
     n = n/2
172Rounding gives (m
Trăn trừ nhị phân
x)x =
Trăn trừ nhị phân
kx + w
Trăn trừ nhị phân
k - r
Trăn trừ nhị phân
k, where r = 1 if
     n = n/2
149 > or if
     n = n/2
149 =and b0 = 1. 32 Finally,(m
Trăn trừ nhị phân
x) - (m
Trăn trừ nhị phân
xx) = mx - x mod(
Trăn trừ nhị phân
k) + w
Trăn trừ nhị phân
k - (
Trăn trừ nhị phân
kx + w
Trăn trừ nhị phân
k - r
Trăn trừ nhị phân
k)
= x - x mod(
Trăn trừ nhị phân
k) + r
Trăn trừ nhị phân
k.
And once again, r = 1 exactly when rounding
     n = n/2
175 to p - k places involves rounding up. Thus Theorem 14 is proven in all cases. z

Theorem 8 (Kahan Summation Formula)

Suppose thatis computed using the following algorithm
     n = n/2
2
     n = n/2
08
     n = n/2
4
     n = n/2
5
     n = n/2
6
     n = n/2
7
     n = n/2
8
     n = n/2
9
Then the computed sum S is equal to S =
Trăn trừ nhị phân
xj (1 +
Trăn trừ nhị phân
j) + O(N
Trăn trừ nhị phân
2)
Trăn trừ nhị phân
. xj. , where .
Trăn trừ nhị phân
j.
Trăn trừ nhị phân
2
Trăn trừ nhị phân
.

Proof

First recall how the error estimate for the simple formula
Trăn trừ nhị phân
xi went. Introduce s1 = x1, si = (1 +
Trăn trừ nhị phân
i) (si - 1 + xi). Then the computed sum is sn, which is a sum of terms, each of which is an xi multiplied by an expression involving
Trăn trừ nhị phân
j's. The exact coefficient of x1 is (1 +
Trăn trừ nhị phân
2)(1 +
Trăn trừ nhị phân
3) . (1 +
Trăn trừ nhị phân
n), and so by renumbering, the coefficient of x2 must be (1 + 
Trăn trừ nhị phân
3)(1 +
Trăn trừ nhị phân
4) . (1 +
Trăn trừ nhị phân
n), v.v. The proof of Theorem 8 runs along exactly the same lines, only the coefficient of x1 is more complicated. In detail s0 = c0 = 0 andyk = xkck - 1 = (xk - ck - 1) (1 +
Trăn trừ nhị phân
k)sk = sk - 1
Trăn trừ nhị phân
Trăn trừ nhị phân
yk = (sk-1 + yk) (1 +
Trăn trừ nhị phân
k)ck = (sksk - 1)yk= [(sk - sk - 1) (1 +
Trăn trừ nhị phân
k) - yk] (1 +
Trăn trừ nhị phân
k)where all the Greek letters are bounded by
Trăn trừ nhị phân
. Although the coefficient of x1 in sk is the ultimate expression of interest, in turns out to be easier to compute the coefficient of x1 in sk - ck and ck. When k = 1,c1 = (s1(1 +
Trăn trừ nhị phân
1) - y1) (1 + d1)= y1((1 + s1) (1 +
Trăn trừ nhị phân
1) - 1) (1 + d1)= x1(s1 +
Trăn trừ nhị phân
1 + s1g1) (1 + d1) (1 + h1)s1 - c1 = x1[(1 + s1) - (s1 + g1 + s1g1) (1 + d1)](1 + h1)= x1[1 - g1 - s1d1 - s1g1 - d1g1 - s1g1d1](1 + h1)Calling the coefficients of x1 in these expressions Ck and Sk respectively, thenC1 = 2
Trăn trừ nhị phân
+ O(
Trăn trừ nhị phân
2)
S1 = +
Trăn trừ nhị phân
1 -
Trăn trừ nhị phân
1 + 4
Trăn trừ nhị phân
2 + O(
Trăn trừ nhị phân
3)
To get the general formula for Sk and Ck, expand the definitions of sk and ck, ignoring all terms involving xi with i > 1 to getsk = (sk - 1 + yk)(1 +
Trăn trừ nhị phân
k)= [sk - 1 + (xk - ck - 1) (1 +
Trăn trừ nhị phân
k)](1 +
Trăn trừ nhị phân
k)= [(sk - 1 - ck - 1) -
Trăn trừ nhị phân
kck - 1](1+
Trăn trừ nhị phân
k)ck = [{sk - sk - 1}(1 +
Trăn trừ nhị phân
k) - yk](1 +
Trăn trừ nhị phân
k)= [{((sk - 1 - ck - 1) -
Trăn trừ nhị phân
kck - 1)(1 +
Trăn trừ nhị phân
k) - sk - 1}(1 +
Trăn trừ nhị phân
k) + ck - 1(1 +
Trăn trừ nhị phân
k)](1 + 
Trăn trừ nhị phân
k)= [{(sk - 1 - ck - 1)
Trăn trừ nhị phân
k -
Trăn trừ nhị phân
kck-1(1 +
Trăn trừ nhị phân
k) - ck - 1}(1 +
Trăn trừ nhị phân
k) + ck - 1(1 +
Trăn trừ nhị phân
k)](1 +
Trăn trừ nhị phân
k)= [(sk - 1 - ck - 1)
Trăn trừ nhị phân
k(1 +
Trăn trừ nhị phân
k) - ck - 1(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k
Trăn trừ nhị phân
k))](1 +
Trăn trừ nhị phân
k),sk - ck = ((sk - 1 - ck - 1) -
Trăn trừ nhị phân
kck - 1) (1 +
Trăn trừ nhị phân
k)- [(sk - 1 - ck - 1)
Trăn trừ nhị phân
k(1 +
Trăn trừ nhị phân
k) - ck - 1(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k
Trăn trừ nhị phân
k)](1 +
Trăn trừ nhị phân
k)= (sk- 1 - ck - 1)((1 +
Trăn trừ nhị phân
k) -
Trăn trừ nhị phân
k(1 +
Trăn trừ nhị phân
k)(1 +
Trăn trừ nhị phân
k))+ ck - 1(-
Trăn trừ nhị phân
k(1 +
Trăn trừ nhị phân
k) + (
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k
Trăn trừ nhị phân
k)) (1 +
Trăn trừ nhị phân
k))= (s- 1 - ck - 1) (1 -
Trăn trừ nhị phân
k(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k
Trăn trừ nhị phân
k))+ ck - 1 - [
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k
Trăn trừ nhị phân
k) + (
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k(
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k +
Trăn trừ nhị phân
k
Trăn trừ nhị phân
k))
Trăn trừ nhị phân
k]Since Sk and Ck are only being computed up to order
Trăn trừ nhị phân
2, these formulas can be simplified toCk= (
Trăn trừ nhị phân
k + O(
Trăn trừ nhị phân
2))Sk - 1 + (-
Trăn trừ nhị phân
k + O(
Trăn trừ nhị phân
2))Ck - 1Sk= ((1 + 2
Trăn trừ nhị phân
2 + O(
Trăn trừ nhị phân
3))Sk - 1 + (2
Trăn trừ nhị phân
+
Trăn trừ nhị phân
(
Trăn trừ nhị phân
2))Ck - 1Using these formulas givesC2 =
Trăn trừ nhị phân
2 + O(
Trăn trừ nhị phân
2)
S2 = 1 +
Trăn trừ nhị phân
1 -
Trăn trừ nhị phân
1 + 10
Trăn trừ nhị phân
2 + O(
Trăn trừ nhị phân
3)
and in general it is easy to check by induction thatCk =
Trăn trừ nhị phân
k + O(
Trăn trừ nhị phân
2)
Sk = 1 +
Trăn trừ nhị phân
1 -
Trăn trừ nhị phân
1 + (4k+2)
Trăn trừ nhị phân
2 + O(
Trăn trừ nhị phân
3)

Finally, what is wanted is the coefficient of x1 in sk. To get this value, let xn + 1 = 0, let all the Greek letters with subscripts of n + 1 equal 0, and compute sn + 1. Then sn + 1 = sn - cn, and the coefficient of x1 in sn is less than the coefficient in sn + 1, which is Sn = 1 +

Trăn trừ nhị phân
1 -
Trăn trừ nhị phân
1 + (4n + 2)
Trăn trừ nhị phân
2 = (1 + 2
Trăn trừ nhị phân
+
Trăn trừ nhị phân
(n
Trăn trừ nhị phân
2)). z

Differences Among IEEE 754 Implementations


Note – This section is not part of the published paper. It has been added to clarify certain points and correct possible misconceptions about the IEEE standard that the reader might infer from the paper. Tài liệu này không phải do David Goldberg viết, nhưng nó xuất hiện ở đây với sự cho phép của ông

Bài báo trước đã chỉ ra rằng số học dấu phẩy động phải được triển khai cẩn thận, vì các lập trình viên có thể phụ thuộc vào các thuộc tính của nó để đảm bảo tính đúng đắn và chính xác của chương trình của họ. Đặc biệt, tiêu chuẩn IEEE yêu cầu triển khai cẩn thận và chỉ có thể viết các chương trình hữu ích hoạt động chính xác và mang lại kết quả chính xác trên các hệ thống tuân thủ tiêu chuẩn. Người đọc có thể muốn kết luận rằng các chương trình như vậy sẽ có thể di chuyển được đến tất cả các hệ thống của IEEE. Thật vậy, phần mềm di động sẽ dễ viết hơn nếu nhận xét "Khi một chương trình được di chuyển giữa hai máy và cả hai đều hỗ trợ số học IEEE, thì nếu bất kỳ kết quả trung gian nào khác đi, thì đó phải là do lỗi phần mềm, không phải do sự khác biệt về số học," là

Thật không may, tiêu chuẩn IEEE không đảm bảo rằng cùng một chương trình sẽ mang lại kết quả giống hệt nhau trên tất cả các hệ thống phù hợp. Hầu hết các chương trình sẽ thực sự tạo ra các kết quả khác nhau trên các hệ thống khác nhau vì nhiều lý do. Thứ nhất, hầu hết các chương trình liên quan đến việc chuyển đổi số giữa định dạng thập phân và nhị phân và tiêu chuẩn IEEE không chỉ định đầy đủ độ chính xác mà các chuyển đổi đó phải được thực hiện. Mặt khác, nhiều chương trình sử dụng các chức năng cơ bản do thư viện hệ thống cung cấp và tiêu chuẩn hoàn toàn không chỉ định các chức năng này. Tất nhiên, hầu hết các lập trình viên đều biết rằng các tính năng này nằm ngoài phạm vi của tiêu chuẩn IEEE.

Nhiều lập trình viên có thể không nhận ra rằng ngay cả một chương trình chỉ sử dụng các định dạng số và hoạt động theo tiêu chuẩn IEEE có thể tính toán các kết quả khác nhau trên các hệ thống khác nhau. Trên thực tế, các tác giả của tiêu chuẩn dự định cho phép các triển khai khác nhau thu được các kết quả khác nhau. Ý định của họ thể hiện rõ trong định nghĩa về thuật ngữ đích trong tiêu chuẩn IEEE 754. "Đích có thể được chỉ định rõ ràng bởi người dùng hoặc được cung cấp ngầm bởi hệ thống (ví dụ: kết quả trung gian trong các biểu thức con hoặc đối số cho các thủ tục). Một số ngôn ngữ đặt kết quả tính toán trung gian ở đích ngoài tầm kiểm soát của người dùng. Tuy nhiên, tiêu chuẩn này xác định kết quả của một hoạt động theo định dạng của đích đó và các giá trị của toán hạng. " (IEEE 754-1985, p. 7) Nói cách khác, tiêu chuẩn IEEE yêu cầu mỗi kết quả được làm tròn chính xác theo độ chính xác của đích mà nó sẽ được đặt vào, nhưng tiêu chuẩn không yêu cầu độ chính xác của đích đó được xác định bởi chương trình của người dùng. Do đó, các hệ thống khác nhau có thể đưa kết quả của chúng đến các đích với độ chính xác khác nhau, khiến cùng một chương trình tạo ra các kết quả khác nhau (đôi khi rất đáng kể), mặc dù các hệ thống đó đều tuân theo tiêu chuẩn

Một số ví dụ trong bài báo trước phụ thuộc vào một số kiến ​​thức về cách làm tròn số học dấu phẩy động. Để dựa vào các ví dụ như thế này, một lập trình viên phải có khả năng dự đoán cách một chương trình sẽ được diễn giải và đặc biệt, trên hệ thống IEEE, độ chính xác của đích đến của mỗi phép toán số học có thể là bao nhiêu. Than ôi, lỗ hổng trong định nghĩa đích của tiêu chuẩn IEEE làm suy yếu khả năng của lập trình viên để biết chương trình sẽ được diễn giải như thế nào. Do đó, một số ví dụ nêu trên, khi được triển khai dưới dạng các chương trình di động rõ ràng bằng ngôn ngữ cấp cao, có thể không hoạt động chính xác trên các hệ thống IEEE thường cung cấp kết quả tới đích với độ chính xác khác với mong đợi của lập trình viên. Các ví dụ khác có thể hoạt động, nhưng việc chứng minh rằng chúng hoạt động có thể nằm ngoài khả năng của một lập trình viên bình thường

Trong phần này, chúng tôi phân loại các triển khai số học IEEE 754 hiện có dựa trên độ chính xác của các định dạng đích mà chúng thường sử dụng. Sau đó, chúng tôi xem xét một số ví dụ từ bài báo để chỉ ra rằng việc cung cấp kết quả với độ chính xác rộng hơn mong đợi của chương trình có thể khiến chương trình tính toán kết quả sai mặc dù nó có thể đúng khi sử dụng độ chính xác dự kiến. Chúng tôi cũng xem lại một trong những bằng chứng trong bài báo để minh họa nỗ lực trí tuệ cần thiết để đối phó với độ chính xác ngoài dự kiến ​​ngay cả khi nó không làm mất hiệu lực chương trình của chúng tôi. Những ví dụ này cho thấy rằng bất chấp tất cả những gì tiêu chuẩn IEEE quy định, sự khác biệt mà tiêu chuẩn này cho phép giữa các triển khai khác nhau có thể ngăn chúng ta viết phần mềm số hiệu quả, di động có hành vi mà chúng ta có thể dự đoán chính xác. Sau đó, để phát triển phần mềm như vậy, trước tiên chúng ta phải tạo ra các ngôn ngữ lập trình và môi trường hạn chế tính biến đổi mà tiêu chuẩn IEEE cho phép và cho phép các lập trình viên thể hiện ngữ nghĩa dấu phẩy động mà chương trình của họ phụ thuộc vào.

Triển khai IEEE 754 hiện tại

Việc triển khai số học IEEE 754 hiện tại có thể được chia thành hai nhóm được phân biệt theo mức độ chúng hỗ trợ các định dạng dấu phẩy động khác nhau trong phần cứng. Các hệ thống dựa trên mở rộng, được minh họa bởi dòng bộ xử lý Intel x86, cung cấp hỗ trợ đầy đủ cho định dạng độ chính xác kép mở rộng nhưng chỉ hỗ trợ một phần cho độ chính xác đơn và kép. chúng cung cấp các hướng dẫn để tải hoặc lưu trữ dữ liệu ở độ chính xác đơn và kép, chuyển đổi dữ liệu nhanh chóng sang hoặc từ định dạng kép mở rộng và chúng cung cấp các chế độ đặc biệt (không phải mặc định) trong đó kết quả của các phép toán số học được làm tròn thành đơn . (Bộ xử lý sê-ri Motorola 68000 làm tròn kết quả cho cả độ chính xác và phạm vi của định dạng đơn hoặc kép trong các chế độ này. Intel x86 và bộ xử lý tương thích làm tròn kết quả thành độ chính xác của định dạng đơn hoặc kép nhưng vẫn giữ nguyên phạm vi như định dạng kép mở rộng. ) Các hệ thống đơn/kép, bao gồm hầu hết các bộ xử lý RISC, cung cấp hỗ trợ đầy đủ cho các định dạng chính xác đơn và kép nhưng không hỗ trợ định dạng chính xác kép mở rộng tuân thủ IEEE. (Kiến trúc IBM POWER chỉ cung cấp hỗ trợ một phần cho độ chính xác đơn, nhưng với mục đích của phần này, chúng tôi phân loại nó là hệ thống đơn/kép. )

Để xem cách tính toán có thể hoạt động khác nhau trên hệ thống dựa trên mở rộng so với trên hệ thống đơn/kép, hãy xem xét phiên bản C của ví dụ từ phần Các khía cạnh hệ thống

     n = n/2
30
     n = n/2
31
     n = n/2
32
     n = n/2
33
     n = n/2
34
     n = n/2
35
     n = n/2
9

Ở đây các hằng số 3. 0 and 7. 0 được hiểu là số dấu phẩy động có độ chính xác kép và biểu thức 3. 0/7. 0 kế thừa kiểu dữ liệu

     n = n/2
176. Trên một hệ thống đơn/kép, biểu thức sẽ được đánh giá với độ chính xác kép vì đó là định dạng hiệu quả nhất để sử dụng. Do đó,
     n = n/2
053 sẽ được gán giá trị 3. 0/7. 0 được làm tròn chính xác để tăng gấp đôi độ chính xác. Trong dòng tiếp theo, biểu thức 3. 0/7. 0 một lần nữa sẽ được đánh giá với độ chính xác kép và tất nhiên kết quả sẽ bằng với giá trị vừa được gán cho
     n = n/2
053, vì vậy chương trình sẽ in "Bằng" như mong đợi

Trên một hệ thống dựa trên mở rộng, mặc dù biểu thức 3. 0/7. 0 có loại

     n = n/2
176, thương số sẽ được tính trong một thanh ghi ở định dạng kép mở rộng và do đó ở chế độ mặc định, nó sẽ được làm tròn thành độ chính xác kép mở rộng. Tuy nhiên, khi giá trị kết quả được gán cho biến
     n = n/2
053, thì nó có thể được lưu trữ trong bộ nhớ và vì
     n = n/2
053 được khai báo là
     n = n/2
176, nên giá trị sẽ được làm tròn thành độ chính xác gấp đôi. Trong dòng tiếp theo, biểu thức 3. 0/7. 0 may again be evaluated in extended precision yielding a result that differs from the double precision value stored in
     n = n/2
053, causing the program to print "Not equal". Tất nhiên, các kết quả khác cũng có thể xảy ra. trình biên dịch có thể quyết định lưu trữ và do đó làm tròn giá trị của biểu thức 3. 0/7. 0 ở dòng thứ hai trước khi so sánh nó với
     n = n/2
053 hoặc nó có thể giữ
     n = n/2
053 trong sổ đăng ký với độ chính xác mở rộng mà không cần lưu trữ. Trình biên dịch tối ưu hóa có thể đánh giá biểu thức 3. 0/7. 0 tại thời điểm biên dịch, có lẽ ở độ chính xác kép hoặc có thể ở độ chính xác kép mở rộng. (Với một trình biên dịch x86, chương trình sẽ in "Bằng" khi được biên dịch với tối ưu hóa và "Không bằng" khi được biên dịch để gỡ lỗi. ) Cuối cùng, một số trình biên dịch cho các hệ thống dựa trên mở rộng tự động thay đổi chế độ làm tròn độ chính xác để khiến các hoạt động tạo ra kết quả trong các thanh ghi để làm tròn các kết quả đó thành độ chính xác đơn hoặc kép, mặc dù có thể với phạm vi rộng hơn. Do đó, trên các hệ thống này, chúng tôi không thể dự đoán hành vi của chương trình chỉ bằng cách đọc mã nguồn của nó và áp dụng hiểu biết cơ bản về số học IEEE 754. Chúng tôi cũng không thể buộc tội phần cứng hoặc trình biên dịch không cung cấp môi trường tuân thủ IEEE 754;

Cạm bẫy trong tính toán trên các hệ thống dựa trên mở rộng

Sự khôn ngoan thông thường cho rằng các hệ thống dựa trên mở rộng phải tạo ra kết quả ít nhất là chính xác, nếu không muốn nói là chính xác hơn kết quả được cung cấp trên hệ thống đơn/kép, vì hệ thống trước luôn cung cấp độ chính xác ít nhất bằng và thường cao hơn hệ thống sau. Các ví dụ tầm thường như chương trình C ở trên cũng như các chương trình phức tạp hơn dựa trên các ví dụ được thảo luận bên dưới cho thấy sự khôn ngoan này tốt nhất là ngây thơ. một số chương trình rõ ràng là di động, thực sự là di động trên các hệ thống đơn/kép, cung cấp kết quả không chính xác trên các hệ thống dựa trên mở rộng chính xác do trình biên dịch và phần cứng đôi khi cung cấp độ chính xác cao hơn mong đợi của chương trình

Các ngôn ngữ lập trình hiện tại khiến chương trình khó xác định độ chính xác mà nó mong đợi. Như phần Ngôn ngữ và Trình biên dịch đã đề cập, nhiều ngôn ngữ lập trình không chỉ định rằng mỗi lần xuất hiện của một biểu thức như

     n = n/2
035 trong cùng một ngữ cảnh sẽ đánh giá cùng một giá trị. Some languages, such as Ada, were influenced in this respect by variations among different arithmetics prior to the IEEE standard. Gần đây hơn, các ngôn ngữ như ANSI C đã bị ảnh hưởng bởi các hệ thống dựa trên mở rộng phù hợp với tiêu chuẩn. Trên thực tế, tiêu chuẩn ANSI C rõ ràng cho phép trình biên dịch đánh giá biểu thức dấu phẩy động với độ chính xác rộng hơn mức thường được kết hợp với loại của nó. Do đó, giá trị của biểu thức
     n = n/2
035 có thể khác nhau tùy thuộc vào nhiều yếu tố. liệu biểu thức được gán ngay cho một biến hay xuất hiện dưới dạng biểu thức con trong biểu thức lớn hơn;

Các tiêu chuẩn ngôn ngữ không hoàn toàn đổ lỗi cho sự thất thường của việc đánh giá biểu thức. Các hệ thống dựa trên mở rộng chạy hiệu quả nhất khi các biểu thức được đánh giá trong các thanh ghi độ chính xác mở rộng bất cứ khi nào có thể, tuy nhiên các giá trị phải được lưu trữ được lưu trữ ở độ chính xác hẹp nhất được yêu cầu. Hạn chế một ngôn ngữ để yêu cầu đánh giá

     n = n/2
035 với cùng một giá trị ở mọi nơi sẽ áp dụng hình phạt hiệu suất đối với các hệ thống đó. Thật không may, việc cho phép các hệ thống đó đánh giá
     n = n/2
035 khác nhau trong các ngữ cảnh tương đương về mặt cú pháp sẽ tự áp đặt hình phạt đối với các lập trình viên của phần mềm số chính xác bằng cách ngăn họ dựa vào cú pháp của chương trình để thể hiện ngữ nghĩa dự định của họ

Các chương trình thực tế có phụ thuộc vào giả định rằng một biểu thức đã cho luôn đánh giá cùng một giá trị không?

_______537____538____539____540
     n = n/2
41
     n = n/2
42____543
     n = n/2
44

Trên hệ thống dựa trên mở rộng, trình biên dịch có thể đánh giá biểu thức

     n = n/2
190 
     n = n/2
044 
     n = n/2
006 ở dòng thứ ba với độ chính xác mở rộng và so sánh kết quả với
     n = n/2
190. Tuy nhiên, khi biểu thức tương tự được chuyển đến hàm nhật ký ở dòng thứ sáu, trình biên dịch có thể lưu trữ giá trị của nó trong bộ nhớ, làm tròn nó thành độ chính xác đơn. Do đó, nếu
     n = n/2
006 không nhỏ đến mức
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 làm tròn thành
     n = n/2
190 với độ chính xác mở rộng nhưng đủ nhỏ để
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 làm tròn thành
     n = n/2
190 với độ chính xác đơn, thì giá trị trả về bởi
     n = n/2
03 sẽ bằng 0 thay vì
     n = n/2
006 và lỗi tương đối sẽ . Tương tự, giả sử phần còn lại của biểu thức ở dòng thứ sáu, bao gồm cả sự xuất hiện lại của biểu thức con
     n = n/2
190 
     n = n/2
044 
     n = n/2
006, được đánh giá theo độ chính xác mở rộng. Trong trường hợp đó, nếu
     n = n/2
006 nhỏ nhưng không đủ nhỏ để
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 làm tròn thành
     n = n/2
190 với độ chính xác duy nhất, thì giá trị do
     n = n/2
03 trả về có thể vượt quá giá trị chính xác gần như bằng với
     n = n/2
006 và một lần nữa, lỗi tương đối có thể đạt đến một. Ví dụ cụ thể, lấy
     n = n/2
006 là 2-24 + 2-47, vì vậy,
     n = n/2
006 là số chính xác đơn lẻ nhỏ nhất sao cho
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 làm tròn thành số lớn hơn tiếp theo, 1 + 2-23. Khi đó
     n = n/2
20 
     n = n/2
044 
     n = n/2
22 xấp xỉ 2-23. Bởi vì mẫu số trong biểu thức ở dòng thứ sáu được đánh giá ở độ chính xác mở rộng, nó được tính toán chính xác và mang lại giá trị
     n = n/2
006, do đó,
     n = n/2
03 trả về khoảng 2-23, lớn gần gấp đôi giá trị chính xác. (Điều này thực sự xảy ra với ít nhất một trình biên dịch. When the preceding code is compiled by the Sun WorkShop Compilers 4. 2. 1 Trình biên dịch Fortran 77 dành cho các hệ thống x86 sử dụng cờ tối ưu hóa
     n = n/2
25, mã được tạo sẽ tính toán
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 chính xác như mô tả. Kết quả là, hàm trả về 0 cho
     n = n/2
29 và
     n = n/2
30 cho
     n = n/2
31. )
Trăn trừ nhị phân
. Similarly, suppose the rest of the expression in the sixth line, including the reoccurrence of the subexpression
     n = n/2
190 
     n = n/2
044 
     n = n/2
006, is evaluated in extended precision. In that case, if
     n = n/2
006 is small but not quite small enough that
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 rounds to
     n = n/2
190 in single precision, then the value returned by
     n = n/2
03 can exceed the correct value by nearly as much as
     n = n/2
006, and again the relative error can approach one. For a concrete example, take
     n = n/2
006 to be 2-24 + 2-47, so
     n = n/2
006 is the smallest single precision number such that
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 rounds up to the next larger number, 1 + 2-23. Then
     n = n/2
20 
     n = n/2
044 
     n = n/2
22 is approximately 2-23. Because the denominator in the expression in the sixth line is evaluated in extended precision, it is computed exactly and delivers
     n = n/2
006, so
     n = n/2
03 returns approximately 2-23, which is nearly twice as large as the exact value. (This actually happens with at least one compiler. When the preceding code is compiled by the Sun WorkShop Compilers 4.2.1 Fortran 77 compiler for x86 systems using the
     n = n/2
25 optimization flag, the generated code computes
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 exactly as described. As a result, the function delivers zero for
     n = n/2
29 and
     n = n/2
30 for
     n = n/2
31.)

Để thuật toán của Định lý 4 hoạt động chính xác, biểu thức

     n = n/2
190 
     n = n/2
044 
     n = n/2
006 phải được đánh giá theo cùng một cách mỗi khi nó xuất hiện; . Tất nhiên, vì
     n = n/2
38 là một hàm nội tại chung trong Fortran, trình biên dịch có thể đánh giá biểu thức
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 với độ chính xác mở rộng xuyên suốt, tính logarit của nó với cùng độ chính xác, nhưng rõ ràng chúng ta không thể cho rằng trình biên dịch sẽ làm như vậy. (Người ta cũng có thể tưởng tượng một ví dụ tương tự liên quan đến hàm do người dùng định nghĩa. Trong trường hợp đó, một trình biên dịch vẫn có thể giữ đối số ở độ chính xác mở rộng mặc dù hàm trả về một kết quả chính xác duy nhất, nhưng rất ít nếu có bất kỳ trình biên dịch Fortran hiện có nào làm được điều này. ) Do đó, chúng tôi có thể cố gắng đảm bảo rằng
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 được đánh giá một cách nhất quán bằng cách gán nó cho một biến. Thật không may, nếu chúng ta khai báo biến
     n = n/2
45 đó, thì chúng ta vẫn có thể bị trình biên dịch làm hỏng khi thay thế một giá trị được lưu trong sổ đăng ký với độ chính xác mở rộng cho một lần xuất hiện của biến và một giá trị được lưu trong bộ nhớ với độ chính xác duy nhất cho một lần xuất hiện khác. Thay vào đó, chúng ta cần khai báo biến có kiểu tương ứng với định dạng chính xác mở rộng. FORTRAN 77 tiêu chuẩn không cung cấp cách để thực hiện việc này và trong khi Fortran 95 cung cấp cơ chế
     n = n/2
46 để mô tả các định dạng khác nhau, nó không yêu cầu rõ ràng việc triển khai đánh giá các biểu thức với độ chính xác mở rộng để cho phép các biến được khai báo với độ chính xác đó. Nói tóm lại, không có cách di động nào để viết chương trình này bằng Fortran tiêu chuẩn được đảm bảo để ngăn biểu thức
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 bị đánh giá theo cách làm mất hiệu lực bằng chứng của chúng tôi

Có những ví dụ khác có thể gặp trục trặc trên các hệ thống dựa trên mở rộng ngay cả khi mỗi biểu thức con được lưu trữ và do đó được làm tròn với cùng độ chính xác. Nguyên nhân là làm tròn kép. Ở chế độ chính xác mặc định, một hệ thống dựa trên mở rộng ban đầu sẽ làm tròn từng kết quả thành độ chính xác kép mở rộng. Nếu kết quả đó sau đó được lưu trữ với độ chính xác gấp đôi, nó sẽ được làm tròn lại. Sự kết hợp của hai cách làm tròn này có thể mang lại một giá trị khác với giá trị có được bằng cách làm tròn kết quả đầu tiên một cách chính xác để tăng gấp đôi độ chính xác. Điều này có thể xảy ra khi kết quả được làm tròn thành độ chính xác kép mở rộng là "trường hợp giữa chừng", tôi. e. , nó nằm chính xác ở giữa hai số chính xác kép, do đó, lần làm tròn thứ hai được xác định theo quy tắc làm tròn số chẵn. Nếu lần làm tròn thứ hai này quay cùng hướng với lần đầu tiên, thì sai số làm tròn ròng sẽ vượt quá nửa đơn vị ở vị trí cuối cùng. (Tuy nhiên, xin lưu ý rằng việc làm tròn hai lần chỉ ảnh hưởng đến các phép tính có độ chính xác kép. Người ta có thể chứng minh rằng tổng, hiệu, tích hoặc thương của hai số p-bit hoặc căn bậc hai của một số p-bit, trước tiên được làm tròn thành q bit và sau đó thành p bit cho cùng một giá trị như thể kết quả là . Do đó, độ chính xác kép mở rộng đủ rộng để các tính toán chính xác đơn lẻ không bị làm tròn kép. )

Trăn trừ nhị phân
2p + 2. Thus, extended double precision is wide enough that single precision computations don't suffer double-rounding.)

Một số thuật toán phụ thuộc vào làm tròn chính xác có thể thất bại với làm tròn kép. Trên thực tế, ngay cả một số thuật toán không yêu cầu làm tròn chính xác và hoạt động chính xác trên nhiều loại máy không tuân thủ IEEE 754 cũng có thể thất bại khi làm tròn hai lần. Hữu ích nhất trong số này là các thuật toán di động để thực hiện nhiều phép tính số học chính xác mô phỏng được đề cập trong phần Phép toán được làm tròn chính xác. Ví dụ: quy trình được mô tả trong Định lý 6 để tách một số dấu phẩy động thành các phần cao và thấp không hoạt động chính xác trong số học làm tròn kép. cố gắng chia số chính xác kép 252 + 3 × 226 - 1 thành hai phần, mỗi phần có tối đa 26 bit. Khi mỗi thao tác được làm tròn chính xác thành độ chính xác gấp đôi, phần có thứ tự cao là 252 + 227 và phần có thứ tự thấp là 226 - 1, nhưng khi mỗi thao tác được làm tròn trước thành độ chính xác gấp đôi mở rộng rồi đến độ chính xác gấp đôi, quy trình sẽ tạo ra giá trị cao . Số thứ hai chiếm 27 bit, vì vậy bình phương của nó không thể được tính chính xác với độ chính xác gấp đôi. Tất nhiên, vẫn có thể tính bình phương của số này với độ chính xác gấp đôi mở rộng, nhưng thuật toán kết quả sẽ không còn khả dụng cho các hệ thống đơn/kép. Ngoài ra, các bước sau trong thuật toán nhân chính xác bội giả định rằng tất cả các tích từng phần đã được tính toán với độ chính xác kép. Xử lý chính xác hỗn hợp các biến kép kép và mở rộng sẽ khiến việc triển khai trở nên đắt đỏ hơn đáng kể

Tương tự như vậy, các thuật toán di động để thêm nhiều số chính xác được biểu thị dưới dạng mảng các số chính xác kép có thể không thành công trong phép tính làm tròn hai lần. Các thuật toán này thường dựa trên một kỹ thuật tương tự như công thức tính tổng của Kahan. Như lời giải thích không chính thức về công thức tính tổng được đưa ra trên Errors In Summation gợi ý, nếu

     n = n/2
50 và
     n = n/2
042 là các biến dấu phẩy động với.
     n = n/2
50.
Trăn trừ nhị phân
.
     n = n/2
042. và chúng tôi tính toán.

     n = n/2
45
     n = n/2
46

sau đó trong hầu hết các phép tính,

     n = n/2
54 khôi phục chính xác lỗi làm tròn đã xảy ra trong máy tính
     n = n/2
55. Tuy nhiên, kỹ thuật này không hoạt động trong số học làm tròn hai lần. nếu
     n = n/2
50 = 252 + 1 và
     n = n/2
042 = 1/2 - 2-54, thì
     n = n/2
50 
     n = n/2
044 
     n = n/2
042 làm tròn đầu tiên thành 252 + 3/2 với độ chính xác kép mở rộng và giá trị này làm tròn thành 252 + 2 với độ chính xác kép theo các mối quan hệ . Một lần nữa, ở đây, có thể khôi phục lỗi làm tròn bằng cách tính tổng ở độ chính xác kép mở rộng, nhưng sau đó chương trình sẽ phải thực hiện thêm công việc để giảm kết quả cuối cùng trở lại độ chính xác gấp đôi và việc làm tròn kép có thể ảnh hưởng đến quá trình này, . Vì lý do này, mặc dù các chương trình di động để mô phỏng nhiều phép tính chính xác bằng các phương pháp này hoạt động chính xác và hiệu quả trên nhiều loại máy, nhưng chúng không hoạt động như quảng cáo trên các hệ thống dựa trên mở rộng

Cuối cùng, một số thuật toán thoạt nhìn có vẻ phụ thuộc vào việc làm tròn chính xác trên thực tế có thể hoạt động chính xác với phép làm tròn hai lần. Trong những trường hợp này, chi phí đối phó với việc làm tròn hai lần không nằm ở việc triển khai mà nằm ở việc xác minh rằng thuật toán có hoạt động như quảng cáo hay không. Để minh họa, ta chứng minh biến thể sau của Định lý 7

Định lý 7'

Nếu m và n là các số nguyên có thể biểu diễn theo độ chính xác kép của IEEE 754 với. m. < 252 và n có dạng đặc biệt n = 2i + 2j, sau đó (mn)
Trăn trừ nhị phân
n = m, miễn là cả hai phép toán dấu phẩy động đều được làm tròn chính xác thành gấp đôi độ chính xác hoặc làm tròn trước thành gấp đôi mở rộng .

Proof

Giả sử không mất mát rằng m > 0. Đặt q = mn. Chia tỷ lệ theo lũy thừa hai, chúng ta có thể xem xét một cài đặt tương đương trong đó 252

Trăn trừ nhị phân
m < 253 và tương tự như vậy đối với q, sao cho cả m và q đều là số nguyên có bit ít quan trọng nhất chiếm vị trí đơn vị (i. e. , ulp(m) = ulp(q) = 1). Trước khi chia tỷ lệ, chúng tôi giả sử m < 252, vì vậy sau khi chia tỷ lệ, m là một số nguyên chẵn. Also, because the scaled values of m and q satisfy m/2 < q < 2m, the corresponding value of n must have one of two forms depending on which of m or q is larger. nếu q 

Gọi e là lỗi làm tròn trong phép tính q, sao cho q = m/n + e và giá trị được tính toán q

Trăn trừ nhị phân
n sẽ là giá trị được làm tròn (một hoặc hai lần) . Trước tiên, hãy xem xét trường hợp trong đó mỗi thao tác dấu phẩy động được làm tròn chính xác để tăng gấp đôi độ chính xác. Trong trường hợp này,. e. < 1/2. Nếu n có dạng 1/2 + 2-(k + 1) thì ne = nq - m là bội số nguyên của 2-(k + 1) và. ne. < 1/4 + 2-(k + 2). Điều này ngụ ý rằng. ne.
Trăn trừ nhị phân
1/4. Nhớ lại rằng hiệu giữa m và số có thể biểu thị lớn hơn tiếp theo là 1 và hiệu giữa m và số có thể biểu thị nhỏ hơn tiếp theo là 1 nếu m > 252 hoặc 1/2 nếu m = 252. Vì vậy, như. ne.
Trăn trừ nhị phân
1/4, m + ne sẽ làm tròn thành m. (Ngay cả khi m = 252 và ne = -1/4, tích sẽ làm tròn thành m theo quy tắc làm tròn số chẵn. ) Tương tự, nếu n có dạng 1 + 2-k, thì ne là bội số nguyên của 2-k và. ne.
Trăn trừ nhị phân
1/2. Chúng ta không thể có m = 252 trong trường hợp này vì m hoàn toàn lớn hơn q, vì vậy m khác với các lân cận có thể biểu diễn gần nhất của nó bằng ±1. Vì vậy, như. ne.
Trăn trừ nhị phân
1/2, lại m + ne sẽ làm tròn thành m. (Thậm chí nếu. ne. = 1/2, tích sẽ làm tròn thành m theo quy tắc làm tròn số chẵn vì m chẵn. ) Điều này hoàn thành bằng chứng cho số học được làm tròn chính xác.

Trong số học làm tròn hai lần, vẫn có thể xảy ra trường hợp q là thương được làm tròn chính xác (mặc dù nó thực sự đã được làm tròn hai lần), vì vậy. e. < 1/2 như trên. Trong trường hợp này, chúng tôi có thể kháng cáo các lập luận của đoạn trước miễn là chúng tôi xem xét thực tế rằng q

Trăn trừ nhị phân
n sẽ được làm tròn hai lần. Để giải thích điều này, hãy lưu ý rằng tiêu chuẩn IEEE yêu cầu định dạng kép mở rộng mang ít nhất 64 bit có nghĩa, sao cho các số m ± 1/2 và m ± 1/4 có thể biểu diễn chính xác ở độ chính xác kép mở rộng. Do đó, nếu n có dạng 1/2 + 2-(k + 1), sao cho. ne.
Trăn trừ nhị phân
1/4, sau đó làm tròn m + ne thành độ chính xác kép mở rộng phải tạo ra kết quả khác với m nhiều nhất là 1/4 và như đã lưu ý ở trên, giá trị này sẽ làm tròn thành m gấp đôi . Tương tự, nếu n có dạng 1 + 2-k, sao cho. ne.
Trăn trừ nhị phân
1/2, sau đó làm tròn m + ne thành độ chính xác kép mở rộng phải tạo ra kết quả khác với m nhiều nhất là 1/2 và giá trị này sẽ làm tròn thành m ở độ chính xác kép. (Nhắc lại m > 252 trong trường hợp này. )

Cuối cùng, chúng ta còn lại để xem xét các trường hợp trong đó q không phải là thương được làm tròn chính xác do làm tròn hai lần. Trong những trường hợp này, chúng ta có. e. < 1/2 + 2-(d + 1) trong trường hợp xấu nhất, trong đó d là số bit thừa ở định dạng kép mở rộng. (Tất cả các hệ thống dựa trên mở rộng hiện có đều hỗ trợ định dạng kép mở rộng với chính xác 64 bit có nghĩa; đối với định dạng này, d = 64 - 53 = 11. ) Vì làm tròn hai lần chỉ tạo ra kết quả làm tròn không chính xác khi làm tròn thứ hai được xác định theo quy tắc làm tròn về số chẵn, nên q phải là một số nguyên chẵn. Do đó, nếu n có dạng 1/2 + 2-(k + 1), thì ne = nq - m là bội số nguyên của 2-k và

ne. < (1/2 + 2-(k + 1))(1/2 + 2-(d + 1)) = 1/4 + 2-(k + 2) + 2-(d + 2) + 2-

Nếu k

Trăn trừ nhị phân
d, điều này ngụ ý. ne.
Trăn trừ nhị phân
1/4. Nếu k > d, ta có. ne.
Trăn trừ nhị phân
1/4 + 2-(d + 2). Trong cả hai trường hợp, lần làm tròn đầu tiên của tích sẽ mang lại kết quả khác với m nhiều nhất là 1/4 và theo các đối số trước đó, lần làm tròn thứ hai sẽ làm tròn thành m. Tương tự, nếu n có dạng 1 + 2-k thì ne là bội số nguyên của 2-(k - 1) và

ne. < 1/2 + 2-(k + 1) + 2-(d + 1) + 2-(k + d + 1)

Nếu k

Trăn trừ nhị phân
d, điều này ngụ ý. ne.
Trăn trừ nhị phân
1/2. Nếu k > d, ta có. ne.
Trăn trừ nhị phân
1/2 + 2-(d + 1). Trong cả hai trường hợp, lần làm tròn đầu tiên của tích sẽ mang lại kết quả khác với m nhiều nhất là 1/2 và một lần nữa theo các đối số trước đó, lần làm tròn thứ hai sẽ làm tròn thành m. z

Bằng chứng trước cho thấy tích chỉ có thể làm tròn hai lần nếu thương thực hiện, và thậm chí sau đó, nó làm tròn thành kết quả đúng. Bằng chứng cũng chỉ ra rằng việc mở rộng lập luận của chúng ta để bao gồm khả năng làm tròn hai lần có thể là một thách thức ngay cả đối với một chương trình chỉ có hai phép tính dấu phẩy động. Đối với một chương trình phức tạp hơn, có thể không thể tính toán một cách có hệ thống các tác động của phép làm tròn hai lần, chưa kể đến các kết hợp tổng quát hơn của phép tính chính xác kép mở rộng và kép

Hỗ trợ ngôn ngữ lập trình cho độ chính xác mở rộng

Không nên dùng các ví dụ trước để gợi ý rằng bản thân độ chính xác mở rộng là có hại. Nhiều chương trình có thể được hưởng lợi từ độ chính xác mở rộng khi lập trình viên có thể sử dụng nó một cách có chọn lọc. Thật không may, các ngôn ngữ lập trình hiện tại không cung cấp đủ phương tiện để lập trình viên chỉ định thời điểm và cách sử dụng độ chính xác mở rộng. Để chỉ ra những hỗ trợ nào là cần thiết, chúng tôi xem xét các cách mà chúng tôi có thể muốn quản lý việc sử dụng độ chính xác mở rộng

Trong một chương trình di động sử dụng độ chính xác kép làm độ chính xác làm việc danh nghĩa của nó, có năm cách chúng ta có thể muốn kiểm soát việc sử dụng độ chính xác rộng hơn

  1. Biên dịch để tạo mã nhanh nhất, sử dụng độ chính xác mở rộng nếu có thể trên các hệ thống dựa trên mở rộng. Rõ ràng hầu hết các phần mềm số không yêu cầu nhiều phép tính hơn là lỗi tương đối trong mỗi thao tác được giới hạn bởi "máy epsilon". Khi dữ liệu trong bộ nhớ được lưu trữ ở độ chính xác kép, epsilon của máy thường được coi là lỗi làm tròn tương đối lớn nhất ở độ chính xác đó, vì dữ liệu đầu vào (đúng hoặc sai) được cho là đã được làm tròn khi chúng được nhập và kết quả sẽ . Do đó, trong khi tính toán một số kết quả trung gian ở độ chính xác mở rộng có thể mang lại kết quả chính xác hơn, thì độ chính xác mở rộng là không cần thiết. Trong trường hợp này, chúng tôi có thể muốn trình biên dịch chỉ sử dụng độ chính xác mở rộng khi nó sẽ không làm chậm chương trình một cách đáng kể và nếu không thì sử dụng độ chính xác kép
  2. Sử dụng định dạng rộng hơn gấp đôi nếu định dạng đó đủ nhanh và đủ rộng, nếu không thì sử dụng định dạng khác. Một số tính toán có thể được thực hiện dễ dàng hơn khi có độ chính xác mở rộng, nhưng chúng cũng có thể được thực hiện với độ chính xác kép chỉ với nỗ lực lớn hơn một chút. Xem xét tính toán định mức Euclide của một vectơ số chính xác kép. Bằng cách tính bình phương của các phần tử và tích lũy tổng của chúng ở định dạng kép mở rộng IEEE 754 với phạm vi số mũ rộng hơn, chúng ta có thể tránh được tình trạng tràn hoặc tràn quá sớm đối với các vectơ có độ dài thực tế. Trên các hệ thống dựa trên mở rộng, đây là cách nhanh nhất để tính định mức. Trên các hệ thống đơn/kép, một định dạng kép mở rộng sẽ phải được mô phỏng trong phần mềm (nếu một phần mềm được hỗ trợ) và việc mô phỏng như vậy sẽ chậm hơn nhiều so với việc chỉ sử dụng độ chính xác kép, kiểm tra các cờ ngoại lệ để xác định xem có xảy ra tràn hoặc tràn . Lưu ý rằng để hỗ trợ việc sử dụng độ chính xác mở rộng này, một ngôn ngữ phải cung cấp cả chỉ báo về định dạng có sẵn rộng nhất, tốc độ hợp lý, để chương trình có thể chọn phương pháp sử dụng và các tham số môi trường cho biết độ chính xác và phạm vi của từng định dạng . g. , rằng nó có phạm vi rộng hơn gấp đôi)
  3. Sử dụng định dạng rộng hơn gấp đôi ngay cả khi nó phải được mô phỏng trong phần mềm. Đối với các chương trình phức tạp hơn ví dụ chuẩn Euclide, lập trình viên có thể chỉ muốn tránh viết hai phiên bản của chương trình và thay vào đó dựa vào độ chính xác mở rộng ngay cả khi nó chậm. Một lần nữa, ngôn ngữ phải cung cấp các tham số môi trường để chương trình có thể xác định phạm vi và độ chính xác của định dạng có sẵn rộng nhất
  4. Không sử dụng độ chính xác rộng hơn; . Đối với các chương trình được viết dễ dàng nhất phụ thuộc vào số học chính xác kép được làm tròn chính xác, bao gồm một số ví dụ được đề cập ở trên, ngôn ngữ phải cung cấp cách để lập trình viên chỉ ra rằng không được sử dụng độ chính xác mở rộng, mặc dù có thể tính toán kết quả trung gian . (Các kết quả trung gian được tính toán theo cách này vẫn có thể phải làm tròn hai lần nếu chúng bị tràn khi lưu vào bộ nhớ. nếu kết quả của một phép toán số học trước tiên được làm tròn thành 53 bit có nghĩa, sau đó được làm tròn lại thành ít bit có nghĩa hơn khi nó phải được không chuẩn hóa, thì kết quả cuối cùng có thể khác với kết quả thu được bằng cách làm tròn chỉ một lần thành số không được chuẩn hóa. Tất nhiên, hình thức làm tròn kép này rất khó có thể ảnh hưởng xấu đến bất kỳ chương trình thực tế nào. )
  5. Làm tròn kết quả chính xác cho cả độ chính xác và phạm vi của định dạng kép. Việc thực thi nghiêm ngặt độ chính xác kép này sẽ hữu ích nhất cho các chương trình kiểm tra phần mềm số hoặc bản thân số học gần giới hạn của cả phạm vi và độ chính xác của định dạng kép. Các chương trình kiểm tra cẩn thận như vậy có xu hướng khó viết theo cách di động; . Do đó, một lập trình viên sử dụng hệ thống dựa trên mở rộng để phát triển phần mềm mạnh mẽ phải có khả năng di động đối với tất cả các triển khai IEEE 754 sẽ nhanh chóng đánh giá cao khả năng mô phỏng số học của các hệ thống đơn/kép mà không cần nỗ lực phi thường

Không có ngôn ngữ hiện tại nào hỗ trợ tất cả năm tùy chọn này. Trên thực tế, một số ngôn ngữ đã cố gắng cung cấp cho lập trình viên khả năng kiểm soát việc sử dụng độ chính xác mở rộng. Một ngoại lệ đáng chú ý là ISO/IEC 9899. Ngôn ngữ lập trình 1999 - Tiêu chuẩn C, phiên bản mới nhất của ngôn ngữ C, hiện đang ở giai đoạn chuẩn hóa cuối cùng

Tiêu chuẩn C99 cho phép triển khai đánh giá các biểu thức ở định dạng rộng hơn định dạng thường được liên kết với loại của chúng, nhưng tiêu chuẩn C99 khuyến nghị sử dụng một trong ba phương pháp đánh giá biểu thức duy nhất. Ba phương pháp được đề xuất được đặc trưng bởi mức độ mà các biểu thức được "thăng cấp" thành các định dạng rộng hơn và việc triển khai được khuyến khích xác định phương pháp nào nó sử dụng bằng cách xác định macro tiền xử lý

     n = n/2
62. nếu
     n = n/2
62 là 0, mỗi biểu thức được đánh giá theo định dạng tương ứng với loại của nó; . (Việc triển khai được phép đặt
     n = n/2
62 thành -1 để chỉ ra rằng phương pháp đánh giá biểu thức là không thể xác định được. ) Tiêu chuẩn C99 cũng yêu cầu tệp tiêu đề
     n = n/2
72 xác định các loại
     n = n/2
73 và
     n = n/2
74, có độ rộng tối thiểu tương ứng là
     n = n/2
65 và
     n = n/2
176 và nhằm khớp với các loại được sử dụng để đánh giá các biểu thức
     n = n/2
65 và
     n = n/2
176. Ví dụ: nếu
     n = n/2
62 là 2, cả
     n = n/2
73 và
     n = n/2
74 đều là
     n = n/2
70. Cuối cùng, tiêu chuẩn C99 yêu cầu tệp tiêu đề
     n = n/2
72 xác định các macro tiền xử lý chỉ định phạm vi và độ chính xác của các định dạng tương ứng với từng loại dấu phẩy động

Sự kết hợp các tính năng được yêu cầu hoặc khuyến nghị bởi tiêu chuẩn C99 hỗ trợ một số trong năm tùy chọn được liệt kê ở trên nhưng không phải tất cả. Ví dụ: nếu một triển khai ánh xạ loại

     n = n/2
70 sang định dạng kép mở rộng và định nghĩa
     n = n/2
62 là 2, thì lập trình viên có thể giả định một cách hợp lý rằng độ chính xác mở rộng tương đối nhanh, vì vậy các chương trình như ví dụ định mức Euclide có thể chỉ cần sử dụng các biến trung gian loại
     n = n/2
70 ( . Mặt khác, việc triển khai tương tự phải giữ các biểu thức ẩn danh ở độ chính xác mở rộng ngay cả khi chúng được lưu trữ trong bộ nhớ (e. g. , khi trình biên dịch phải tràn các thanh ghi dấu phẩy động) và nó phải lưu trữ kết quả của các biểu thức được gán cho các biến được khai báo
     n = n/2
176 để chuyển đổi chúng thành độ chính xác kép ngay cả khi chúng có thể được giữ trong các thanh ghi. Do đó, cả loại
     n = n/2
176 và loại
     n = n/2
74 đều không thể được biên dịch để tạo mã nhanh nhất trên phần cứng dựa trên mở rộng hiện tại

Tương tự như vậy, tiêu chuẩn C99 cung cấp các giải pháp cho một số vấn đề được minh họa bằng các ví dụ trong phần này nhưng không phải tất cả. Phiên bản tiêu chuẩn C99 của hàm

     n = n/2
91 được đảm bảo hoạt động chính xác nếu biểu thức
     n = n/2
190 
     n = n/2
044 
     n = n/2
006 được gán cho một biến (thuộc bất kỳ loại nào) và biến đó được sử dụng xuyên suốt. Tuy nhiên, một chương trình tiêu chuẩn C99 di động, hiệu quả để chia một số có độ chính xác kép thành các phần cao và thấp, tuy nhiên, khó khăn hơn. làm cách nào chúng tôi có thể phân tách ở đúng vị trí và tránh làm tròn hai lần nếu chúng tôi không thể đảm bảo rằng các biểu thức
     n = n/2
176 được làm tròn chính xác để đạt được độ chính xác gấp đôi? . Định lý 14 nói rằng chúng ta có thể tách ở bất kỳ vị trí bit nào miễn là chúng ta biết độ chính xác của phép tính cơ bản, và
     n = n/2
62 và macro tham số môi trường sẽ cung cấp cho chúng ta thông tin này

Đoạn sau đây cho thấy một triển khai có thể

     n = n/2
47
     n = n/2
47
     n = n/2
49
     n = n/2
50
     n = n/2
51
     n = n/2
52
     n = n/2
53
     n = n/2
54
     n = n/2
55
     n = n/2
56
     n = n/2
57
     n = n/2
58
     n = n/2
59
     n = n/2
0
     n = n/2
1

Tất nhiên, để tìm ra giải pháp này, lập trình viên phải biết rằng biểu thức

     n = n/2
176 có thể được đánh giá với độ chính xác mở rộng, rằng vấn đề làm tròn hai lần sau đó có thể khiến thuật toán gặp trục trặc và độ chính xác mở rộng đó có thể được sử dụng thay thế theo Định lý 14. Một giải pháp rõ ràng hơn chỉ đơn giản là xác định rằng mỗi biểu thức được làm tròn chính xác để tăng gấp đôi độ chính xác. Trên các hệ thống dựa trên mở rộng, điều này chỉ yêu cầu thay đổi chế độ chính xác làm tròn, nhưng thật không may, tiêu chuẩn C99 không cung cấp cách di động để thực hiện việc này. (Bản nháp ban đầu của Chỉnh sửa dấu phẩy động C, tài liệu làm việc chỉ định các thay đổi được thực hiện đối với tiêu chuẩn C90 để hỗ trợ dấu phẩy động, khuyến nghị rằng việc triển khai trên các hệ thống có chế độ chính xác làm tròn cung cấp các hàm
     n = n/2
99 và
     n = n/2
300 để lấy và đặt . Đề xuất này đã bị xóa trước khi các thay đổi được thực hiện đối với tiêu chuẩn C99. )

Thật trùng hợp, cách tiếp cận của tiêu chuẩn C99 để hỗ trợ tính di động giữa các hệ thống có khả năng số học số nguyên khác nhau gợi ý một cách tốt hơn để hỗ trợ các kiến ​​trúc dấu phẩy động khác nhau. Mỗi triển khai tiêu chuẩn C99 cung cấp tệp tiêu đề

     n = n/2
72 xác định các loại số nguyên mà triển khai hỗ trợ, được đặt tên theo kích thước và hiệu quả của chúng. ví dụ:
     n = n/2
304 là loại số nguyên có độ rộng chính xác là 32 bit,
     n = n/2
305 là loại số nguyên nhanh nhất có độ rộng tối thiểu 16 bit của triển khai và
     n = n/2
306 là loại số nguyên rộng nhất được hỗ trợ. Người ta có thể tưởng tượng một sơ đồ tương tự cho các loại dấu phẩy động. ví dụ:
     n = n/2
307 có thể đặt tên cho loại dấu phẩy động với độ chính xác chính xác là 53 bit nhưng có thể có phạm vi rộng hơn,
     n = n/2
308 có thể đặt tên cho loại triển khai nhanh nhất với độ chính xác ít nhất là 24 bit và
     n = n/2
309 có thể đặt tên cho loại nhanh hợp lý rộng nhất được hỗ trợ. Các loại nhanh có thể cho phép trình biên dịch trên các hệ thống dựa trên mở rộng tạo mã nhanh nhất có thể chỉ tuân theo ràng buộc rằng giá trị của các biến được đặt tên không được thay đổi do tràn thanh ghi. Các loại chiều rộng chính xác sẽ khiến trình biên dịch trên các hệ thống dựa trên mở rộng đặt chế độ làm tròn độ chính xác thành làm tròn theo độ chính xác đã chỉ định, cho phép phạm vi rộng hơn chịu cùng một ràng buộc. Cuối cùng,
     n = n/2
74 có thể đặt tên cho một loại có cả độ chính xác và phạm vi của định dạng kép IEEE 754, cung cấp đánh giá kép nghiêm ngặt. Cùng với các macro tham số môi trường được đặt tên phù hợp, sơ đồ như vậy sẽ sẵn sàng hỗ trợ tất cả năm tùy chọn được mô tả ở trên và cho phép người lập trình chỉ ra ngữ nghĩa dấu phẩy động mà chương trình của họ yêu cầu một cách dễ dàng và rõ ràng.

Hỗ trợ ngôn ngữ cho độ chính xác mở rộng có phức tạp như vậy không? . Tuy nhiên, các hệ thống dựa trên mở rộng đặt ra những lựa chọn khó khăn. chúng không hỗ trợ tính toán chính xác kép thuần túy cũng như tính toán chính xác mở rộng thuần túy hiệu quả như hỗn hợp của cả hai và các chương trình khác nhau yêu cầu các hỗn hợp khác nhau. Hơn nữa, không nên để việc lựa chọn thời điểm sử dụng độ chính xác mở rộng cho những người viết trình biên dịch, những người thường bị cám dỗ bởi các điểm chuẩn (và đôi khi được các nhà phân tích số nói thẳng ra) coi số học dấu phẩy động là "không chính xác vốn có" và do đó không xứng đáng cũng như không có khả năng. . Thay vào đó, sự lựa chọn phải được trình bày cho các lập trình viên và họ sẽ yêu cầu các ngôn ngữ có khả năng thể hiện sự lựa chọn của họ.

Sự kết luận

Những nhận xét ở trên không nhằm mục đích chê bai các hệ thống dựa trên mở rộng mà để vạch trần một số sai lầm, điều đầu tiên là tất cả các hệ thống IEEE 754 phải cung cấp kết quả giống hệt nhau cho cùng một chương trình. Chúng tôi đã tập trung vào sự khác biệt giữa các hệ thống dựa trên mở rộng và hệ thống đơn/kép, nhưng còn có sự khác biệt nữa giữa các hệ thống trong mỗi họ này. Ví dụ: một số hệ thống đơn/đôi cung cấp một lệnh duy nhất để nhân hai số và cộng số thứ ba chỉ với một lần làm tròn cuối cùng. Hoạt động này, được gọi là phép cộng-thêm kết hợp, có thể khiến cùng một chương trình tạo ra các kết quả khác nhau trên các hệ thống đơn/kép khác nhau và giống như độ chính xác mở rộng, nó thậm chí có thể khiến cùng một chương trình tạo ra các kết quả khác nhau trên cùng một hệ thống tùy thuộc vào việc . (Một phép cộng nhân hợp nhất cũng có thể làm hỏng quá trình tách của Định lý 6, mặc dù nó có thể được sử dụng theo cách không di động để thực hiện nhiều phép nhân chính xác mà không cần tách. ) Mặc dù tiêu chuẩn IEEE không lường trước được hoạt động như vậy, nhưng nó vẫn tuân thủ. sản phẩm trung gian được chuyển đến một "điểm đến" ngoài tầm kiểm soát của người dùng đủ rộng để giữ chính xác sản phẩm đó và tổng cuối cùng được làm tròn chính xác để khớp với điểm đến chính xác đơn hoặc kép

Ý tưởng rằng IEEE 754 quy định chính xác kết quả mà một chương trình nhất định phải cung cấp dù sao cũng rất hấp dẫn. Nhiều lập trình viên muốn tin rằng họ có thể hiểu hành vi của một chương trình và chứng minh rằng nó sẽ hoạt động chính xác mà không cần tham khảo trình biên dịch biên dịch nó hoặc máy tính chạy nó. Theo nhiều cách, hỗ trợ niềm tin này là một mục tiêu đáng giá cho các nhà thiết kế hệ thống máy tính và ngôn ngữ lập trình. Thật không may, khi nói đến số học dấu phẩy động, mục tiêu hầu như không thể đạt được. Các tác giả của các tiêu chuẩn IEEE biết điều đó và họ đã không cố gắng đạt được nó. Kết quả là, mặc dù gần như tuân thủ (hầu hết) tiêu chuẩn IEEE 754 trong toàn ngành công nghiệp máy tính, các lập trình viên của phần mềm di động vẫn phải tiếp tục đối phó với số học dấu chấm động không thể đoán trước.

Nếu các lập trình viên khai thác các tính năng của IEEE 754, họ sẽ cần các ngôn ngữ lập trình giúp dự đoán số học dấu phẩy động. Tiêu chuẩn C99 cải thiện khả năng dự đoán ở một mức độ nào đó với chi phí yêu cầu các lập trình viên viết nhiều phiên bản chương trình của họ, mỗi phiên bản cho một

     n = n/2
62. Liệu các ngôn ngữ trong tương lai có chọn thay vào đó để cho phép các lập trình viên viết một chương trình duy nhất với cú pháp thể hiện rõ ràng mức độ phụ thuộc vào ngữ nghĩa của IEEE 754 hay không vẫn còn phải xem. Các hệ thống dựa trên mở rộng hiện tại đe dọa triển vọng đó bằng cách cám dỗ chúng ta giả định rằng trình biên dịch và phần cứng có thể biết rõ hơn lập trình viên về cách thực hiện tính toán trên một hệ thống nhất định. Giả định đó là sai lầm thứ hai. độ chính xác được yêu cầu trong kết quả tính toán không phụ thuộc vào máy tạo ra nó mà chỉ phụ thuộc vào kết luận rút ra từ nó, và của người lập trình, trình biên dịch và phần cứng, tốt nhất chỉ có người lập trình mới có thể biết những kết luận đó có thể là gì