Morgan - Numerical Methods (523161), страница 38
Текст из файла (страница 38)
Two such limitations are Ymax and Eps. Ymax is the maximum allowableargument for the function that will produce accurate results with minimum error, andEps is the smallest allowable argument. The values for these are chosen based on thesize of the data types and the functions being approximated. They’re important in thecalculation of a number of elementary functions, notably flsin (discussed later).Square RootsThe first function presented here, flsqr, computes the square roots of floatingpoint numbers.
Simply, this function finds the square root of the mantissa portion ofthe float and then the root of 2exponent. It then reconstructs the float and returns.To begin with, the function frxp is called to constrain the radicand to a small,relatively linear region, .5 x < 1 (this represents an exponent of -1).
Within this269NUMERICAL METHODSregion, all square roots adhere to the relationship, nput_raidcand < root < 1,precisely, all roots must exist from about .7071067 to 1.0. This makes it much easierto come up with an initial estimate that is very close. Just taking the mid-range valuefor the first estimate would improve it considerably. Recall that Newton’s Methoddelivers about about twice the number of accurate bits for each iteration; that is, ifthe initial estimate is accurate to x bits, after the first iteration, will have about 2*x + 1accurate bits.
But even this can require an unknown number of iterations to converge,so the estimate must be improved.The most popular solution is the formula for a straight line, y = m*x+b.Calculating the values for m and b that provide the best fit to the square-root curveyields slightly different values depending on the approach you take. Cody and Waiteuse the values .59016 for m and .41731 for b, which will always produce an initialestimate that’s less than one percent in error. Solving for y in the equation for astraight line yields the first estimate, and only two passes through Newton’s Methodproduces a result for a 24-bit mantissa.Finding the root of 2exponent is simple if the exponent is even: divide by 2, just aswith logarithms. If the exponent is odd, however, it cannot be divided evenly by two,so it must first be incremented by one.
To compensate for this adjustment, we dividethe root of the input mantissa by sqrt(2). In other words, the exponent represents log,of the input number; to find its root, simply divide by two. If the exponent must beincremented before the division, the root of that additional power must be removedfrom the mantissa to keep the result correct.It’s then a simple matter of reassembling the float using the new mantissa andexponent.Flsqr: Algorithm1. Test input to see whether it is greater than zero.If it's equal to zero, or less, exit with error through step 7.2. Use frxp to get the exponent from the input float, and to set itsexponent to zero, constraining it to .51.
Multiply thisnumber, f, by .59016 and add .41731 for our first approximation.3. Make two passes through r=(x/r+r)/2.270THE ELEMENTARY FUNCTIONS4. Inspect the exponent, n, derived earlier with frxp.If it's odd, multiply our best estimate from Heron's formula by thesquare root of .5 and increment n by 1.Even or odd, divide n by two.5. Add n back into the exponent of the float.6. Write the root to the output.7.
Leave.Flsqr: Listing; ******; flsqrflsqr proclocaluses bx cx dx si di, fp0:qword, fp1:wordresult:qword, temp0:qword, temp1:qword, exp:byte,xn:qword, f:qword, y0:qword,m:bytepushfcldrepleasubmovstoswdi, word ptr xnax, axcx, 4invokeflcomp, fp0, zerocmpjecmpjemovsubmovstoswrepnotandmovjmpgot-result:;error, entry value too;largeax, 1okax, 0got-resultdi, word ptr fplax, axcx, 4axax, 7f80hword ptr result[4], axflsqr_exit;make it plus infinity271NUMERICAL METHODSrepmovsubmovstoswjmpdi, word ptr fplax, axcx, 4invokefrxp, fp0, addr f, addr expinvokeinvokeflmul, f, y0b, addr temp0fladd, temp0, y0a, addr y0invokeinvokemovshlsubshrmovfldiv, f, y0, addr temp0fladd, y0, temp0, addr temp0ax, word ptr temp0[4]ax, 1ah, 1ax, 1word ptr temp0[4], axinvokeinvokemovshlsubshrmovmovfldiv, f, temp0, addr temp1fladd, temp0, temp1, addr temp0ax, word ptr temp0[4]ax, 1ah, 1ax, 1word ptr y0[4], axax, word ptr temp0[2]movmovmovsubmovword ptrax, wordword ptrax, axword ptrmovmovsarjncal, byte ptr expcl, alal, 1evnflsqr_exitok:heron:;get exponent;two passes through;(x/r+r)/2 is all we need;should always be safe;subtracts one half;by decrementing the;exponent one;should always be safe;subtracts one half;by decrementing the;exponent oney0[2], axptr temp0y0, axy0[6], axchk_n:272;arithmetic shift, pleaseTHE ELEMENTARY FUNCTIONSodd:invokeflmul, y0, sqrt_half, addr y0movincsaral, clalal, 1;bump exponent on odd;divide by twomovcl, al;n/2->mmovshladdax, word ptr y0[4]ax, 1ah, cl;adjustment for uneven;exponentevn:power:write_result:shrmovleamovmovmovswrepax, 1word ptr y0[4], axsi, word ptr y0di, word ptr fplcx, 4flsqr_exit:popfretflsqr endpSines and CosinesThe final routine implements the sine function using a minimax polynomialapproximation.
A minimax approximation seeks to minimize the maximum errorinstead of the average square of the error, which can allow isolated error spikes. Theminimax method keeps the extreme errors low but can result in a higher averagesquare error. Ultimately, what this means is that the function is resolved using apower series whose coefficients have been specially derived to keep the maximumerror to a minimum value.This routine defines the input argument as some integer times π plus a fractionequal to or less than π/2. It expects to reduce the argument to the fraction f, byremoving any multiplies of π It then approximates the sine (f) based on the273NUMERICAL METHODSevaluation of a small interval symmetric about the origin, f, and puts the numberback together as our result.
It solves for the cosine by adding π /2 to the argument andproceeds as with the sine (see Figure 6-3).In this function we again encounter Ymax and Eps. These limitations depend onthe precision available to the arithmetic in the particular machine and help guaranteethe accuracy of your results. According to Cody and Waite, Ymax should be nogreater than π*2(t/2) and Eps, no less than 2(-t/2), where t is the number of bits availableto present the significand9.
In this example, t is 11 bits, but that doesn’t take theextended precision into account.The algorithm is a fairly straightforward implementation. If the input argumentis in range, this function initially reduced it to xn initially by multiplying by l/ π(floating-point multiplication is generally faster than division) and calling intrnd toget the closest integer. Multiplying xn by π and subtracting the result from theabsolute value of the input argument extracts a fraction, f, which is the actual angleto be evaluated with Cody and Waite’s minimax approximation.The polynomial R(g) is evaluated using a table of precomputed coefficients andHorner’s rule, except that in this implementation, the usual loop (see Polyeval in thelast section) was unrolled.R(g) = (((r4*g+r3)*g+r2)*g+rl)*gwhere g = f*f. The values r4 through r1 are coefficients stored in the table sincos.After the evaluation, R(g) is multiplied by f and f is added to it.
The only thingleft to do is adjust the sign of the result according to the quadrant. The pseudocodefor this implementation of flsin is as follows.flsin: Algorithm1. See that the input argument is no greater than Ymax.If it is, exit with error through step 8.2. Take the absolute value of the input argument.Multiply by 1/π to remove multiple components of πUse intrnd to round to the closest integer, xn.274THE ELEMENTARY FUNCTIONSTest xn to see whether it's odd or even. If it's odd, there is a signreversal; complement sign.3.
Reduce the argument to f through (|x|-xn*c1)-xn*c2 (in other words,subtract the rounded value xn multiplied by p from the input argument).4. Compare f with Eps.If it is less, we have our result, exit through step 8.5. Square f, (f*f->g) and evaluate r(g)6. Multiply f by R(g), then add f.7. Correct for sign; if sign is set, negate result.8. Write the result to the output and leave.Flsin: Listing; *****.datasincos qword404900000000h, 3a7daa20968bh, 0be2aaaa8fdbeh, 3c088739cb85h,0b94fb2227f1ah, 362e9c5a91d8h.code;; ******; flsinflsin proclocaluses bx cx dx si di, fp0:qword, fp1:word, sign:byteresult:qword, temp0:qword, temp1:qword,y:qword, u:qwordpushfcldinvokeflcomp, fp0, ymaxcmpjlax, 1absx;error, entry value too;largeerror-exit:275NUMERICAL METHODSrepleasubmovstoswjmpdi, word ptr resultax, axcx, 4writeoutabsx :movorjnsandmovax, word ptr fp0 [4]ax, axdeconstruct_exponentax, 7fffhword ptr fp0[4], ax;make absolutedeconstruct_exponent:flmul, fp0, one_over_pi, addr resultinvoke;(x/pi)invokeintrnd, result, addr temp0;intrnd(x/pi)movax, word ptr temp0[2]movmovanddx, word ptr temp0[4]cx, dxcx, 7f80hshlmovsubsubexponentjsincorjeextract_int:shlrclrcllooptestjenotnot_odd:276cx, 1cl, chch, chcl, 7fh;determine if integer;has odd or even;number of bits;get rid of sign and;mantissa portion;subtract bias (-1) fromnot-oddclcl, clnot-oddax, 1dx, 1word ptr bx, 1extract_intdh, 1not_oddbyte ptr sign;position as fixed pointTHE ELEMENTARY FUNCTIONSxpi:invokeinvoke;extended-precision;multiply by piflmul, sincos[8*0], temp0, addr result;intrnd(x/pi)*c1flsub, fp0, result, addr result;|x|-intrnd(x/pi)invokeflmul, temp0, sincos[8*1], addr temp1;intrnd(x/pi)*c2invokeflsub, result, temp1, addr y;Ychk_eps:invokerepflabs, y, addr temp0;is the argument less;than eps?invokeorjnsleasubmovstoswjmpflcomp, temp0, epsax, axr_gdi, word ptr resultax, axcx, 4invokeflmul, y, y, addr uinvokeinvoke;evaluate r(g);((r4*g+r3)*g+r2)*g+rl)*gflmul, u, sincos[8*5], addr resultfladd, sincos[8*4], result, addr resultinvokeinvokeflmul, u, result, addr resultfladd, sincos[8*3], result, addr resultinvokeinvokeflmul, u, result, addr resultfladd, sincos[8*2], result, addr resultinvokeflmul, u, result, addr resultwriteoutr_g277NUMERICAL METHODS;result== zfxr:invokeflmul, result, y, addr resultinvokefladd, result, y, addr result:r*r+fhandle_sign:cmpjnexorbyte ptr sign, -1writeoutword ptr result[4], 8000h;result * signwriteout:movleamovmovswrepdi, word ptr fplsi, word ptr resultcx, 4flsin_exit:popfretflsin endpDeriving the elementary functions is both a science and an art.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
