SOF'I'WARE-PRACTICE AND EXPERIENCE, VOL. 8, 171-176 (1978) Rational Arithmetic for Minicomputers BERTHOLD K. P. HORN Artificial Intelligence Laboratory, Massachusetts Institute of Technology, 545 Technology Square, Cambridge, Mass. 02139 U.S.A. SUMMARY A representation for numbers using two computer words is discussed, where the value represented is the ratio of the corresponding integers. This allows for better dynamic range and relative accuracy than single-precision fixed point, yet is less costly than floating point arithmetic. The scheme is easy to implement and particularly well suited for minicomputer applications that call for a great deal of numerical computation. The techniques described have been used to implement a mathematical function subroutine package for a minicomputer as well as a number of applications programs in the machine vision and machine manipulation area. KEY WORDS Rational arithmetic Minicomputers Number representation INTRODUCTION Typically, minicomputers come equipped with arithmetic instructions that operate on single words. Frequently double-word addition and subtraction is also provided for. Unfortunately, integers are inconvenient for many calculations and require the programmer (or at least the compiler) to remember scale factors relating the number represented to the integer actually stored. Even with scale factors one frequently runs out of dynamic range or has difficulty because of the need to estimate the approximate magnitude of variables ahead of time. Floating-point representation is the obvious, but costly answer. Software to perform the common floating-point operations requires substantial amounts of memory and is usually very slow, while hardware that carries out the same operations-if availableis expensive. RATIONAL A R I T H M E T I C An oft-overlooked alternative that lies between fixed-point and floating-point arithmetic in complexity, cost and capability is rational arithmetic (see Figure 1). Variables are stored in two words and the number represented is the ratio of the integers in these words. This allows a fair dynamic range and possibly higher precision than single-word integer arithmetic, without the complexity of extracting exponents and mantissas as required for floating-point arithmetic. I t is also easy to provide software that performs the usual arithmetic operations at relatively high speed. Since the number of instructions required is relatively small, one may elect to compile them in line instead of providing them as subroutine calls. Similarly, it should be easy to microprogram a machine to provide these operations directly as part of its instruction repertoire. Notice that we are not discussing exact rational arithmetic which requires the 0038-6644/78/02084171 SO1.OO @ 1978 by John Wiley & Sons, Ltd. Received I5 July 1977 171 172 BERTHOLD I<.P. HORN COMPUTER REPRESENTATION MATHEMATICAL ABSTRACTION FIXED POINT 'NO POINT' FLOATING POINT INTEGERS RATIONALS REAL NUMBERS Figure 1. Analogy between mathematical nrrmber systems and computer representations, The proposed 'no point' representation corresponds to the rational numbers ability to represent two integers of arbitrary size.lp2 Our two-word representations for numbers will be approximations in the same way that floating-point numbers are approximations to real numbers. T h e minimal capabilities required to implement rational arithmetic in this sense are single-word multiplications (with double-word results), double-word addition, subtraction and shifting. I t is helpful to have four registers available to hold partial results. Addition, subtraction, multiplication or division of rational numbers produces a pair of doubleword results which have to be 'compressed' into two single words before they can be used in further calculations. NORMAL1 ZATION T h e main decision to be made when implementing rational arithmetic is the choice of an appropriate 'normalization' procedure, which will find a rational number that can be represented in the two-word format, and which approximates the double-size (four-word) result of the elementary arithmetic operations. T h e fastest method is certainly one where both double-word results are simply shifted right until each fits into one word. I t must be pointed out that this does not wually result in the best approximation possible in the two-word format. In fact, it is fairly easy to see that the relative accuracy of this method is no better than that obtained in single-precision integer arithmetic. Its simplicity however makes this scheme very attractive in situations where the extended dynamic range and not the higher accuracy is of most interest. Next, we notice that one can divide both terms of the result by any common factors they might have (provided the computer has a divide instruction of course). With luck this will lead to sufficient reduction in the size of the results so that they will now fit into the twoword format. I n this case the answer is exact, but some work has to be done to find the greatest common divisor in the first place (see Figure 2). T h e number of iterations in the GCD algorithm is approximately equal to the logarithm of the n ~ m b e r s One . ~ certainly cannot rely on this method alone always to provide the required compression of the two double-word results. This is the appropriate method however for exact rational arithmetic.1~2 GCD(N, M): I F NI = 0,THEN N ELSE GCD(M, REMAINDER(N, M)) Figure 2. A simple algorithm for finding the greatest common divisor of two numbers. This algorithm is a slkht modification of Euclid's, and more elaborate, faster algorithms exist. The function R E M A I N D E R determines the remainder of the indicated division. That is, R E M A I N D E R ( N , M ) = N - LN/~VIJ*iVI. The GCD algorithm is useful for normalization after arithmetic operations or for simplifying terms before they enter into arithmetic. The tail-end recursion can of course be turned into a more eflcient iteration T h e Euclidian algorithm for finding the greatest common divisor suggests a related method for finding optimal approsimants using number theorctic ideas. T h e method is RATIONAL ARITHMETIC FOR MINICOMPUTERS 173 based on an algorithm for finding successive convergents of a continued fraction expansion (see Figure 3). T h e method is described in detail in Reference 3 and reproduced in FRACTS(S,N): L E T P l + O ; Q ~ + l ; P g + l ; Q ~ + O ; I + O DO UNTIL P 2 > N O R Q,>NOR S = I LETPo+P1; Qo+Ql; Pl+P2; Ql+Q2 LET I <- INTEGER-PART@) LET P2+-Po+I *PI; Qz+Qo+I *Q1 IF S # I, THEN LET S + 1/(S - I), ELSE LET PI + P,; Q1+ Qz END END Figure 3. A simple algorithnr for deternrining good rational approximations to a given nunrber S. The successive approxinrants (P2/Q2)are alternately larger than S and smaller than S . The algorithm terminates wlren either numerator or denominator exceeds the given bound N, or when the ratio exactly equals the number S. The algorithm is based on a method for finding continued fraction expansions Reference 4. This technique involves more computation, but since a result of the form p/q is accurate to within l/q" we are now dealing with a method that has similar relative accuracy as one might expect from double-precision integer arithmetic (see Figure 4). Figure 4. Illustration of the theorem about cunvergents. The grid points correspond to rational numbers. Those falling zoithin the shaded triangle represent approximations to the number S. These rational numbers have a denominator less than N and lie between S - E and S E in value. If the area of the triangle, E N 2 exceeds one, we expect to find one or more grid-points within its boundaries, since there is one grid-point per tinit area. This suggests that a number S can be approximated with accuracy I/N2 with rational nunzbers that have denominator less than or equal to N + T h e largest approximant that fits into the two-word format has been called the mediant conversion value.4 Fortunately, the algorithm produces the approximants in order of increasing size, and so one simply keeps going until one of the two factors no longer fits into a single word. See Figure 5 for examples of small approximants for common mathematical constants. 174 BERTHOLD K. P. HORN RATIONAL APPROXIMATION MATHEMATICAL CONSTANT Figure 5. Rational approximants for some common mathematical constants, accurate to at least 16 bits or about 5 decimal digits I We have presented here two extreme case algorithms for the normalization, the one exceedingly fast, but of limited accuracy, the other complicated, but with exceptional accuracy. Other points on the spectrum of possibilities can be explored. For example, a compromise acceptable for many real applications is a technique which involves simply dividing both double-word result terms by the first (high-order) word of the larger plus one. This ensures that the two resulting numbers are as large as possible. Clearly the method is not applied if the results already occupy only a single word each. All of the above methods involve a certain amount of bookkeeping to deal properly with negative numbers and the idiosyncrasies of the particular scheme used to represent negative numbers. I ALGORITHMS Let n[a, b] stand for the normalization operation just discussed and let us represent a pair of words containing the numbers a and b as (alb). Evidently, the four usual arithmetic operations produce the following results: + (alb) (cld) +n[ad +be, bd] (alb) - (cld) +n[ad- be, bd] (alb) x (cld) +n[ac, bd] (alb) + (cld) +n[ad, be] Curiously, with this representation, multiplication and division involve a little less effort than addition and subtraction. Also notice that division by zero does not lead to immediate disaster. For reasons of efficiency one may want to retain integers represented by single words. T h e conversion between such single words and the double-word rational representation is straightforward of course and may be used to implement mixed arithmetic using a form of 'rational contagion'. That is, integers about to enter into operations with rationals are converted to rational form first. Alternatively, one can easily code the mixed arithmetic operations separately in order to gain speed by avoiding futile multiplications by one. Conversion between an external decimal representation and the internal rational format for input and output purposes requires division, but is also simple. MATHEMATICAL F U N C T I O N SUBROUTINES Most common methods for computing elementary transcendental functions involve argument reduction and approximation by a polynomial or a ratio of two low-order polynomials in the independent variable. Such techniques can easily be implemented for I * RATIONAL ARITHMETIC FOR MINICOMPUTERS 175 minicomputers using rational arithmetic. This is sensible even when the rest of the program uses other representations since the resultant subroutines are comparatively fast and short. In this case one can often simplvy the normalization operation after arithmetic steps since the magnitudes of the various coefficients and the range reduced input are known. An additional idea worth exploring when using trigonometric functions is a representation for angles which lets 360" correspond to the largest integer that can be represented by a single word, plus one. This provides for automatic wrap-around of angles larger than 360" or less than zero on most machines. I t also makes for the best possible angular resolution for a given word size. Further, it is a good idea to let arc-tangent be a function of two variables. T h e result is the arc-tangent of the ratio, with the quadrant picked according to the signs of the two arguments. This avoids the usual two-way ambiguity in the result of inverse trigonometric functions. Such a scheme also avoids difficulties with angles near 90" and 270" for which the argument would otherwise become excessively large. COMPARISON W I T H OTHER REPRESENTATIONS The rational arithmetic method has already been compared with the usual single precision fixed-point and floating-point methods. A comparison with two other simple schemes, using two words to represent a number, may be called for. T h e first is a double-word representation with an imaginary binary point between the two words. That is, the first word is the integer part and the second the fractional part. This representation has the same dynamic range and comparable accuracy, but while addition and subtraction are simpler, the other operations, and particularly division, are much more difficult to is less implement. Secondly, we could consider a simple floating point form, with one word used for the exponent, the other for the mantissa. While this avoids the difficulty of extracting substrings from words inherent in the usual floating-point representations, it is also wasteful. T h e dynamic range of course is larger than it is for the rational representation, but the accuracy is less and this technique is also more difficult to implement. CASE S T U D Y A N D SUMMARY T h e techniques described here have been used to implement a mathematical function package for a P D P l l minicomputer. T h e approximations shown in Figures 6, 7 and 8 were used, the results are accurate up to the limitations of the word-length of the machine and the execution times are short (160 ps for ATAN, 185 ps for SIN/COS and 115 ps for S Q R T on a PDP11/40). I n this case, the normalization problem was simple since the expected ranges of the numbers involved were known in advance. Similar techniques have also been employed in a number of programs involving applications of machine vision and machine manipulation techniques." FOR ( X 1 < 1 / 2 Figure 6 . Rational function approJcimations for A T A N ( X ) for range-reduced X . The first approximation is obtained from an application of Aitken series acceleration, the second from continued fraction expansion. All coeficients are rational numbers and the results have suficient accuracy for 16-bit mim'combuters. The results are in radians, but can easily be converted into f r a & n a l revoluti& representation by multiplying by 3551113 176 BERTHOLD K. P. HORN X ( l - X 2/6 (1 - X 2/20)) E S I N (X) FOR I X I <m/4 1-X2/2(1-X2/12(1-X2/30))~COS(X) FORIXI<m/4 Figure 7. Truncated Taylor series approximations for sine and cosine for range-redrrced A ' (in radians). The ~esrrltsare of nrmient accuracy for 16-bit conzputers. A single function producing both values acts as a convenient inverse to the proposed arc-tangent function which takes tzco arguments SQRT(Y): L E T Xl t 0 ; Xo t 1 D O U N T I L Xl = Xo OR XI = Y/Xo L E T Xot X1; X1t (Xo Y/Xo)/2 END RETURN X, + END Figure 8. Simple Newton-Raphson double-precision sqrrare-root algorithm suitable for miniconrputer rise. Divisions are intended to produce integer quotients. Division by two can obviously be implemented rising a right shift. Other initial guesses, such as Y, Y/2 or (Y+ 1)/2 corrld be used instead of 1 T h e 'no-point' or rational representation for numbers has been shown to make for easy-to-implement arithmetic operations, yet is effective in providing capabilities somewhere between those available with the more traditional fixed- and floating-point schemes. T h c low complexity and cost in terms of memory and execution time recommend this method for applications involving minicomputers. Users of minicomputers need no longer fear matrix arithmetic, coordinate transformations, evaluations of polynomials or calculation of transcendental functions since rational arithmetic will allow them to do a fair job of tackling these tasks, while requiring only a modest effort. I t is the rational choice! REFERENCES 1. P. Henrici, 'A subroutine for computations with rational numbers',J. A C M , 3, 10-15 (1956). 2. D. E. Knuth, The Art of Computer Progranmnzing, Volume 2, Semi-numerical Algorithms, AddisonWesley, Reading, Mass., 1969, pp, 290-292. 3. G. H. Hardy and E. M. Wright, A n Introdrrction to the Theory of Numbers, Clarendon Press, Oxford and London, 1934, pp. 137-138. 4. D. W. Matula, 'Fixed-slash and floating-slash rational arithmetic', Proc. 3rd Symp. on Comprrter Arithmetic, I E E E , 90-91 (1975). 5. B. I<. P. Horn, 'A problem in computer vision: orienting integrated circuit chips', Comprrter Graphics and Image Processing, 4, 294-303 (1975). View publication stats L e