Morgan - Numerical Methods (523161), страница 18
Текст из файла (страница 18)
Be sure to push ax, cx, dx,; and bx before destroying them; they are used here. The value for the x; coordinate is in bx, and the value for the y coordinate is in dx.orjnsaddaddaddjmpdpositive:addaddaddchk_loop:loopretlineendpax, axdpositivebx, word ptr xstepdx, word ptr ystepax, incrshort chk loop;calculate new position and;update the decision variablebx, word ptr xstep_diagdx, word ptr ystep_diagax, word ptr diag incrline loopWhen fixed-point operands become very large, it sometimes becomes necessaryto find alternate ways to perfom arithmetic.
Multiplication isn’t such a problem; ifit exists as a hardware instruction on the processor you are using, it is usually fasterthan division and is easily extended.Division is not so straightforward. When the divisors grow larger than the sizeallowed by any hardware instructions available, the programmer must resort to other104REAL NUMBERSmethods of division, such as CDIV (described in Chapter 2), linear polation (used inthe floating-point routines), or one of the two routines presented in the followingpages involving Newton-Raphson approximation and iterative multiplication.
Thefirst two will produce a quotient and a remainder, the last two return with a fixed pointreal number. Choose the one that best fits the application.3, 4Division by InversionA root of an equation exists whenever f(x)=0. Rewriting an equation so thatf(x)=0 makes it possible to find the value of an unknown by a process known asNewton-Raphson Iteration. This isn’t the only method for finding roots of equationsand isn’t perfect, but, given a reasonably close estimate and a well behaved function,the results are predictable and correct for a prescribed error.Look at Figure 3-3. The concept is simple enough: Given an initial estimate ofa point on the x axis where a function crosses, you can arrive at a better estimate byevaluating the function at that point, f(x,), and its first derivative of f(x0) for the slopeof the curve at that point.
Following a line tangent to f(x0) to the x axis produces animproved estimate. This process can be iterated until the estimate is within theallowed error.The slope of a line is determined by dividing the change in the y axis by thecorresponding change in the x axis: dy/dx.
In the figure, dy is given by f(x0), thedistance from the x axis at x0 to the curve, and dx by (x1- x0), which results inf´(x0) = f(x0)/(xl-x0)Solving for x, givesx1= x0 - f (x0)/f’(x0)Substituting the new x, for x0 each time the equation is solved will cause theestimate to close in on the root very quickly, doubling the number of correctsignificant bits with each pass.To using this method to find the inverse of a number, divisor, requires that the105NUMERICAL METHODSFigure 3-3.Newton-Raphsoniteration.equation be formulated as a root. A simple equation for such a purpose isf(x) = 1/x - divisorFrom there, it’s a short step tox = 1/ divisorNow the equation for finding the root becomes an equation for finding theinverse of a number:x1,=((1/x)- divisor)/(-1/ divisor 2 )106REAL NUMBERSwhich simplifies to:xnew = xold * (2 -divisor(x)old)In his book Digital Computer Arithmetic, Joseph Cavanagh suggests that thisequation be simplified even further to eliminate the divisor from the iterative processand use two equations instead of one.
He makes the term divisor(x) equal to onecalled unity (because it will close in on one) in this routine, which reduces theequation to:xnew= Xold * (2 - unity)Multiplying both sides of this equation by the divisor, divisor, and substitutingagain for his equality, divisor(x) = unity, he generates another equation to producenew values for unity without referring to the original divisor:( divisor (x))new = ( divisor (x))old * (2 - unity )=unity new= unity old* (2 - unity)This breaks the process down to just two multiplies and a two’s complement.When these two equations are used together in each iteration, the algorithm willcompute the inverse to an input argument very quickly.To begin, there must be an initial estimate of the reciprocal.
For speed, this canbe computed with a hardware instruction or derived from a table if no suchinstruction exists on your processor. Multiply the initial estimate by the divisor to getthe first unity. Then, the two equations are evaluated as a unit, first generating a newdivisor and then a new unity, until the desired precision is reached.The next routine expects the input divisor to be a value with the radix pointbetween the doublewords of a quadword, fixed-point number.
The routine finds themost significant word of the divisor, then computes and saves the number of shiftsrequired to normalize the divisor at that point-that is, position the divisor so that itsmost significant one is the bit just to the right of the implied radix point: .1XXX...107NUMERICAL METHODSFor example, the number 5 is101.000Bradix pointNormalized, it is:000.101Bradix pointAfter that, the actual divisor is normalized within the divisor variable as if theradix point were between the third and fourth words.
Since the greatest numberproportion or divisor will see is two or less, there is no danger of losing significantbits. Placing the radix point there also allows for greater precision.Instead of subtracting the new proportion from two, as in the equation, we two’scomplementproportion and the most significant word is ANDed with 1H to simulatea subtraction from two. This removes the sign bits generated by the two’s complement and leaves an integer value of one plus the fraction bits.Finally, the reciprocal is realigned based on a radix point between the doublewordsas the fixed-point format dictates, and multiplied by the dividend.divnewt: Algorithm1.Set the loop counter,lp, for three passes.
This is a reasonable numbersince the first estimate is 16-bits. Check the dividend and the divisorfor zero.If no such error conditions exist, continue with step 2,Otherwise, go to step 10.2. Find the most significant word of the divisor.Determine whether it is above or below the fixed-point radix point.In this case, the radix point is between the doublewords.Test to see if it is already normalized.If so, go to step 5.3.Shift a copy of the most significant word of the divisor left or rightuntil it is normalized, counting the shifts as you proceed.4.Shift the actual divisor until its most significant one is the MSB of108REAL NUMBERSthe third word of the divisor.
This is to provide maximum precision forthe operation.5. Divide 10000000H by the most significant word of the divisor for a firstapproximation of the reciprocal. The greater the precision of this firstestimate, the fewer passes will be required in the algorithm (the resultof this division will be between one and two.)Shift the result of this division left one hex digit, four bits, to alignit as an integer with a three-word fraction part. This initial estimateis stored in the divisor variable.Divide this first estimate by two to obtain the proportionality variable,proportion.6.Perform a two's complement on proportion to simulate subtracting it fromtwo. Multiply proportion by divisor. Leave the result in divisor.7.Multiply proportion by the estimate, storing the results in bothproportion and estimate.
Decrement lp.If it isn't zero, continue with step 6.Otherwise, go to step 8.8.Using the shift counter, shift, reposition divisor for the finalmultiplication.9. Multiply divisor, now the reciprocal of the input argument, by thedividend to obtain the quotient. Write the proper number of words to theouput and exit.10. An attempt was made to divide zero or divide by zero; exit with error.divnewt: Listing; *****;divnewt- division by Raphson-Newton zeros approximation;;;divnewtproc uses bx cx dx di si, dividend:qword, divisor:qword,quotient:wordlocaltemp[8]:word, proportion:qword, shift:byte, qtnt_adjust:byte,lp:byte, tmp:qword, unity:qword;upwardcldsubcx, cx109NUMERICAL METHODSbyte ptr lp, 3qtnt_adjust, clcx, word ptr dividend[O]cx, word ptr dividend[2]cx, word ptr dividend[4]cx, word ptr dividend[6]ovrflw;zero dividendsubororororjecx, cxcxr wordcx, wordcx, wordcx, wordovrflw;zero divisorsubmovfind_msb:leadecdeccmpjemovmovsubcmpjbjatestjneshift_left:decshltestjnejmpnormalizeshift_right:incshror110;should only take three passesmovmovororororjeptrptrptrptrdivisor [O]divisor [2]divisor [4]divisor [6]ax, axbx, 8;look for MSB of divisordi, word ptr divisorbxbxword ptr [di] [bx], axfind_msbbyte ptr gtnt_adjust, blax, word ptr [di][bx]cx, cxbx, 2hshift_leftshift_rightword ptr [di][bx], 8000hnorm dvsrcxax, 1ah, 80hsave_shiftshift_ leftcxax, 1ax, ax;di is pointing at divisor;get most significant word;save shifts here;see if already normalized;normalized?;it's already there;count the number of shifts toREAL NUMBERSjejmpsave shift:movsubshift back:cmpjeshrrcrrcrrcrjmpsave-shiftshift right;count the number of shifts to;normalizebyte ptr shift, clax, axword ptr [di][6], axnorm_dvsrwordptr [di][6],word ptr [di][4],word ptr [di][2],word ptr [di] [0],shift back1111norm dvsr:testjneshlrclrcljmpword ptr [di][4],make_firstword ptr [di][O],word ptr [di][2],word ptr [di][4],norm_dvsrmake first:movsubmovdx, l000hax, axbx, word ptr [di][4]divsubmovcorrect dvsr:shlrclloopmovmovsubmovmovshr;we will put radix point;at word three8000h111;the divisor;truly normalized;for maximum;this should normalize divisor;first approximation;could come from a tablebxdx, dxcx, 4;keep only the four least bitsax, 1dx, 1correct_dvsrword ptr divisor[4],word ptr divisor[6],cx, cxword ptr divisor[2],word ptr divisor[0],dx 1;don't want to waste time with;a big shiftaxdxcxcx;don't want to waste time111NUMERICAL METHODS;with a big shiftrcrmu1shlax, 1bxax, 1rclmovsubmovmovmovdx, 1word ptrcx, cxword ptrword ptrword ptr;reconstruct for first attempt;don't want to waste time;with a big shiftunity[4], dxunity[6], cxunity[2], cxunity, cx;this could be done with;a tablemakeproportion:movsubmovmovmovinvert_proportion:notnotnotnegword ptrax, axword ptrword ptrword ptrproportion[4], dxwordwordwordwordproportion[6]proportion[4]proportion[2]proportionptrptrptrptrproportion[6], axproportion[2], axproportion, axjcaddadcadcmloopword ptr proportion[2], 1word ptr proportion[4], 0word ptr proportion[6], 0andword ptr proportion[6], 1;attempt to develop with;2's complementmloop:invokerep112;make it look like it was;subtracted from 2mu164, proportion, divisor, addr templealeamovmovswsi, word ptr temp[6]di, word ptr divisorcx, 4invokemu164, proportion, unity, addr templeasi, word ptr temp[6]REAL NUMBERSreprepleamovmovswlealeamovmovswdi, word ptr unitycx, 4decjejmpbyte ptr lpdivnewt_shiftinvert_proportionsubnotmovmovstoswjmpax, axaxcx, 4di, word ptr quotientsi, word ptr temp[6]di, word ptr proportioncx, 4;six passes for 64 bitsovrflw:repdivnewt_shift:leamovorjsqtnt_right:movsubmovsubjmpqtnt_left:negsubaddqtlft:shlrclrclrclloopdivnewt_mult:;make infinite answerdivnewt_exitdi, word ptr divisorcl, byte ptr shiftcl, clqtnt_left;get shift count;positive, shift leftch, 10hch, clcl, chch, chqtlftClch, chcl, 10hword ptrword ptrword ptrword ptrqtlft;we Want to take it to the MSB[dil[O],[di][2],[di1[4],[di][6],1111;multiply reciprocal by dividend113NUMERICAL METHODSreprep;see that temp is clearsubmovleastoswax, axcx, 8di, word ptr tempinvokemovaddmovleaaddcmpjaemovmovswjmpmul64, dividend, divisor, addr tempbx, 4;adjust for magnitude of resultbl, byte ptr qtnt_adjustdi, word ptr quotientsi, word ptr tempsi, bxbl, 0ahwrite zerocx, 4write_zero:movrepmovswsubstoswdivnewt_exit:popfretdivnewt endpdivnewt exitcx, 3ax, axDivision by MultiplicationIf the denominator and numerator of a fraction are multiplied by the same factor,the ratio does not change.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
