Web And Email Hosting, Jsp Web Hosting Blog

4.2.1 SINGLE-PRECISION CALCULATIONS 207 4) Unforeseen rounding overfIow. When a number like .999999997 is rounded to 8 digits, a carry will occur to the left of the decimal point, and the result must be scaled to the right. Many people have mistakenly concluded that rounding overflow is impossible during multiplication, since they look at the maximum value of lfJV 1, which is 1 - 2b-p + b- p; and this cannot round up to 1. The fallacy in this reasoning is exhibited in exercise 11. Curiously, it turns out that the phenomenon of rounding overflow is impossible during floating point division (see exercise 12). There is a school of thought that says it is harmless to round a value like .999999997 to .99999999 instead of to 1.0000000, since this does not increase the worst-case bounds on relative error. The floating point number 1.0000000 may be said to represent all real values in the interval [l.OOOOOOO-5 X 10P8, 1.0000000+5 X lo-s], while .99999999 represents all values in the much smaller interval (.99999999 -5 X lo- , .99999999 + 5 x lo- ). Even though the latter interval does not contain the original value .999999997, each number of the second interval is contained in the first, so subsequent calculations with the second interval are no less accurate than with the first. This ingenious argument is, however, incompatible with the mathematical philosophy of floating point arithmetic expressed in Section 4.2.2. 5) Rounding before normalizing. Inaccuracies are caused by premature round- ing in the wrong digit position. This error is obvious when rounding is being done to the left of the appropriate position; but it is also dangerous in the less obvious cases where rounding is first done too far to the right, followed by rounding in the true position. For this reason it is a mistake to round during the scaling-right operation in step A5, except as prescribed in exercise 5. (The special case of rounding in step N5, then rounding again after rounding overflow has occurred, is harmless, however, because rounding overflow always yields ~1.0000000 and this is unaffected by the subsequent rounding process.) 6) Failure to retain enough precision in intermediate calculations. Detailed analyses of the accuracy of floating point arithmetic, made in the next section, suggest strongly that normalizing floating point routines should always deliver a properly rounded result to the maximum possible accuracy. There should be no exceptions to this dictum, even in cases that occur with extremely low probability; the appropriate number of significant digits should be retained throughout the computations, as stated in Algorithms A and M. C. Floating point hardware. Nearly every large computer intended for scientific calculations includes floating point arithmetic as part of its repertoire of built-in operations. Unfortunately, the design of such hardware usually includes some anomalies that result in dismally poor behavior in certain circumstances, and we hope that future computer designers will pay more attention to providing the proper behavior than they have in the past. It costs only a little more to build the machine right, and considerations in the following section show that substantial benefits will be gained. Yesterday s compromises are no longer appropriate for modern machines, based on what we know now.

206 ARITHMETIC 4.2.1 2) Failure to treat exponent underfIow or overfIow properly. The size of e, should not be checked until after the rounding and normalization, because preliminary tests may give an erroneous indication. Exponent underflow and overflow can occur on floating point addition and subtraction, not only during multiplication and division; and even though this is a rather rare occurrence, it must be tested each time. Enough information should be retained so that mean- ingful corrective actions are possible after overflow or underflow has occurred. It has unfortunately become customary in many instances to ignore exponent underflow and simply to set underflowed results to zero with no indication of error. This causes a serious loss of accuracy in most cases (indeed, it is the loss of all the significant digits), and the assumptions underlying floating point arithmetic have broken down, so the programmer really must be told when underflow has occurred. Setting the result to zero is appropriate only in certain cases when the result is later to be added to a significantly larger quantity. When exponent underflow is not detected, we find mysterious situations in which (U @ V) @ w is zero, but u @ (V @ w) is not, since u @ ZJ results in exponent underflow but u @ (U @J w) can be calculated without any exponents falling out of range. Similarly, we can find positive numbers a, b, c, d, and y such that (~69~ CD b)0(cC3~ $ 4 = 3, (11) (a CJ3 b0~)0(c Cl3 d0y) = 1 if exponent underflow is not detected. (See exercise 9.) Even though floating point routines are not precisely accurate, such a disparity as (11) is certainly unexpected when a, b, c, d, and y are all positive! Exponent underflow is usually not anticipated by a programmer, so he needs to be told about it.* 3) Inserted garbage. When scaling to the left it is important to keep from introducing anything but zeros at the right. For example, note the “ENTX 0″ instruction in line 21 of Program A, and the all-too-easily-forgotten ENTX 0 instruction in line 04 of the FLOT subroutine (10). (But it would be a mistake to clear register X after line 27 in the division subroutine.) *On the other hand, it must be admitted that today s high-level programming languages give the programmer little or no satisfactory way to make use of the information that a floating point routine wants to tell him; and the MIX programs in this section, which simply HLT when errors are detected, are even worse. There are numerous important applications in which exponent underflow is relatively harmless, and it is desirable to find a way for programmers to cope with such situations easily and safely. The practice of silently replacing underflows by zero has been thoroughly discredited, but there is another alternative that has recently been gaining much favor, namely to modify the definition we have given for floating point numbers, allowing an unnormalized fraction part when the exponent has its smallest possible value. This idea of gradual underflow, which was first embodied in the hardware of the Electrologica X8 computer, adds only a small amount of complexity to the algorithms, and it makes exponent underflow impossible during addition or subtraction. The simple formulas for relative error in Section 4.2.2 no longer hold in the presence of gradual underflow, so the topic is beyond the scope of this book. However, by using formulas like round(z) = ~(1 -6) + c, where 161 < ib -P and 1~1 < ib–P–4, one can show that gradual underflow succeeds in many important cases. See W. M. Kahan and J. Palmer, ACM SIGNUM Newsletter (Oct. 1979), 13-21.

4.2.1 SINGLE-PRECISION CALCULATIONS 205 13 FDIV STJ EXITF Floating point division subroutine: 14 JOV OFLO Ensure overflow is off. 15 STA TEMP TEMP + v. 16 STA FV (0 : 4) FV+-ItffffO. 17 LDI TEMP (EXP) 18 LD2 ACC(EXP) 19 DEC2 -Q,l r12 e e, -e, + q. Z?O ENTX 0 .21 LDA ACC 22 SLA 1 rA+ fu. 23 CMPA FV(I : 5) 24 JL *+3 Jump if lful < IfJ. .25 SRA I Otherwise, scale f,, right 26 INC2 1 and increase r12 by 1. .27 DIV FV Divide. .28 JNOV NORM Normalize, round, and exit. 29 DVZRO HLT 3 Unnormalized or zero divisor 1 The most noteworthy feature of this program is the provision for division in lines 23-26, which is made in order to ensure enough accuracy to round the answer. If Iful < Ifv I, straightforward application of Algorithm M would leave a result of the form & 0 f f f f in register A, and this would not allow a proper rounding without a careful analysis of the remainder (which appears in register X). So the program computes fw + fu/fv in this case, ensuring that fw is either zero or normalized in all cases; rounding can proceed with five significant bytes, possibly testing whether the remainder is zero. We occasionally need to convert values between fixed and floating point representations. A fix-to-float routine is easily obtained with the help of the normalization algorithm above; for example, in MIX, the following subroutine converts an integer to floating point form: 01 FLOT STJ FXITF Assume that rA = u, an integer. 02 JOV OFLO Ensure overflow is off. 03 ENT2 Q+5 Set raw exponent. (10) 04 ENTX 0 05 JMP NORM Normalize, round, and exit. 1 A float-to-fix subroutine is the subject of exercise 14. The debugging of floating point subroutines is usually a difficult job, since there are so many cases to consider. Here is a list of common pitfalls that often trap a programmer or machine designer who is preparing floating point routines: 1) Losing the sign. On many machines (not MIX), shift instructions between registers will affect the sign, and the shifting operations used in normalizing and scaling numbers must be carefully analyzed. The sign is also lost frequently when minus zero is present. (For example, Program A is careful to retain the sign of register A in lines 30-34. See also exercise 6.)

204 ARITHMETIC 4.2.1 The rather long section of code from lines 25 to 37 is needed because MIX has only a 5-byte accumulator for adding signed numbers while in general 2p+ 1 = 9 places of accuracy are required by Algorithm A. The program could be shortened to about half its present length if we were willing to sacrifice a little bit of its accuracy, but we shall see in the next section that full accuracy is important. Line 55 uses a nonstandard MIX instruction defined in Section 4.5.2. The running time for floating point addition and subtraction depends on several factors that are analyzed in Section 4.2.4. Now let us consider multiplication and division, which are simpler than addition, and which are somewhat similar to each other. Algorithm M (Floating point multiplication or division). Given base b, excess 9, pdigit, normalized floating point numbers u = (eu, fiL) and v = (e,, fv), this algorithm forms the product w = u @ v or the quotient w = u @ v. Ml. [Unpack.] Separate the exponent and fraction parts of the representations of LL and V. (Sometimes it is convenient, but not necessary, to test the operands for zero during this step.) M2. [Operate.] Set . . . e, + e, + ev - 4, fw + fu fv for multiplication; (g) ew +- e, -e, + q + 1, fW c (F1fu)/fV for division. (Since the input numbers are assumed to be normalized, it follows that either fW = 0, or l/b* 5 If,,,] < 1, or a division-by-zero error has occurred.) If necessary, the representation of fW may be reduced to p + 2 or p + 3 digits at this point, as in exercise 5. M3. [Normalize.] Perform Algorithm N on (e,, fW) to normalize, round, and pack the result. (Note: Normalization is simpler in this case, since scaling left occurs at most once, and since rounding overflow cannot occur after division.) 1 The following MIX subroutines, which are intended to be used in connection with Program A, illustrate the machine considerations necessary in connection with Algorithm M. Program M (Floating point multiplication and division). 01 Q EQU BYTE/2 q is half the byte size 02 FMUL STJ EXITF Floating point multiplication subroutine: 03 JOV OFLO Ensure overflow is off. 04 STA TEMP TEMP + w. 0.5 LDX ACC rx +-u. 06 STX FU(O:4) FU+ztffffO. 07 LDI TEMP (EXP) 08 LD2 ACC@XP) 09 INC2 -&,I r12 +- e, + ev -q. 10 SLA 1 11 MUL FU Multiply fu times fv. 12 JMP NORM Normalize, round, and exit.

4.2.1 SINGLE-PRECISION CALCULATIONS 203 19 4H INCI 0,2 20 5H LDA FV a ENTX 0 22 SRAX 0,l 23 6H ADD FU 24 JOV N4 25 JXZ NORM 26 CMPA =O=(l:l) 27 JNE N5 28 SRC 5 29 DECX 1 30 STA TEMP 31 STA HALF(O:O) 32 LDAN TEMP ADD HALF ADD HALF SRC 4 JMP N3A CON l//2 CON 0 CON 0 JAZ ZRO CMPA =O=(l:i) JNE N5 SLAX 1 DEC2 1 JMP N2 ENTX 1 SRC 1 INC2 1 CMPA =BYTE/2=(5:5) JL N6 JG 5F JXNZ 5F STA TEMP LDX TEMP(4:4) JXO N6 STA *+l(O:O) INCA BYTE JOV N4 J2N EXPUN ENTX 0,2 SRC I DEC2 BYTE STA ACC J2N * 33 34 35 36 37 HALF 38 FU 39 FV 40 NORM 41 N2 42 43 N3 44 N3A 45 46 N4 47 48 49 N5 50 51 52 53 54 55 56 5H 57 58 59 N6 60 N7 61 62 ZRO 63 8H 64 EXITF rI1 + e, -e,. (Step A4 unnecessary.) A5 Scale rig& Clear rX. Shift right e, -e, places. A6. Add. A7. Normdize. Jump if fraction overflow. Easy case? Is f normalized? If so, round it. IrXl ++ (rA[. (rX is positive.) (Operands had opposite signs, registers must be adjusted before rounding and normalization.) Complement least significant portion. Jump into normalization routine. One half the word size (Sign varies) Fraction part fU Fraction part fv Nl. Test fL N2. Is f normalized? To N5 if leading byte nonzero. N3. Scale left. Decrease e by 1. Return to N2. N4. Scale rig& Shift right, insert 1″ with proper sign. Increase e by 1. N5. Round. Is [tail1 < ib? Is ltaill > +b? Jtaill = 4 b; round to odd. To N6 if rX is odd. Store sign of rA. Add bd4 to IfI. (Sign varies) Check for rounding overflow. N6. Check e. Underflow if e < 0. N7. Pack. rX + e. r12 + e -b. Exit, unless e 2 b. 65 EXPOV HLT 2 Exponent overflow detected 66 EXPUN HLT 1 Exponent underflow detected 67 ACC CON 0 Floating point accumulator 1

202 ARITHMETIC 4.2.1 be expressed as a normalized floating point number in the required range, special action is necessary.) N7. [Pack.] Put e and f together into the desired output representation. 1 Some simple examples of floating point addition are given in exercise 4. The following MIX subroutines, for addition and subtraction of numbers having the form (4), show how Algorithms A and N can be expressed as computer programs. The subroutines below are designed to take one input u from symbolic location ACC, and the other input v comes from register A upon entrance to the subroutine. The output w appears both in register A and location ACC. Thus, a fixed point coding sequence LDA A; ADD B; SUB C; STA D (7) would correspond to the floating point coding sequence LDA A, STA ACC; LDA B, JMP FADD; LDA C, JMP FSUEl; STA D. (8) Program A (Addition, subtraction, and normalization). The following program is a subroutine for Algorithm A, and it is also designed so that the normalization portion can be used by other subroutines that appear later in this section. In this program and in many others throughout this chapter, OFLO stands for a subroutine that prints out a message to the effect that MIX s overflow toggle was unexpectedly found to be on. The byte size b is assumed to be a multiple of 4. The normalization routine NORM assumes that r12 = e and rAX = f, where rA = 0 implies rX = 0 and r12 < b. 00 BYTE EQU 1(4:4) Byte size b 01 EXP E&U I:1 Definition of exponent field 0.2 FSUB STA TEMP Floating point subtraction subroutine: 03 LDAN TEMP Change sign of operand. 04 FADD STJ EXITF Floating point addition subroutine: 05 JOV OFLO Ensure overflow is off. 06 STA TEMP TEMP +-v. 07 LDX ACC rX + 21. 08 CMPA ACC @XI ) Steps Al. A2, A3 are combined here: 09 JGE IF Jump if e, 2 e, . 10 STX FU (0 : 4) FU+ fffffo. 11 LD2 ACC@XP) r12 + e, . 12 STA FV (0 : 4) 13 LDIN TEMP @XI ) rI1 + -e, . 1-4 JMP 4F 15 IH STA FU (0 : 4) FU t f f f f f 0 (u, v interchanged) . 16 LD2 TEMP (EXP) r12 + e, . 17 STX FV (0 : 4) 18 LDlN ACC (!ZXF ) rI1 t -e, .

4.2.1 SINGLE-PRECISION CALCULATIONS 201 N4. Scale right Nl. underflow N7. Pack Fig. 3. Normalization of (e, f). base-b digits to the right of the radix point. If such a large accumulator is not available, it is possible to shorten the requirement to p + 2 or p + 3 places if proper precautions are taken; the details are given in exercise 5.1 A6. [Add.] Set fw +- fU + fV. A7. [Normalize.] (At this point (G,,, fiu) represents the sum of u and 21, but lfwj may have more than p digits, and it may be greater than unity or less than l/b.) Perform Algorithm N below, to normalize and round (ezu, jiu) into the final answer. 1 Algorithm N (Normalization). A raw exponent e and a raw fraction f are converted to normalized form, rounding if necessary to p digits. This algorithm assumes that ]f 1 < b. Nl. [Test f .] If If] 2 1 ( fraction overflow ), go to step N4. If f = 0, set e to its lowest possible value and go to step N7. N2. [Is f normalized?] If ] f ] 2 l/b, go to step N5. N3. [Scale left.] Shift f to the left by one digit position (i.e., multiply it by b), and decrease e by 1. Return to step N2. N4. [Scale right.] Shift f to the right by one digit position (i.e., divide it by b), and increase e by 1. N5. [Round.] Round f to p places. (We take this to mean that f is changed to the nearest multiple of b-p. It is possible that (bpf)modl = f so that there are two nearest multiples; if b is even, we choose the one that makes bPf + $b odd. Further discussion of rounding appears in Section 4.2.2.) It is important to note that this rounding operation can make If] = 1 ( rounding overflow ); in such a case, return to step N4. N6. [Check e.] If e is too large, i.e., larger than its allowed range, an exponent overflow condition is sensed. If e is too small, an exponent underflow condition is sensed. (See the discussion below; since the result cannot

200 ARITHMETIC 4.2.1 Al. Unpack i-; A7. Normalize I Fig. 2. Floating point addition. Note: Since floating point arithmetic is inherently approximate, not exact, we will use round symbols $7 8, 8, 0 to stand for ffoating point addition, subtraction, multiplication, and division, respectively, in order to distinguish approximate operations from the true ones. The basic idea involved in floating point addition is fairly simple: Assuming that eu 2 e,, we take e, = eu, fW = fu. + fV/beu+~ (thereby aligning the radix points for a meaningful addition), and normalize the result. Several situations can arise that make this process nontrivial, and the following algorithm explains the method more precisely. Algorithm A (Floating point addition). Given base b, excess 4, p-digit, normalized floating point numbers u = (ezL, fu) and Y = (e,, fV), this algorithm forms the sum w = u$ w. The same procedure may be used for floating point subtraction, if -v is substituted for V. Al. [Unpack.] Separate the exponent and fraction parts of the representations of u and v. A2. [Assume e, 2 e,.] If e, < e,, interchange u and w. (In many cases, it is best to combine step A2 with step Al or with some of the later steps.) A3. [Set e,.] Set e, c ezL . A4. [Test e, -e, .] If e, -e, 2 p+ 2 (large difference in exponents), set fW c fu. and go to step A7. (Actually, since we are assuming that ZL is normalized, we could terminate the algorithm; but it is occasionally useful to be able to normalize a possibly unnormalized number by adding zero to it.) A5. [Scale right.] Shift f,, to the right e, -e, places; i.e., divide it by Vu- u. [Note: This will be a shift of up to p + 1 places, and the next step (which adds fiL to fV) thereby requires an accumulator capable of holding 2p + 1

4.2.1 SINGLE-PRECISION CALCULATIONS 199 used for this purpose, notably characteristic and mantissa ; but it is an abuse of terminology to call the fraction part a mantissa, since this concept has quite a different meaning in connection with logarithms. Furthermore the English word mantissa means a worthless addition. ) The MIX computer assumes that its floating point numbers have the form (4 Here we have base b, excess q, floating point notation with four bytes of precision, where b is the byte size (e.g., b = 64 or b = loo), and q is equal to [$b]. The fraction part is f j j j j, and e is the exponent, which lies in the range 0 2 e < b. This internal representation is typical of the conventions in most existing computers, although b is a much larger base than usual. B. Normalized calculations. A floating point number (e, j) is normalized if the most significant digit of the representation of f is nonzero, so that l/b I VI < 1; (5) or if j = 0 and e has its smallest possible value. It is possible to tell which of two normalized floating point numbers has a greater magnitude by comparing the exponent parts first, and then testing the fraction parts only if the exponents are equal. Most floating point routines now in use deal almost entirely with normalized numbers: inputs to the routines are assumed to be normalized, and the outputs are always normalized. Under these conventions we lose the ability to represent a few numbers of very small magnitude-for example, the value (0, .OOOOOOOl) can t be normalized without producing a negative exponent-but we gain in speed, uniformity, and the ability to give relatively simple bounds on the relative error in our computations. (Unnormalized floating point arithmetic is discussed in Section 4.2.2.) Let us now study the normalized floating point operations in detail. At the same time we can consider the construction of subroutines for these operations, assuming that we have a computer without built-in floating point hardware. Machine-language subroutines for floating point arithmetic are usually writ- ten in a very machine-dependent manner, using many of the wildest idiosyn- crasies of the computer at hand; so floating point addition subroutines for two different machines usually bear little superficial resemblance to each other. Yet a careful study of numerous subroutines for both binary and decimal computers reveals that these programs actually have quite a lot in common, and it is possible to discuss the topics in a machine-independent way. The first (and by far the most difficult!) algorithm we shall discuss in this section is a procedure for floating point addition, (eu, .fd $ (euj .fv) = (k, full. (6)

198 ARITHMETIC 4.2 4.2. FLOATING POINT ARITHMETIC IN THIS SECTION, we shall study the basic principles of doing arithmetic on floating point numbers, by analyzing the internal mechanisms underlying such calculations. Perhaps many readers will have little interest in this subject, since their computers either have built-in floating point instructions or their computer manufacturer has supplied suitable subroutines. But, in fact, the material of this section should not merely be the concern of computer-design engineers or of a small clique of people who write library subroutines for new machines; every well-rounded programmer ought to have a knowledge of what goes on during the elementary steps of floating point arithmetic. This subject is not at all as trivial as most people think; it involves a surprising amount of interesting information. 4.2.1. Single-Precision Calculations A. Floating point notation. We have discussed fixed point notation for numbers in Section 4.1; in such a case the programmer knows where the radix point is assumed to lie in the numbers he manipulates. For many purposes it is considerably more convenient to let the position of the radix point be dynamically variable or floating as a program is running, and to carry with each number an indication of its current radix point position. This idea has been used for many years in scientific calculations, especially for expressing very large numbers like Avogadro s number N = 6.02252 x 1023, or very small numbers like Planck s constant h = 1.0545 X 1O-27 erg sec. In this section we shall work with base b, excess q, Aoating point numbers with p digits: Such numbers will be represented by pairs of values (e, f), denoting (e, f) = f x bepq. 0) Here e is an integer having a specified range, and f is a signed fraction. We will adopt the convention that Ifl < 1; in other words, the radix point appears at the left of the positional representation of f. More precisely, the stipulation that we have pdigit numbers means that bpf is an integer, and that -bP < bPf < bP. (2) The term floating binary implies that b = 2, floating decimal implies b = 10, etc. Using excess-50 floating decimal numbers with 8 digits, we can write, for example, Avogadro s number N = (74, +.60225200) ; Planck s constant h = (24, +.10545000). (3) The two components e and f of a floating point number are called the exponent and the fraction parts, respectively. (Other names are occasionally