Digitized by the Internet Archive in 2013 http://archive.org/details/continuedproduct920chow 7?^ltx-•<' M' 1^' 7 UIUCDCS-R-78-920 Continued-Product Algorithms for Division, Logarithm, Exponential and Square Root in Radix Ten by Catherine Yuk-Fun Chow UILU-ENG 78 1713 ^ May, 1978 m DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN The Library of the JUN 1 5 1978 liversity or Illinois hana-Cha' URBANA, ILUNOIS CONTINUED-PRODUCT ALGORITHMS FOR DIVISION, LOGARITHM, EXPONENTIAL AND SQUARE ROOT IN RADIX TEN BY CATHERINE YUK-FUN CHOW B.S., University of Michigan, 1974 THESIS Submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science in the Graduate College of the University of Illinois at Urbana-Champaign, 1978 Urbana, Illinois m Acknowledgment I wish to thank my advisor, Professor James E. Robertson, for his advice and guidance. I also wish to thank Mr. John L. Larson for his help in proofreading and Ms. Cathy Gallion for her excellent typing- IV Table of Contents 1. Introduction 1 2. Division 2 2.1 Finding the Selection Rules 3 2.2 The Division Algorithm 12 3. Logarithm 15 3.1 Natural Logarithm 15 3.2 Logarithm to the Base 10 20 4. Exponential to the Base e 26 4.1 Finding the Selection Rules 27 4.2 Algorithm for Exponential . 34 5. Exponential to the Base 10 38 5.1 Finding the Selection Rules 39 5.2 The Algorithm for Exponential to the Base 10 ... 45 6. Square Root 50 6.1 The Normalization Procedure 51 6.2 The Algorithm for Square Root 58 7. Summary 61 References 62 1. Introduction This paper presents radix 10 continued product (sum) algorithms for evaluating quotients, logarithms, exponentials and square roots. These algorithms generate the results digit by digit beginning with the most significant digit. The number of evaluation steps depends on the number of significant digits required in the results. The redundant digit set used in the continued products (sums) directly influences the rules for digit selection in the individual steps [3]. This work is based on the binary algorithms given by De Lugish [1] and the radix 16 algorithms by Ercegovac [2]. The particular form of con- tinued products (sums) used in these algorithms transforms to recursive evaluation steps, each of which typically sums 2 or 3 numbers together. While the radix 2 and the radix 16 algorithms are suitable for hardware im- plementation, the radix 10 algorithms are more suitable in a paper-and-pencil environment. With the necessary tables and an adding machine, the functions mentioned above can be easily evaluated to m significant digits in about (m+1) steps. This is an advantage over evaluating the functions by conven- tional methods which usually require full word-length multiplication or exponentiation and whose convergence rates are usually dependent on the operands. 2. Division This section describes how division in radix 10 can be performed through finding the factors of a convergent infinite product (continued product). Given any two real numbers Y and X, they can be normalized such that Y = yo • 10" and X = Xq • 10^ where ly^l, |xq| e [0.1, 1) and a , /3 are integers. Then Y/X = Y^/Xq • 10^""^^ where Yq = (sign of Xq) • y^ and Xq = |xq|. Since the difference (a-jS) can be easily found, the following discussion concentrates on finding Q = Y^/X^. Suppose Y^ and X^ are each multiplied by a sequence of factors M|^ (k = 0, 1, ..., m). Then Q = If we can find the appropriate factors such that the product X^ • n. _q M. approximates 1, Q would be approximated by the product Y^ • ^^^q M. . Let M. be defined as (1 + S. • 10 ) where S. is a digit constant in the set {7, 6, ..., 6,7} and k is an {0,1, ..., m}. It will be later shown that given any Xq in the interval [0.1, 1), the constants S. can be chosen such that the product XqITj^^q M. approximates 1 with an error less than 7/9 • 10"^. ^0 •"l •"2- . ^0 ^0 • ":=o \ '^o . M, •"2- m ■ K-o \ If IXq -n^^o^- 111 7/9 •10-'" then Yq . n;^.o ^k = Vh' h • K-0 \ = Q • (1+p) where the relative error p is less than 7/9 • 10'"' in magnitude. This means that Q can be evaluated by forming the continued product Yq • n._Q (1 + S. -10" ) in which the constants are so chosen that m — k X^ . n._Q (1 + S. -10 ) approximates 1 to m digits. 2.1 Finding the Selection Rules This section develops the rules for finding the constants S, such that the continued product Xq • n'|"_ (1 + S. -10" ) approximates 1 with an error less than 7/9 • 10" . Let x,,i = Xg-n^.o(i+s..io-^') then X|^^., = ''k • (1 ■^ ^k'''"''^' and X^^l = Xo-n^^jjd .S^.IO-'^). The constants S, are to be chosen in such a way that S. will make X.^, approximate 1 to k digits. \\+^ - 1| 1 7/9 . 10■^ (2.1) for k=0,l,...,m. In other words, at step k, X. already agrees with 1 to (k-1) A.L. digits, the selection of S. aims at forcing the k digit of Xj^ to agree with the k*^ digit of 1. Let the quantities (X.^-. - 1 ) be scaled by 10 so that the scaled remainders are bounded by the same quantity and the selection of S. depends on the same digital positions for each k. Let Rk+i = 10^(Xk,i - 1) (2.2) then R,^^^ = lo'' • [(10"'''^^-R^ + 1)(1 + S^'IO"'') - 1] = lOR^ + S^ + R^-S^-IO"''^^ (2.3) For the recursion (2.3), the criterion for selection S. becomes -7/9 < Rfc+i 1 7/9 for k = 0, 1 , .. ., m. Consider the following rule for selecting S. : Sj^l = the integer part of (|lOR. | + 0.5) (sign of Sj^) = -(sign of R|^) (2.4) This rule is applicable only if the last term in Equation (2.3) satisfies the constraint IR^-S^-IO""^"^^! 1 7/9-1/2 = 5/18. (2.5) It is found that IR^.S^-IO'"^"^^! < 0.06 (lO'"^"*""^). Therefore for k > 3, [R^.S^.IO -k+1 < 0.06 < 5/18 which satisfies (2.5). This means the constants S^, S^, ...» S can be selected by rounding lOR. to one digit as defined in (2.4). -k+1 Furthermore, the term (R. -S. -10 ) does not affect the first (k-3) digits of R.^-, since R,-S,.10- ■^■^^1 < 0.06 (lO"'"'^^). Given R. with k >^ 3, the digits S. , S. ^, , ...» Spi^,^ are actually available (after receding the most significant digits of R|^). When k ^ (m+3)/2, the -k+1 term R. -S. -10 will have no effect on the choice of the remaining digits The recursion (2.3) for evaluating R.^-, simplifies to ^k+1 = '^'h'h (2.6) for k > (m+3)/2. To complete the normalization, the rules for selecting Sq, S, and $2 still have to be found. It is more convenient to work backwards from S^ to S-j. The reason for doing this will become apparent as the reader follows along. The interval of R2 over which a certain Sp can be chosen is R3 - $2 R 3 - S2 ^10(l+S2l0"^) ' 10(l+S2l0"^)_ R. and R. denote the lower and upper bounds of R. respectively. The following table lists the intervals of IOR2 for $2 in {7, 6^, ...» 6, 7}. a 1 lORg 1 b 7 6 5 4 3 2 1 T 2 3 4 5 6 7 •7.2 •6.3 •5.4 ■4.5 •3.5 •2.6 ■1.7 ■0.7 0.3 1.3 2.4 3.4 4.5 5.6 6.8 -5.9 -5.0 -4.1 -3.2 -2.2 -1.3 -0.3 0.7 1.7 2.8 3.8 4.9 6.0 7.1 8.3 Table 1. Selection Intervals for {S^} All neighboring selection intervals for the S^'s overlap. How- ever, the rounding-to-one-digit selection rule defined in (2.4) can be applied only to R^ in the interval [-0.54, 0.55). The approach taken here is to define the selection rule of $« to be (2.4), and to restrict R^ to the interval [-0.54, 0.55). The selection interval of R, for S-. in the given digit set is ^2 -h ^2-h 10{l+S,10"b ' 10{1+S^10"M. S a < lOR^ < b 7 -4.4 -3.8 6 -4.0 -3.5 5 -3.6 -3.0 4 -3.2 -2.5 3 -2.7 -1.9 2 -2.1 -1.3 1 -1.4 -0.5 -0.5 0.55 T 0.53 1.7 2 1.9 3.1 3 3.6 5.0 4 5.8 7.5 5 9.0 11.1 6 13.7 16.3 7 21.6 25.1 Table 2. Selection Intervals for S^ The selection intervals below the clotted line in Table 2 cannot be used because they are not contiguous. Therefore, S, can only take the values 1, 0, ..., 7. If the selection rule (2.4) is used to select S, , R, has to be restricted to the narrow interval -.14 £ R-, < 0.15. This interval does not result in contiguous intervals for the selection of S^. Hence, the selection of S, has to be different from (2.4). Figure 1 presents a way of selecting S, . Finally, the selection intervals for Sq can be studied. Table 2 indicates that R, is restricted to the interval (-0.44, 0.17). Since R^ = Xq(1+Sq) - 1, the interval of X,. for a digit Sq is R^ + 1 1 + S, R^ + 1 1 + S J Since Xq is in the interval [0.1, 1), the values of S^ which are listed in Table 3 are sufficient. a £ Xq £ b 7 6 5 4 3 2 1 ,07 ,08 ,10 ,12 ,14 ,19 28 ,56 .14 .16 .19 .23 .29 .39 .58 1.17 Table 3. Selection Intervals for S 2.0 1 1.7 1.0 - )S,=1 0.53 -^ 0.0 - S^=0 -1.0 - -2.0 - -3.0 - -4.0 - -5.0 0.5 -1.3 -2.0 -2.5 -3.0 -3.5 -4.0 -4.4 Intervals of lOR^ S^ = l S^ = 2 S^=3 S^=4 S^=5 S^=6 S^ = 7 1.0 -I — 1.0 ^ 0.8 - 0.6 0.4 - 0.2 - 0.0 J ) V^ 0.58 ^ )so=i 0.35 -^ 0.19 ^0=2 0.1 j V^ Intervals of X, ^0 Figure 1, Selection of S, Figure 2. Selection of S. 10 As shown in Table 3, Sq is defined over the entire interval of X« ^ [0.1, 1). Figure 2 defines the selection of Sq over the interval 0.1 1 Xq < 1. m - k It has been shown that to normalize X^ • n (1 + S. 10 ) to 1 , two initial special steps are necessary. For k >_ 2, S. can be selected by rounding |10R. | to one digit and taking the sign opposite to that of R. . The entire normalization procedure is summarized in the following algorithm. The Normalization Algorithm N Objective: Given X^ in [0.1, 1), find the digit constants S. , k = 0, 1, ..., m, such that Xq • n^^^ (1 + Sj^-lO'"^) approximates 1 to m digits of accuracy. Basic Selection Rule: S|^| = the integer part of (|10R|^| + 0.5) (sign of S^) = -(sign of Rj^) Step : (i) Choose S^ according to Figure 2. (ii) R^ = Xq(1 + Sq) - 1 Step 1 : (i) Choose S, according to Figure 1. (ii) Form R^ using the recursion (N.2) (N.l) Basic Recursion: \+l ^ ""^^^^ ^ \'"'°'''^ ^ \ ^^'^^ Constants: k, = the smallest integer k such that k >^ (m+3)/2 The Normalization Procedure: n For k = 2, 3, . . . , m do Step k : (i) Choose S. using the rule (N.l) (ii) If k < k^ then form Rj^^^ with (N.2) else form R k+1 by: Example 1. Normalization Procedure Given X^ = 0.39, find the constants S. , k = 0, 1 , . . . , m, such m — k that Xq • n._Q (1 + S. -10 ) approximates 1 to m digits. Xq = 0.39, m = 9, k^ = 6 The normalization steps are listed below: k \ Selection Rules Used •^k+l Formula Used 1 Figure 2 -.22 ^1 = h^^'^0^ - 1 1 3 Figure 1 .14 N.2 2 T Rule N.1 .3860000001 N.2 3 4 Rule N.l -.1554399990 N.2 4 2 Rule N.l .4452891299 N.2 5 4 Rule N.l .4527131830 * N.2 6 4 Rule N.l .527... ^+1 = lOR, + S^ 7 5 Rule N.l .271... II 8 3 Rule N.l -.28... if 9 3 Rule N.l .1... II A- J. .-11.. J. U_ 1 1^ JJ_.*J._ _^ / n_\ _^j 1 £ .- Ac^tua^lly, the leading digits of (-Rs), after being recoded into the digit set {5, 4, ..., 4, 5}, would yield tjie last four selection digits. In this case. -Rg = -.45271 = .4533T b/- 4, by - 0, oq ~ J, bq = o. 12 2.2 The Division Algorithm A division algorithm utilizing the above multiplicative normaliza- tion can now be defined. The definition of the normalization procedure ensures that where the error p has a magnitude less than 10" . Then the continued product ^0 "k=0 ^^ "*" ^k'^^"*^^ ^^^^^ ^^^ quotient of Y^/X to m correct digits: = Yq/Xq • (1 .p). However, this does not take rounding errors into consideration. Let the number of extra digits necessary to effect a final result of m correct digits be Am, then (m+1) . 1/2 . 10"^'"^^"^^ < 1/2 • lO"'" (2.6) ^ 1 logiQ(m+l). Therefore, for m digits of accuracy, m+1 steps of normalization as well as a precision of (m + Am) digits for the partial results are necessary. The division algorithm is outlined below. 13 The Division Algorithm D Objective: To find (Y/X), where X ^ 0. (I) Initialization: (i) Normalize Y and X such that Y = Yq • lo"" and X = Xq • 10^ with IYqI, |xq| e [0.1, 1) (1i) Xq = |xq| Qq = (sign of Xq) • Yq (II) Evaluating the continued products: For k = 0, 1 , . . . , m do Step k: (i) Find S. and R.^, according to Step k of the normalization algorithm N. (i^") Vi -\' ^^ ^s^-10'^). (Ill) Final Step: (Y/X) = Q^^^ • 10^"-^^ ■ While it is true that the selection digits S. , S. ^-j , ...» S^^ can be found bY receding R. (see the Normalization Algorithm in Section 2.1) ■^1 one can find Q , di recti y from Q. and R. . n where R. is R. rounded to (m+l-k,) places. This is true because 14 n^(i.s.-io-^) = 1.2^ s.-io-^'.p 10" . The term (2.. S. -10" ) is not changed by the receding procedui I K^ 1 where p consists of terms with powers less than -{m+4) and |p| is less than ire, If the algorithm is not hardware implemented, (2.7) can be easily used. Although the continued product method of division is more complicated than conventional division, the computational procedures of the continued product can be applied to other functions such as logarithm, exponential and square root. Example 2. Division Given X = 0.03125 Y = 0.0009765625 (= X^) , find Q = Y/X with m=9. (I) Xq = .3125 a = -1 (II) ^0 = .9765625 /3 = -3 %- .9765625 k \ '^k+l "k 2 -.0625 2.9296875 1 1 .3125 3.22265625 2 3 .03125 3.125976563 3 .3125 3.125976563 4 3 .1240625 3.12503877 5 T .240624 3.12500752 ^10 = Qg(l - Q = Q^Q^ lo"'' * i. 12500752 (1 - .00002406) = 3.125 X 10 = 3. 125000001 = 0.03125 15 3. Logarithm 3.1 Natural Logarithm The normalization algorithm N given in Chapter 2 can be applied to evaluate the natural logarithm of a number. A table of constants (-ln(l + S, -10' )) is required for about half of the evaluation steps. Given X = Y^. • 10 with an integer exponent and X^ in the interval [0.1, 1), ln(X) = In(XQ) + a-ln(lO). Knowing a and In(lO), the last term can be easily evaluated. The problem is then reduced to finding ln(Xj.), where X^ ^ [0.1, 1). From the identity ° "k=o(^ ^ ^'^""^^ IhCXq) = ln[XQ • n^^^d + S^-IO""^)] + SJJ^Q -ln(l + S^-IO-'^) With the normalization procedure N, 1^0 ' "k=0^^ "^ \''^^'^^ ■ "•! - ^/^ * ^^'^• Then, for an accuracy of m digits, ln(XQ • nj^^^ci + S^-IO""")) = 0, and In(XQ) = S^^q -ln(l + S^-IO"""). The constants -ln(l + S.10~ ) are necessary for the evaluation of InCXg). 16 Note that ln(l + Sj^-lO"'^) = S^-IO""^ - (S,^^/2)-10"^'^+ (higher order terms). When k satisfies k 1 ^-^ ^^ m/2 + 0.9. (3.1) ln(l + Si^lO'"^) = S|^•^0"^ to m digits. Therefore, for k >_ m/2 + 0.9, the constants ln(l + S. -10 ) need not be precomputed. To provide for rounding errors, the table of constants and the intermediate results require an extended precision of (m +^m) digits, where ^m >• log,Q(m+l) as shown in (2.6). Algorithm for Natural Logarithm Objective: To find ln(X), where X > 0. (t) Initialization: (i) Normalize X such that X = X^ • l(f with X 6 [0.1, 1) (ii) k, = the smallest integer not less than (m/2 + 0.9) (iii) Lq = (II) Finding In(XQ): . . For k = 0, 1, . . . , m do Step k (i) Find S. and R.^, by performing Step k of the algorithm N. (ii) If k < k^ then L,^^^ = "-i^ + (-Ind+S^-lO"*^)), '^'' 4+1 = Lk - Sk-10'^- 17 (III) Final Step: ln(X) = L^^^ +a.ln(10) > A further simplification is possible. Once L. and R. in the above algorithm is obtained, L , can be found without carrying out the rest of the steps: ^m.l = \ ' \ ' 10 -k^+1 (3.2) n where R. is R. rounded to (m+l-k, ) places, n "^l ' Table 4 provides the necessary constants for evaluating In(XQ) to an accuracy of 9 digits (i.e., m = 9). The following example illustrates the algorithm with the simplification stated in (3.2). Exampl e 3. Natural logari thm Find ln(0. 03125) with m=9. (I) X = = 0.03125 ^0 = 0.3125 a = -1 ^0 = In(lO) = 2.302585093 (II) k \ 2 1 1 2 3 3 4 3 5 1 [From Table 4] ■^k+l ln(l + S^IO""^) 4.1 = 4-^"(i ' h .0625 -1.0986122887 -1.0986122887 .3125 -0.095310179 -1.1939224685 .03125 0.0304592075 -1.1634632610 .3125 0-0 -1.1634632610 .1240625 0.0003000450 -1.1631532160 .240625 0.000100000 -1.1631532160 18 (Iin ho = Lg^fgXlO-5 = -1.163153216 + 0.000002406 = -1.1631508100 a.ln(lO) = -2.3025850930 ln(0. 03125) = L^Q + a-ln(lO) = -3.4657359030 ■ Example 4. Natural logarithm Given Y = 0.0009765625, find ln(Y) with m = 9- (0.0009765625 = 0.03125 , compare with Example 3) (I) Xj,= .9765625 a - -3 ^0 = (II) k \ \+^ -ln(l + S^-IO'*^) -.0234375 0.0 1 1 .7421875 -.0953101798 2 7 -.09765625 .0725706928 3 1 .0224609375 -.0009995003 4 .224609375 0.0 5 2 .2460488281 .0000200002 ho = Lg . fgXlO-5 = -.0237189871 + .000002460 = -0.0237165271 a-ln(lO) = -6.907755279 4.1 = 4-i"(i ' V'^ 0.0 -.0953101798 -.022739487 -.0237389873 -.0237389873 -.0237189871 19 (III) ln(y) = L^Q + a-ln(lO) = -6.931471806 k = k = 1 2 .2231435513 1 - .6931471806 T .1053605157 2 -1.0986122887 1 -.0953101798 5 -1.7917594692 2 -.1823215568 3 -.2623642645 4 -.3364722366 5 -.4054951081 6 -.4700036292 7 -.5306282511 k = 2 k = 3 k = 4 k = 5 7 .0725706928 .0070245249 .0007002451 .0000700025 6 .0618754037 .0060180723 .0006001801 .0000600018 5" .0512932944 .0050125418 .0005001250 .0000500013 4 .0408219945 .0040080214 .0004000800 .0000400008 3 .0304592075 .0030045090 .0003000450 .0000300004 2 .0202027073 .0020020027 .0002000200 .0000200002 T .0100503359 .0010005003 .0001000050 .0000100000 1 -.0099503309 -.0009995003 -.0000999950 -.0000099999 2 -.0198026273 -.0019980027 -.0001999800 -.0000199998 3 -.0295588022 -.0029955090 -.0002999550 -.0000299995 4 -.0392207132 -.0039920213 -.0003999200 -.0000399992 5 -.0487901642 -.0049875415 -.0004998750 -.0000499987 6 -.0582689081 -.0059820717 -.0005998201 -.0000599982 7 -.0676586485 -.0069756137 -.0006997551 -.0000699975 Table 4. The Constants -ln(l + S lo''^) for m = 9 20 3.2 Logarithm to the Base 10 The algorithm for evaluating logi^(X) is Mery similar to that for evaluating ln(X). In fact, the same continued product normalization is used. Again, let X = Xq • 10*" where a is an integer and Xq^ [0.1, 1). Then log^Q(X) = log^Q(XQ) +a The value of 1o9iq(Xq) can, of course, be obtained from ln(XQ)/ln(10). But it can also be evaluated through a continued product. Since and |Xq • ^\^^{^ + S^IO"'') - 1| < 7/9 • lO"'", to m digits of accuracy Again, a table of the constants -log,Q(l + S. 10 ) is necessary. This time, the constants for k >_ {{m/2) + 0.9) are not as simple as those in the natural logarithm. For k >^ ((m/2) + 0.9), log^Qd + S^IO'"") = log^Q(e)-S^-10"''. However, given one complete set of constants, -log,^(l + S. 10 ), where 21 k, = m/2 + 0.9, the other sets of constants for k > k, can be easily formed by shifting the radix point (k - k, ) places to the left. Table 5 gives the necessary constants for m=9. Since the number of extra digits required to avoid rounding errors (^) is 1 for m=9, the constants in Table 5 have a 10-digit precision. The following algorithm for evaluating log-|o(X) is basically the same as the natural logarithm algorithm. Algorithm for the logarithm to the base 10 Objective: To find log,Q(X), where X > 0, to m digits. (I) Initialization: (i) Normalize X such that X = Xq • lo" with Xq G [0.1, 1). (ii) Lq = (II) For k = 0, 1, 2, ..., m do Step k (i) Find S. and R.^, by Step 1 of Algorithm N. (III) Final step: ^"aio'x' -^^^ '" Example 5. Logarithm to the base 10 Find log^Q{0. 03125) with m=9 (I) Xq = .3125 a = -1 ^0 = 22 k = k = 1 1 .3010299957 1 .0457574906 2 .4771212547 1 -.0413926852 5 .7781512504 2 3 4 5 6 7 -.0791812461 -.1139433523 -.1461280357 -.1760912592 -.2041199827 -.2304489214 ^k k = 2 k = 3 k = 4 k = 5 7 .0315170514 .0030507515 .0003041126 .0000304017 6 .0268721464 .0026136156 .0002606549 .0000260585 5 .0222763947 .0021769193 .0002172015 .0000217153 J .0177228767 .0017406616 .0001737525 .0000173721 3 .0132282657 .0013048417 .0001303079 .0000130290 2 .0087739243 .0008694587 .0000868676 .0000086860 T .0043648054 .0004345118 .0000434316 .0000043430 1 -.0043213728 -.0004340776 -.0000434273 .0000043429 2 -.0086001718 -.0008677215 -.0000868520 .0000086858 3 -.0128372247 -.0013000933 -.0001302688 .0000130286 4 -.0170333393 -.0017337128 -.0001736831 .0000173714 5 -.0211892991 -.0021660618 -.0002170930 .0000217142 6 -.0253058653 -.0026979807 -.0002604985 .0000260560 7 -.0293837777 -.0030294706 -.0003038998 .0000303995 Tabl -k. e 5, The constants -109,^(1 + S, -10 ) for m=9 23 k = 6 k = 7 k = 8 k = 9 5 .0000021715 .0000002171 .0000000217 .0000000022 4 .0000017372 .0000001737 .0000000173 .0000000017 3 .0000013029 .0000001303 .0000000130 .0000000013 2 .0000008686 .0000000869 .0000000087 .0000000009 T .0000004343 .0000000434 .0000000043 .0000000004 1 -.0000004343 -.0000000434 -.0000000043 -.0000000004 2 -.0000008686 - . 0000000869 -.0000000087 -.0000000009 3 -.0000013029 -.0000001303 -.0000000130 -.0000000013 4 -.0000017372 -.0000001737 -.0000000174 -.0000000017 5 -.0000021715 -.0000002171 -.0000000217 -.0000000022 Tabl -k e 5. The constants -log,^(l + S. -10' ) for m=9 24 (II) k k ^k •^k+l 4.1 2 -.0625 -.4771212547 1 -.625 -.4771212547 2 6 -.2875 -.5024271200 3 6 .1241375 -.5050251007 4 3 * .2413625863 -.5051553695 5 T -.5051510265 6 2 -.5051501579 7 4 -.5051499842 8 T -.5051499799 9 4 -.5051499782 *Recoding the leading digits of i-^c) into {5, ..., 5} gives fi' 7' R' ^ in* -Rg = .024136... = -0.24144 = 0.24T4 (III) log^Q(0. 03125) = L^q + a = -1 .505149978 Example 6. Logarithm to the base 10 Find log^Q(0. 0009765625) with m=9. (0.0009765625 = .03125^, compare with Example 5) (I) Xq = 0.9765625 a = -3 "-0 = ° 25 (II) k ^k+1 k+l -.0234375 1 -.234375 2 2 -.390625 -.0086001718 3 4 .078125 -.0103338846 4 T -.218828125 . -.0102904530 5 2 *-. 18832501 56 -.0102991388 6 2 -.0103000074 7 T -.0102999640 8 2 -.0102999553 9 3 -.0102999566 M-Rg) = = 0." 18832... = 0.21232..., The leading digits yield Sg to S,q. (Ill) log^Q(. 0009765625) = L^q + a = -3.010299957 26 4. Exponential to the Base e X The exponential e can be evaluated digit by digit with a scheme similar to that used in division and logarithm. The radix 2 algorithm defined by De Lugish [1] is adapted to radix 10 in the following discussion. Since log,^(e) • In(lO) =1, we have the identity: Let I = the integer part of X-log,^(e) and f = the fractional part of X-log,^(e) Equation (4.1) becomes then I + f = X-log^Q(e) gX ^ g(I+f)-ln(10) = gMn(lO) . gf-ln(lO) = 10^ . e'o X where X^ = f • In(lO). After the above manipulation, finding e is X mainly reduced to finding e , where X^ < In(lO) * 2.30. The bound for X ^ e becomes |e °| < 10. ^0 To find e , the following identity is used: 27 where Mj^ = 1 + S|^'10"'^ and Sj^ e {7, 6, ..., 6, 7}. Normalizing the quantity (Xq - 2._q ln(M. )) to will result in ® ^ k=0 "^k Equation (4.2) means that e can be evaluated in the form of m ~ k a continued product n (1 + S.-IO ), where the constants S. are chosen such that Xq - S^^q ln(M|^) - 0. 4.1 Finding the Selection Rules Let X^^^ = Xq - S^^Q ln(l + S.-lO'""), then X^^^ = X|^ - ln(l + Si^-lO"*^). Again, let one digit of X^ be normalized to at each step k. That is, \\+-^\ 1 7/9 . 10■^ (4.3) Let R, . 1 be the scaled remainder such that k+1 R^+i = lo'^EX^ - ln(l + S^-10"'')]. Then R,^^^ = lo'^ElO''^"^^ -R^ - ln(l + Si^-lO""^)] = lOR^ - lo'^-lnd + S|^-10"'^) (4.4) k -k To evaluate R^+i . the constants 10 • ln(l + ^k'^O ) are required. The constants are essentially the same as those used in the algorithm for natural logarithm. Table 6 provides the constants for m=9. k = 1 k = 2 28 k = 3 1 -12.0397280433 -7.2570692835 -7.0246149370 6 - 9.1629073187 -6.1875403718 -6.0180723256 5 - 6.9314718056 -5.1293294388 -5.0125418235 4 - 5.1082562377 -4.0821994520 -4.0080213975 3 - 3.5667494394 -3.0450207485 -3.0045090203 2 - 2.2314355131 -2.0202707318 -2.0020026707 l" - 1.0536051566 -1.0050335854 -1.0005003336 1 .9531017980 .9950330853 .9995003331 2 1.8232155679 1.9802627296 1.9980026627 3 2.6346426447 2.9558802242 2.9955089798 4 3.3647223662 3.9220713153 3.9920212695 5 4.0546510811 4.8790164169 4.9875415110 6 4.7000362925 5.8268908124 5.9820716775 7 5.3062825106 6.7658648474 6.9756137364 k = ^ 7.0024511440 6.0018007204 5.0012504168 4.0008002135 3.0004500900 2.0002000267 1.0000500034 .9999500033 1.9998000267 2.9995500899 3.9992002133 4.9987604165 5.9982007197 6.9976611427 k = 5 -7.0002450115 -6.0001800077 -5.0001250043 ■4.0000800027 ■3.0000450011 -2.0000200009 ■1.0000050004 .9999949996 1.9999800002 2.9999550004 3.9999200020 4.9998750036 5.9998200070 6.9997550108 k = 6 -7.0000245038 -6.0000180062 -5.0000125016 -4.0000080041 -3.0000045066 -2.0000020021 -1.0000005046 .9999994975 1.9999979950 2.9999954995 3.9999919970 4.9999874945 5.9999819991 6.9999754967 Table 6. lo"^ ln(l + Sj^lO""^) for m=9 29 From ln(l + Si^'lO""^) = Sj^lO'"^ - {S^^/2)'}Q'^^ + (S^'^/3)-10"^'^, ..., where e^ = lo'^((S^2/2).10"2'' - (S^^/3)-10"^'^ + ...). Since e. decreases as k increases, Equation (4.5) strongly suggests choosing S. by rounding lOR. to one digit: S|^| = the integer part of (|lORj + 0.5) (4.6) (sign of Sj^) = (sign of Rj^) Since \\+-\\ i. 7/9, e. in Equation (4.5) must satisfy |e,^| 1 7/9 - 1/2 = 5/18 for the selection rule (4.6) to be legally applied. For k >_ 3, |e. | is found to be less than 0.25 for S. e {T, 6^, ,.., 6, 7}. Therefore, S^, S., ..., S can be selected by rounding lOR. to a single digit as defined in (4.6). Again for k >^ ((m/2) + 0.9), ln(l + Sj^-lO"*^) = Sj^-lO"*^ to m digits accuracy. The recursion (4.4) simplifies to for k >_ ((m/2) + 0.9). Let k, be the smallest integer not less than ((m/2) + 0.9), then once R. is obtained, all the remaining digits S. , S. ^, , ..., S are available. 30 To find out the rules for selecting S^, S, , and Mq, the selec- tion intervals have to be studied. For a digit S. , the selection interval IS 10R,^e [R^^^ + lo'^-lnd + S^IO"''), R^^^ + lo'^-lnd + S^IO"*^)]. (4.8) where R.^, and R.^, are the lower and upper bounds of R.^, respectively, Table 7 shows the selection intervals for (S^}. The values used for R- and RT are -.749 and -749 respectively. ^2 7 6 5 4 3 2 1 T 2 3 4 5 6 7 Table 7. a < IOR2 1 b 6.1 7.5 5.1 6.5 4.2 5.6 3.2 4.6 2.3 3.7 1.3 2.7 0.3 1.7 -0.7 0.7 -1.7 -0.3 -2.7 -1.3 -3.7 -2.3 -4.8 -3.4 -5.8 -4.4 -6.9 -5.5 -8.0 -6.6 Selection of $2 31 All neighboring intervals overlap. However, the simple rounding rule (4.6) cannot be applied to 7. In order that S^ can be selected by rounding, Rp is restricted to the interval (-6.5, 7.5) and 7 is excluded from {Sp}. Using (4.8), the selection intervals for S, are found. Table 8 shows that S, cannot have the values 7, 6, 5, 4, and 3 because the corresponding intervals are not contiguous. Furthermore, S, can be selected by rounding lOR. if R, is restricted to the interval (-2.5, 3.3]. Restricting R, to this interval, the selection intervals for the initial step can be arranged as shown in Table 9. ^1 a < lOR^ < : b 4.7 6.0 4.2 5.4 3.4 4.8 2.8 4.1 2.0 3.3 1.2 2.5 0.4 1.7 -0.6 0.7 -1.7 -0.4 ;2.8 ___:]. 5. -4.2 -2.9 -5.7 -4.4 -7.5 -6.2 -9.8 -8.5 12.7 -11.3 Table 8. Selection Intervals for S, 32 In(MQ) M„ for m=9 1.8 < ''o < 2.31 2.0 7.3890560989 1.3 < "o < 1.8 ' 1.5 4.4816890700 0.8 < "^o < 1.3 1.0 2.7182818285 0.3 < "o < 0.8 0.5 1.6487212707 0.2 < "o < 0.3 0. 1.0000000000 0.7 < ''o < -0.2 -0.5 0.6065306597 1.2 < "o < -0.7 -1.0 0.3678794412 1.7 < ^0 < -1.2 -1.5 0.2231301601 2.2 < ''o < -1.7 -2.0 0.1353352832 2.31 < ''o < -2.2 -2.5 0.0820849986 Table 9. The initial step of choosing In(M^) With the initial step defined in Table 9, all S. for k >^ 1 can be selected by rounding lOR. to one digit as defined in (4.6). The normal i: satisfy normalization of Xq - ^^^^ ln(l + Sj^lO" ) to has been defined to Xq - S^^Q ln(l + S.-IO""") < lO""^ for k = 0, 1 , . . . , m. To m digits of accuracy, The initial manipulation of the operand and the evaluation of e IS summarized in the following. 33 4.2 Algorithm for Exponential Objective: Given X, find e to m digits. (I) Initialization: (i) N = X . log^Q(e) (ii) I = the integer part of N. (iii) f = the fractional part of N. (iv) k, = the smallest integer not less than ((m/2) + 0.9) (v) Xq = f . In(lO) Note that if |X| < InlO ~ 2.3, Xq = X. (II) Finding E^^^ (= e "): Step : (i) Choose In(M^) according to Table 9. (ii) Ri = Xq - In(MQ) (iii) E^ = Mq For k = 1, 2, . . . , m do Step k : (i) Select S. by: S|^ = the integer part of (|10R|^| + 0.5) (sign of Sj^) = (sign of Rj^) (ii) If k < k-j, then R|^+T = lOR,^ - lo'^-lnd + S^-10"'') else 3h (III) Final step: As explained in (2.7), R. and E. can yield E , directly to: Em+l = E.-d^R^;^ -10-^^1) (4.9) n where R. is R. rounded to m+l-k, places. K^ K^ I To provide for rounding errors, the intermediate results and the constants lOK • ln(l + S. -10 ) should have (m +Am) digits, where ^m >^ log,Q(m+l). The constants in Table 6 are used in the following examples, Example 8. Exponential (base e) Given X = -3.965735903 (from Example 3), X Find e with m=9. (I) N = -3.465735903 x log^Q(e) = -1.505149978 I = -1 f = -0.505149978 Xq = f • In(lO) = -1.16315081 k^ = 6 (II) Step In(MQ) = -1.0 Mq = -0.3678794412 R^ = Xq - In(MQ) = -0.16315081 E^ = 0.3678794412 35 k \ ■^k+l ^k.l 1 2 .5999274131 .294303553 2 6 .172383321 .3119617662 3 2 -.27415959 .3125856897 4 3 .2587551 .312491914 5 3 -.412399 .3125012888 E^Q = Egd + Rg X 10"^) = .3125012888 x .999995876 = 0.3125 (III) e^ = E .1 10 10' = 0.03125 Example 7. Exponential (base e) Given X = -6.931471806, V find e with m=9. (Compare with Example 4) Theconstants In(lO) and log,Q(e) are required, In(lO) = 2.302585093 log^Q(e) = 0.4342944819 (I) N = X • log^Q(e) = -3.010299956 I = -3 f = -0.010299957 Xq = f X In(lO) = -0.023165274 k^ = 6 (II) Step In(MQ) = 0, Mq = 1 R^ = -0.0237165274 E^ = 1.0 36 k h '^k+1 ^k.l 1 -.237165274 1.0 2 2 -.3513820083 .98 3 4 .4942013145 .97608 4 5 -.056737272 .9765804 5 T .43263728 .9765582743 E^Q = 1.000004326 x .9765582743 = .9765624989 e = E 10 10^ = .0009765625 37 5. Exponential to the Base 10 The general logic used in arriving at the continued product X X evaluation of e can be used for the case of 10 . As one might expect, the constants log,Q(l + S. -10 ) are required. It is also found that the digit selection is slightly more complicated in the algorithm for X finding 10 . Besides the special initial step, the constant S. is selected by a scaled rounding procedure for a 'normal' step. Let I = the integer part of the given operand X, and Xq = the fractional part of X. Then 10^ = 10^ + Xq I ^0 = 10 • 10 ^. ^0 > . Again, the problem is essentially reduced to finding 10 , where |X^| < 1. From the identity ^0 ^^0-^k=0 l°9io(Mk)) logioK=0^k^ , X 10 ^ = 10 ^ K U lU K . iQ lU K U K ^ (5 Ij Where M,^ = 1 + Sj^-lO""^, and S^ e {7, 6, ..., 6, 7}. If the expression ['] .m is normalized to such that "m+l = ^0 - ^r=0 l°9ioK) i 7/9 • 10-"' 38 Assuming lO" << 1 , 10 = e ^ (1 + io-'")ln(10) « (1 + In(lO) • lO"'") < (1 + 2.31 • 10'^) From (5.1) rm 10 ^ =10 . 10 '^ "^ '-^ "^ = lo^'"'^^ . n"^ M then rr_„ M. approximates 10 to (m-1) digits. For m digits of accuracy, one additional step would be necessary. That is, X k=0 10 ° = n'J;^!! (1 + S. -lO'"^) is m digits. 5.1 Finding the Selection Rules Let X^^^ = Xq - S^^Q log^Qd + S.-IO"'*), .(5.2) and R^^^ = lo"" • X^^^ (5.3) The constants S. are selected in such a way that Vll ^ 7/9 (5.4) for k = 0, 1 , .. . , m. From (5.2) and (5.3), Rk+i = lOR^ - lo"" • log^gd + S^-10"'') (5.5) 39 Since S^ . 10"'^ < 1, ln(l + S^-IO'"^) = S^-IO"'^ - (S,^^/2) . lO"^*^ + (S^^/3) • 10"^'^ (5.6) - higher order terms Let e,^ = 10^(-S^2/2) • lO'^"" + (S^^/3) • lO"^"^ - ...) • log^^Ce) (5.7) then R^^^ = lOR^ - log^Q(e) • ln(l + S^-10"^) = lOR^ - log^Q(e) • S,^ - e^ (5.8) From (5.8), it seems that a scaled rounding selection rule for S. would be appropriate. Usually, S. is selected if lOR. lies in the interval (S. - 0.5, S. + 0.5). If the selection interval is scaled by log,->(e), S. is selected when lOR. is the interval [log,^(e) • (S. - 0.5), logi^(e) • (S. + 0.5)]. Figure 3 defines a selection rule for S. . This rule is based on the scaled intervals. Before discussing the applicability of the selection rules in Figure 3, let us look at the properties of this rule. S. can be selected if and only if R. lies in the interval [-.325, .325]. Furthermore, the quantity (lOR. - log-,^(e)-S. ) is always less than 0.27. This can be easily verified by comparing the value of log,^(e)-S. to the bounds of the selection intervals for S. e {7, 6, ..., 6, 7}. Returning to the recursion (5.8), where e. may be expressed as ^k = ^^^ ' ""o^io^e) • (ln(l + S^.IO"*^) - S^-IO""^). i+0 ,325 ,285 ,24 ,20 ,15 ,11 ,07 ,02 ,00 lO.R.I 3.25-r 2.85-- 2.4-- 2.0-- 1.5-- 1.1-- 0.7-- 0.2-1- 0.0 ) } } } ) } ) I) 7 6 5 4 3 2 1 Figure 3. Intervals for Selecting S. (i) (sign of S,^) = (sign of R,^) (ii) Select |S, | according to |R. | For k >^ 3 and the chosen digit set |S. | , the magnitude of e. is found to be less than 0.011. This means for k > 3 'k+1 lOR,^ - log^Q(e)-S|^ - e,^| < .27 + 0.011 = 0.287 Since the selection rule in Figure 3 only requires Rj^^^l <_ 0.325. (5.9) it can be applied to select S. for k >_ 3. For the selection of S„, the allowable selection interval for a particular digit S^ is lORp in [R3 + lO^.log^^d + S^-IO"''). R3 + lO^.log^Qd + S^IO""")], 41 (R. and R. denote the lower and upper bounds of R. respectively.) The following table lists the allowable intervals for the selection of each of the digits in {S^}. From (5.9), the values for R- and RT are -.325 and .325 respectively. ^2 7 6 5 4 3 2 1 T 2 3 4 5 6 7 a £ IOR2 1 b 2.7 3.25 2.3 2.85 1.8 2.4 1.4 2.0 1.0 -1.6 0.6 1.1 0.2 0.7 -0.3 0.3 -0.7 -0.2 -1.2 -0.6 -1.6 -1.0 -2.0 -1.5 -2.5 -2.0 -3.0 -2.4 -3.4 -2.83 Table 10. Selection of Sp k2 S„ is defined for |Rp| £0.325, and the selection intervals of all digits in {S^} enclose the corresponding intervals defined in Figure 3. Therefore, the selection rule in Figure 3 can also be used to select Sp. For iRpI £0.325, the selection intervals of {S,} are those listed in Table 11. 7 6 5 4 3 2 1 T 2 3 4 5 6 7 Table 11. Selection Intervals for {S,} a < L lo-Ri 1 b 2.0 2.6 1.8 2.3 1.5 2.0 1.2 1.7 0.9 1.4 0.5 1.1 0.1 0.7 0.3 0.3 0.7 -0.2 1.29 -0.7 1.8 -1.23 2.5 -2.0 3.3 -2.7 4.3 -3.7 5.5 -5.0 i+3 The digits 4, 5, 6 and 7 cannot be used for S, because the corresponding intervals are not contiguous. Comparing the intervals shown in Table 11 with the intervals shown in Figure 3, the scaled rounding rule in Figure 3 can only be applied to R, in the interval [-0.11, 0.14]. It is found that given |Xq| < 1 , X^^ can be easily mapped into R, in this interval. Therefore, S, can be selected by the rule defined in Figure 3. Since R-i = ^^ - log^QlMg), the following choice of logio(MQ) will result in R, falling in the interval [0.1, 0.1], which is the sub- interval of the desired interval. 1.0 - 0.9 - 0.7 -- 0.5 -- 0.3 " ■0.1 -- ■0.1 -- ■0.3 -- •0.5 -- •0.7 -- •0.9 -' •1.0 - log,,(M(,) "o Mq for m=9 } 0.9 eO-9 7.9432823472 } 0.8 eO.8 6.3095734448 ) 0.6 eO.5 3.9810717055 } 0.4 eO.4 2.5118864315 } 0.2 3O.2 1.5848931925 } 0.0 1 1.0000000000 ) -0.2 e-0.2 0.6309573445 } -0.4 e-0.4 0.3981071706 } -0.6 e-0.6 0.2511886432 } -0.8 e-0.8 0.1584893192 } -0.9 e-0.9 0.1258925412 Figure 4. Selecting log,Q(MQ) Log,Q(MQ) is selected according to the given X-. kh The selection rules for the m+2 steps of the normalization have been found. There is one special initial step and (m+1 ) uniform steps. The normalization of Xq - ^^ ^^^lO^'^k^ ^"^ ^^^ evaluation of Y 10 are summarized next. 5.2 The Algorithm for Exponential to the Base 10 Y Objective: To evaluate 10 to m digits. (I) Initialization (i) I = the integer part of X. (ii) X^ = the fractional part of X. (II) Evaluating 10 ^ (= E^^2) Step (i) Select "IoQiqCMo) according to Figure 4. (ii) E^ = Mq (iii) Ri = Xq - log^Q(MQ) For k = 1 , 2, ...» m+1 do Step k (i) Select S. as indicated in Figure 3. (ii) R^^^ = lOR^ - lo'^.log^Qd + S^.IO"''). (iii) E^+l =^k-(^ ^S^-IO"^). (Ill) Final step 10^ = E -10^ lU t^^^ lU Since the last few steps are not simplified as in the previous algorithms this algorithm is not as efficient as any of them. The following k -k examples illustrate this. The constants 10 • log,Q(l + S. -10 ) for m=9 are given in Table 12. i+5 k = 1 k = 2 k = 3 7 -5.2287874528 -3.1517051446 -3.0507515046 6 -3.9794000867 -2.6872146400 -2.6136156027 5 -3.0102999566 -2.2276394711 -2.1769192543 4 -2.2184874962 -1.7728766960 -1.7406615763 3 -1.5490195999 -1.3228265734 -1.3048416883 2 - .9691001301 - .8773924308 - .8694587126 T - .4575749056 - .4364805402 - .4345117740 1 .4139268516 .4321373783 .4340774793 2 .7918124605 .8600171762 .8677215312 3 1.1394335231 1.2837224705 1.3009330204 4 1.4612803568 1.7033339299 1.7337128090 5 1.7609125906 2.1189299070 2.1660617565 6 2.0411998266 2.5305865265 2.5979807199 7 2.3044892138 2.9383777685 3.0294705536 ^k k = 4 k = 5 k = 6 7 -3.0411258916 -3.0401677805 -3.0400720152 6 -2.6065489343 -2.6058450678 -2.6057747114 5 -2.1720154587 -2.1715266982 -2.1714778389 4 -1.7375254559 -1.7372126723 -1.7371814038 3 -1.3030789173 -1.3029029895 -1.3028854029 2 - .8686758343 - .8685976501 - .8685898333 T - .4343161981 - .4342966535 - .4342947010 1 .4342727686 .4342923103 .4342942637 2 .8685021165 .8685802780 .8685880930 3 1.3026880523 1.3028639026 1.3028814912 4 1.7368305846 1.7371431849 1.7371744520 5 2.1709297223 2.1714181243 2.1714559785 6 2.6049854739 2.6056887214 2.6057590737 7 3.0389978481 3.0399549759 3.0400507317 Table 12. lo"^ loQio^^ ^ ^k'^^'"^^ ^°^ ^"'^ he k = 7 k = 8 k = 9 7 -3.0400624420 -3.0400617313 -3.0400631784 6 -2.6057677036 -2.6057671852 -2.6057671148 5 -2.1714729778 -2.1714726436 -2.1714741375 4 -1.7371782953 -1.7371781061 -1.7371780748 3 -1.3028836564 -1.3028835731 -1.3028850984 2 - .8685890608 - .8685890444 - .8685890365 T - .4342945087 - .4342945200 - .4342960610 1 .4342944344 .4342942071 .4342929747 2 .8685888562 .8685887184 .8685859490 3 1.3028832346 1.3028832254 1 . 3028820085 4 1.7371775696 1.7371777280 1.7371749820 5 2.1714718612 2.1714722263 2.1714710408 6 2.6057661093 2.6057667203 2.6057640134 7 3.0400502831 3.0400612099 3.0400500713 Table 12. lo"^ ^OQio^^ "^ Si^-lO'*^) for m=9 Example 9. Exponential (Base 10) Find 10-^-505149978 ^.,,^^9 X = -1.505149978 (I) I = -1 Xq = -.505149978 h7 (Compare with Example 5.) (II) step log^C ,(«o) = -.6 ^1 = ■ "o = .2511885432 ^1 = ■■ h - log^o(MQ) = .094850022 k \ '^k+l ^k.l = ¥' ' \'''''^ 1 2 - .1566877595 .3014263718 2 4 - .136456334 .3134834267 3 3 - .059721652 .3125429764 4 T - .162900322 .3125117221 5 4 0.108209453 .3124992217 6 2 0.213506437 .3124998467 7 5 - .036507491 .3125000029 8 T .07021961 .3124999998 9 2 - .1664129365 .3125000004 10 4 .3125000003 (III) lO'' = E 11 = .03125 10 -1 kQ Example 10. Exponential (Base 10) Given X = -3.0102999566 find 10 . (Compare with Example 6.) (I) X = -3.0102999566 I = -3 Xq = -0.0102999566 (II) Mq = 1. logTo(MQ) = ^1 = = 1 «1 = = -.0102999566 k \ "^k+l 1 - .102999566 2 2 - .1526032292 3 4 .214629284 4 5 - .024636882 5 T .1879278335 6 4 .142103883 7 3 .118155596 8 3 - .121327275 9 3 0.0896106957 2 1 "k+1 .98 .97608 .97656804 .9765582743 .9765621806 .976524735 .9765625028 .9765624999 .976525001 (III) 10^ = E^Q • 10""^ = 0.0009765625 h9 6. Square Root One way to use a continued product to find the square root of a positive X is as follows: Let X = Xq • lo", where a is an even integer and Xq is in the interval [0.01, 1). Xq is then normalized to 1 by multiplying to it a sequence of factors (J^ ), where i = 0, 1, . . . , m. If then n^.o(Ji') = ^^h and 'l^oJi ' 'f^o- Therefore, The factor J. used in (6.1) is defined to be (1 + (l/2)-S.-10'"') with S. in the digit set {9, 8, ..., 8, 9}. The next section will develop a normalization procedure so that IXq n^^^Q J.^ - 1| < lO'"^ for k = 1, 2, ... . 50 For k=m. That is. "_ 3, it is less than 1/5. Therefore, S^, S-, ..., S can be selected by rounding lOR. to one digit as defined in (6.8). Furthermore, the last three terms in the recursion (6.7) would not affect the most significant (k-3) digits of R^+i- This means for k >_ (m+3)/2, the recursion (6.7) becomes simply Like most of the algorithms presented before, about half of the constants are found by simple recoding. Since the rules for selecting S. for k >^ 3 are settled, the selection intervals for S^ will be studied. The selection interval for a constant S. is ^.1 - \ - (^/^)\' • ^0"' Vi ■ h - ^'^'^\' ' 10'' < lOR. < (1 + 1/2 . S|^ . 10""^)^ ^ (1 + 1/2 • S^ . lO"*^)^ where R.^-, and R.^, are the lower and upper bounds of R.^, respectively. Using -0.95 and 0.95 for R- and rT, the selection intervals for (Sp) are computed and listed in Table 13. 53 From Table 13, all digits of {9, 8, ..., 8, 9} can be used for $2, but the simple rounding rule (6.8) cannot be applied to the digit 9^, which is thus excluded from {Sp}. Therefore, R« is restricted to the interval [-0.929, 0.85). The selection intervals for S, are listed in Table 14. All intervals are continguous, but the digits 6, 5, ...» 8, 9 are sufficient for the possible values of R, . From the table, R-, should be in the interval [-0.56, 1). Figure 5 defines one way of selecting S, . ^2 9 8 7 6 5 4 3 2 1 T 2 3 4 5 6 7 8 9 Table 13. Selection Intervals for {S^} a J i 10R2 < b -9.2 -7.6 -8.4 -6.7 -7.5 -5.8 -6.6 -4.9 -5.7 -4.0 -4.7 -3.0 -3.8 -2.1 -2.8 -1.1 -1.9 -0.1 -0.9 0.9 0.1 1.9 1.1 2.9 2.1 4.0 3.2 5.1 4.2 6.1 5.3 7.2 6.4 8.4 7.5 9.5 8.7 10.6 5h h 9 8 7 6 5 4 3 2 1 T 2 3 4 5 6 7 8 9 a 1 lOR^ < b -5.6 -4.9 -5.3 -4.5 -5.0 -4.1 -4.6 -3.6 -4.1 -3.1 -3.7 -2.5 -3,1 -1.8 -2.5 -1.1 -1.7 -0.2 -0.9 0.8 0.1 2.0 1.2 3.3 2.6 5.0 4.2 6.9 6.2 9.2 8.6 12.1 11.5 15.6 15.2 20.1 20.0 25.8 Table 14. Selection Intervals for {S,} To complete the normalization procedure, the initial step is 2 defined in Figure 6. Since ^i = ^^ • Jq - 1, the selection of Jq as indicated in Figure 6 will ensure that R-, will stay in the desired interval [-0.56, 1). 55 «1 Intervals lOR^ 0.10 1.0 0.87 8.7 0.66 6.6 0.47 4.7 0.30 3.0 0.16 1.6 0.05 0.5 0.05 -0.5 0.15 -1.5 0.20 -2.0 0.25 -2.5 0.33 -3.3 0.39 -3.9 0.43 -4.3 0.47 -4.7 0.50 -5.0 0.57 -5.7 of _6_ S_ .?_ I_ 2 'l _0_ 1 Figure 5. Selection of S Choose S, according to R, Intervals of Xq 1.00 r 56 0.50 0.20 0.10 0.04 0.01 (Jq)' 1 1 2 4 3 9 4 16 !_____ _._49 Figure 6. The Initial Step Select Jq according to X... The entire normalization of Xq to 1 through the continued product Xq • n'l^^^d + (l/2)-S^-10"'^) has been defined. The initial step is specialized and the rule for selecting S-, is also different from the selection rules of other constants S, for k ^ 2. For k ^ 2, S, is selected by rounding lOR. to one digit. The recursion formula for R.^-, simplifies to R.^, = lOR. + S, for k >; {m+3)/2. The normalization results in |Xq •n'J=o^^ ■" (V2)-s^.io-'^)2 - i| < 10""' (6.10) in m+1 steps. From (6.3), ^ = Xon^=o(^ "■ (V2)s^-io-^) to m digits. 57 6.2 The Algorithm for Square Root Objective: For a given X > 0, find •Y to m digits. (I) Initialization Normalize X such that X = X^-IO where a is an even integer and X^ lies in the interval [0.01, 1). (II) Finding /X^ (= T^^^ ) Step (i) For the given X^, find J^ from Figure 6. (ii) R^ = Xq-Jq^ - 1. (iii) T^ = Xq-Jq Step 1 (i) Choose S, according to Figure 5. (ii) R2 = 10R^-(1 + (l/2)-S^-10"b^ + S^ + (l/4)-S^^-10'^ (iii) T2 = T^-(l + (l/2)-S^-10"'') For k=2, 3, ...,mdo Step k (i) Select S. by |S|^| = the integer part of (|10R. | + 0.5) (sign of Sj^) = -(sign of R,^) (ii) T^+1 = T^-(l + (l/2)-S^-10"'') (iii) If k > (m+3)/2 then R^^^ = lORj^d + (1/2)S,^-10"'')^ + S^ + (^/4).s^2•lo■^ else R^^^ = lOR^ + S^. 58 (III) Final step Actually, the algorithm can be further simplified as in division. Let k, be the smallest integer not less than (m+3)/2. Once R. and S. are found, there is no need to carry out the rest of the digit selection. J _^_■, can be found from the following formula: T..1 = \ ■ (1 - iy2)\^o-^"^) (6.11) n where R. is R. rounded to (m+l-k) significant digits. The following examples illustrate the algorithm for evaluating square roots and the simplification stated above. Example 11. Square Root Given X = 0.0009765625, find yx (I) a = -2 Xq = 0.09765625 (II) Step Jq = 4, j/ = 16 ^1 " ^O'^O^ ■ ^ " 0-5625 T^ = 0.390625 Step 1 ^1 = ~' R^ = 0.0 T = T^(l + (l/2)S^.10'b = 0.3125 Since R^ = 0.0, all remaining S., k = 3, ...» 9 are zero. .-. T^Q = 0.3125 .-. ^/l = T^o ^ ^°"^^ " 0.03125 59 Example 12. Square Root Find the square root of 75.41926 to 9 digits (I) m = 9 a = 2 X = 0.7541916 (II) Step ^0 ^ 1, R^ step 1 1 ^1 = 3, R2 Steps 2 to 5 = -.2458085 = -.0258160900 S^ = 0, R3 = -.2581609000 S3 - 3, R^ h - ^' ^5 h - ^' ^6 .4128903644 .1276522977 .2765127116 ^-5> From (6.11), "^10 = "^6^^ ^ (l/2)-(0.2765)-10"^) = Tg(l - 0.0000013825) = .8684420534 (III) 75.41926 = T^Q-10^ = 8.68442053 7541926 ,86732034 ,86732034 8686213205 ,8684475962 ,8684432540 60 7. Summary Continued product (sum) algorithms for radix 10 division, logarithm, exponential and square root have been studied. Each algorithm consists of two parts: the normalization of a quantity to a constant and the evaluation of the value of the function. The normalization process generates a sequence of digits S^, S, , Sp, ..., one digit in each step. The function evaluation process uses these digits in the same sequence to generate the result, yielding one more digit of accuracy in each step. One can either perform both processes in parallel or perform the entire normalization process first. In most of the algorithms, about half of the steps can be simplified. The total number of steps is 1 or 2 more than the number of significant digits desired. The use of redundant representation simplifies digit-selection which then depends only on the most significant positions of the partially normalized quantities. All of the algorithms presented can be carried out in a paper- and-pencil environment as long as the necessary tables of constants are available. Even in the absence of an adding apparatus, the algorithm for natural logarithm is still practical. 61 References [1] B. G. De Lugish, "A class of algorithms for automatic evaluation of certain elementary functions in a binary computer," Report 399, Department of Computer Science, University of Illinois, Urbana, June 1970. [2] M. D. Ercegovac, "Radix 16 division, multiplication, logarithms, and exponential algorithms based on continued product representation," Report 541, Department of Computer Science, University of Illinois, Urbana, August 1976. [3] J. E. Robertson, Lecture notes for computer science course 364, Fall 1976, Department of Computer Science, University of Illinois, Urbana. [4] J. E. Robertson, "A method of decimal division based on the use of continued products," private communication, December 1973. BIBLIOGRAPHIC DATA SHEET 1. Report No. UIUCDCS-R-78-920 3. Recipient's Accession No. 4. Title and Subtitle CONTINUED-PRODUCT ALGORITHMS FOR DIVISION, LOGARITHM, EXPONENTIAL AND SQUARE ROOT IN RADIX TEN 5. Report Date May 1978 6. 7. Author(s) Catherine Yuk-Fun Chow 8< Performing Organization Rept. No. 9. Performing Organization Name and Address Department of Computer Science University of Illinois at Urbana-Champaign Urbana. IL. 61801 10. Project/Task/Work Unit No. 11. Contract/Grant No. NSF MCS 77-18642 12. Sponsoring Organization Name and Address National Science Foundation Washington, D.C. 13. Type of Report & Period Covered 14. 15. Supplementary Notes 16. Abstracts This paper studies radix 10 continued product/sum algorithms for division, logarithm, exponentiation, and square root. A convergent continued product/stim series transforms the evaluation of a function to a set of recursive steps. An evaluation step typically is a summation of 2 or 3 numbers. A redundant number representation (signed-digit) simplifies the digit selection in each recursion. Convergence of the result is linear: one additional digit of accuracy is obtained in each recursion. 17. Key Words and Document Analysis. 17a. Descriptors continued products/sums, radix 10, division, logarithm, exponential, redundant number representation, linear convergence, signed digits, digit selection intervals, normalization, scaled remainder 17b. Identif icrs/Open-Ended Terms 17c. COSATI Field/Group 18. Availability Statement 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 22. Price FORM NTIS-35 (10-70) USCOMM-DC 40329-P7 1 M 2 1978 m :3