Division and square root for mobile and scientific computing markets

UNIVERSIDADE DE SANTIAGO DE COMPOSTELA DEPARTAMENTO DE ELECTRONICA E COMPUTACION DOCTOR OF PHILOSOPHY DISSERTATION Division and Square root for Mobile and Scientific computing markets VIJAYKUMAR HOLIMATH October 2007 Dr. Javier D´ ıaz Bruguera at the Department of Electronics and Computer Engineering in the University of Santiago De Compostela. Hereby Certifies: Research work entitled ”Division and Square Root for Mobile and Scientific Computing Markets” has been developed by Vijaykumar Holimath under my supervision in the Department of Electronics and Computer Engineering of the University of Santiago de Compostela, qualifies him for the PhD degree. Santiago de Compostela, October 2007 ————- Certifica: Que el trabajo de investigaci´on titulado ”Division y Ra´ız Cuadrada para Aplicaciones M´oviles y Cient´ıficas” ha sido realizado por D. Vijaykumar Holimath bajo mi supervisi´on en el Departamento de Electr´onica y Computaci´on de la Universidad de Santiago de Compostela, y constituye la Tesis que presenta para optar al grado de Doctor. Santiago de Compostela, October 2007 Prof. Javier D´ ıaz Bruguera, Director of this PhD dissertation. Javier D´ ıaz Bruguera, Chair of Department of Electronics and Computer Engineering. Vijaykumar Holimath, PhD Candidate. Acknowledgment First of all I would like to thank my supervisor, Prof. Javier D´ıaz Bruguera. I have learned many techniques working with him and his advice, ideas and guidelines have been invaluable. I would also like to thank the people I met during my visits, for their assistance in doing the research supported by my supervisor: Prof. Naofumi Takagi (University of Nagoya, Japan), Prof. David Matula (Southern Methodist University, USA) and Prof. Jean-Michel Muller (ENS-Lyon, France). I would also like to thank all the people in the Department of Electronics and Computer Engineering and the Computer Architecture Group for their help and cordiality. The time spent in Santiago De Compostela has been very productive, both from a professional and personal perspective. Finally, I would like to thank my parents and my sisters (Beena and Rashmi) whose support enabled me to stay in Spain. Santiago De Compostela, October 2007 Contents Objectives 1 1 Introduction 7 1.1 Number Representation . . ..................... 7 1.2 Rounding . . ............................ 9 1.3 Special values ............................ 13 1.4 Exceptions . . ............................ 16 1.5 Importance of Division and Square root Functions . . ....... 16 2 Classification of Hardware Methods 19 2.1 Non-iterative Methods . . ..................... 20 2.1.1 Table-based methods . . . . . . .............. 20 2.2 Iterative methods . . . . . ..................... 25 2.2.1 Goldschimdt algorithm . . . . . .............. 25 2.2.2 Newton-Raphson Algorithm . . .............. 28 2.3 Digit-recurrence Algorithm ..................... 30 2.3.1 Digit-recurrence division . . . . .............. 30 2.3.2 Accurate Quotient Approximation method . ....... 34 i ii CONTENTS 3 Modified N-R Algorithm without a LUT 39 3.1 N-R Algorithm . .......................... 40 3.1.1 Algorithm .......................... 41 3.2 CS Booth digit representation ................... 43 3.3 Error Computation . . . . . . ................... 49 3.4 Architecture and Performance ................... 51 3.5 Modified N-R square root algorithm without a LUT ........ 56 3.6 Comparisons . . .......................... 58 3.7 Conclusion . . . .......................... 60 4 Modified AQA Method Using a LUT 63 4.1 Proposed Method .......................... 64 4.2 Initial Approximation . . . . ................... 70 4.2.1 Error Computation . . ................... 73 4.2.2 Goldschimdt Algorithm . . . . . . ............ 76 4.2.3 Error computation . . ................... 76 4.2.4 Architecture of Minimax Approximation . . . . ..... 78 4.2.5 Delay Estimates of Initial Minimax Approximation . . . . 81 4.2.6 Architecture and Performance of initial approximation . . 83 4.3 Algorithm .............................. 88 4.4 CS to Redundant Booth digit representation . . . . . . ..... 90 4.5 Error Computation . . . . . . ................... 95 4.6 Architecture and Performance ................... 97 4.7 Comparisons . . ..........................104 4.8 Conclusions . . . ..........................106 CONTENTS iii Conclusions 109 Future Work 115 Resumen en Espa˜ nol 117 Bibliography 123 x SD Sign Digit SRT Sweeney Robertson Tocher UPC Ultra Personal Computer VLSI Very Large Scale Integration Abstract Division and square root functions are the basic functions in many floating point applications. The uses of these functions are segmented. For example, they are frequent in 3D graphics animation, virtual reality (especially lighting effects), view port transformation, motion synthesis, signal processing, scientific computing applications etc. These types of applications are suitable for the high performance computing market. Division and square root functions are used less in internet applications, business applications such as audio and video streaming, charting, presentation, snap shot viewing, day to day accounting etc. Most of these applications are suitable for the mobile computing market such as PDA,UPC, tablet PC. In both ofthese markets, the siliconarea is crucial and inthe highperformance market speed is the challenge. In state-of-the-art processors, the silicon area may be regarded as expensive and the speed is to some extent acceptable. Most designers are willing to sacrifice the speed in favor of low cost and complexity. We therefore aim to develop efficient algorithms both in terms of area and the speed. We propose a modified Newton-Raphson reciprocation algorithm for single precision computation, without using a look-up table. We regard this algorithm as the best silicon efficient algorithm and the speed is average compared to the state- of-the-art processors. This algorithm could be suitable for the mobile computing market. xi xii It is a variable latency algorithm, requiring a maximum of 22 iterations, which provide a linear convergence. Multiplying with the dividend, results in division. Our methodrequires a cycle for each iteration, performingmultiply addand Booth recoding in one cycle. Initial approximation is a two’s complement of the divisor, which can be performed during partial products summation. Multiples of the divisor are constant throughout iteration and are generated only once, in the second iteration. To evaluate the proposed method, we used a radix-4 multiplier and a radix-8 multiplier. We then compared it with conventional processors. The comparison shows that it offers a good trade-off between performance and area, making it suitable for mobile computing applications such as PDA,UPC, mobile phones, tablet PC, etc. We also propose a division algorithm using a modified Accurate Quotient Approximation Method for double precision computation and this is extended to a square root algorithm,with a look-uptable. It is a better silicon efficient algorithm with a significantly higher performance than the state-of-the-art processors. This algorithm can be suitable for both mobile and scientific computing markets. The proposed method employs a second degree minimax approximation to obtain an initial 15 bits estimate of reciprocal and inverse square root values. This is then passed through an iteration of Goldschimdt reciprocal and inverse square root algorithms, to enhance the initial approximation to 30 bits and to make the initial approximation more silicon efficient. Following this, we perform an iteration of the Accurate Quotient Approximation method to obtain double precision results. To speed-up these functions, a high seed 30 bits is used and intermediate results are computed in carry-save form. To evaluate the proposed method, we used a radix-4 multiplier and a radix-8 multiplier and then compared it with conventional processors. This comparison shows that there is a significant improvement in performance, with better silicon efficiency. It offers some advantages over N-R method, GLD algorithm, SRT method and Cyrix algorithm, if an iteration algorithm is considered. It is xiii expected that, with the proposed method, processor architects will have an option for their next generation processors. xiv List of Publications International Conferences •V. Holimath and J. Bruguera, A Linear Convergent Functional Iterative Division withouta Look-up Table, Proc 9th EUROMICRO Conference on Digital System Design (DSD), pages 236-239, Dubrovink (Croatia), Sept 2006. Journals and Transactions •V. Holimath and J. Bruguera, Division for Mobile Computing Applications, submitted to Transactions on High-Performance Embedded Architectures and Compilers (HiPEAC), 2007. •V. Holimath and J. Bruguera, Divisionand Square root using Accurate Quotient Method, submitted to IEEE Transactions on Computers, 2007. •V. Holimath and J. Bruguera, Division and Square root for Mobile and Scientific computing Markets, submitted to Journal of the ACM, 2007. xv xvi Objectives Modern applications comprise several floating point operations, among them addition, multiplication, division and square root. The division and square root functions are less frequent than the two basic arithmetic operations (+/-,×) - only somewhat less than 0.8% of the total number of instructions executed [?] in general purpose applications. The poor performance of many processors when computing these operations can cause an overall execution time that is comparable to the time spent performing addition and multiplication. While the methodology for designing efficient high-performance adders and multipliers is well understood, the design of division and square root still remains a serious challenge, often viewed as a Black Art among system designers [?]. They are frequent in 3D graphics animation, virtual reality [?,?,?,?,?,?,?, ?,?,?] (especially lighting effects), view port transformation, motion synthesis, digitalsignalprocessingapplicationsandscientific computations. Thecost intime of the division and square root function depends on how smartly the compiler can schedule, in order to maximize the distance between the production of the results of other functions and the dependency with respect to division/square root. These functions can be expressed as the product of the dividend and reciprocal of the divisor, Q=X×(1/Y )and Q=Y×(1/√Y). There are several techniques described in the literature [?] and carried out in state-of-the-art processors such as N-R method, GLD algorithm [?,?,?], SRT algorithm [?,?]. The performance of these functions depends on the precision; the precision depends on the market targeted. For example, in 3D graphics single precision is used and 1 2 in scientific computing applications double precision is used. To speed-up the division and the square root functions, somewhat more than ten bits of initial approximation of reciprocal and inverse square root are used. The values of the reciprocal and the inverse square root are stored in ROM also called LUT in terms of bits. Another option to obtain the initial approximation is, storing the coefficient values of linear/quadratic approximation algorithms of reciprocation and inverse square root in the LUT’s. Then, the values of the coefficients are added together or the coefficient values are multiplied by the respective input operands of the approximation algorithms and then added. The linear/quadratic approximation algorithms reduce the area for the LUT size. To obtain the required precision, the LUT’s or linear/quadratic approximation algorithms are coupled with the N-R method, GLD algorithm or SRT method. The N-R and GLD methods are quadratic convergence in character and the SRT method linear convergent. On average, the speed of these operations is more than 20 cycles in either single precision or double precision computation in state-of-the art processors. To enhance performance according to market demands, there is a trade-off between the performance and silicon area. The area for these functions accumulates from LUT’s. The more bits of initial approximation that are used, the more silicon area is required to speed-up these functions. This means that the silicon area is directly proportional to the market demands. Most designers are willing to sacrifice the speed in favor of low cost and design complexity. In the state-of-the-art processors such as Intel, AMD, IBM, the silicon area for the initial approximation may be regarded as expensive and the performance of these functions is to some extent acceptable. These processors are targeted at various markets such as the mobile computing market, the desktop computing market, the scientific computing market etc. In this thesis, we therefore aim to develop an algorithm that is: a) a best silicon efficiency algorithm and b) a better silicon efficient and best performance 3 algorithm: 1. A modified N-R reciprocation/division algorithm without using a LUT for single precision computation is proposed. It is a best silicon efficient algorithm and the performance is about the same as the average performance of the state-of-the-art processors. 2. The double precision computation of AQA method for both division and square root, using a LUT is proposed. It is a better silicon efficient algorithm and the performance is significant better compared to the state-of-the- art processors. We have chosen this approach to provide an option for processor architects that suits their target market. The modified N-R algorithm is suitable for the mobile computing market. In this market area is crucial. Most architects look for a division algorithm with a low cost, a reduced area with acceptable performance that can run internet applications, business applications (such as audio and video streaming), charting, presentations, snapshot viewing, day to day accounting etc. The modifiedAQA method is suitable for the high performance market. In this market the performance is the key. Most architects look for high performance division and square root algorithms with acceptable cost and complexity that can be applied to scientific computingapplications, 3D graphics, signal processing applications etc. The modified AQA method is also suitable for the mobile computing market. In the modified N-R algorithm, we achieve our objective by taking an initial approximation as a two’s complement of the divisor and modifying some steps of the N-R algorithm [?]. This initial approximation can be performed during partial products summation of the multiplicand (divisor) and it is silicon efficient. The average performance is maintained by executing a multiply add operation and Booth recodingin thesame cycleand keeping the divisormultiples constantin 10 Chapter 1. Introduction point number. On the other hand, the Xis inside the range of the floating point number, which is 1≤X<2for the IEEE standard. We decompose Xinto two parts Xfand Xdsuch that X=Xf+Xd×r−f(1.6) with 0≤Xd<1.Xfhas the precision of the significand in the floating point system and Xdrepresents the rest. IEEE standard-754 defines four different rounding modes, which must be supported: Rounding to Nearest (RNE)(Unbiased, Tie to even) value and to the even number when the result is exactly halfway between two machine numbers. In the RNE mode the value represented is the closest possible to the exact value (infinite precision represented by X), it produces the smallest absolute error. In terms of the operation on the infinite precision significand, RNE can be described as follows: RNE(X)=⎧ ⎨ ⎩ Xf+r−fif Xd≥1 2 Xfif Xd<1 2 (1.7) The RNE consists of adding r−f/2to the infinite precision significand and keeping the resulting ffractional digits. i.e. RNE(X)=((X+r−f 2)rf)r−f(1.8) For Xd≥1 2, the addition of r−fcan produce a significand that cannot be represented (significand overflow). In such a case, the resulting significand is multiplied by β−1and the exponent 1.2. Rounding 11 incremented by 1. The absolute error is ABRE[RNE]=⎧ ⎨ ⎩−Xdr−f×βEif Xd<1 2 (1 −Xd)r−f×βEif Xd≥1 2 The maximum absolute error occurs when Xd=1 2, resulting in MABRE[RNE]=r−f 2×βEmax (1.9) We nowconsider the bias. As indicated above, theabsolute error for Xd=aand for Xd=1−a(for a<1 2) have the same magnitude, but different sign. Consequently, with respect to the bias, these errors cancel each other out. The only remaining case is for a=1 2, which produces a positive error. To have a bias equal to 0, is a case to be treated in a special manner. The IEEE standard specifies this case, as rounding is done to even1. i.e, the unbiased round to nearest is RNE(X)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Xfif Xd<1 2 Xf+r−fif Xd>1 2 Xfif Xd=1 2and Xf=even Xf+r−fif Xd=1 2and Xf=odd Consequently, for this mode RB[RNE]=0 (1.10) The bias (RB). This is defined as the average absolute error 1Rounding to odd in the tie case also has a bias of zero. However, round to even is preferable because it leads to less error when the result is divided by 2–a common computation. 12 Chapter 1. Introduction considering an unsigned significand (if the signed significand is used, the bias wouldbe zero for most roundingmodes) and measures the tendency toward error of a particular sign. To compute this average, it is necessary to consider a frequency distribution of the values of unsigned significand. The usual assumption is a uniform frequency distribution (this may not occur in typical applications), in which case RB = lim t→∞ X∈{Xm+t}(Rmode(X)−X) #X(1.11) where {Xm+t}is the set of all unsigned significands with m+t bits, and #Xis the number of significands in the set. The RNE (unbiased) produces the smallest possible absolute error and has zero bias. Rounding towards 0(RZ) (truncation), in which the result should be the value closestto and no greater in magnitude thanthe exact value. In terms of the operation on the infinite precision significand, the rounding significand is obtained by discarding Xd. I.e. RZ(X)=(X×rf)r−f=Xf The absolute error is ABRE[RZ]=−Xdr−f×βE Since Xd<1the maximum absolute error is MABRE[RZ]≈r−f×βEmax This absolute error is larger than the RNE. Furthermore, for an 1.3. Special values 13 unsigned significand, the absolute error is always negative and the bias is significant. Its value is RB[RZ]≈−1 2r−f Rounding towards+∞RP or upwards, in which the result should be the value closest to and no less than the exact value. Roundingtowards-∞RMor downwards,in whichthe resultshould be the value closest to and no greater than the exact value. These two directed modes are useful in interval arithmetic, in which the operands and the result of an operation are intervals. This permits the monitoring the accuracy of the result. In terms of the infinite precision siginificand and the sign, RP(X)=⎧ ⎨ ⎩ Xf+r−fif Xd>0and s=0 Xfif Xd=0or s=1 RM(X)=⎧ ⎨ ⎩ Xf+r−fif Xd>0and s=1 Xfif Xd=0or s=0 By default RNE is the active rounding mode for IEEE standard. More information on the rounding can be obtained from [?]. 1.3 Special values The exponents in the floating point formats have a biased representation. The bias representation is preferred because it simplifies the comparison of floating point numbers by making it a fixed point comparison and the minimum exponent is represented by 0, so that the representation of the floating point value 0is all zeros (0sign, 0exponent, 0significand). 14 Chapter 1. Introduction In a biased representation with bias B, the signed integer Exis represented by the positive integer denoted by EBsuch that EB=Ex+B(1.12) To represent the minimum exponent by EB=0, we obtain B=−Emin (1.13) For a symmetric range −B≤Ex≤B(1.14) resulting in −0≤EB≤2B(1.15) If eis the number of bits of the binary representation of EB, then 2B≤2e−1(1.16) Consequently, for Binteger we obtain B≤1 2(2e−2) (1.17) for example e=8(single precision computation) we can make B= 127 and EB=Ex+ 127 (1.18) for a symmetric exponent range of −127 ≤Ex≤127. Note that the maximum value of EBis 255, so that this value can be used to represent Ex= 128 1.3. Special values 15 (non-symmetric range) or as a singularity condition. Two values are not used for real numbers, these special values correspond to the smallest possible exponent, Emin which is used for +0,−0and denormalized numbers, and the largest possible exponent, Emax which is reserved for −∞,+∞ and NaNs. The interpretations of a floating point value are the following: •A maximal exponent and a non-zero significand represents NaN. The value in the siginificand may be set according to the machine representation allowing different NaNs. •A maximal exponent and a zero significant represents −∞ or +∞according to the sign bit. •A minimal exponent and a non-zero signinficand represent a denormalized number. •A minimal exponent and a zero significand represents +0 or −0according to the sign bit. •The remaining exponent values and significands represents normal numbers. The uses of the special values are: the overflow operation or division of a real number by zero is −∞ or +∞according to the sign of the intermediate result or the dividend. The IEEE standard 754 calls for two NaNs:signaling and quiet. The signaling NaN is intended to provide for uninitialized variables and arithmetic-like enhancements not included in the standard. The quiet NaN is intended to provide retrospective diagnostic information for invalid or unavailable operands and results. The result of an under flowing operation is +0 or −0, according to the sign of the intermediate result. Denormalized numbers are intended to support gradual underflow [?,?]. 16 Chapter 1. Introduction 1.4 Exceptions There are five types of exception that must be signaled when detected (i.e. a status flag must be set or a trap routine is called): •Invalid operation, when an operand is not valid for the operation to be performed. The result is a quiet NaN. •Division by zero, if the dividend is a finite non-zero number and the divisor is zero. The result is a +∞or −∞. •Overflow, if the magnitude of a result exceeds the largest finite number representable in the floating point format of the operation. The delivered result may be +∞or −∞ or plus or minus the largest representable number in the floating point format, depending on the rounding mode. •Inexactresult, if the roundedresult of anoperation isnot exact or if overflow occurs but there is no overflow trap. 1.5 Importance of Division and Square root Functions Use of the division and square root functions are segmented. They are frequent in 3D graphicsanimation, virtualreality (especially lightingeffects), viewport transformation, motion synthesis, signal processing, scientific computing applications etc. These types of applications one may call high computing applications. These functions are less intense in internet applications, business applications (such as audio and video streaming), charting, presentation, snap shot viewing, day to day accounting etc. One may distinguish these types of applications as low computing applications. In early days, computing division and square root functions were carried out in software emulations. However, most recent processors [?,?,?,?] complying 1.5. Importance of Division and Square root Functions 17 with IEEE-754 standard are equipped with FPU’s performing FP addition, subtraction, multiplication, divisionand square root functionsin VLSI. We encourage readers to obtain FP addition, subtraction, and multiplicationalgorithms from [?]. In the next chapter we discuss division and square root algorithms. A survey performed on FPU’s [?] reveals that while the majority of the microprocessors surveyed carry out both FP addition and multiplication in 3 cycles, FP division consumes about 60 cycles and usually FP square root slower than FP division. The survey shows that the designers have put their emphasis on the development of faster FP adders and multipliers. This negligence intentionally widens the performance gap by downplaying FP division and square root functions. Software developers take advantage of FP addition and multiplication algorithms to avoid these complex operations. Another investigation performed by Oberman [?] reveals the relationship between the latency of FP division/square root and the system performance. The study shows that FP adders and multipliers are the consumers for 27% and 18% of FP divider/square root results respectively. This means that if an inefficient FP divider and square root is used in the FPU, the processor interlock period generally increases, because the FP divider/square root result consumers have to wait longer for the data. Now days, there is a trend towards using addition operation in a multiplication. This is also called multiply add fused operation for division and square root functions, such as is found in IBM 603eTM microprocessor FPU [?]. The multiply-add fused operation can increase the performance of division and square root functions. Dealing with the division and square root more seriously and the better balancing of performance among other FP units, is more reasonable than compromising the overall performance of the whole processor. 18 Chapter 1. Introduction Chapter 2 Classification of Hardware Methods In this chapter we discuss hardware methods for the division and square root functions. They can be classified into two groups: non-iterative methods and iterative methods, we briefly discuss these methods. The non-iterative methods are limited to single precision computation, due to prohibitive area requirements, leadingto itsuse asan initialapproximatorfor theiterativemethods. The methods we discuss in this chapter are incorporated in most state-of-the-art processors. The main types of techniques used for approximating division and square root in hardware can be classified in two groups: Non-iterative Methods [?,?,?,?,?,?,?,?] Direct table look-up. Polynomial or rational approximations. Table-based methods. Iterative Methods [?,?,?,?,?,?] Functional iterative methods. Digit-recurrence and on-line algorithms. In the first group are direct table look-up, polynomial and rational approxima- 19 26 Chapter 2. Classification of Hardware Methods finding a sequence K1,K 2,K 3,...such that ri=YK 1K2...K i−→ 1(2.9) and therefore qi=XK1K2...K i−→ X Y(2.10) The reciprocal (Q=1/Y) can be computed as specific case of the division: qi=K1K2...K i−→ 1 Y(2.11) The first factor K1is usually obtained as a low precision approximation of the reciprocal. If ri=YiKihas the form ri=1−α, then it is possible to obtain Ki+1 =1+α, by a simple two’s complementing ri, leading to ri+1 =1−α2, which guarantees the convergence of the algorithm. In summary, the steps to be performed are the following: •The first factor K1is a low-accuracy approximation of the reciprocal 1/Y , and can be obtained by a table-based method, such as direct table look-up, bipartite table algorithms or polynomial approximation. Let this value have m-bit accuracy, i.e.: 1−2−m<K 1Y<1+2 −m(2.12) •Define v=1−K1Y.|v|<2−m. -Computation of r1=YK 1=1−v. -Computation of q1=XK1(this multiplication is not necessary for the reciprocal computation). •By 2’s complementing r1,K2=1+vis obtained. 2.2. Iterative methods 27 -Computation of r2=r1K2=1−v2. -Computation of q2=q1K2. •By 2’s complementing r2,K3=1+v2is obtained. -Computation of q3=q2K3. At this point, q3<X/Y <q 3(1 + 2−8m). The number of iterations required to obtain a result, accurate to a certain precision, is a function of the accuracy of an initial approximation (seed) value. By using a more accurate seed value, the total number of iterations can be reduced. Square root/inverse square root computation-Let us assume an n-bit input operand Ysatisfying 1≤Y<2. In this case a sequence K1,K 2,K 3,...should be found such that: K1,K 2,...K i−→ 1 √Y(2.13) and therefore qi=YK 1,K 2,...K i−→ √Y(2.14) The stepstobe performedin thesquare rootor inversesquare root computation are: Let Y=1.d1,d 2...d n−1and define  Y=1.d1,d 2...d mwhere m<<n. Two look-up tables can be employed to obtain both K1r=1/ Yand K1=1/ Ywith m-bit accuracy. These tables must be designed so that each table look-up value K1rcorre- sponds at full target accuracy to the square pf the table look-up value K1. Another way of computing K1rand K1is obtaining K1with a table driven method (which leads to an important area reductionregardingdirect tablelook-up) andthen computingthe squaring operation K1r=(K1)2. The low precision approximations K1rand K1are then employed 28 Chapter 2. Classification of Hardware Methods in some calculations: -Computation of r1=YK 1r. -Computation of q1=YK 1(only for the square root computation; for inverse square root computation, q1=K1). v 1=1−r1and compute the independent multiplications: -Computation of q2=(1+v1 2)q1. -Computation of (1 + v1 2)2=(1+v1 2)(1 + v1 2). -Computation of r2=(1+v1 2)2r1. Define v2=1−r2and compute the independent multiplications: -Computation of q3=(1+v2 2)q2. At this point, q3=√Y(1+α), when computingthe squareroot function, andq3=(1+α)/√Y, when computing inverse square root, with |α|<28m−2. The total number of iterations to be performed again depends on both the target precision and on the accuracy of the initial approximation. One of the main drawbacks of using the functional iterative method is the difficulty in obtaining a correctly rounded result. 2.2.2 Newton-Raphson Algorithm In the N-R algorithm [?], a primitive function is chosen that has a root at the reciprocal or inverse square root, depending on the function to be computed. The selection of this primitive function is based on the convenience of the resulting iterative form. Division/Reciprocal computation, the most widely used target root, is the divisor reciprocal 1/Y, which is the root of the primitive function. f(q)=1 q−Y=0 (2.15) 2.2. Iterative methods 29 The well known quadratically converging N-R equation is qi+1 =qi−f(qi) f(qi)(2.16) Applying this equation to the primitive function 2.15, gives the resulting iterative formula used to compute the reciprocal: qi+1 =qi×(2 −Y×qi)(2.17) with quadratically converging error: i+1 =Y 2 i(2.18) Each iteration involves two multiplications and a subtraction, which is commonly replaced by a two’s complement operation. The computation of the division requires a final multiplication by the dividend (X/Y =X×(1/Y )). The multiplications to be performed are dependent operations and therefore cannot be computed in parallel, which leads to performance degradation. Square root/inverse square root computation- Applying the Newton-Raphson equation to the primitive function: f(q)= 1 q2−Y=0 (2.19) the following iteration formula is obtained 2qi+1 =qi×(3 −Y×q2 i)(2.20) with quadratically converging error: 2i+1 =(3 √Y 2 i+Y 3 i)(2.21) 30 Chapter 2. Classification of Hardware Methods Each iteration involves three multiplicationsand the computation of the square root requires a multiplication by the operand Y, since the N-R algorithm converges to an inverse square root. 2.3 Digit-recurrence Algorithm The algorithmbelongs to thesame type of initialapproximationof the functions in the hardware, usually known as digit-by-digit iterative methods due to their linear convergence. This means that a fixed number of bits of the result are obtained in each of the iterations. Implementation of this type of algorithm, typically of low complexity, utilizes a small area and relatively large latencies. The fundamental choices in the design of a digit-by-digit algorithm are the radix, the allowed coefficients or digits and the representation of the partial remainder (residual). 2.3.1 Digit-recurrence division Digit-recurrence algorithms use subtractive methods to calculate quotients, one digit per iteration. SRT divisionis the name of the most commondigit-recurrence divisionalgorithm[?,?]. The quotientisdefined tocomprise Nradix-rdigits,with r=2 b(2.22) and N=n/b(2.23) Such a division algorithm requires N iterations to compute the final n-bit result, since b bits of quotient are retired per iteration. The following recurrence is used in every one of the iterations of the SRT algorithm: W[i+1]=rW[i]−Yq i+1,(2.24) 2.3. Digit-recurrence Algorithm 31 where qi+1 is the (i+1) th digit of quotient, q[i]is the quotient in step i, and W[j] is the partial remainder in step i, with an initial value of W[0] = X. In each of the iterations, one digitof the quotientis determinedby the selection function qi+1 =SEL(W[i],Y)(2.25) The final quotient after N iterations is then Q= N  i=1 qir−i(2.26) Choice of radix-The fundamental method for decreasing the overall latency of the algorithm is to increase the radix r, typically chosen as a power of 2 in order to perform the products by ras shiftsof the operands. To achieve the same precision, the number of iterations required to computethe final result will be reduced. However, the reduction has some side effects, such as increasing the complexity of the selection function, increasing the area and, in most cases, also increasing the cycle time. Therefore, an analysis of trade-offs between area and speed is necessary for determining the values of the radix r, which results in efficient implementation. Choice of digitset-In thedigit-by-digitalgorithm, somerange of digitsmustbe chosen forthe allowed values ofthe quotientin each of theiterations. The simplest case is where, for radix r, there are exactly rallowed values (non-redundant digitset). However, to increase the performance of the algorithm, a redundant digit-set is usually employed. Such a digit-set can be composed of symmetric signed-digit consecutive integers, where the maximum digit is a: qi∈{−a, −a+1,...,−1,0,1...,a−1,a}(2.27) Thus, to make a digit-set redundant, it must contain more than rconsecutive 32 Chapter 2. Classification of Hardware Methods integer values including zero and amust satisfy a≥r/2(2.28) The redundancy of a digit set is determined by the value of the redundancy factor ρwhich is defined as ρ=a r−1,ρ> 1 2(2.29) When a=r/2the representation is called minimally redundant, while that with a=r−1is called maximally redundant, with ρ=1. A representation is known as non-redundant if a=(r−1)/2, while a representation a>r−1is called over-redundant. For the nextresidual W[i+1] to be bounded when a redundant digit-set is used, the value of the quotient digit must be selected so that: |W[i+1]|<ρ×Y(2.30) By using a large number of allowed digits a(a large value of ρ), the complexity and latency of the selection function can be reduced. However, choosing the smaller number of allowed digits simplifies the generation of multiples of the divisor Y. Multiples that are a power of two can be generated by simply shifting. If a multiple is required that is not a power of two (e.g. three), an extra operation, such as addition, may be required, adding to the complexity and latency of generating the divisor multiple. The complexityof digit selection and that of generating multiples of the divisor Y must be balanced. After the redundancy factor ρis chosen, it is possible to derive the selection function. A containment condition determines the selection intervals. A selection interval is the region in which a particular quotient digit can be chosen. These 2.3. Digit-recurrence Algorithm 33 expressions are given by Uk=(ρ+k)X Lk=(−ρ+k)X, (2.31) where Uk(Lk)is the largest (smallest) value of residual W[i], such that it is possible for qi+1 =kto be chosen while still keeping the next partial remainder bounded. However, when a high radix is used, the only practical method is selection by rounding [?,?,?,?]. Two alternatives have been used to allow selection by rounding: performing a scaling of recurrence [?,?] and performing selection by table in the first iteration until the convergence conditions are met [?]. Choice of residual representation-The residual can be represented in two different forms, either redundant or non-redundant. A conventional two’s complement is an example of the non-redundant form, while CS two’s complement and SD binary representation are examples of the redundant form. Additions and subtractions are involved in digit-recurrence algorithms. If these operations are performed in non-redundant form, they require a full width adder with carry propagation, increasing the cycle time. If the residual is computed in redundant form, a carry free adder, such as a CSA, can be used in recurrence, minimizing the cycle time. However, the complexity of the selection function increases and so, additionally, twice as many registers are required to store the residual between iterations. It is possible to convert the quotient digits as they are produced, in order to avoid the extra cycle required to perform addition. This scheme is known as onthe-fly conversion [?]. It can also be extended to perform on-the-fly rounding of the result [?]. Combined Division and Square root Computation: for computing the square root, the employed recurrence is 34 Chapter 2. Classification of Hardware Methods W[i+1]=rW[i]−f[i]qi+1 (2.32) By defining f[i]=2q[i]+qi+1r−(i+1) for the square root and f[i]=Yfor the division, a single unit for both recurrences can be obtained. A common set of selection constraints must be valid for both operations, and f[i]must be generated without adding extra delay to the critical path of the division recurrence. 2.3.2 Accurate Quotient Approximation method Digit recurrence algorithms are applicable to low radix division and square root implementation, as radix increases the quotient digit selection hardware becomes complex, increasing the cycle time, area or both. AQA is a variant of a digit recurrence algorithm, for achieving a very high radix division/square root with an acceptable increase in cycle time, area and precise rounding with simpler quotient digit selection hardware. The term ”very high radix” is applied when large numbers of bits are retired in the iteration. The high radix algorithm was independently proposed by: Wong and Flynn [?,?], and Matula [?,?]-it is carried out in Cyrix processor. Both these algorithms are identical. The algorithm is as follows: 1. Initially, set the quotient Q=0,P=Xfor the division (Q i=X/Y ) and P=Yfor the square root (Q i=Y/√Y). Then, obtain an approximation of 1/Yh,1/√Yhfrom a look-up table, using the top mbits of Y, returning an hbit approximation. Only m−1bits are actually required to index into the table, as the guaranteed leading bit can be assumed. 2. Scale the truncated divisor by the reciprocal approximation and inverse square root approximation and obtain the partial quotient for the square 2.3. Digit-recurrence Algorithm 35 root. Similarly, scaling the dividend by the reciprocal approximation partial quotient can be used for the division. Y=(1/Yh)×Yq 1=P×1/Yh qi=P×1/Yh...(i≥1)qi=P/2×1/Yh...(i≥2) (2.33) ...division ...square root The Yis invariant across the division iterations, therefore, only needs to be performed once. Subsequent iterations use only one multiplication: Y×P, The product P×1/Yhand P×1/√Yhcan be viewed as the next partial quotient bits. 3. Perform the general recurrence to obtain the next partial remainder: P i=P−P×Y...for division (2.34) P i=P−qi×(2Q+qi)...for square root (2.35) 4. Compute the new quotient as Q i=Q+qi...for division (2.36) Q i=Q+qi...for square root (2.37) The new quotientis then developed by forming the product P×(1/Yh),P× (1/√Yh)and adding the result to the old quotient Q. 5. The new partial remainder P iis normalized by left shifting to remove any leading 0s. 6. Variables are adjusted such that Q=Q i,and P=P i. 42 Chapter 3. Modified N-R Algorithm without a LUT If each tripel t 2i+1,t  2i,t  2i−1substituted by a tripel (rt3i,rt 2i,t y1i), the value of the multiplier changes to t= m−1  i=0 (−2.rt3i+rt2i+rt1i).4i(3.15) Therefore, one can define the set of tripels(rt3i,rt 2i,rt 1i), with 0≤i<m  to be a redundant Booth-2 digit representation of t,ifft=t. Where m is a PP’s reduction factor for multiples of the multiplicand with bbit multiplier (b=29bit in our case discussed in the next section) having t m+1 =t m=t −1=0 then m=(b+1)/2. Thus, it does notchange the product if fed to the multiplier or the set of tripels t 2i+1,t  2i,t  2i−1based on the non-redundant representation of tor a redundant Booth-digit representation (rt3i,rt 2i,rt 1i)of t. Initial approximation is a two’s complement of the divisor Y. To obtain t[1], first making a one’s complement of the divisor and thereafter performing a one’s complemented multiple generation of t[0]. Then, performing a PP’s summation of t[0] using Yas the multiplier, during the PP’s summation, the constant 1can be added. During 21 iterations (i.e. to generate t[i+1]), the divisor Yis used as a multiplicand and t[i]as a multiplier. Truncating the intermediate result, t[i], at the 29th position, deletes the fraction portion t30,t 31 ..., keeping a value in the fraction range (-1/2,1/2)ulp. We assume to have a carry-save representation of t+const, that already includes the additive constant const =m−1 i=0 2.4i. The const is required to obtain tin eqn.3.15 form. Looking at bit windows with a width of 2 in this carry-save representation of the number t+const, each window contains 4 bits (see Figure 3.1); two with a weight of one and two with a weight of two. The binary value wjof the part of the number within a window iis in the range wi∈{0,...6}. The number t+const can then be written by: 3.2. CS Booth digit representation 43 S’ S S’ S S’ S S’ S S’ S S’ S S’ S 2i+7 2i+6 2i+5 2i+4 2i+3 2i+2 2i+1 2i 2i-1 2i-2 2i-3 2i-4 2i-5 2i-6 C2i-6 C2i C2i-2 C2i-4 C2i+2 C2i+4 C2i+6 2 bit adder 2 bit adder 2 bit adder 2 bit adder 2 bit adder adder 2 bit 2 bit windows carrysave number representation of <t>+ const 2-bit adder row redundant Booth digit representation of <t> cs (t , t ) Figure 3.1: Compression from carry-save to redundant Booth-digit representation (a) 2 bit addition in CSA (b) 3 bit addition in CSA Figure 3.2: CS-Booth digit representation in a CSA cell 44 Chapter 3. Modified N-R Algorithm without a LUT t+const = m−1  i=0 wi.4i(3.16) If we input the 4 bits of a window iinto a 2-bit adder, we get three output bits c2i+2,s 2i+1 and s2i, that represent the value of the window by wi=4.c2i+2 +2.s2i+1 +s2i(3.17) t+const = m−1  i=0 (4.c2i+2 +2.s2i+1 +s2i)4i,(3.18) Subtract the additive constant const on both sides, and obtain the value of t t+const − m−1  i=0 2.4i= m−1  i=0 (4.c2i+2 +2(s2i+1 −1) + s2i)4i,(3.19) As x−1≡− xfor x∈{0,1}, then one can substitute s2i+1 −1by −s2i+1 and we can get t= m−1  i=0 (4.c2i+2 −2s2i+1 +s2i)4i(3.20) In radix-4 recoding [?,?,?,?,?] the multiplier is recoded with s2m+1 = s2m=c0=0, then we have 3.2. CS Booth digit representation 45 t+const = m  i=0 (2.s2i+1 +s2i+c2i)4i(3.21) Subtracting the additive constant const on both sides, and obtain the value of t t+const − m−1  i=0 2.4i= m  i=0 (2.(s2i+1 −1) + s2i+c2i)4i(3.22) t= m  i=0 (−2s2i+1 +s2i+c2i)4i(3.23) Therefore, the set of triples (s2i+1,s 2i,c 2i)is a redundant Booth-2 digit representation of tand the sequence of operations described before is a partial compression from a carry-save representation of t+const to a redundant Booth 2 digit representation of t. The implementation of this partial compression is depicted in Figure 3.1. In the radix-4 multiplier [?], the const can be introduced by replacing the fourth level of 3:2CSA with 5:2CSA as shown in Figure 3.5 (on the right hand side). In the proposed method, performing CS-Booth-2 digit representation is just wire crossing in adjacent CSA cells. Figure 3.2 (a) depicts the CS–Booth-2 digit representation in a CSA for radix-4 multiplier. Where tsi,t ci are the intermediate sums and carries of t. We assumed a 3:2CSA consisting of a combination of X-or’s, OR’s and AND gates, however it can be built with NAND or NOR gates. In the radix-8 multiplier, the constant is const =m−1 i=0 4.(8)i, the addition of const to t(i.e. t+const) can be performed by replacing the fifth level 46 Chapter 3. Modified N-R Algorithm without a LUT CS-Redundant Booth digit representation S2i+2 S2i+1 S2i C2i Select: 4M 2M 3M M S Figure 3.3: CS-Booth recoding of 3:2CSA with 4:2CSA (the 3:2CSA requires two sums and a carry from the previous level, by replacing 3:2CSA with 4:2CSA we introduce another input which is const). A row of 3 bit adders are required for partial compression from CS representation of t+const to redundant Booth-3 digit representation of t, which is a wire crossing in an adjacent CSA cell, Figure 3.2(b) depicts the redundant Booth-3 digit representation in the final sixth level of 3:2CSA.An extra radix-8 recoder is required at the end of the CSA tree, in order to recode in the same pipe line stage of PP’s summation. In the radix-8 multiplier, the recoder is coupled with a 3Mmultiple generation, which is another pipe line stage [?]. Figure 3.3 depicts CS-radix-8 Booth recoder. The results of this recoder can be directly applied to the multiplexer to select the multiplicands. In the radix-4 multiplier, an extra radix-4 recoder [?] is not required, because the radix-4 recoder of the multiplier is placed in the same pipe line stage of the PP’s summation. More details on the proposed method of architecture can be obtained in Section 3.4. Note that, the multiplexersin Figure 3.2 (a)& (b) may not add additional delay in the cycle time, because they can be selected before the start of the PP’s summation. In the case of the radix-4 multiplier, it can be selected at the start of the first stage and in the case of the radix-8 multiplier; it can be selected at the start of 3.3. Error Computation 47 the second stage. This partial compression can easily be extended to the carry-save representation to higher radix Booth recoding multipliers. For compressions to redundant radix-2kBooth digits, const has to change to const =m−1 i=0 2k−1.(2k)iand a row of k-bit adders has to be used for partial compression. 3.3 Error Computation For single precision computation, an estimate accurate to n+3bits before rounding must be obtained. In principle, a result accurate to n+1bits before rounding would suffice to guarantee exact rounding. However, when a multiplier and adder are used in the iterations, performing either truncation or rounding at (n+1) th position of the intermediate result, returns less than nbits of precision. In this case, an accuracy of n+3bits of the estimate is required instead of n+1. In order to guarantee an exactly rounded result, we performed exhaustive simulation of the reciprocation algorithm for all 224 arguments. We obtained an accuracy of ±0.5ulp (for 22 iterations). Under the described conditions, the result before rounding tsatisfies: −2−n−3<t<2−n−3(3.24) −2−n−3<Q<2−n−3(3.25) Rounding is performed by adding 2−n−3to tand then truncating the resulting value to n+2. Therefore, we must guarantee ti+1n+2 ≤2−26, since the range of results is 0.5<t≤1. Error computation for reciprocation can be set as: t[0] =t[0] + 0(3.26) 48 Chapter 3. Modified N-R Algorithm without a LUT s= (t[i]+i).Y +sr (3.27) t[i+1]n+2 =t[i]+i+t[i].Y +i.Y +sr +tr (3.28) where sr and tr are the error introduced by S and the t[i+1]operation due to finite word length at n=29, sr, tr <2−29.0is the error introduced by the initial approximation at n=4. 0<2−4, therefore for iiterations i<2−26. Therefore, the error bound for S (for iiterations), with respect to truncation, is s≤2−26. The error in the final result must be bounded by: t[i+1]n+2 =i+i.Y +sr +tr ≤2−26 +2 −26 +2 −29 +2 −29 (3.29) The result before rounding Qmust satisfy eqn.3.25. Rounding can be performed by computing a corresponding remainder: rem =X−Y×Q. By observing the sign and magnitude of rem and the value of n+1, all IEEE rounding modes can be implemented by selecting either Q, Q+2−nor Q−2−n. The action table for correctly rounding Q is shown in Table 3.2 [?]. Guard Bit Remainder RN RP(+/-) RM(+/-) RZ 0=0 trunc trunc trunc trunc 0 - trunc trunc/dec dec/trunc dec 0 + trunc inc/trunc trunc/inc trunc 1=0 RNE inc/trunc trunc/inc trunc 1 - trunc inc/trunc trunc/inc trunc 1 + inc inc/trunc trunc/inc trunc Table 3.2: Action table for rounding For directed rounding modes RP and RM, the action depends on the sign of estimate: those entries that contain two operations such as pos/neg correspond to the final result being positive or negative respectively. Most industrial multipliers 3.4. Architecture and Performance 49 incorporate these modes. 3.4 Architecture and Performance The architecture of the proposed method is, simply to introduce a radix-8 Booth recoder at the end of CSAtree in the radix-8 multiplierand an extra 3:2CSA isrequired for CS–Booth-2 representation in the radix-4 multiplier. We used a radix-4 multiplier of IBM Power PC 603eTM FPU [?], which is targeted at the mobile computing market and a radix-8 multiplier of IBM G5 FPU [?] which is targeted at high performance applications. The radix-8 multiplier (64 ×56 bit) is pipelined in three stages: in the first stage, generating 3Xmultiplicand and Booth recoding, in the second stage, 19 PP’s summation in 6 levels of 3:2CSA, and the third and final stage is for carry propagate addition (see Figure 3.4). The proposed method requires 24(Y)×29(t[i]) bit multiplicationand 10 PP’s are accumulated outof 19PP’s. It requires 4+1levelsof CSAtree forthe multiply add operation: The first 3 levels of 3:2CSA’s for 9 PP’s summation. In the 4th level, 4:2CSA summation for 10th PP, a carry, a sum and a carry from previous 3:2CSA’s. That plus 1 level of 3:2CSA for adding t[i]+S(not shown in Figure 3.4). As mentioned, in the sixth level, CS–Booth-3 digit representation is performed with an extra radix-8 recoder at this sixth level. The results for this recoder can be directly applied to the 4:1 multiplexer to select the multiplicands. The two’s complement for the initial approximation can be done by first making a one’s complement of the multiplicand and then adding 1 during the PP’s summation. Note that, in the case of the radix-4 MAF (i.e. 24(B) + 24(A)×24(C)multiplied result added to another 24 bit operand say, B(B+A×C))[?], commonly 50 Chapter 3. Modified N-R Algorithm without a LUT multiplier reg.multiplicand reg. 3X reg. X reg. Booth select reg. Booth decode 4:1 multiplexer carry and sum reg. Product reg. 3X 120−bit adder Radix−8 recoding 19−to−2, 6 level 3:2 CSA tree 1M 4M partial products 1M 4M Figure 3.4: Proposed architecture in the radix-8 multiplier implemented in industry due to its simplicity and ease in creating multiples of the multiplicand, a similar architecture and performance can be obtained to that of the radix-8 multiplier. This is an architecture with a three stage pipeline. In the first stage, Booth recoding and 14 PP’s summation in three levels of 4:2CSA tree. The CS result is added to the aligned B operand in the final fourth level 3:2CSA as shown in the first stage of Figure 3.5. The CS result is passed to the second stage, to obtain the result in non redundant form. The third stage is for normalization and rounding. The Cis 28 bits, it can be made by concatenating four 0’s for a single precision multiplier. It can also support double precision MAF . More details of this architecture are discussed in the next chapter. 3.4. Architecture and Performance 51 FRC C exp 0 Booth recoding aligned Bexp−bexp−a sumcarry LZD Normalization,Rounding,Tag and Exception tmp exp exp adj Multiply Add 161−bit Adder Norm/Round 53x 28 exp adj 161−bit mantissashift cnt Alignment shift count C mantA mantB mantA expB exp FRAFRB St2 St3 St1 CC 3:2 CSA(106b) 3 levels of 4:2 CSA Multiplier Array Feed back CS t[i] sum carry sum carry t[i] S 3:2 CSA To Booth re c t[i] CS Alignment Shifter 4 to 1 Mux Const 5:2 CSA Figure 3.5: Proposed architecture in the MAF With the 24(Y)×29(t[i]) bit multiplication in the proposed method, the first stage is for Booth recoding and 15 PP’s summation in three levels of 4:2CSA tree and in the final fourth level t[i]+Sis performed (see Figure 3.5 on the right hand dashed lines) by replacing 3:2CSAwith 5:2CSA. To obtain the t[i]in CS–Booth- 2 representation, an extra 5th level of 3:2CSA is required at the end of 5:2CSA of the fourth level in the first stage, as shown in Figure 3.5 with dashed lines. An extra 3 bit Booth-2 recoder to 14 bit Booth-2 recoder and an extra row of 5:1 multiplexer is required during multiplexing of the multiplicands to generate 15 PP’s. The advantages of the proposed method with respect to N-R method are that 58 Chapter 3. Modified N-R Algorithm without a LUT the multiplier and, also, a large LUT is not available i.e., the performance is about the same as the average performance of state-of-the-art processors. 3.7 Conclusion Single precision computation of a modified Newton-Raphson reciprocation algorithm is proposed. It is a variable latency algorithm and only requires a maximum of 22 iterations providing a linear convergence in character. A variable latency algorithm is useful in speed independent design, also called a self-timed divider. The self-timed divider provides processing capability according to the task at hand. In this way, it can reduce power consumption without sacrificing performance. The advantages of our method are that it does not require a look-up table and each of the iterations is a cycle, unlike the pairs of cycles found in the state-of- the-art algorithms. The performance of the proposed method is about the same as the average performance of state-of-the-art processors. To maintain the average performance, intermediate results are computed in CS form. The multiplyadd operation and Booth recoding can be performed in the same cycle. Initial approximation in our case is a two’s complement of the divisor, which can be performed during partial products summation. The multiples of the divisor are constant throughouteach of the iterations. Therefore we generate them only once, in the second cycle with a radix-4 multiplier and in the third cycle with a radix-8 multiplier, i.e. in the second iteration. This leads to a reduction of a cycle in each iteration, unlike the N-R method. The eqn.3.10-eqn.3.13 could be considered as a multiplication in hardware terms. In the N-R method the number operation required for reciprocation are: a multiplication and multiply-add operation (see eqn.3.2-eqn.3.4) and each of the iterations requires two cycles. The multiply-add operation requires the multiplicand (S) to be in non-redundant form. It is possible to use the redundant form, but the data width required is twice that of the non-redundant form, unlike the pro- 3.7. Conclusion 59 posed method. Without the use of a LUT for the N-R algorithm, the performance is about the same. This is because each iteration requires two multiplications, which is about 6 cycles (we assumed 3 cycles for a multiplication), it requires 4 iterations to obtain the single precision result. With respect to the SRT method, the proposed method does not require a quotient digit selectionhardware to select the next quotient digits, which increases the cycle time with respect to the radix. We used both a radix-4 and a radix-8 multiplier. In a radix-8 multiplier, there may be a delay of an OR gate in cycle time, due to the extra radix-8 Booth recoder. In a radix-4 multiplier, it is 1τslower in cycle time, due to another level of 3:2CSA for CS–Booth-2 digit representation. We set-out a reciprocation/division algorithm with better silicon efficiency than conventional methods as shown in Table 3.4. The proposed method can fit into a radix-4 multiplier with extra hardware of about 31fa and about 6fa for a radix-8 multiplier, which is less than a 1.5%increase in area with respect to industrial multipliers. There is about 16%improvementin performance with respect to IBM and Itanium processors, but 4%decrease with respect to IBM 603eand AMD processors, due to the fact that 1.09 bits retired per iteration, the multipliers are shared and a large LUT is not available i.e., performance is about the same as the average performance of state-of-the-art processors. The proposed method can be further extended to double precision computation, which requires about 50 cycles. This can be reduced by employing an appropriate look-up table with market targeting. It offers a good trade-off between performance and area making it suitable for the mobile computing market. In the mobile computing market area is crucial. Most architects look for a divisionalgorithm with a low cost and reduced area, and an acceptableperformance thatcan runinternetapplications, businessapplications such as audio and video streaming, charting, presentations, snapshot viewing, day to day accounting etc. It is expected that, with the proposed method, processor 60 Chapter 3. Modified N-R Algorithm without a LUT architects will have an option for their next generation processors. Chapter 4 Modified AQA Method Using a LUT In this chapter we propose a division and square root using the Accurate Quotient Approximation Method for double precision computation. The proposed method employs a second degree minimax approximation to obtain an initial 15 bits estimate of reciprocal and inverse square root values. This is then passed through an iteration of Goldschimdt reciprocal and inverse square root algorithms, to enhance the initial approximation and to make initial approximation more silicon efficient. Thereafter we perform an iterationof the Accurate Quotient Approximation method to obtain a double precision result. To evaluate the proposed method, we used a radix-4 MAF of IBM 603eTM and a radix-8 multiplier of G5 FPU. We then compared it with conventional processors and it shows that there is a significant improvement in performance with better silicon efficiency. To obtain the best performance, intermediate results are computed in CS form and some computations are overlapped. We used the AQA method because there are some advantages over the N-R method, GLD algorithm, Cyrix algorithm and SRT algorithm, it is detailed in this chapter. 61 62 Chapter 4. Modified AQA Method Using a LUT Minimax approx 15 15 30 Goldschimidt itn Modified AQA method Rounding 30 m=4&3 55 15 30 3 0 1/ Y , 1/ Y 30 1/ Y , 1/ Y X/Y , Y/ Y Figure 4.1: Proposed Concept 4.1 Proposed Method Some steps of the AQA method can be modified for a single iteration algorithm, to obtain the double precision results, using high seed, 30 bits, for the reciprocal and inverse square root. Figure 4.1 exhibits the concept of the proposed method. We interface an initial 15 bit minimax approximation with GLD reciprocation and inverse square root algorithms to make the initial approximation silicon efficient. Interfacingthishighseed30bitsinitialapproximationwithmodifiedAQA method boosts the performance, allows the overlapping of partial remainder computation and, furthermore, offers more parallelism. A brief aspect of these features is depicted in Section 4.6. 4.1. Proposed Method 63 Recalling eqn.2.34 - eqn.2.37 and algorithm steps 1-7 of subsection 2.3.2 for i=1and 2Q=0 P 1=P−P×YP 1=P−q1×(2Q+q1) P 1=X−X×Y×(1/Y30)P 1=Y−Y2×(1/Y30)2(4.1) Q 1=Q+q1Q 1=Q+q1(4.2) ...division (Q i=X/Y )...square root (Q i=Y/√Y) Where Pis current partial remainder, Qis current partial quotient, P 1is next partial remainder, Q 1is next partial quotient and q1is next partial quotient bits. One would substitute Qas X/Y30 for the division and Y/√Y30 for the square root, assuming that these operations guarantee the first 28 bits of the quotient. Multiplying the next partial remainder P 1with respective initial approximations i.e. P 1×(1/Y30)for the division and P 1/2×(1/√Y30)for the square root can be viewed as a current partial quotient bits q1in our case, adding this partial quotient to the respective Qone can obtain a final Q 1. Substituting these assumptions results in Q 1=(X)×(1/Y30)+(P 1)×(1/Y30)...division Q 1=(Y)×(1/Y30)+(P 1/2) ×(1/Y30)...square root (4.3) We require these assumptions because AQA is a linear convergent algorithm. I.e. h−2bits are generated in each iteration (28 bits in our case). By performing X/Y30 and Y/√Y30 we can guarantee the first 28 bits of the quotient. To obtain the double precision result, we need to pass an iteration to obtain the lower 28 bits of the quotient, which is P 1×(1/Y30)and P 1/2×(1/√Y30)for the respective functions. Therefore, the upper 28 bits of Q 1and lower 28 bits of Q 1can be performed in parallel. Simplifying the eqn.4.1 and the eqn.4.3 both the next partial remainder and 64 Chapter 4. Modified AQA Method Using a LUT the quotient can be viewed as P =−Y×X×(1/Y30)P =−Y2×(1/Y30)2(4.4) Q=(2X+P)×(1/Y30)Q=(S+P/2) ×(1/Y30)(4.5) ...division (Q=X/Y30)where S=Y+Y/2...square root (Q=Y/√Y30) To allow a fair comparison, we assumed a 30 bit initial approximation for the Cyrix processor algorithm. The concept comparisons are explained using Figure 4.2, Table 4.1 and Table 4.2. In the Cyrix square root algorithm, after obtaining the squared initial approximation(1/√Y30)2and Y2(initiallyQ=0), they are multipliedfirstly for iteration i=1. Follow the algorithm steps 1 to 7 of subsection 2.3.2 and Figure 4.2. Then, this result is subtracted from Y(P=Y) for i=1and the result is stored in a register P i, this is a next partial remainder. During the multiplication of Y2with (1/√Y30)2, an independent operation q1=P×1/√Y30 is performed for i=1, thereafter it is added to Q(Q=0for i=1) and the result is stored in a register Q i, this is a next partial quotient. The next partial remainder and quotient are updated as a current partial remainder and quotient i.e. P=P 1,Q=Q 1in the first iteration to allow iteration for the second iteration i=2. In the second iteration, the current partialremainder Pis right shifted and then multiplied with 1/√Y30 to obtain q2. Afterwards, the q2is added to a left shifted result of Qi.e. 2Q+q2. The result of this addition is then multiplied with q2(see eqn.2.35). These operations are shown in a dotted square box of Figure 4.2. The multiplied result is then subtracted from P. The result of this operation is the next partial remainder P 2. The next partial quotient Q 2is obtained by adding q2to Q. An example of this algorithm using radix-10 is depicted in Table 4.1. In the proposed method for the square root, after obtaining the squared initial approximation(1/√Y30)2and Y2, they are multiplied. The result of thisoperation we call current partial remainder P. Follow the eqn.4.4&eqn.4.5 and Figure 4.1. Proposed Method 65 Proposed Method Cyrix Processor Square Root Algorithm Y=1.546949120382167. 1/√Y30 =0.804010989. (1/√Y30)2= 0.646433670432758121. Y2=2.393051581051160208. S=Y+Y/2= 2.3204236805732505. P =−Y2×(1/√Y30)2, Q=(S+P/2) ×(1/√Y30). Y=1.546949120382167. 1/√Y30 =0.804010989. Q=0,P =Y,q1= P×1/√Y30,(1/√Y30)2= 0.646433670432758121, P2=2.393051581051160208 for i=1.qi=P/2×1/√Y30 for i≥2.P i=P−qi×(2Q+qi), Q i=Q+qi.P=P i,Q=Q i i=1 P = −1.546949117073816457, Q=1.2437640935411212 q1=1.243764092211146147. P 1=3.30835050 ×10−9,Q 1= 1.243764092211146147 P=P 1,Q=Q 1 i=2 q2=1.32997500 ×10−9. P 2=3.74683 ×10−17, Q 2=1.2437640935411212 Table 4.1: Comparison of Proposed method with Cyrix algorithm 4.2. During the multiplication of Y2with (1/√Y30)2, an independent operation Y+Y/2(S) is performed. We then add the result of this addition to a right shifted P i.e. S+P/2. The result of this addition is then multiplied by 1/√Y30, which in our case is a quotient. An example of our algorithm is depicted in Table 4.1 (for radix-10). In the Cyrix division algorithm, after obtaining Y, which is a product of 1/Y30 ×Y(follow the algorithm steps 1 to 7 of subsection 2.3.2 and Figure 4.2), it is multiplied by P(P=X) for the first iteration i=1and is then subtracted from P. This result is a next partial remainder and is stored in a register P i. During the multiplication of Pby Y, an independent multiplication Pwith 1/Y30 is performed to obtain q1. This result is added to Q(Q=0for i=1). The result is a next partial quotient and is stored in a register Q i. Both the next partial remainder and quotient are updated as a current partial remainder and quotient i.e. P=P 1, Q=Q 1in the first iteration to allow to iterate for second iteration i=2. In the second iteration, current partial remainder Pis multiplied by Yand the 66 Chapter 4. Modified AQA Method Using a LUT Y’ Q Q ... Q Q.. ... .. P P P .. .. .. ... (i=1)(i=1) Y 30 Y 2 (i=1) 1/Y , Y (1 P=X, P=Y (i=1) 30 P P/2 30 Q=0 (i=1) (i=1) (i>=2) P=X, 2Q (i>=2) h 1h ih h 1h ih (i>=1) (i>=2) (i>1) (i>=2) (i>=2)(i>1) (i=1) Sqrt P, P P=X,P=Y P Q (i>1)(i>=2) (i>=2) (i=1) P, P (i>=2) (i>=2) i=i+1)( (i>= 2) (i>1)(i=i+1) (i>=1) (i>1) 12 3 q ,q , q ... .. q 123 p , p , p ....p h h 2 / Y ) 30 1/Y , 1/ Y 30 P/ Y , P/ Y 30 Q’ i P’ i (a) AQA method 12 3 P’’ 30 30 1/ Y , ( 1/ Y ) 2 X, Y Y Y Y/2 2X 30 1/ Y , 1/ Y 30 p , p , p ...... .... ... ..... ....p 59 Q 28 Q 28 123 q , q , q .. .. .. ... q 28 Q’ (b) Modified AQA method Figure 4.2: Concept comparison of Cyrix algorithm with the proposed method 4.1. Proposed Method 67 result of this product is subtracted from P, the result is a next partial remainder P 2. Note that the Yis generated only once, in the first iteration. During the multiplication Yby P, another independent multiplication Pby 1/Y30 is performed to obtain q2and the result of this multiplication is added to Qto obtain the next partial quotient Q 2. An example of this algorithm using radix-10 is shown in Table 4.2. In the proposed method (for division), after obtaining X×Yproduct and 1/Y30, they are multiplied. The result of this multiplication we call current partial remainder P, follow the eqn.4.4&4.5 and Figure 4.2. The result of this multiplication is subtracted from a left shifted Xi.e. 2X+P and then this result is multiplied by 1/Y30 to obtain, in our case, the quotient Q. An example of the proposed method is depicted in Table 4.2 (for radix-10). In summary, the number of dependent operations in Cyrix algorithms is (for both division and square root): an addition (in square root), multiplication, subtraction and an addition–follow right to left of the eqn.2.34-eqn.2.37. It requires two iterations to obtain double precision results. The first 28 bits and the last 28 bits of the quotient are generated serially. In the proposed method, the number of dependent operations for the division and square root are the same: a multiplication, an addition and multiplication– follow the eqn.4.4&4.5. Note that X×Y,Y2and Y+Y/2are independent operations. It requires an iteration to obtain double precision results. The first 28 bits and the last 28 bits of Qare generated in parallel and it is added during multiplication, to obtain final Q(eqn.4.5). The proposed method can be operated in fixed point mode, unlike Cyrix algorithms. The 30 bit initial approximation for the reciprocal and the inverse square root can be obtained by indexing the top mbits of Yto the look-uptable, which returns more than mbits of precision for these functions. Somewhat more than mbits of precision can be obtained from [?,?,?,?]. To obtain a 30 bit result, one would consider 9≤m≤20. 74 Chapter 4. Modified AQA Method Using a LUT at position 2−r−1before such truncation, the rounding error would be bounded by 2−r−1instead. In the minimax method round to the nearest is preformed, by adding one to the C0in advance at position 2−r−1i.e. C 0=C0+2−16 this is stored in the LUT instead of C0. We then truncate the intermediate result at position 2−16 to guarantee 15 bit results. Summarizing the total error in the final result can be expressed as: total ≤ approx +|C2|Y2 2<2−16 (4.14) 4.2.2 Goldschimdt Algorithm To obtain the double precision resultfor the divisionand the square root i.e. 53 bit, we need 30 bit initial approximation, passing this 15 bit approximation through an iteration of GLD reciprocation and inverse square root algorithms [?], we can obtain 30 bit initial approximation. Following are the GLD equations 1/Y30 =1/Y15 ×(2 −(1/Y15)×Y30)(4.15) 2(1/Y30)=1/Y15 ×(3 −(1/Y15)2×Y30)(4.16) 4.2.3 Error computation The error computation for the reciprocal and inverse square root can be set as |1/Y30 |=(1/Y15 +|1/Y15 |)×(2 −[(1/Y15 +|1/Y15 |)×(Y30 +Y30 )] + vt) +t(4.17) 2|1/√Y30 |=(1/√Y15 +|1/√Y15 |)×(3 −[(1/√Y15 +|1/√Y15 |)2×(Y30 +Y30 )] + vt)+t(4.18) 4.2. Initial Approximation 75 The maximum absolute error when computing the reciprocal and the inverse square root, |1/Y30 |,|1/√Y30 |are proportional to 2 1/Y15 and 2 1/√Y15 respectively. The |1/Y30 |and |1/√Y30 |are also depend on the the use of finite arithmetic, this is an intermediate word length truncating at t, introducing error vt and tat the multiplierand thefinal result. Furthermore, it alsodepends onthe absoluteerror of the initial approximations |1/Y15 |,|1/√Y15 |and the error of the truncated divisors, Y30 . Therefore, the error in the final result must be bounded by: |1/Y30 |=21/Y15 .Y30 Y15 +Y30 Y2 15 +vt Y15 +2 1/Y15 .Y30 +21/Y15 .Y30 Y15 +vt. 1/Y15 +2 1/Y15 .Y30 +t<2−30 (4.19) |1/√Y30 |=|1/√Y15 |+3|1/√Y15 |.Y30 2.Y15 +Y30 2.Y15.√Y15 +vt 2.√Y15 +52 1/√Y15 .Y30 2.√Y15 + |3 1/√Y15 |.Y30 +3|1/√Y15 |.Y30 2+vt.|1/√Y15 | 2+2 1/√Y15 .Y30 2√Y15 +t<2−30 (4.20) The target 30 bit precision can be met Y30 , vt <2−32 and t<2−30 by keeping |1/Y15 |,|1/√Y15 |<2−16: |1/Y30 |<2−15 +2 −32 +2 −32 +2 −32 +2 −47 +2 −48 +2 −64 +2 −30 <2−29 (4.21) |1/√Y30 |<2−16 +2 −17 +2 −33 +2 −33 +2 −33 +2 −48 +2 −49 +2 −49 +2 −65 +2−30 <2−29 (4.22) 76 Chapter 4. Modified AQA Method Using a LUT The conservative bounds on Y30 , vt <2−32 and t<2−31 and |1/Y15 |, |1/√Y15 |<2−16 can be met by employing a truncation on the multipliers. 4.2.4 Architecture of Minimax Approximation When m=3for the inverse square root, the lower part Y2has 12 bits and Y2 2 has 6 leading zeros with a word length of 12 bits, if truncation at position 2−18 is performed for Y2 2. Similarly when m=4for the reciprocal, the lower part Y2has 11 bits and Y2 2has 8 leading zeros with a word length of 10 bits, if truncation at position 2−18 is performed for Y2 2. This is because Y2has mleading zeros, which means that Y2 2will have 2mleading zeros. We used a specialized squaring unit for squaring operation Y2 2. More information on the squaring operation can be obtained from [?,?,?]. We use Y2and Y2 2as multipliers. To use Y2 2(Y2is in non redundant form) as a multiplier, it is converted from CS to Booth digit representation (discussed in Section 4.4). The Figure 4.6 shows an initial 15 bit approximation in the standard radix-4 and radix-8 multipliers. The size of the LUT’s for the inverse square root are 2×23×(15 + 10 + 8) = 0.52Kbits=0.064KB. In the architecture we included the LUT for the reciprocal, which is 24×(15+8+7)=0.47Kbits=0.059KB. Therefore, the total size of the LUT’s to be employed is 0.99Kbits which is 0.123KB. In Figure 4.6 (a), is a two stage pipeline architecture. In the first stage, the tables are addressed using the operand Y1to look-up co-efficient values C0,C 1,C 2, in parallel computation of Y2 2. After that, multiplexing for multiplicand and multiplier–to select either for division/square root or general multiplication and then, recoding of Y2and Y2 2. Thereafter, PP’s generation and 10 PP’s summation (6 PP’s of C1+4PP’s of C2) plus C0in three levels of 4:2CSA tree. In the case of the reciprocal, we generate the same 10 PP’s. Figure 4.5 shows the PP’s accumulation of the initial 15 bit approximation. 4.2. Initial Approximation 77 C0 1 C Y2 1 C Y 2 2 4 pp’s of C 2 6 pp’s of C 1 2 2 2 2 2 2 0 −5 −10 −15 −20 −25 Figure 4.5: Minimax accumulation of PP’s We have chosen a two stage pipeline architecture because it can be easily incorporated in the standard radix-4 multiplier (53 ×53 bit), the multiplier accumulates 27 PP’s in four levels of 4:2CSA tree. The last stage is a CPA of the accumulation tree that performs assimilation from the CS-non redundant form of the result. A similar architecture for a radix-8 multiplier (53×53 bit ) is shown in Figure 4.6 (b). In the case of the radix-8 multiplier, it is a three stage pipeline. In the first stage it looks-up values of C0,C 1,C 2from look-up tables, in parallel with squaring of Y2. Then uses multiplexing for the multiplicand and multiplier–to select, either for division/square root or general multiplication. Thereafter 3M multiple generation, which is in parallel with recoding of the multiplier Y2and Y2 2. In the second stage, multiplexing of multiplicands and 7 PP’s summation for inverse square root and reciprocal (4 PP’s of C1and 3PP’s of C2) plus C0in four levels of 3:2CSA tree. The standard radix-8 multiplier requires six levels of 3:2CSA tree for 18 PP’s summation. Finally the third stage is a CPA, to obtain 78 Chapter 4. Modified AQA Method Using a LUT C C C 01 2 PP_gen CSA tree CPA 15 1 0 Mux Unit Squaring Y 2 Y 2 2 53 Y 1 C 0 12 (C +C ) 8 53 Y 2 Y X Y Mux 8 CS feed−bac k Booth recode 00 10 11 12 15 25 25 2525 15&18 4&3 15 1/ Y , 1/ Y 15 (a) Initial approx in a radix-4 multiplier C C C 01 2 1 0 0 1 Unit Squaring Y 2 Y 2 2 53 Y 1 8 12 (C +C ) 8 53 Y 2 Y X Y 15 CPA C 0 1 0 4 to 1 Mux’s 3M Booth recode CSA tree CS−Booth recoding Mux Mux Mux 12 15 25 25 25 25 4&3 15&18 15 1/ Y , 1/ Y 15 (b) Initial approx in a radix-8 multiplier Figure 4.6: Minimax approximation in the Industrial Multipliers 4.2. Initial Approximation 79 the result in the non-redundant form. The access time for LUT’s are 1/4th of cycle time and 3Mmultiplegeneration for (C1+C2)is half the cycle time of standard multipliers. 4.2.5 Delay Estimates of Initial Minimax Approximation In this section we analyze the cycle time for conventional multipliers and the initial approximation. We will show that the cycle time of initial approximation can be matched with the cycle time of the conventional multipliers. The cycle time we analyzed in terms τ-the delay of one full adder [?,?]. These estimates are reliable approximations. The actual speed-up depends on the technology employed and its implementation. We do not consider inter-connection delay because it is technology dependent. In principle, the cycle time depends upon the scheduling and market targeting. The pipelining of our architecture is flexible and can be performed differently, according to the technology to be employed. We used a radix-4 multiplier of IBM Power PC 603eTMFPU [?] and a radix-8 multiplier of IBM G5 FPU [?]. The estimated cycle time of the proposed initial approximation is compared with these industrial multipliers. The Power PC FPU uses a radix-4 multiply-add fused (53 + 53 ×53 bit) operation i.e. 53(A)×53(C)multiplied result added to another 53 bit operand say, B (B+A×C). It isa four stage pipelinefor doubleprecision MAF (Multiply- Add Fused) and is a three stage pipeline for single precision MAF. In the case of double precision MAF,53(A)×28(C)-the lower 28 bits of C multiplied and added to aligned B operand in the final 3:2 CSA as shown in the first stage of Figure 4.13 of Section 4.6, the result (in CS form) is fed back to the second pass. In the second pass, 53(A)×25(C)-the first 25 bits multiplied1, the PP’s of this pass assimilates with the first pass CS result so, 16 PP’s were added (14 PP’s + carry and sum of the first pass) with three levels of 4:2CSA’s. In both 1In second pass the first 25 bits of C made to 28 bits with a concatenating three 0s to generate 14 PP’s. 80 Chapter 4. Modified AQA Method Using a LUT Radix-4MAF Proposed method in radix-4 MAF first pass & second pass (trecoding +t5:1mux)+3× t4:2CSA +1×t3:2CSA +treg (trecodingY2andY 2 2+t5:1mux)+ 3×t4:2CSA +treg 2τ+3×1.5τ+1τ+1τ=8.5τ2τ+3×1.5τ+1τ=7.5τ Third stage tcpa texp adj +tLZD +treg tcpa +treg 7τ+1τ+1τ=9τ7τ+1τ=8τ Fourth stage tnorm texp adj +tround +treg – 5τ+2τ+1τ=8τ– Table 4.3: Delay estimates in radix-4 multiplier passes 14 PP’s are generated. In the second stage, CPA which is in parallel with exponent aligning (exp adj)) and leading zero detection (LZD). And in the third stage, normalization and rounding. The total latency of MAF operation is 4 cycles for double precision: first two cycles for dual pass multiplication, in the third cycle, carry propagate addition, and the final fourth cycle for normalization and rounding. In the case of single precision it is 3 cycles: in the first cycle single pass multiplication, in the second cycle carry propagate addition and the final third cycle is for normalization and rounding. The cycle time is set by the delay of the slowestpath in the architecture, which may be in the third stage. Table 4.3 shows the delay estimates for MAF and the initial approximation in the MAF. Initial approximation in the MAF requires single pass multiplication. The cycle time of MAF can be estimated as 9τ. Note that the LUT’s and squaring unit can be placed either in the first stage or before the first pipeline stage. However, we assumed that it was placed before the first pipeline stage and is coupled with 6:1 multiplexers (see Figure 4.13 on the top right). The radix-8 multiplier (64 ×56 bit) [?] consists of three stages: in the first stage, 3Mmultiple generation and recoding of the multiplier. In the second stage, multiplexing of the multiplicands and 19 PP’s summation in six levels of 3:2CSA’s and, finally, the third stage is for CPA. There is a selection of two pos- 4.2. Initial Approximation 81 Radix-8 multiplier Proposed method in radix-8 multiplier First stage t3Mtrecoding +treg tlut tY2squaring +t2−to−1mux + t3C1,3C2,3Y2tY2andY 2 2recoding +treg 7τ+1τ=8τ1.5τ+0.5τ+3τ+1τ=6τ Second stage t4−to−1mux +6×t3:2CSA+treg t4−to−1mux +4×t3:2CSA + 1×t3:2CSA for cs−booth digit rep + tradix−8rec +treg 1.5τ+6×1τ+1τ=8.5τ1.5τ+4×1τ+1τ+0.3τ+1τ=7.8τ Third stage tcpa +treg tcpa +treg 7τ+1τ=8τ7τ+1τ=8τ Table 4.4: Delay estimates in radix-8 multiplier sible normalization results in the third cycle. This selection signal is built into a custom designed data flow for speed. Table 4.4 shows the delay estimates of the radix-8 multiplier and the initial approximation in the radix-8 multiplier. The cycle time of the radix-8 multiplier can be estimated as 8.5τ. Therefore the cycle time of the initial approximation may not affect the cycle time of the industrial multipliers. The 15 bit approximation can also be obtained from linear approximation algorithms, bipartite method [?,?,?]. In these cases the table size is more than 0.45KB, the quadratic Minimaxapproximationrequires 0.123KB. Otherquadratic methods such as [?,?] can be used, but they may add delay to the critical path for the industrial multipliers. The Minimax method may not add critical path delay to the industrial multipliers. However, the use of novel efficient algorithms in terms of speed and area is welcome. 4.2.6 Architecture and Performance of initial approximation To evaluate the proposed 30 bit initial approximation, we used a radix-4 MAF of IBM 603eTM FPU [?] and a radix-8 multiplier of G5 FPU [?,?]. As discussed 82 Chapter 4. Modified AQA Method Using a LUT before, the latency of radix-4 MAF operation is 3 cycles (single pass multiplication) with a throughput of 1 and the latency of the radix-8 multiplier is 3 cycles with a throughput of 1. In the radix-4 multiplier, after obtaining the 15 bit initial inverse square root approximation in the non redundant form as shown in Figure 4.6 (a), it has to be squared (1/√Y15)2to perform an iteration of GLD inverse square root algorithm– follow the eqn 4.16. For squaring 1/√Y15, using the 1/√Y15 as a multiplier and a multiplicand, we generate 8 PP’s of 1/√Y15 and they are accumulated in two levels of 4:2CSA tree. The product is obtained in CS form to convert to Booth- 2 digit representation (discussed in Section 4.4) for use as a multiplier and Y30 as a multiplicand. Using (1/√Y15)2as a multiplier, we generate 16 PP’s of Y30 and they are accumulated in three levels of 4:2CSA tree, the constant 3can be introducedin the CSAtree (replacing last 4:2CSAof first levelwith 5:2CSA2[?]). The result is obtainedin CS–Booth-2 digit form to use as a multiplier and 1/√Y15 as a multiplicand. We generate 16 PP’s of 1/√Y15 and they are accumulated in a three levels of 4:2CSA tree, thereafter CPA addition of the product to obtain 1/√Y30 in the non-redundant form. Performance: The performance of the initial 1/√Y30 is followed by above explanation, it is shown in the Figure 4.7. First couple of cycles for obtaining the 1/√Y15 in the non redundant form. In the third cycle, obtaining (1/√Y15)2CS–Booth-2 digit representation. The fourth cycle is for performing 3−(1/√Y15)2.Y30, obtain the result in CS–Booth-2 digit representation. And finally the fifth and sixth cycles are to obtain 1/√Y30 in the non redundant form. 2Replacing 4:2CSA with 5:2CSA in the first level of last 4:2CSA the delay in cycle time may be compensated, as this 5:2CSA is a late arriving signal, first preference could be given to this 5:2CSA for summing in a routing algorithm of PP’s summation. 4.2. Initial Approximation 83 115 / Y CS 1/ Y CS non rd St1 St2 non rd St2 St1 15 2 CS St1 3− 15 2Y30 CS St1 1/ Y 30 non rd St1 St2 1 2 3 4 5 6 CS St1 1/Y30 CS St1 non rd St2 Reciprocal Inverse square root (1/ Y ) (1/ Y ) 15 15 2−1/Y Y 30 Figure 4.7: Timing diagram of initial 30 bit approximationin the radix-4 multiplier In the radix-8 multiplier, after obtainingthe 1/√Y15 in thenon-redundant form as shown in Figure 4.6 (b), it is squared. Using the 1/√Y15 as a multiplier and a multiplicand, we generate 6 PP’s of 1/√Y15 and these are accumulated in two levelsof a 3:2CSAtree. The result is obtainedin CS–Booth-3 digit representation (discussed in section 4.4) and recoded using an extra radix-8 recoder. The result of this recoder can be directly applied to the multiplexer to multiplex the multiplicands for the next cycles. During the PP’s summation of 1/√Y15, an independent instruction 3Y30 multiple can be performed, because the multiplier has a throughput of 1. Using (1/√Y15)2as a multiplier and Y30 as a multiplicand, we generate 11 PP’s of Y30 and these are accumulated in five levels of a 3:2CSA tree and the constant 3can be introduced in this tree. The result is obtained in CS–Booth-3 digit representation with recoding. During this accumulation, the 3(1/√Y15)multiple instruction can be performed. Using 3−(1/√Y15)2.Y30 as a multiplier and 1/√Y15 as a multiplicand we generate 11PP’s of 1/√Y15 and these are accumulated in five levels of a 3:2CSA tree. The result is allowed to propagate the carry to obtain 1/√Y30 in the non-redundant form. 90 Chapter 4. Modified AQA Method Using a LUT (a) 2 bit addition in CSA (b) 3 bit addition in CSA Figure 4.11: CS-Booth digit representation in a CSA cell B+const = m−1  i=0 (4.c2i+2 +2.s2i+1 +s2i)4i,(4.28) Subtract the additive constant const on both sides, and obtain the value of B B+const − m−1  i=0 2.4i= m−1  i=0 (4.c2i+2 +2(s2i+1 −1) + s2i)4i,(4.29) As x−1≡− xfor x∈{0,1}, then one can substitute s2i+1 −1by −s2i+1 and we can get B= m−1  i=0 (4.c2i+2 −2s2i+1 +s2i)4i(4.30) 4.4. CS to Redundant Booth digit representation 91 In radix-4 recoding [?,?,?,?,?] the multiplier is recoded with s2m+1 = s2m=c0=0, then we have B+const = m  i=0 (2.s2i+1 +s2i+c2i)4i(4.31) Subtracting the additive constant const on both sides, and obtain the value of B B+const − m−1  i=0 2.4i= m  i=0 (2.(s2i+1 −1) + s2i+c2i)4i(4.32) B= m  i=0 (−2s2i+1 +s2i+c2i)4i(4.33) The set of triples (s2i+1,s 2i,c 2i)is a redundant Booth 2 digit representation of B. The implementation of this partial compression is depicted in Figure 4.10. In the radix-4 MAF [?], the const is introduced by replacing the third level of 4:2CSA with 5:2CSA (the 4:2CSA requires two carries and sums of previous level, by replacing this 4:2CSA with 5:2CSA we introduce another input which is const). Performing CS–Booth-2 digit representation is a just wire crossing in a adjacent CSA cells, it can be performed in the final 3:2CSA level. Figure 4.11 (a) depicts the CS–Booth-2 digit representation in a CSA for radix-4 multiplier. Where Bsi,B ci are the intermediate sums and carries of B. We assumed a 3:2CSA consists of a combination of X-or’s, OR’s and AND gates, however it can be built with NAND or NOR gates. In the radix-8 multiplier [?], the constant is const =m−1 i=0 4.(8)i, it can be introduced by replacing the fifth level of 3:2CSA with 4:2CSA (the 3:2CSA re- 92 Chapter 4. Modified AQA Method Using a LUT CS-Redundant Booth digit representation S2i+2 S2i+1 S2i C2i Select: 4M 2M 3M M S Figure 4.12: CS-Booth recoding quires two sums and a carry of previous level, by replacing this 3:2CSA with 4:2CSA we introduce another input which is const). Performing CS–Booth-3 digit representation is a just wire crossing in an adjacent CSA cells, it can be performed in the final sixth level of 3:2CSA. An extra radix-8 recoder is required at the end of the CSA tree. Because of the need to recode in the pipe line stage of PP’s summation. In the radix-8 multiplier, the recoder is coupled with 3Mmul- tiple generation which is another pipe line stage [?]. Figure 4.11 (b) and Figure 4.12 depicts the 3 bit addition and CS–radix-8 Booth recoding. The results of this recoding can be directly applied to the multiplexer to select the multiplicands. In the radix-4 multiplier, an extra radix-4 recoder [?] is not required because the radix-4 recoder of the multiplier is placed in the same pipeline stage of PP’s summation. The multiplexers in Figures 4.11 (a)& (b) may not add additional delay in cycle time because they can be selected before start time of the PP’s summation. In the case of the radix-4 multiplier, it can be selected at the start of the first stage and in the radix-8 multiplier it can be selected at the start of the second stage. The work presented in [?] requires an additional hardware for 2 bit adder3,in 3Theauthorproposeda 30bitunfoldedarchitectureforredundantreciprocationapproximation. 4.5. Error Computation 93 our design, as mentioned, this can be performed by just wire crossing in the CSA. We suggest that you obtain more information on CS-Booth kdigits from [?]. In [?,?] other partial compressions of redundant number representation were described. Mostof these recoders require at leastfour logiclevels. Each of the two logic levels is a delay of mux & xor gate. From this point of view, the described technique is simpler, faster and the hardware can be shared with the CSA tree. 4.5 Error Computation An accuracy of n+2bits of the final estimate is required, in order to guarantee exact rounding. The result before rounding Q 1must satisfy: −2−n−2<[Q 1]<2−n−2(4.34) Rounding is performed by adding 2−n−2to Q 1and then truncating the result to a value to n+1bits to form Q 1, which has an error of strictly ±0.5ulp. Therefore, we must guarantee |q|<2−55, since the range of output result 0.5≤Q 1<2and 1≤Q 1<2for division and square root. Error computation for the division and square root operation using eqn.4.1& eqn.4.3 can be set as : |p 1|=X−(X.Y +X.Y )×(1/Y30 +|1/Y30 |)+pt(4.35) |q 1|=(X+|p 1|)×(1/Y30 +|1/Y30 |)+qt(4.36) |p 1|=Y−(Y2+Y2)×(1/Y30 +|1/√Y30 |)2+pt(4.37) |q 1|=(Y+|p 1|/2) ×(1/Y30 +|1/√Y30 |)+qt(4.38) error accumulatesfrom theinitialapproximation(Minimax methodand GLDalgorithm) and error due to finite word length of the intermediate result at b= 94 Chapter 4. Modified AQA Method Using a LUT 59,pt, X.Y , Y2and qt≤2−59. The target precision can be met keeping |1/Y30 |,|1/√Y30 |<2−29 (truncating at position 2−30), solving the above equations results |q 1|<−X.Y 1/Y30 Y30 −X.Y 2 1/Y30 +pt Y30 −X.Y Y2 30 −2X.Y 1/Y30 Y30 +pt|1/Y30 | −X.Y 2 1/Y30 +qt<2−55 (4.39) |q 1|<2−29 +2 −58 +2 −59 +2 −59 +2 −87 +2 −88 +2 −117 +2 −59 <2−55 (4.40) |q 1|<Y30Y30|1/√Y30 |/2−Y30Y 2 1/√Y30 +pt/2Y30 −Y2|3 1/√Y30 |/2+ pt|1/√Y30 |/2+qt<2−55 (4.41) |q 1|<2−30 +2 −58 +2 −60 +2 −88 +2 −89 +2 −59 <2−55 (4.42) These conservative bounds on pt, X.Y , Y2and qtcan be met by employing truncation in the multiplier, which carry out P 1and Q 1at b=59. The result before rounding Q 1must satisfy eqn.4.34. Rounding can be performed by computing a corresponding remainder. By observing the sign and magnitude of remainder and the value of n+1, all IEEE rounding modes can be implemented by selecting either Q 1,Q  1+2 −nor Q 1−2−n. The action table for correctly rounding Q 1is shown in Table 4.5 [?]. For directed rounding modes RP and RM, the action depends on the sign of estimate: those entries that contain two operations such as pos/neg correspond to the final result being positive or negative respectively. Most industrial multipliers incorporate these modes. 4.6. Architecture and Performance 95 Guard Bit Remainder RN RP(+/-) RM(+/-) RZ 0=0 trunc trunc trunc trunc 0 - trunc trunc/dec dec/trunc dec 0 + trunc inc/trunc trunc/inc trunc 1=0 RNE inc/trunc trunc/inc trunc 1 - trunc inc/trunc trunc/inc trunc 1 + inc inc/trunc trunc/inc trunc Table 4.5: Action table for rounding 4.6 Architecture and Performance The architecture of the proposed method is, simply to introduce LUT’s and a squaring unit during input operand multiplexing in a multiplier. We used the IBM Power PC 603eTM microprocessor FPU [?,?], which is targeted at portable computers and an IBM G5 FPU [?,?] that is targeted at scientific computing applications. The proposed method requires minimum hardware requirements in these multipliers. The Power PC FPU uses radix-4 MAF operation. As discussed previously, the latency of MAF operation is 4 cycles: the first two cycles for dual pass multiplication, in the third cycle, carry propagate addition, and the final fourth cycle for normalization and rounding. In a radix-4 multiplier [?], after obtaining the squared 30 bit initial approximation (1/√Y30)2,Y2and Sin the square root, we need to multiply i.e. an iteration of −Y2×(1/√Y30)2is required to obtain double precision result (eqn.4.23). By using Y2as a multiplicand and (1/√Y30)2as a multiplier, we generate 15 PP’s (P) in each pass and compressed in three levels of 4:2CSA, the S(algorithm step 2 of section 4.3) can be subtracted during the PP’s summation of the second pass (using the space of addition feed back of CS result i.e. replacing the last 4:2CSA of the first level with 5:2CSA see [?]) as shown in Figure 4.13 on the right dashed arrow. And then, with this result (acting as a multiplier) multiplied by 1/√Y30 (59 ×30 bit), we generate 30 PP’s of 1/√Y30-a dual pass multiplica- 96 Chapter 4. Modified AQA Method Using a LUT tion, thereafter the result in CS form is allowed to propagate the carry to obtain Q. Similarly, for the division, after obtaining the 30 bit initial approximation 1/Y30 and X×Yproduct, an iteration of −X×Y×(1/Y30) is required to obtain the double precision result. By using the result of X×Yas a multiplicand and 1/Y30 as a multiplier, we generate 15 PP’s and this is summed in three levels of a 4:2CSA. The operand 2Xcan be introduced during this PP’s summation (using the space of addition feed back of CS result i.e. replacing one of 2:1 multiplexer for CS feed back with 3:1 multiplexer, which is in parallel with the 5:1 multiplexer for multiplicand selection see [?]) as shown in Figure 4.12 on the right dashed arrow. And then, with this result multiplied by 1/Y30 (59 ×30 bit), we generate 30 PP’s of 1/Y30-a dual pass multiplication and the result Qis obtained in CS form. Note that, 15 PP’s can be generated, by adding an extra 3 bit Booth recoder to the 14 bit Booth recoder and adding an extra 5:1 multiplexer to the row of 5:1 multiplexers (for selecting the multiplicands) without additional delay. By replacing 4:2CSA with 5:2CSA in the first level of last 4:2CSA the delay in cycle time may be compensated. As this 5:2CSA is a late arriving signal, first preference could be given to this 5:2CSA for summing in routing algorithm of PP’s summation. As mentioned, the results of (1/√Y15)2and 1/Y15,(1/√Y30)2and 1/Y30, 3−(1/√Y15)2Y30 and 2−1/Y15Y30,S+P/2,2X+P and Q(for division) are obtained in CS form to enhance the performance of the division and square root, to use as a multipliers. The results of these operations are represented in Booth digit representation by just wire crossing in the final level of 3:2CSA (see Figure 4.10& 4.11 (a) followed by eqn’s. 4.25-4.33). To use these operands as a CS multipliers, one has to replace third level of 4:2CSA with 5:2CSA to introduce the const =m−1 i=0 2.4i. The third level of 4:2CSA requires two sums and carries from the previous level, by replacing 4:2CSA with 5:2CSA we introduce another input which is const (not shown in the Figure 4.13). 4.6. Architecture and Performance 97 FRC C exp 0 Booth recoding aligned Bexp−bexp−a sumcarry LZD Normalization,Rounding,Tag and Exception tmp exp exp adj Multiply Add 161−bit Adder Norm/Round 53x 28 Multiplier Array Alignment Shifter exp adj 161−bit mantissashift cnt Alignment shift count C mantA mantB mantA expB exp FRAFRB St2 St3 And Squaring unit 2 Y and Y 2 2 12 (C + C ) C 0 St1 CS feed−back CC S 2X 3:2 CSA(106b) 3 levels of 4:2 CSA Sqrt LUT’s. Rec and Inverse Figure 4.13: Proposed architecture Performance: The square root performance shown in Figure 4.14. In the first two cycles, generating a 15 bit minimax approximation, a single pass multiplication with carry propagate addition is obtained. During two to five cycles Y2, a dual pass multiplication is performed. Four to seven cycles for the generation of 1/√Y30 approximation through GLD algorithm (eqn.4.16). During four and five cycles some instructionsof 1/√Y30 can be performed because MAF[?] has a throughput of 1. 98 Chapter 4. Modified AQA Method Using a LUT 1 2 3 4 5 6 7 8 9 10 11 1 15 / Y non rd Q’ .... .... .... 13 Rounding 14.... ..... ..... ... 17 non rd Y 2 (1/ Y ) 15 3− 15 2 Y 30 2 1( St1 CS St1 St3 1/ Y CS non rd 15 30 1/Y non rd CS Q’ CS Rounding X.Y St1 St2 St1St1 St2 St3 St2 St1 St1 St1 St1 St1 St2 St3 CS non rd St2 St1 CS St1 St1 CS CS St1 St1 St1 CS St1 CS St1 St1 St1 St2 Square root Division St2 St3 1/ Y 30 / Y ) 30 15 2−1/Y Y 30 S S+P’’/2 2X+P’’ 2 (1/ Y ) St3 St1 St2 Figure 4.14: Performance in the MAF multiplier Eight to ten cycles to obtain partial part of Q(eqn.4.24). Note the nine and the ten cycles for dual pass multiplication to obtain S+P/2in CS–Booth digit form. The PP’s summation of these cycles assimilates with added Y+Y/2(S), which was performed in the eighth cycle (as mentioned, the MAF has a throughput of 1). The S+P/2operation acts as a multiplier in the next cycles. Thereafter, elevento thirteencyclesfor Qacting1/√Y30 as a multiplicand and the final, fourteen to seventeen cycles for rounding. This can be performed by computing a remainder: rem = Y−Q2. Thus the square root is performed in 17 cycles. Similarly, for the division, use the timing diagram shown in Figure 4.14. It can be expected to perform in 12 cycles. In the radix-8 multiplier [?,?], similar architecture and performance can be obtained with a three stage pipeline: in the first stage booth recoding and 3M multiple generation, in the second stage PP’s summation and in the third stage carry propagate addition. There is a selection signal built into the customdesigned 4.6. Architecture and Performance 99 1/Y non rd CS St1 Square root Division 30 Y2 St1 115 St3 non rd non rd Y+Y/2 St3 12 Booth rec St1 15 2−Y Y 30 CS St2 non rd CS X.YBooth rec 15 St1 St2 St3 3X.Y St1 2X+P’’ CS St2 St1 Q’ CS Rounding St1 St2 St3 St1 St2 1/Y15 CS St2 non rd St3 Booth rec CS St3 St2 3− 15 non rd Y (1 30 CS CS St2 St2 15 St1 St2 CS 3Y2 Y ( 2St2 CS S+P’’/2 3(1/ St1 CS St1 St2 Q’ non rd Rounding 3 St2St1 St3 St1 3Y30 St2 CS St3 St1 St1 30 3(1/Y ) 3Y30 St1.....St2 / Y (1/ Y )2 (1/ Y )2 3(1/ Y ) 15 15 (1/ Y ) 1/ Y 30 30 3(1/ Y ) Y ) 30 30 1/ Y )2 30 / Y )2 1 2 3 4 5 6 7 8 9 10 11 13..........14 15...... .....17 3(1/Y ) Figure 4.15: Performance in the radix-8 multiplier data flow for the speed. This selection signal selects two possible normalization results in the third cycle. The square root operation requires 59(Y2)×60(1/√Y30)2bit multiplication and 20 PP’s are accumulated out of 19 PP’s (64 ×56 bit), all six levels of 3:2 CSA’s are required to compress the result into sum and carry. An extra 3:2CSA addition i.e. another cycle is required to subtract a compressed result P/2to S, this is because the (1/√Y30)2operation accumulates all levels of 3:2CSA’s. In the final sixth level, CS–Booth-3 representation is performed as shown in Figure 4.11 (b). Thereafter it is recoded using an extra radix-8 recoder (see Figure 4.12). The performance of the square root is similar to that of the radix-4 multiplier as shown in Figure 4.15. Similarly, for the division, with the 59(X.Y )×30(1/Y30)bit multiplication, 10 PP’s were accumulated out of 19 PP’s, the operand 2Xcan be introduced in a 3:2CSA tree. It requires five levels of 3:2 CSA’s out of six and plus 1 level for CS-Booth digit representation along with radix-8 recoder. Multiplying this result with the 1/Y30 results in Q. The division can be expected in 11 cycles, a cycle is reduced because it is a three stage pipeline. Figure 4.15 depicts the division performance in the radix-8 multiplier. 106 Chapter 4. Modified AQA Method Using a LUT with respect to state-of-the-art processor LUT’s. There is a significant improvementin performance withrespect to the state-of- the-art processors. This is because we used high accuracy of initial approximation (30 bit) and intermediate results are computed in carry save form. Firstly we obtain a 15 bit initial approximation from the Minimax method, which is silicon efficient and may not add critical path delay to industrial multipliers. We then interface this approximation with a GLD algorithm to obtain a 30 bit initial approximation. However, the use of state-of-the-art efficient initial approximation algorithms in terms of speed and area is encouraged. An iteration is required after the 30 bit initial approximation to obtain double precision results. It offers the best trade-off between performance and area, making it suitable for both the mobile and scientific computing market. In both the mobile computing and high performance computing market area is crucial. Moreover, the speed is the key in the high computing market. Most architects look for division and square root algorithmsthat can run graphics applications, scientific computingapplications, signal processing applications, audio and video streaming applications, charting etc. It is expected that, with the proposed method, processor architects will have an option for their next generation processors. Conclusions Summary We discussed the division and square root computation for the single precision and double precision computation. Basics of various floating point precisions, rounding modes, special values and exceptions are explained in Chapter 1. The use of division and square root functions depends on the target market. For example, in mobile computing markets such as PDA,UPC, tablet PC , where it can run internet applications, day to day accounting, snapshot viewing, presentations etc. In these markets, the use of division and square root is less intense and so singleprecision can suffice. Also, in these marketsthe processor architects look for low cost and low complexity, but with acceptable performance. The modified N-R reciprocation/division algorithm could be suitable for these markets. The square root function without a LUT, by modifying the N-R algorithm, is proposed. However, this adversely affects performance, because it has to be squared in each iteration and each iteration requires at least three cycles. Somewhat more than 50 cycles are required to compute the square root. More research on the “without a LUT for modified N-R square root algorithm” would be very welcome. These functions are frequent in 3D graphics, scientific computing applications, mathematical computing applications etc. Double precision computation is required in many applications such as the previously mentioned scientific and mathematical computing applications. In these markets, the processor architects 107 108 Conclusions look for higher performance with acceptable levels of cost and complexity. The modified AQA method can be applied for these markets. Furthermore, the modified AQA method can be used in the mobile computing market for future processors. To obtain a higher performance for these functions, the initial approximation in the form of LUT or LUT with linear/quadratic approximation algorithms can be used to reduce the area required for these functions. The initial approximations are usually coupled with the N-R method, the GLD algorithm and the Digit recurrence algorithm, also called the SRT method. A discussion of the initial approximation algorithms and the N-R method, the GLD algorithm and the Digit recurrence algorithm are detailed in Chapter 2. These functions are slower because they require sequences of addition or subtraction, multiplication and multiply-add operations, blocking the FPU throughout the computation process. Some processors implemented these functions in the software, such as in Itanium processors and others implemented them in the hardware, such as in IBM and AMD processors. Implementing these functions in the software provides a high degree of parallelism and does not block the FPU during the computation, but performance is lower. By implementing these functions in the hardware, performance can be enhanced, but this is complex and the cost is higher. Most designers are willing to sacrifice the speed in favor of low cost and low complexity. Performance and complexity depends on the precision and the precision is dependent on the target market or on market demands. We proposed the modified AQA method, which offers significant improvements in performance, with better silicon efficiency than the conventional processors for double precision computation. We also proposed the average performance modified N-R reciprocation algorithm, which offers the best silicon efficiency when compared to the conventional processors used for single precision computation. In Chapter 3, a single precision computation of a modified N-R reciprocation algorithmwithout a LUT isproposed, which is the best siliconefficient algorithm. Conclusions 109 It is a variable latency algorithmand only requires a maximum of 22 iterations with a linear convergence. A variable latency algorithm is useful in a self-timed divider, because the processing capability depends on the task at hand. In this way, it can reduce power consumption without sacrificing performance. To evaluate the proposed method, we used a radix-4 and radix-8 multiplier. In a radix-8 multiplier, there may be a delay of an OR gate in cycle time, due to the extra radix-8 Booth recoder. In a radix-4 multiplier, it is 1τslower in cycle time, due to another level of 3:2CSA for CS–Booth-2 digit representation. We compared the proposed method with the state-of-the-art processors in Table 3.4. The proposed method can fit into a radix-4 multiplier with extra hardware of about 31fa and of about 6fa for a radix-8 multiplier, which is less than a 1.5% increase in area with respect to industrial multipliers. There is about 16%improvement in performance with respect to IBM and Itanium processors, but a decrease of 4%with respect to IBM 603eTM and AMD processors. This is due to the fact that 1.09 bits are retired per iteration; it is shared with the multiplier and a large LUT is not available. Therefore, the performance is about the same as the average performance of the state-of-the-art processors. The significance of the proposed N-R algorithm is discussed in the next section. In Chapter 4, we proposeddivisionand square root functionsby modifyingthe AQA method, using a LUT for high performance computing applications. This requires iteration after the 30 bit initial approximation. The proposed modified AQA method offers some attractive strategies additional to the conventional N-R, GLD,SRT and Cyrix algorithms, which are discussed in the following section. To evaluate the proposed modified AQA method, we used a radix-4 multiplier of IBM 603eTM FPU , which is targeted at the mobile computing market and a radix-8 multiplier of IBM G5 FPU it targeted at high performance applications. We then compared it with the conventional processors for double precision computation in Table 4.7. Comparison shows that better silicon efficiency with significant improvement in performance can be expected. It can be fitted into a radix-4 multiplier with extra hardware of about 60fa and 110 Conclusions with about 72fa in a radix-8 multiplier. This, on average, is somewhat more than a70% decrease in the silicon area, with respect to the state-of-the-art processor LUT’s. Performance is significantly better, because we used high seed, 30 bits, and the intermediate resultsare computed in CSform, i.e. the intermediate results are converted to CS to Booth digit representation for use as multipliers. This is wire crossing in our design. Contributions The advantages/contributions of the proposed modified N-R reciprocation algorithm are that it does not require a look-up-table and that performance is about the same as the average performance of the state-of-the-art processors. The average performance can be maintained by computing intermediate results in CS form. The initial approximation is a two’s complement of the divisor which can be performed during partial products summation. Each of the iterations is a cycle, unlike the pairs of cycles found in the state-of-the-art processors. The multiples of the divisorare constant throughout each of the iterations. Therefore we generate them only once, in the third cycle in the case of the radix-8 multiplier and in the second cycle in the case of the radix-4 multiplier, i.e. in the second iteration. This leads to a reduction of a cycle for each of the iterations. Moreover, the multiply-add operation and Booth digit representation can be performed in the same cycle of partial products summation. It is a variable latency algorithm and requires only 22 iterations. Variable latency algorithms are useful in speed independent designs that provide power according to the task at hand. The eqn.3.10-eqn.3.13 could be considered to be a multiplication in hardware terms. In the N-R method, to obtain the reciprocation, the required operations are: a multiplication and multiply add operation (see eqn.3.2-eqn.3.4) and each of the iterations require two cycles. In the N-R method, the multiply-add operation requires (eqn.3.4) the multiplicand Sto be in the non-redundant form. It is possible to use the multiplicand in the redundant form, but this doubles the data width. Conclusions 111 These features are different in the proposed method (eqn.3.10-eqn.3.13–could be considered to be a multiplication in hardware terms). It is possible to use the N-R method as a without LUT algorithm, as in the proposed method, but performance is about the same. This is because each iteration requires two multiplications i.e. about 6 cycles (we assumed 3 cycles for a multiplication) and requires 4 iterations to obtain the single precision result. With respect to the SRT method, the proposed method does not require a quotient digit selection hardware to select the quotient digits which increases the cycle time with respect to the radix. The proposed modified N-R algorithm can also be extended to the square root computation and it only requires a maximum of 22 iterations. But it does not offer the average performance of the conventional processors. This is due to the squaring of the quotient in each of the iterations. Each of the iterations requires at least three cycles. Therefore, more than 50 cycles are needed to obtain the result. Breakthrough research on this algorithm is substantial. The advantages/contributionsof the proposedAQA methodare shownin Table 4.6. •Unlike the N-R method, only the 30 bit initial approximation is required to obtain the final quotient. •With respect to the GLD algorithm the proposed method is similar if the 30 bit initial approximation was used, but the method of computation in the hardware is different. We perform intermediate results in CS form used them as multipliers and, furthermore, we use the overlapping of some computations. These features increase performance and the comparisons made can be seen in Figures 4.15 and 4.16. •The cycle time might not be increased, although large numbers of bits are retired in the iteration, unlike the digit recurrence method which increases the cycle time with respect to the radix. 112 Conclusions •With respect to the Cyrix algorithm, the modified AQA requires an iteration to obtain double precision results. It requires a multiplication, an addition and a multiplication dependent operation. The first 28 bits and the last 28 bits of the quotient generated in parallel, is added during multiplication. To obtain double precision results, the Cyrix algorithm requires two iterations andanaddition(in squareroot)and multiplication,subtractionand additiondependent operationsare required. The first 28bits and the last28 bitsof the quotientare generated sequentially, if the same 30 bits initialapproximation is assumed. Furthermore, the proposed method can be operated in fixed point mode, unlike the Cyrix algorithm. •The look-up table size for the initial approximation is less than half the table size of the conventional processors. Booth recoding can be performed in the same cycle of partial products summation and overlapping some of the computations of the partial remainder (see Figures 4.14 & 4.15) which increases performance. •To make the proposed method silicon efficient, we first obtain a 15 bit initial approximation using the Minimax method, which is silicon efficient and may not add critical path delay to industrial multipliers. We then interface this approximation with the GLD algorithm to obtain a 30 bit initial approximation. However, the use of novel initial approximation algorithms that are efficient in terms of speed and area, is encouraged. •An iteration is required after a 30 bit initial approximationto achieve double precision results. The modified AQA method can also be extended to the mobile computing market, as it can support high performance applications for future processors. The proposed modified AQA method can fit into industrial multipliers with extra hardware of less than 70fa. The choice is up to the architects. Future Work Performing manual calculations is a part of human nature. Computers readily support binary arithmetic. A recent processor, the IBM z900 series, is the only one capable of performing decimal instructions in the hardware [?,?]. However, its decimal computation capability is limited to integer operands. Recently, decimal arithmetic has become more attractive in the financial and commercial world such as banking, tax calculation, currency conversion, insurance and accounting, etc. The following facts may explain this interest: •Asurveyof commercialdatabases[?] showsthat more than90% ofnumbers are stored in decimal or integer form, while more than half of them are represented in a purely decimal format. •It is well understood that, when converting decimal and binary formats, most fractional decimal numbers are approximately represented in binary floating point representation and, therefore, may loose precision [?,?]. This means that using binary floating point numbers in financial applications, where errors cannot be tolerated, does not necessarily mean that correct results will be obtained. •Regulations, such as the European Directorate General II [?], specify decimal digits for currency calculations. It is likely that in the near future, most high-end processors will perform decimal operations in decimal floating point format using dedicated decimal floating point units. 113 114 Future Work In decimal arithmetic there are some complex functions such as sequential multiplication, division and square root. The fundamental operations of division and square root are based on multiplication. Computing partial products of the multiplication and adding each newly computed product to the previous partial product, cannot be accomplished in a reasonable time [?] without a fast circuit for decimal addition. This can adversely affect speed if division and square root were also taken into consideration. One method for computing partial products is by taking advantage of carry free addition. The building of fast circuits for decimal carry free addition may be a possibility. Designing the computation of division and square root, based on fast addition and multiplication, is challenging. It may be interesting to focus the Accurate Quotient Approximation method on decimal arithmetic. The decimal divider and square root can demand a larger area and, thus, may obtain a reasonable speed-up. Achieving the highest speed in the smallest chip area is one of the most important goals. However, there may be some other areas for designers to consider. Among these is the design of methods using low power consumption [?], motivated by battery operated devices demanding intensive computation in portable environments, and for the proposed use of divider and square root computations, it might be interesting to study whether they can be redesigned for low power use, with or without performance loss in small calculations. Resumen en Espa˜ nol En esta Tesis se analizan los algoritmos para el c´alculo de la divisi´onylara´ız cuadrada para una representaci´on punto flotante con simple y doble precisi´on y se proponen varias implementaciones alternativas. Los aspectos b´asicos relacionados con la representaci´on punto flotante, modos de redondeo, valores especiales y excepciones, se resumen en el cap´ıtulo 1. El uso que se hace de las funciones de divisi´on y ra´ız cuadrada depende, en gran medida, del mercado al que se dirige la implementaci´on. Por ejemplo, en mercados de computaci´on m´ovil, tales como PDAs, tablet PC, d´onde se requiere la ejecuci´on de aplicaciones de internet, visi´on de fotogramas, presentaciones, etc., no se realiza un uso masivo de estas operaciones; por este motivo, una representaci´on de punto flotante con precisi´on simple es suficiente. Por otra parte, en estos segmentos de mercado, se imponen procesadores con una arquitectura de bajo coste y baja complejidad, pero con rendimiento aceptable. El algoritmo Newton–Raphson (N–R) modificado para c´alculo del rec´ıproco y divisi´on, que se propone en esta memoria, puede ser adecuado para este mercado. Otra de las implementaciones propuestas, ra´ız cuadrada sin tabla de aproximaci´on inicial (Look–Up Table, LUT) est´a basado una modificaci´on del algoritmo N–R. Sin embargo, su rendimiento es algo deficiente, porque cada iteraci´on requiere tres ciclos y son necesarias varias iteraciones. El n´umero total de ciclos que se necesitan para el c´alculo de la ra´ız cuadrada con este algoritmo es, aproximadamente, 50. Creemos que es necesario profundizar todav´ıa m´as en este algoritmo para logar una implementaci´on eficiente. 115 122 BIBLIOGRAPHY [9] A. D. Booth. A signed binary multiplicationtechnique. Quart. J. Mech. Appl Math, 4(2):236–240, 1951. [10] W. S. Briggs, T. B.Brightman, and D.W. Matula. Methodand ApparatusFor Performing Division Function Using A Rectangular Aspect Ratio Multiplier. US Patent 5 046 038, US Patent Office, September 1991. [11] W. S. Briggs, T. B.Brightman, and D.W. Matula. Methodand ApparatusFor Performing The Square Root Function Using A Rectangular Aspect Ratio Multiplier. US Patent 5 060 182, US Patent Office, October 1991. [12] F. Y. Busaba, C. A. Krygowski, W.H. Li, E. M.Schwarz, and S. R. Carlough. The ibm z900 decimal arithmetic unit. In Proc 35th Asilomar Conference on Signals, Systems and Computers, pages 1335–1339, November 2001. [13] J. Cao and B. Wei. High-performance hardware for function generation. In Proc 13th Int’l symp.on Computer arithmetic, pages 184–188, 1997. [14] J. Cao, B. Wei, and J. Cheng. High-performance architectures for elementary function generation. In Proc 15th Int’l symp.on Computer arithmetic, pages 136–144, 2001. [15] T. C. Chen. A binary multiplication scheme based on squaring. IEEE Trans. Computers, 20:678–680, 1971. [16] W. D. Clinger. How to read floating point numbers accurately. In Proc of the Conference on Programming Language Design and Implementation, pages 92–101, June 1990. [17] J. N. Coleman, E. I. Chester, C. I. Softely, and J. Kaldec. Arithmetic on the european logarithmic microprocessor. IEEE Trans. Computers, 49(7):702– 715, July 2000. [18] M. Daumas and D. Matula. Recoders for Partial Compression and Rounding. Laboratoire de I’Informatique du Parallelisme, report no.97-01, January 1997. BIBLIOGRAPHY 123 [19] K. Diefendorff, P. K. Dubey, R. Hochprung, and H. Scales. Altivec extension to powerpc accelerates media processing. IEEE Micro, pages 85–95, March/April 2000. [20] E.Antelo, T.Lang, and J.D.Bruguera. Computation of x/d in a veryhigh radix combined division/square-root unit with scaling and selection by rounding. IEEE Trans. Computers, 47(2):15–161, February 1998. [21] E.Antelo, T.Lang, and J.D.Bruguera. Very-high radix cordic vectoring with scaling and selection by rounding. In Proc 14th IEEE symp.on Computer arithmetic, pages 204–213, April 1999. [22] E.Antelo, T.Lang, and J.D.Bruguera. High-radix cordic rotation based on selection by rounding. Journal of VLSI Signal Processing, 25(2):141–153, 2000. [23] M. D. Ercegovac, L. Imbert, D. Matula, J. M. Muller, and G. Wei. Improving goldschimdt division, square root and square root reciprocal. IEEE Trans. Computers, 49(7):759–763, 2000. [24] M. D. Ercegovac and T. Lang. On-the-fly conversion of redundant into conventional representations. IEEE Trans. Computers, 36(7):895–897, July 1987. [25] M. D. Ercegovacand T. Lang. On-the-fly rounding. IEEE Trans.Computers, 41(12):1497–1503, December 1992. [26] M. D. Ercegovac and T. Lang. Division and Square Root: Digit Recurrence Algorithms and Implementations. Kluwer Academic Publishers, 1994. [27] M. D. Ercegovac and T. Lang. Digital Arithmetic. Kluwer Academic Publishers, 2004. [28] M. D. Ercegovac, T. Lang, and P. Montuschi. Very-high radix division with selection by rounding and prescaling. IEEE Trans. Computers, 43(8):909– 918, 1994. 124 BIBLIOGRAPHY [29] M. D. Ercegovac, T. Lang, J.-M. Muller, and A. Tisserand. Reciprocation, square root, inverse square root and some elementary functions using small multipliers. IEEE Trans. Computers, 49(7):628–637, 2000. [30] M. J. Flynn. On division by function iteration. IEEE Trans. Computers, 19:702–706, 1970. [31] D. M. Gay. Correctly Rounded Binary-Decimal and Decimal-Binary Conversions. Numerical Analysis Manuscript 90-10. AT & T Bell Laboratories, 1990. [32] D. Goldberg. What every computer scientist should know about floatingpoint arithmetic. ACM Computing Surveys, 23(1):5–47, March 1991. [33] R. E. Goldschimdt. Applicationsof Divisionby Convergence. Massachussets Institute of Technology Master dissertation, June 1964. [34] D. Harris. A powering unit for an opengl lighting engine. In Proc 35th AsilomarConference on Signals, Systemsand Computers, pages 1641–1645, 2001. [35] V. K. Jain, S. A. Wadecar, and L. Lin. A universal nonlinear component and its application towsi. IEEE Trans. Components, Hybridsand Manufacturing Technology, 16(7):656–664, 1993. [36] T. Jayarshee and D. Basu. On binary multiplication using the quarter square algorithm. In Proc. Spring Joint Computer Conf., pages 957–960, 1974. [37] R. Jessani and C. Olson. The floating-point unit of the powerpc 603e. IBM J.Res.Develop, 40(5):559–566, September 1996. [38] R. Jessani and M. Putrino. Comparison of single-and dual-pass multiply-add fused floating-point units. IEEE Trans. Computers, 47(9):927–937, September 1998. [39] I. Koren. Evaluating elementary functions in a numerical co-processor based on rational approximations. IEEE Trans. Computers, 40:1030–1037, 1994. BIBLIOGRAPHY 125 [40] T. Lang and P. Montuschi. Very-high radix square root with prescaling and rounding and a combined division/square root unit. IEEE Trans. Computers, pages 827–841, August 1999. [41] D. M. Lewis. 114 mflops logarithmic number system arithmetic unit for dsp applications. IEEE Journal of Solid-State Circuits, 30(12):1547–1553, 1995. [42] P. E. Madrid, B. Millar, and E. E. Swartzlander. Modified booth algorithm for high radix multiplication. In IEEE Computer design conference, pages 118–121, 1992. [43] P. Markstein. IA-64 and Elementary Functions. Hewlett-Packard Professional Books, 2000. [44] C. May, E. Silha, R. Simpson, and H. Warren. The PowerPC Architecture: a specification for a New Family of RISC Processors. Morgan Kaufman Publishers Inc., San Francisco CA, 1994. [45] J. M. Muller. Elementary Functions. Algorithms and Implementation. Birkhauser, 1997. [46] J. M. Muller. A few results on table based methods. Reliable Computing, 5(3), 1999. [47] A. Naini, A. Dhablania, W. James, and D. Sarma. 1-ghz hal sparc64 dual floating point unit with ras features. In Proc 15th IEEE symp.on Computer arithmetic, pages 173–183, June 2001. [48] S. Oberman, G. Favor, and F. Weber. Amd-3d now! technology: Architecture and implementations. IEEE Micro, 19(2):37–48, 1999. [49] S. F. Oberman. Floating point division and square root algorithms and implementation in the amd-k7 microprocessor. In Proc 14th IEEE symp.on Computer arithmetic, pages 106–115, April 1999. 126 BIBLIOGRAPHY [50] S. F. Oberman and M. J. Flynn. Design issues in division and other floatingpoint operations. IEEE Trans. Computers, 46(2):154–161, February 1997. [51] S. F. Oberman and M.J.Flynn. Division algorithms and implementations. IEEE Trans. Computers, 46(8):833–854, August 1997. [52] B. Parhami. Computer Arithmetic: Algorithms and Hardware Designs. Oxford University Press, 2000. [53] W. J. Paul and P. M. Seidel. On the complexity of booth recoding. In 3rd Conference on Real Numbers and Computers, RNC3, pages 199–218, April 1998. [54] J. A. Pi˜neiro and J. D. Bruguera. High-speed double precision computation of reciprocal, division, square root and inverse square root. IEEE Trans. Computers, 51(12):1377–1388, December 2002. [55] J. A. Pi˜neiro, S. F. Oberman, J. M. Muller, and J. D. Bruguera. High-speed function approximation using a minimax quadratic interpolator. IEEE Trans. Computers, 54(3):304–318, March 2005. [56] S. Playstation2. 2.44-gflops 300-mhz floating-point vector-processing unit for high performance 3d graphics computing. IEEE Journal of Solid-State Circuits, 35(7):1025–1033, July 2000. [57] S. Playstation2. Vector unit architecture for emotion synthesis. IEEE Micro, 20(2):40–47, 2000. [58] J. Rabaey, A. Chandrakasan, and B. Nikolic. Digital Integrated Circuits. Prentice Hall Publishers, 2002. [59] D. Sarma and D. Matula. Faithful interpolation in reciprocal tables. In 13th Symp.on Computer arithmetic, volume 13, pages 82–91, 1997. [60] D. Sarma and D. W. Matula. Faithful bipartite rom reciprocal tables. IEEE Trans. Computers, 47(11):1216–1222, November 1998. BIBLIOGRAPHY 127 [61] M. J. Schulte and J. E. Stine. Approximating elementary functions with symmetric bipartite tables. IEEE Trans.Computers, 48(8):842–847, 1999. [62] M. J. Schulte and J. E. Stine. The symmetric addition method for accurate function approximation. Journal of VLSI signal processing, 21(2):167–177, 1999. [63] M. J. Schulte and K. E. Wires. High-speed inverse square roots. In Proc 14th Int’l symp.on Computer arithmetic, pages 124–131, April 1999. [64] E. Schwarz, L. Sigal, and T. McPherson. Cmos floating-point unit for the s/390 parallel enterprise server g4. IBM J. Res. Develop, 41(4/5):475–488, July/September 1997. [65] E. M. Schwarz, R. M. AverillIII, and L. J. Sigal. A radix-8 cmos s/390 multiplier. In 13th Symp.on Computer arithmetic, Asilomar, CA, pages 2–9, July 1997. [66] E. M. Schwarz, M. A. Check, C. L. K. Shum, T. Koehler, S. B. Swaney, J. D. MacDougall, and C. A. Krygowyski. The micro architecture of the ibm eserver z900 processor. IBM Journal Research and Development, 46(4/5):381–394, July/September 2002. [67] E. M. Schwarz and M. J. Flynn. Hardware starting approximation for the square root operation. In Proc 11th Symp.on Computer arithmetic, pages 103–111, 1993. [68] E. M. Schwarz and C. A. Krygowski. The s/390 g5 floating-point unit. IBM J. Res. Develop, 43(5/6):707–721, September/ November 1999. [69] P.-M. Seidel. High-speed redundant reciprocal approximation. Integration the VLSI journal, 28:1–12, 1999. [70] H. C. Shin, J. A. Lee, and L. S. Kim. A minimized hardware architecture of fast phong shader using taylor series approximation in 3d graphics. In Proc 128 BIBLIOGRAPHY Int’l Conference on Computer Design, VLSI in Computers and Processors, pages 286–291, 1998. [71] P. Soderquist and M. Leeser. Area and performance tradeoffs in floating point divide and square root implementations. ACM Computing Surveys, 28(3):518–564, 1996. [72] N. Takagi. Powering by a table look-up and a multiplication with operand modification. IEEE Trans.Computers, 46(11):1216–1222, 1998. [73] P. T. P. Tang. Table-driven implementation of the logarithm function in ieee floating-point arithmetic. ACM Trans. Mathematical Software, 4(16):378– 400, December 1990. [74] P. T. P. Tang. Table Look-up Algorithms for Elementary Functions and their Error analysis. Argonne National Lab. Report, MCS-P194-1190, January 1991. [75] A. T. Tsang and M. Olschanwosky. A Study of Database 2 Customer Queries. Technical Report TR-03.413. IBM Santa Teresa Laboratory, April 1991. [76] D. Wong and M.Flynn. Fastdivisionusingaccurate quotientapproximations to reduce the number of iterations. IEEE Trans. Computers, 41(8):981–995, August 1992. [77] W. F. Wong and E. Goto. Fast hardware-based algorithms for elementary function computations. IEEE Trans.Computers, 43(3):278–294, March 1994. [78] K. C. Yeager. The mips r10000 superscalar microprocessor. IEEE Micro, 16(2):28–40, 1996. [79] M. Zhang, S. Vassiliadis, and J. Delgado-Frias. Sigmoid generators for neural computing using piecewise approximations. IEEE Trans.Computers, 45:1045–1050, 1996.