Morgan - Numerical Methods (523161), страница 33
Текст из файла (страница 33)
We also need a way to determine the quadrant thefunction is in that fixes the sign (see Figure 6-3).Dcsin will produce either a sine or cosine depending upon a switch, cs_flag.Dcsin: Algorithm1. Clear sign, clear the output variable, and check the input argumentfor zero.Figure 6-2. Sine and cosine are the same function 90 degrees out of phase.224THE ELEMENTARY FUNCTIONSIf it is zero, set the output to 0 for sines and 1 for cosines.Otherwise, continue with step 2.2. Reduce the input argument by dividing by 360 (we are dealing indegrees) and take the remainder as our angle.If the result is negative, add 360 to make it positive.3. Save a copy of the angle in a register, divide the original again by90 to identify which quadrant it is in.
The quotient of this divisionremains in AX.4. Check cs-flag to see whether a sine or cosine is desired.A0h requests sine; continue with step 9.Anything else means a cosine; continue with step 5.5. Compare AX with zero.If greater, go to step 6.Otherwise, continue with step 13.6. Compare AX with one.Figure 6-3. Quadrants for sine/cosine approximations.225NUMERICAL METHODSIf greater, go to step 7.Otherwise, set sign.Two's complement the angle.Add 180 to make it positive again.Continue with step 14.7. Compare AX with two.If greater, go to step 8.Otherwise, set sign.Subtract 180 from the angle to bring it back within 90 degrees.Continue with step 13.8.
Two's complement the angle.Add 360 to point it back into the table.Continue with step 14.9. Compare AX with zero.If greater, go to step 10.Otherwise, 2's complement the angle.Add 90 to point it into the table.Continue with step 14.10. Compare AX with one.If greater, go to step ll.Otherwise, subtract 90 from the angle to bring it back onto thetable.Continue with step 13.11. Compare AX with two.If greater, go to step 12.Otherwise, two's complement the angle,Add 270, so that the angle points at table.Set sign.Continue with step 14.12. Set sign.Subtract 270 from the angle.13.
Use the angle to point into the table.226THE ELEMENTARY FUNCTIONSGet f(x0)from the table in the form of the nominal estimation of thesine.Check to see if any fraction bits require linear interpolation.If not, continue with step 15.Get f(x1) from the table in the form of the next greater approximation.Subtract f(x0) from f(x1) and multiply by the fraction bits.Add the result of this multiplication to f(x0).Continue with step 15.14. Use the angle to point into the table.Get f(x0) from the table in the form of the nominal estimation of thesine.Check to see if any fraction bits requiring linear interpolation.If not, continue with step 15.Get f(x1) from the table in the form of the next smaller approximation.Subtract f(x0) from f(x1) and multiply by the fraction bits.Add the result of this multiplication to f(x0).Continue with step 15.15.
Write the data to the output and check sign.If it's set, two's complement the output and exit.Otherwise, exit.Dcsin: Listing; *****.data;sines(degrees)sine_tblword0ffffh,0fe98h,0fa67h,0f378h,0egdeh,0ddb3h,0cf1bh,0be3eh,0fff6h,0fe17h,0f970h,0f20dh,0e803h,0db6fh,0cc73h,0bb39h,0ffd8h,0fd82h,0f865h,0f08fh,0e617h,0d919h,0cgbbh,0b826h,0ffa6h,0fcdgh,0f746h,0eeffh,0e419h,0d6b3h,0c6f3h,0b504h,0ff60h, 0ff06h,0fclch, 0fb4bh,0f615h, 0f4d0h,0ed5bh, 0eba6h,0e208h, 0dfe7h,0d43bh, 0d1b3h,0c4lbh, 0c134h,0bld5h, 0ae73h227NUMERICAL METHODSword0ab4ch,09679h,08000h,0681fh,04flbh,03539h,01ac2h,Oh0a7f3h,092d5h,07c1ch,06406h,04ad8h,030d8h,0164fh,0a48dh,08f27h,0782fh,05fe6h,04690h,02c74h,011dbh,0a11bh,08b6dh,07438h,05bbeh,04241h,0280ch,00d65h,09d9bh,087a8h,07039h,0578eh,03deeh,023aOh,008efh,09a10h,083d9h,06c30h,05358h,03996h,01f32h,00477h,.code;sines and cosines using a table and linear interpolation;(degrees)dscin proc uses bx cx si di, argument:word, cs_ptr:word, cs_flag:bytelocalpowers_of_two:byte, sign:bytepushfstdrepsubmovmovmovaddstoswaddmovmovaddmovrepejejmpzero_exit:cmpjnejmpcos_0:228;increment downax, axbyte ptr sign, alcx, 4di, wordptr cs_ptrdi, 6di, 8si, didi, word ptr argumentdi, 6cx, 4cmpsw;clear sign flag;clear sin/cos output;first check arguments;for zero;reset pointer;find the first nonzero,or;retumzero_exitprepare_argumentsbyte ptr cs_flag, alcos_0exit;ax is zero;sin(0) = 0THE ELEMENTARY FUNCTIONSincincaxaxadddecmovjmpsi,axaxword ptr [si][4],axexitprepare-arguments:movsi, word ptr argumentax, word ptr [si][4]movsubmovidivdx, dxcx, 360cxorjnsadddx, dxquadrantdx, 360quadrant:movbx, dxmovsubmovdivax, dxdx, dxcx, 90cx;point di at base of;output;make ax a one;cos(0)= 1;one;get integerportion;of angle;modulararithmeticto;reduceangle;we want the remainder;angle has to be;positive for this;to work;we will use this to;compute the value;of the function;put angle in ax;and this to compute;the sign ax holds;an index to the quadrantswitch:cmpbyte ptr cs_flag, 0jedo_sincos_range:cmpjgjmpax, 0cchk_l80walk_up;what do we want?;a zero=sine;anything else=cosine;use incrementing method229NUMERICAL METHODScchk_180:cmpjgnotnegaddjmpax, 1cchk_270byte ptr signbxbx, 180walk_backcchk_270:cmpjgnotsubjmpax, 2clast_90byte ptr signbx, 180walk_upclast_90:negaddjmp;;;do_sin:cmpjgnegaddjmp;set sign flag;use decrementing method;set sign flagbxbx, 360walk_back;find the range of the;argumentax, 0schk_180bxbx, 90walk_back;use decrementing methodax, 1schk_270bx, 90walk_up;use incrementing methodschk_180:cmpjgsubjmpschk_270:cmpjgnotnegaddjmp230ax, 2slast_90byte ptr signbxbx, 270walk_back;set sign flagTHE ELEMENTARY FUNCTIONSslast_90:notsubjmp;;walk_up:shladdmovmovorjeincincmovmovsubjncnegmulnotnegjcincjmppos_res0:mulleave_walk_up:addjmpwalk_back:shladdmovmovorbyte ptr signbx, 270walk_up;set sign flagbx, 1;use angle to point into;the tablebx, offset word ptr sine_tbldx, word ptr [bx]ax, word ptr [si][2]ax, axwrite_resultbxbxcx, dxax, word ptr [bx]ax, dxpos_res0axword ptr [si][2]dxaxleave_walk_updxleave_walk_up;get cos/sine of angle;get fraction bits;linear interpolation;get next approximation;find difference;multiply by fractionbitsword ptr [si] [2]dx, cx;multiply by fraction bits;and add in anglewrite-resultbx,bx,dx,ax,ax,1offset word ptr sine_tblword ptr [bx]word ptr [si][2]ax;point into table;get cos/sine of angle;get fraction bits231NUMERICAL METHODSjewrite_resultdecdecmovmovbxbxcx, dxax, word ptr [bx]subjncnegmulnotnegjcincjmppos_resl:mulleave_walk_back:addax, dxpas_reslaxword ptr [si][2]dxaxleave_walk_backdxleave_walk_backwrite_result:movmovmov232;multiply by fraction bitsdx, cx;multiply by fraction bits;and add in angledi, word ptr cs_ptrword ptr [di], axword ptr [di][2], dxax, axmovmovcmpjenotnotnotnegjcaddadcadcword ptr [di][4], axword ptr [di][6], axbyte ptr sign, alexitword ptr [di] [6]word ptr [di][4]word ptr [di][2]word ptr [di][0]exitword ptr [di][2],1word ptr [di][4],axword ptr [di][6],axpopfretdcsin endp;multiply by fraction bitsword ptr [si][2]subexit:;get next incremental;cos/sine;get difference;stuff result into variable;setup output for qword;fixed point;radix point between the;double wordsTHE ELEMENTARY FUNCTIONSComputing With TablesAlgebra is one of a series of fascinating lectures by the physicist RichardFeynman2.
In it, he discusses the development of algebra and its connection togeometry. He also develops the basis for a number of very interesting and usefulalgorithms for logarithms and powers, as well as the algebraic basis for sines andcosines using imaginary numbers.In algebra, Feynman describes a method of calculating logarithms, and thereforepowers, based on 10 successive square roots of 10. The square root of 10 (10.5) is3.16228, which is the same as saying log10(3.16228) = .5. Since log10(a*c) = log10(a)+ log10(c), we can approximate any power by adding the appropriate logarithms, ormultiplying by the powers they represent. For example, 10.875 = 10(.5+.25+.125) =3.16228 * 1.77828 * 1.33352 = 7.49894207613.As shown in Table 6-1, that taking successive roots of any number is the sameas taking that number through successively more negative powers of two.The following algorithm is based on these ideas and was suggested by DonaldKnuth3.
The purpose of pwrb is to raise a given base to a power x, 0 x < 1. Thisis accomplished in a manner similar to division. We do this by testing the inputargument against successively more negative powers of b, and subtracting those thatdo not drive the input argument negative.
Each power whose logarithm is less thanthe input is added to the output multiplied by that power. If a logarithm of a certainpower can not be subtracted, the power is increased and the algorithm continues. Theprocess continues until x = 0.233NUMERICAL METHODSnumberpower of 10power of 210.0123.16228l/22-l1.77828l/42-21.33352l/82-31.15478l/162-41.074607l/322-51.036633l/642-61.018152l/1282-71.0090350l/2562-81.0045073l/5121.00225111/10242-92-100Table 6-1. Computing with TablesIn pseudocode:Pwrb: Algorithm1. Set the output, y, equal to 1, clear the exponent counter, K.2.
Test our argument, x, to see if it is zero.If so, continue at step 6.If not, go to step 3.3. Use k to point into a table of logarithms of the chosen base tosuccessively more negative powers of two. Test x < logb(1+2-k).If so, continue with step 5.Else, go to step 4.4. Subtract the value pointed to in the table from x.Multiply a copy of the output by the current negative power of twothrough a series of right shifts.Add the result to the output.Go to step 2.234THE ELEMENTARY FUNCTIONS5.
Bump our exponent counter, k, by one,Go to step 2.6. There is nothing left to do, the output is our result,Exit.Pwrb: Listing; *****.datapower10qword4d104d42h, 2d145116h, 18cf1838h, 0d1854ebh,6bd7e4bh, 36bd211h, 1b9476ah,0dd7ea4h, 6ef67ah, 378915h, 1bc802h, 0de4dfh,6f2a7h, 37961h, 1bcb4h, 0de5bh,6f2eh, 3797h, 1bcbh, 0de6h, 6f3h, 379h, 1bdh,0deh, 6fh, 38h, 1ch, 0eh, 7h, 3h, 2h, 1h.code;; ******;pwrb - power to base 2;input argument must be be 1 <= x < 2uses bx cx dx di si, argument:qword, result:wordp w r b procreplocalk:byte, z:qwordmovsubmovstoswincstoswdecstoswmovdi, word ptr resultax, axcx, 2;Y;make y = 1axaxbyte ptr k, al;make k = 0x0:movmovmovsubax, word ptr argumentcx ptr argument [2]dx ptr argument;argument 0 <= x < 1bx, bx235NUMERICAL METHODScmpjnecmpjnecmpjneax, bxnot_done_yetcx, bxnot_doneyetdx, bxnot_doneyetjmppwrb_exitnot_done_yetsubmovcmpjabx, bxbl, byte ptr kbl, 20hpwrb_exitshlshlshlbx, 1bx, 1bx, 1leasi, word ptr power2cmpdx, word ptr [si] [bx] [4]jbjacmpjbjacmpjbincreasereducecx, word ptr [si][bx][2]increasereduceax, word ptr [si] [bx]increasesubsbbsbbmovmovmovax, word ptr [si] [bx]cx, word ptr [si] [bx] [2]dx, word ptr [si][bx][4]word ptr argument, axword ptr argument [2], cxword ptr argument[4], dxsubmovcx, cxcl, byte ptr kmovmovsi, word ptr resultax, word ptr [si];test for 0.0;our pointer and exponent;are we done?;point in to table of qwords;is this log greater than,;equal or less than;x?reduce:236;x<-x-zTHEmovmovcmpjeshiftk:shrrcrrcrloopno_shiftk:addadcadcjmpbx, word ptr [si][2]dx, word ptr [si][4]cl, 0no_shiftkELEMENTARYFUNCTIONS;is this shfit necessary?dx, 1bx, 1ax, 1shiftkword ptr [si], axword ptr [si] [2], bxword ptr [si][4], dx;z<-argument>>kx0increase:incjmppwrb_exit:retpwrb endpbyte ptr kx0;bump the counter to the;next level;and continueThere is another, similar, routine in the TRANS.ASM module dealing withlogarithms.CORDIC AlgorithmsCordic is an acronym meaning COordinate, Rotation Digital Computer4.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
