Power-Programmierung

home *** CD-ROM | disk | FTP | other *** search

/ Power-Programmierung / CD1.mdf / forth / compiler / fpc / doc / sfloat.txt < prev next >

Wrap

Text File | 1988-07-16 | 15.1 KB | 302 lines

Implementation Notes for SFLOAT, a software Floating Point package for 80X8X processors By Robert L. Smith INTRODUCTION: This document contains various notes, thoughts, methods, etc. which relate to the software floating-point package which I have called SFLOAT. My main motivation for writing such a package is to help eliminate a common complaint about FORTH, that it isn't useful because it doesn't have floating point. There are various possible answers to that complaint. One answer is that floating point is available in most commercial FORTH systems. On the other hand, most people initially try FORTH in a public domain package, such as F-83. In that case, they will probably not find a floating point package. (One exception is a hardware floating point package for the PC clones. Another exception is a low precision "Zen" floating point published by Martin Tracy.) Another answer to the complaint is that floating point is not needed. In a pure sense, that is certainly true. However, there are many times when the human effort required to do a problem in fixed point arithmetic becomes unreasonable for the problem at hand. In my particular case, I once wasted a great deal of time trying to approximate, in fixed point, a function involving the probability integral. The problem was trivial in floating point. The floating point routines in this package are written for the most part in assembly language (CODE) routines for a specific processor type, the infamous Intel 8088 family. Why did I write routines in code rather than in high level FORTH? For the usual reason of speed. In addition, detailed arithmetic functions are somewhat easier to write in code rather than high level because of the absence of the carry bit in FORTH. I chose the Intel series because of its popularity. An equivalent package for the Motorola 68XXX series should be much easier to write because of the 32 bit architect- ure. I will leave that to others. It was my intent to produce a set of floating point routines which are both fast and accurate. If you find errors or places where I can improve things, please let me know. REPRESENTATION: I was a member of the IEEE P754 committee on Floating Point. At one point I was strongly inclined to use the IEEE floating point representation. There were two reasons why I rejected it. The first reason was that I tried, and found it very difficult to work with. The second reason is that someone else had built a floating point package around the IEEE format (Philip Koopman, Jr., MVP FORTH). According to reports from Glenn Haydon, the routines were painfully slow. Every representation method has its advantages and disadvantages. I chose a 48 bit representation. The most significant bit is the sign bit. The next 15 bits are used to indicate the power of 2, and the remaining 32 bits are used for the significand (or mantissa, as we used to say), without suppressing the leading bit. The binary point is at the left of the leading significant bit. This representation gives about 9 1/2 decimal digits of accuracy (somewhat greater than the IEEE single precision representation) and a range of 10^-4000 to 10^4000, roughly. That is far more than anyone needs, but greatly reduces the need to test for overflow and underflow. I chose to consider any unnormalized number (with the most significant bit cleared) as zero. Finally, I set the bias on the exponent or power of 2 to 3FFF (Hex representation, obviously). I am not completely satisfied with the last two decisions (representation of zero, and exponent bias), but I don't have the energy to re-do all of the work just to see if the decision should be altered. The current version uses a separate floating point stack. In my initial version, I used the standard parameter stack to hold working floating point numbers. This is allowed in the Forth Vendors Group Standard Floating Point Extension. As of this writing, it appears that the ANSI Standard Forth will specify a separate Floating Point stack, so I re-wrote the routines. ACCURACY: In the IEEE Floating Point Standard, the basic four functions (add, subtract, multiply, and divide) and square root are required to be exact within one-half of the least significant bit. It is my intention to meet that requirement. To understand how that is done, the reader should consult some of the original papers. See, for example, "The INTEL Standard for Floating Point," by John F. Palmer, Proceedings Compsac77, IEEE Catalog No. 77CH1291-4C. The basic idea is that all bits below the least significant bits are to be considered. If these bits represent a number greater than one-half of the least significant bit, then the result is incremented by 1. If the bits represent a number less than one-half, the result is unchanged. If the bits represent exactly one-half, then the least significant bit is examined. If the least significant bit is 0, the result is unchanged. Otherwise the result is incremented by 1. This process is called "rounding to even." The trickiest part of this process occurs when subtracting two numbers which have the same sign (or adding two numbers of opposite sign). I know the casual reader of these notes is going to think that going to all of this effort is gross overkill, but the mathematical literature has numerous cases where this is important. In any case, if you are going off on your own to do floating point arithmetic, please do not simply truncate the results unless you have very good reasons for doing so. Some of the very worst examples come from computers designed by Seymore Cray, whose philosophy seems to be that the results don't have to be computed accurately as long as they can be computed rapidly. DIVISION If you are interested in details, the floating point division routine in this package is fairly interesting. It is my second version. The initial version used a more or less standard technique for dividing by successive subtraction and shifting. I stumbled on an interesting variation for Knuth's method for divsion in the documentation of the Forth-based language called "Abundance" by Roedy Green. In that package, he performs a double precision divide by first shifting the divisor left until the most significant bit is set, obtaining an unsigned quotient and remainder, then adjusting the result based on the initial amount of shifting. In a floating point routine, the operands are almost always normalized, so that the shifting is unnecessary, except possibly for an initial right shift of the numerator by one place. The method is somewhat similiar to the technique that you learned in grammar school, except that instead of working in base ten we use base "s", where s = 2^16 . In other words, each 16 bit word is considered as a digit in this new base. We desire a 2 digit quotient and a two digit remainder. We will first guess the leading digit of the quotient, then determine its exact value, along with an appropriate remainder. Then we will determine the second digit and final remainder the same way. It turns out that our initial quess for each quotient digit must be no more than 2 greater than the correct digit. At each stage the numerator consists of 3 digits, and the denominator consists of 2 digits. We may define the first quotient digit and remainder as follows: (s^2)*A + s*B + C r1 ----------------- = q1 + ------ sD + E sD + E The 3 terms in the numerator account for a possible right shift in the case that the initial numerator was greater than or equal to the denominator. We thus can guarantee that A <= D, and q1 < s. We also require that the remainder r1 < (sD + E). With a small amount of algebraic manipulation we can see that r1 = (s^2)*A + s*B + C - s*q1*D - q1*E That doesn't appear to be of much help, but just hold on. We will make our initial guess for q1 based on a single precision division as follows: s*A + B r0 ------- = g0 + --- D D The key to simplifying this calculation is that we do not ignore the remainder in our trial division. Note that the remainder in this division is r0 = s*A + B - g0*D Assume for the time being that q1 = g0. It may not be, but our equations will still be exact! We can calculate a remainder r1' which is r1' = s*r0 + C - g0*E by making a simple single precision multiply and a double precision subtraction. Our initial guess g0 may have been too large. We can tell that because then r1' would be negative. If r1' is >= 0, our initial guess is the correct guess. If r1' is negative, we may correct for that by simultaneously adding the denominator (sD + E) to r1' and decrementing g0 by 1. We check again to see if the new remainder r1'' is negative. If r1'' is also negative, make a final correction by adding the denominator to r1'' and decrementing g0 again. This final quotient is correct, and the remainder is also correct. The process is then repeated for the second quotient digit, by using the first remainder as the new numerator and then estimating the quotient and remainder as above. The result is that we can calculate the double precision quotient by using 2 single precision divides and two single precision multiplies, and some trivial additions and subtractions. The final part of this process is to round the quotient precisely. This is done by comparing twice the remainder with the divisor. If 2*rem is less than the divisor, the quotient is correct. If 2*rem is greater than the divisor, the quotient must be incremented by 1. It 2*rem is exactly equal to the divisor, we round to even. That is, if the least significant bit in the quotient is 1, then add 1 to the quotient. Otherwise, do not add anything to the quotient. There is an additional detail in the above: it is possible that A = D , even though (sA + B ) < (sD + E) . If we were to just divide, we would encounter an overflow situation. Therefore, consider that the trial quotient g0 = s - 1 . The zero-th remainder is r0 = sA + B - (sD + D) , and the first remainder is r1 = (s^2)*A + s*B + C - (s^2)*D + s*(D - E) + E = s*(B + D - E) + (C + E) Again, the initial guess may be too large, which may be discovered by checking if r1 is negative. If it is, then one or two corrections may need to be made, as in the previous discussion. LOGARITHM These notes are for those few souls who are interested in the details of implementation of the FLN function for software floating point. As noted elsewhere, the basic format for floating point is a 16 bit sign and biased exponent, followed by a 32 bit normalized fraction. The biased exponent occupies bits 0 through 14 with a bias of 3FFF (hexadecimal), and the sign bit is bit 15 of the first 16 bits. The binary point of the fraction is just to the left of the 32 bit fractional part (sometimes called the mantissa). If the most significant bit of the mantissa is zero, the whole number is considered zero. If the argument to the FLN function is negative or zero, it is an error condition. After testing for those conditions, we then break down the problem into two parts, depending on the magnitude of the fraction. If the fraction is greater than or equal to 0.78125 (and, by definition, it must be less than 1.0), then we note that ln(X) = ln((2^I) * 0.f) = I*ln(2) + ln(0.f) where I is the unbiased exponent and 0.f is the fractional part. Calculation of ln(0.f) is based on the approximation: ln(1-x) = -(x + (x^2)/2 + (x^3)/3 + (x^4)/4 + . . . ) It may be noted in passing that there are better approximations than the above, but this one is selected to allow a crucial division to be made using single precision divisors. It is necessary to reduce the range of x in the above equation. (If we did not, it would be necessary to use about 24 terms in the series when x = 0.25). This is accomplished by dividing 0.f into intervals of 2^-5 or 1/32 . The procedure is as follows: 0.f = 1 - f1 , where f1 = 1 - 0.f f1 * 2^5 = J + f2 * 2^5 0.f = 1 - f1 = 1 - (2^-5) * J - f2 = (1 - J*(2^-5)) * [1 - f2/(1 - J*(2^-5))] where J is an integer. The calculation can proceed by using a table of values for the logarithm of (1 - J*(2^-5)), and an approximation for ln(1-x), where x = f2/(1 - J*(2^-5)) In this case, the maximum value for x is 0.0435 . The series approximation for 33 bits of accuracy can terminate at the term (x^7)/7 . In order to maximize the accuracy of the calculation, we proceed using y = x*2^4 . The maximum value of y is 0.696 . We then proceed to calculate: (2^4)*(-ln(1-x)) = y + (y^2)*(2^-5) + (2/3)*(y^3)*(2^-9) + + (y^4)*(2^-14) + (4/5)*(y^5)*(2^-18) + + (4/6)*(y^6)*(2^-22) + (8/7)*(y^7)*(2^-27) The last term is kept only for rounding purposes. The factor (8/7) may be approximated by 1. It is necessary to use double precision in only the first 4 terms (well, maybe only the first 3, with extreme care). The above equation should be factored, and the calculations made from right to left so that the smallest terms are calculated first. (Writing equations in ASCII is tedious, and reading them is probably hellish, at best). The calculation is roughly as follows: -(2^4)*ln(1-x) = (((((y*(2^-5) + (2/3))y*(2^-4) + +(4/5))*(y^2)*(2^-4) + Y)*(2^-5) + (2/3))*(Y^3)*(2^-4) + +(Y^2))*(2^-5) + Y where Y indicates the double precision value and y the single precision value. The case of the initial fraction 0.f being less than 0.78125 proceeds along somewhat similiar lines to the above. We first note that the desired result may be expressed as ln(X) = ln(2^(I-1)*(2*0.f)) = (I-1)*ln(2) + ln(2*0.f) where 1.0 <= (2.0*f) < 1.5625 . In a manner similiar to the first case we divide the interval up into 19 subintervals of 2^-5 or 1/32 each. The proceedure is as follows: 2*0.f = 1 + f1 , where f1 = (2*0.f) - 1 f1 * (2^5) = J + f2 * (2^5) 2*0.f = (1 + J*(2^-5)) * [1 + f2/(1 + J*(2^-5))] where J is an integer. The calculation can proceed by using a table of values for the logarithm of (1 + J*(2^-5)) , and an approximation for ln(1+x), where x = f2/(1 + J*(2^-5)) In this case, the maximum value for x is 0.03125 . The approximation proceeds much as in the first case. The approximation used is: ln(1+x) = x - (x^2)/2 + (x^3)/3 - (x^4)/4 + (x^5)/5 - (x^6)/6 + (x^7)/8 Let x = (2^-5)Y . The calculation proceeds along the lines indicated by: (2^5)*ln(1+x) = Y - (2^-6)(Y + (2^-5)(Y^3)(2/3 - (2^-6)(y - (2^-5)* (y^2)(4/5 - (2^-5)y(2/3 - y*(2^-6))))