Morgan - Numerical Methods (523161), страница 32
Текст из файла (страница 32)
E. Seminumerical Algorithms. Reading, MA: AddisonWesley Publishing Co., 1981, Pages 300-312.216PreviousCHAPTER 6The Elementary FunctionsAccording to the American Heritage Dictionary, elementary means essential,fundamental, or irreducible, but not necessarily simple. The elementary functionscomprise algebraic, trigonometric, and transcendental functions basic to manyembedded applications. This chapter will focus on three ways of dealing with thesefunctions: simple table-driven techniques that use linear interpolation; computational fixed-point, including the use of tables, CORDIC functions, and the TaylorSeries; and finally floating-point approximations.
We’ll cover sines and cosinesalong with higher arithmetic processes, such as logarithms and powers.We’ll begin with a fundamental technique of fixed-point arithmetic—tablelookup—and examine different computational methods, ending with the morecomplex floating-point approximations of the elementary functions.Fixed Point AlgorithmsLookup Tables and Linear InterpolationIn one way or another, many of the techniques in this chapter employ tables. Thealgorithms in this section derive their results almost exclusively through tablelookup. In fact, you could rewrite these routines to do only table lookup, if that is allyou require.Some of the fastest techniques for deriving values involve look-up tables.
Asmentioned in Chapter 5, the main disadvantage to table-driven routines is that thetables are finite. Therefore, the results of any table-driven routine depends upon thetable’s resolution. The routines in this section involve an additional step to helpalleviate these problems: linear interpolation.The idea behind interpolation is to approximate unknown data points based uponinformation provided by known data points.
Linear interpolation attempts to do this217HomeNextNUMERICAL METHODSby bridging two known points with a straight line as shown in Figure 6-1. Using thistechnique, we replace the actual value with the function y=f(x), where y = y0+(xx0)(y1-y0)/(x1-x0). This formula represents the slope of the line between the twoknown data points, with [f(x1)-f(x0)]/(x1-x0) representing an approximation of thefirst derivative (the finite divided difference approximation).
As you can see fromFigure 6-1, a straight line is not likely to represent the shape of a function well andcan result in a very loose approximation if too great a distance lies between each pointof known data.Figure 6-1. Linear interpolation produces an approximate value based on a straight linedrawn between known data. The closer the data points, the better the approximation.Consider the problem of estimating log10(2.22) based on a table of common logsfor integers only. The table indicates that log10(2) = 0.3010299957 and log10(3) =0.4771212547. Plugging these values into the formula above, we get:y=0.3010299957+(2.22-2.0)(0.4771212547-0. 3010299957) / (3-2y=0.3010299957 + (.22) (0.1760182552)/(l)y=0.3397540118.218THE ELEMENTARY FUNCTIONSThe true value is 0.3463529745, and the error is almost 2%.
For more accuracy,we would need more data points on a finer grid.An example of this type of table-driven approximation using linear interpolationis the function lg10, presented later in this section.1 The derivation of the table usedin this routine was suggested by Ira Harden. This function produces log10(X) of avalue based on a table of common logarithms for X/128 and a table of commonlogarithms for the powers of two.
Before looking the numbers up in the table, itnormalizes each input argument (in other words, shifts it left until the mostsignificant one of the number is the MSB of the most significant word) to calculatewhich power of two the number represents. The MSB is then shifted off, leaving anindex that is then used to point into the table.If any fraction bits need to be resolved, linear interpolation is used to calculatea closer approximation of our target value.
The log of the power of two is then addedin, and the process is complete.The function lg10 approximates log10(X) using a logarithm table and fixed-pointarithmetic, as shown in the following pseudocode:lg10: Algorithm1. Clear the output variable. Check the input argument for 0.If zero, exitIf not, check for a negative argument.If so, exitIf all OK, continue with step 2.2 . Determine the number of shifts required to normalize the inputargument, that is so that the MSB is a one. Perform that shift firstby moves and then individual shifts.3 .
Perform linear interpolation.First get the nominal value for the function according to the table.This is the f(x0) from the equation above. It must be guaranteed tobe equal to or less than the value sought.Get the next greater value from the table, f(x1). This isguaranteedto be greater than the desired point.Now multiply by the fraction bits associated with the number we usingto point into the table.
These fraction bits represent the difference219NUMERICAL METHODSbetween the nominal data point, x0, and the desired point.Add the interpolated value to the nominal value and continue withstep 4.4. Point into another table containing powers of logarithms using thenumber of shifts required to normalize the number in step 2.
Add thelogarithm of the power of two required to normalize the number.5. Exit with the result in quadword fixed-point format.lg10: Listing; *****.data; log (x/128);To create binary tables from decimal, multiply the decimal;value you wish to use by one in the data type of your;fixed-point system. For example, we are using a 64-bit fixed;point format, a 32-bit fraction and a 32-bit integer.
In;this system, one is 232, or 4294967296 (decimal), convert;the result of that multiplication to hexadecimal and you are;done. To convert p to our format we would multiply 3.14 by;4294967296 with the result 13493037704 (decimal), which we;then convert to the hexadecimal value 3243f6a89H.log10 tbl word00000h,0036bh,0085dh,00d18h,011a2h,015feh,000ddh,00442h,0092ah,00dddh,0125fh,016b4h,001b9h, 00293h,00517h, 005ebh,009f6h, 00ac1h,00ea0h, 00f63h,0131bh, 013dSh,01769h, 0181ch,00Gbdh,00b8ah,01024h,0148fh,018cfh,0078eh,00c51h,0l0e3h,01547h,01980hword01a30h,01e3bh,02222h,025e7h,01adfh,01ee4h,022c5h,02685h,01b8dh,01f8ch,02367h,02721h,01c3ah,02033h,02409h,027bdh,01ceGh,020d9h,024a9h,02858h,01dg1h,0217eh,02548h,028f3hword0298ch,02d14h,03080h,033d1h,0370ah,02a25h,02da8h,0310fh,0345ch,03792h,02abdh,02e3bh,0319eh,034e7h,03818h,02b54h,02ecdh,0322ch,03571h,0389eh,02beah,02f5fh,032b9h,035fah,03923h,02c7fh,02ff0h,03345h,03682h,039a8hword03a2ch, 03abOh, 03b32h, 03bb5h, 03c36h, 03cb7h,03d38h, 03db8h, 03e37h, 03eb6h, 03f34h, 03fb2h,220THE ELEMENTARY FUNCTIONSword0402fh,04312h,045e3h,048a2h,04b50h,040ach,0438ch,04659h,04915h,04bc0h,04128h,04405h,046cfh,04988h,04c31h,041a3h,0447dh,04744h,049fbh,04ca0h.0421eh,044f5h,047b9h,04a6dh,04d10h04298h,0456ch,0482eh04adeh,;log(2**x)log10_power dword 000000h, 004d1Oh, 009a20h, 00e730h, 013441h, 018151h,01ce61h, 021b72h, 026882h, 02b592h, 0302a3h, 034fb3h,039cc3h, 03e9d3h, 0436e4h, 0483f4h, 04d104h, 051e15h,056b25h, 05b835h, 060546h, 065256h, 069f66h, 06ec76h,073987h, 078697h, 07d3a7h, 0820b8h, 086dc8h.
08bad8h,0907e9h, 0954f9h.code;;Logarithms using a table and linear interpolation.;Logarithms of negative numbers require imaginary numbers.;Natural logarithms can be derived by multiplying result by 2.3025.;Logarithms to any other base are obtained by dividing (or multiplying by the;inverse of) the log10.
of the argument by the log10 of the target base.lgl0 proc uses bx cx si di, argument:word, 1ogptr:wordlocalpowers_of_two:bytepushfstdrep;increment down for zero;check to comesubmovmovaddstoswax,cx,di,di,ax4word ptr logptr6movsi, word ptr logptraddmovaddmovorjssubsi, 6di, word ptr argumentdi, 6ax, word ptr [di]ax, axexitax, ax;clear log output;point at output which is;zero;most significant word;point at input;most significant word;we don't do negatives221NUMERICAL METHODSrepemovcmpswcx, 4jeexit;find the first nonzero,;or return;zeroreposition_argument:movsi,addsi,di,movinccxmovax,subax,repshlsubshlshlshlmovmovswmovmovkeep_shifting:orjsshlrclrclrclincjmpdone_with_shiftmovmovsubmov2 22ax,si,ax,ax,ax,bl,word ptr argument6si4cx1ax111al;shift so MSB is a one;point at input;most significant word;shift the one eight times;make this a one;determinenumberof;emptywords;words to bytes;point to first nonnero word;multiply by eight;shiftsi, word ptr argumentax, word ptr [si][6]ax, axdone_with_shiftword ptr [si][0], 1word ptr [si][2], 1word ptr [si][4], 1ax, 1blshort keep_shiftingword ptr [si][6],axbyte ptr powers_of_two, blbx, bxbl, ahshlbl, 1addbx, offset word ptr log10_tbl;shift until MSB is a one;count shifts as powers;of two;normalize;ax will be a pointer;will point into 127-entry;table;get rid of top bit to form;actual pointerTHE ELEMENTARY FUNCTIONS;linear interpolation;get first approximation;(floor)movax, word ptr [bx]incincmovbxbxbx, word ptr [bx]subxchgbx, axax, bxmulmovSubaddbyte ptr [si][6]al, ahah, ahax, bx;multiply by fraction bits;drop fractional placesbl, 31;need to correct for power;of twoget_power:movsubsubshlshlleaaddsubaddadcmovmovbl, bytebh, bhbx,1bx,1si, wordsi, bxdx, dxax, worddx, worddi, wordword ptrmovsubmovmovword ptr [di][4],dxcx, cxword ptr [di],cxword ptr [di] [6],cx;and following approximation;(ceil);find difference;add interpolated value to;originalptr powers_of_two;point into this tableptr log10_powerptr [si]ptr [si][2]ptr logptr[di] [2],ax;add log of power;write result to qword;fixed pointexit:popfretlg10 endpAn example of linear interpolation appears in the TRANS.ASM modulecalled sqrtt.223NUMERICAL METHODSAnother example using a table and linear interpolation involves sines andcosines.
Here we have a table that describes a quarter of a circle, or 90 degrees, whichthe routine uses to find both sines and cosines. Since the only difference is that oneis 90 degrees out of phase with the other, we can define one in terms of the other (seeFigure 6-2). Using the logic illustrated by this figure, it is possible to calculate sinesand cosines using a table for a single quadrant.To use this method, however, we must have a way to differentiate between thevalues sought, sine or cosine.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
