Morgan - Numerical Methods (523161), страница 35
Текст из файла (страница 35)
These series are often expressed in the followingforms:5sin x = x - x3/3! + x5/5! - x7/7! + x9/9! . . . +(-l)n+1x2n-1/(2n-l)!cos x= 1 - x2/2! + x4/4! - x6/6! + x8/8! . . . +(-l)n+1x2n-2/(2n-2)!tan x = x + x3/3 + 2x5/15 + 17x7/315 + . . .ex = 1 + x = x2/2! + . . . + xn/n! + . . .ln(1 + x) = x - x2/2 + x3/3 - .
. . + (-l)n+1+xn/n + . . .A power-series polynomial of infinite degree could theoretically accommodateevery wrinkle in the shape of a given function within a given domain. But it isn’treasonable to attempt a calculation of a series of infinite degree; instead, somemethod is used to determine when to truncate the series. Usually this is at the pointin the series where the terms fail to contribute significantly to the result. Yourapplication may only require accuracy to 16 bits, such as might be needed for247NUMERICAL METHODSgraphics.
It may be error limited, which means that the result is calculated usingenough precision and to a great enough degree to account for any spikes that mightoccasionally occur in the more distant terms.Since the power series are computed in truncated form, they are prone to an errorfrom that truncation as well as any introduced by the arithmetic. A great deal of efforthas gone into finding the source of those errors and limiting it.6 For most embeddedapplications (such as graphics subsystems, digital filtering and feedback controlloops), the truncated Taylor Series provides adequate results.The quadword fixed-point format used in this section has 32 fraction bits to workwith.
The terms contributing to bits outside that range (aside from guard digits, if youwish) are not computed even if an occasional spike might influence the rest of thecomputation. The 13th term of the sine expansion above rounds up to set the leastsignificant of our fraction bits.An alternative to the doubleword integer and doubleword fraction format couldbe implemented for each of the functions. At most, sine and cosine functions needa l-bit integer, leaving 63 bits for at least 18 decimal digits.
On the other hand, theexponential, ex, will quickly lose any mantissa bits unless x is less than one. Youcould rewrite these routines to maximize the precision of the data types you’re usingand provide greater accuracy; the results could be rounded and realigned for the restof the fixed-point routines.
You can do this without disturbing any de facto formatyou may have in place by doing the conversions and alignment within the callingfunction, as taylorsin below. Such handling is often the case anyway, since aparticular series may require the arguments in a certain format to guaranteeconvergence. The sine and cosine functions presented here are examples of this:Their arguments should be constrained to π/2 for the series to function mostefficiently and accurately.Power-series computations are not necessarily table driven, but the executiontime of the evaluation is so much faster when you precompute the coefficients thatyou need a good reason not to.
If you wish to compute the coefficients at runtime,it’s most efficient if you maintain a copy of the previous powers and factorials andcompute each new one based upon that.Homer’s Rule7 allows us to evaluate an N-degree polynomial with only N-lmultiplications and N additions. To use it, we store the coefficients of the polynomial248THE ELEMENTARY FUNCTIONSin an array. If a degree or series of degrees is missing from the polynomials, theircoefficients automatically become zero. To illustrate, assume a polynomial such asf(x) = 5x4 + 3x3 - 4x2 + 2x + 1We put the coefficients in an array:Poly_arrayword1, 2, -4, 3, 5In the following pseudocode, as in the example, the coefficients of the series (orpolynomial) are assumed to be computed in advance, incorporating the sign of theterm with the value.
They’re stored in a table in reverse order of the polynomialexpression; that is, the first element in the array is the degree zero term. Evaluationis then simply a matter of processing the polynomial. Upon entry to the algorithm,we make the result variable equal to the coefficient of the highest power (here it’s 5).We take a pointer into the array, which is the degree of the polynomial, and use it toselect each succeeding coefficient to add to the result variable after multiplying it bythe value of x.taylorsin: Algorithm1. Set an index to the degree of the polynomial (in this case 4).Use this to retrieve the coefficient of the highest power and set theresult variable equal to that.2. Multiply the value of x by the result variable,Decrement the index.If it goes negative, exit through step 3Retrieve the next coefficient and add it to the result variable,Continue at the beginning of step 2.3.
Horner's Method is complete. Exit.In taylorsin, the sine approximation given above truncated at the 11th degree for ourexample:249NUMERICAL METHODSsin x = x - x3/3! + x5/5! - x7/7! + x9/9! - x11/11!To process this expansion with Homer’s Rule, we need a table of coefficientswith 11 terms in it and zeros for those powers not represented in the expansionindecimal:1, 0, -.16666667, 0, .00833333, 0, -.00019841, 0, .00000275, 0,-.00000003Even this can be avoided if we evaluate the expression x3/3! + x5/5! - x7/7! + x9/9! - x1l/l1! separately with x2 instead of x. This eliminates the necessity of processingall the zero coefficients.With these terms stored in a table, the only thing left to do is evaluate thepolynomial.Actually, two routines are involved: polyeval can be made to work with anypolynomial, while taylorsin is only an entry point.
It tells the subroutine polyevalwhich table to use depending on the function to evaluate, the degree of thepolynomial, and where to put the results and passes the argument. Each functionrequiring polynomial evaluation will require a routine such as taylorsin; this is whereany other fixed-point manipulation-such as scaling, altering the placement of theradix point, or rounding-would be done.taylorsin: Listing;*****;taylorsin - Derives a sin by using a infinite series.
This is in radians.;Expects argument in quadword format; expects to returnthe same.;Input must be/2.taylorsinproc uses bx cx dx di si, argument:qword, sine:wordinvokepolyeval, argument, sine, addr polysin, 10rettaylorsin endpPolyeval does the work and can be made to evaluate any polynomial, given theproper coefficients. Here is how it works:250THE ELEMENTARY FUNCTIONSPolyeval: Algorithm1. Clear an accumulator and see that the output is clear.Set an index equal to the degree of the polynomial.2. Using the index, point into the table of coefficients.Add the value pointed at to the accumulator.3. Multiply the accumulator by the argument, x.4. Decrement the table pointer.If it goes negative, exit through step 5.Otherwise, continue with step 2.5.
Write the accumulator to the output and leave.Polyeval: Listing; *****.datapolysin qword100000000h, 0, 0ffffffffd5555555h, 0, 2222222h, 0,0fffffffffff2ff30h, 0, 2e3ch, 0, 0ffffffffffffff94h.code;; *****;polyeval- Evaluates polynomials according to Horner's rule.;Expects to be passed a pointer to a table of coefficients,;a number to evaluate, and the degree of the polynomial.;The argument conforms to the quadword fixed-point format.polyevallocalreppushfcldsubmovmovleastoswprocuses bx cx dx di si, argument:qword, output:word,coeff:word, n:bytecf:qword, result[8]:wordax, axbyte ptr sign, alcx, 4di, word ptr cf;clear the accumulator251NUMERICAL METHODSrepleamovstoswdi, word ptr resultcx, 8movsubmovsi, word ptr coeffbx, bxbl, byte ptr nshlshlshlbx, 1bx, 1bx, 1addmovmovmovmovleaaddadcsi, bxax, wordbx, wordcx, worddx, worddi, wordword ptrword ptreval:adcadcx_by_y:invokereplealeamovmovswchk_done:decjnspolyeval_exit:movleamovmovswrep252;point at table;point at coefficient of;n-degree;this is the beginning of;our approximation;multiply by eight for the;quadwordptr [si]ptr [si] [2]ptr [si] [4]ptr [si] [6]ptr cf[di], ax[di] [2], bx;add new coefficient to;accumulatorword ptr [di] [4], cxword ptr [di] [6], dx;perform a signed multiplysmul64, argument, cf, addr resultsi, word ptr result [4]di, word ptr cfcx, 4byte ptr neval;decrement pointerdi, word ptr outputsi, word ptr cfcx, 4;write to the outputTHE ELEMENTARY FUNCTIONSpopfretpolyeval endpCalculating Fixed-Point Square RootsFinding square roots can be an art.
This section presents two techniques. Thefirst, and perhaps the most traditional, is Newton’s Method. The other is thetechnique you learned in school adapted, to binary arithmetic. In this section, we’llexamine the square-root approximation in its simplest and most elemental form.Later, in the floating-point section, we’ll combine these with other techniques toimprove the first estimate and speed the overall convergence of the algorithm. Thereis no reason those techniques couldn’t also be made to fit a fixed-point application.Newton’s Method for finding square roots is a favorite among programmersbecause of its speed. It’s given by the equation r´= (x/r +r)/2, with x being ourradicand and r the estimate.
If you are interested, cube roots may be calculated r´=(r +(3 *x)/2)/( r *r +x/(2 *r))/2. It is an iterative approach that eventually finds theroot. There is no guarantee how many iterations it might take-that depends uponthe quality of the initial guess—but it should about double the number of correct bitson each iteration.Formulating that initial best guess is the problem. Resolving the routine canrequire an inordinate number of iterations if the first estimate is very far off. Thisroutine is simple; it only knows that it has a 32-bit input and that the greatest possibleroot of such an input is 16 bits. To improve first estimate, therefore, the routine shiftsthe radicand right until it fits within a 16-bit word. Still, there is no way of telling howmany iterations will be required. A loop counter with a large enough count wouldsuffice but could easily require more iterations than would otherwise be necessary.Instead, a copy of the last estimate is saved and compared with the current estimateafter each iteration.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
