Morgan - Numerical Methods (523161), страница 34
Текст из файла (страница 34)
It wasdevised as a way to derive transcendental functions for real-time airborne navigationand has since been used in Intel math coprocessors and Hewlett-Packard calculators.The CORDIC functions are a group of algorithms capable of computing high—quality approximations of the transcendental functions and require very littlearithmetic power from the processor. Any functions listed in Table 6-2 can becalculated using only shifts, adds, and subtractions. These functions make very goodcandidates for the core of a floating-point library for processors with or withouthardware multiplication and division.237NUMERICAL METHODSinput:output:comments:x = x rectangular unitsl/k(xcos(z)-ysin(z))in general casey = y rectangular unitsl/k(ycos(z)+xsin(z))z = z angle0x=1cos(a) multiplierto computey=00the constantz=00circularkx = circulark(constant)cos(a)obtain siney=0sin(a)and cosine ofz=a0al/k(÷(x2+y2)0anglein general casex=xrectangular unitsy=yz=zrectangular unitsanglel/k(xcosh(z)+ysinh(z))in general casel/k(ycosh(z)+xsinh(z))0circular functionsinverse circular functionsx = x rectangular unitsy = y rectangular unitsz=zz+atan-1(y/x)hyperbolic functionsinverse hyperbolic functionsx=xrectangular unitsy=yz=zrectangular unitsangleTable 6-2.
CORDIC Functions238l/k(÷(x2+y2))in general case0z+tanh-1(y/x)THE ELEMENTARY FUNCTIONSThe CORDIC functions makeup a unified core that can derive many otherfunctions, including circular and hyperbolic, as well as powers and roots (see Table6-2). This discussion will focus on the circular functions; routines for the hyperbolicfunctions and inverses for both circular and hyperbolic are in the moduleTRANS.ASM.Before getting into the specifics of the routine, let’s take some time to understandhow the CORDIC functions work. Notice that this algorithm has some things incommon with the circle algorithm presented earlier in a Chapter 3. That routine useda modified rotation matrix:Ra[x,y] = [cos(a)-sin(a),sin(a)+cos(a)]and very small values for sine and cosine to draw a circle with only shifts, additions,and subtractions.
A similar idea is at work here, but it goes a step farther.See why the rotation matrix might help derive the functions listed above, lookat Figure 6-4. In a Cartesian coordinate system (x,y), you can specify a point on aplane by measuring its position relative to the (x,y) axis. In the figure, point P is atx=20, y=l0. If you draw a line from the origin of the axis to that point, it will forma vector of a certain length offset from the x axis by an angle a. To move this vectorabout the origin by some amount, in this case π/10 radians, you can use the rotationmatrix as shown in Figure 6-4. First solving for x = x*cos(α)-y*sin(α), then y =x*sin(α)+y*cos( α ), you will develop a new set of coordinates for point P.
In thisway, you can move around the origin simply by supplying an angle of rotation andthe current coordinates. With a few small changes, this same mechanism can deliverthe sine and cosine of a desired angle and a number of other functions as well.To make this work, x and y are needed plus a new argument, z, which willrepresent the angle or rotation. Next, simplify the equation by factoring out the cosineusing the fundamental identity tan(a) = sin(α )/cos(α ). This leavesR=a[x,y] =cos( α)[l-tan( α ), tan(a )+1]239NUMERICAL METHODSWriting this out long hand, you havex=cos( α)[x-x(tan( α ) ) ]andY= cos(α )[y(tan( α ))+y]One more step can ease the computing burden even more: replacing the twomultiplications, y(tan(α )) and x(tan(α )), with right shifts if a is made the sum of aseries of smaller a’s and each tan(α) is chosen to be a negative power of two.
If everytan(a) is to become a negative power of two, then the small piece of the angle eachrepresents becomes atan(2-i). This means that we will be breaking the input angle, a,Figure 6-4. Rotation Matrix240THE ELEMENTARY FUNCTIONSinto smaller angles equal to atan(2-i) and subtracting each atan(2-i) from the input aafter each evaluation of the rotation matrix in an effort to close on zero. Thissubtraction may involve positive and negative signs depending upon the quadrant weare in as we hover around zero; as the tangent changes sign, so then must the atan.See Figure 6-3 in the previous section for the progression of signs.
Now, the formulaebecomex = cos (α) [x-x (2-l)]andy = cos ( α) [y(2-l)+y]The cos(α ) remains, but it is a constant (circulark) that has been precomputedand is factored in when needed. Because we are using negative powers of two, eachiteration of the algorithm is responsible for a power of two; the result is 32 iterationsfor 32 fraction bits.For the routine, circular, the table of arctangents was precomputed and storedin the table atan_array, as was the constant cos(α ), circulark.
The same was donefor the hyperbolic functions with the table atanh_array and the constant hyperk. Tosolve for particular functions, see Table 6-2 for the correct inputs and the expectedoutputs.As with so many of the functions covered in this chapter, the input argument forthe angle must be confined to the first quadrant. Circular will solve for both sine andcosine, given X = circulark (the constant), Y = 0, and Z = a, if 0 a < π/2.
Reducingthe argument for these routines can be done in the same manner as in the table drivenroutine, dcsin in an earlier section. Divide the input angle by 2π, to remove unwantedcomponents of π, then divide by π/2. Take the remainder as your input argument andthe quotient as an index to the quadrant the angle is in. See Figure 6-3 for the logic.241NUMERICAL METHODSCircular: Algorithm1. The variables X, Y, and Z serve as both input and output variables.Load x, y, and z into local variables smal lx, smal ly, and smal lz.Set the exponent counter, i, to 0.2.
Multiply x and y by 2-i and store in smal lx and smal ly,respectively. (The multiplication is accomplished by shifting (arithmetically) x and y to the right by the current count in i.)Load z with table entry pointed at by atan_array+i.3. Test zIf true, go to step 5.Else, continue with step 4.4. Add smal ly to x.Subtract smal lx from y.Add smal lz to z.Continue with step 6.5. Subtract smal ly from x.Add smal lx to y.Subtract smal lz from z.Continuewithstep6.6.
Bump the exponent counter, i.Test iIf yes, got to step 7.Otherwise, go to step 2.7. Since we have been using the output variables for intermediatestorage of our results, the output is current and we may exit.Circular: Listing; ******.dataatan_array.code242dword0c90fdaa2h, 76b19c16h, 3eb6ebf2h, 1fd5ba9bh,0ffaaddch, 7ff556fh, 3ffeaabh, 1fffd55h,0ffffabh, 800000h, 3fffffh, 200000h, 100000h,80000h, 40000h, 20000h, l0000h, 8000h, 4000h,2000h, 1000h, 800h, 400h, 200h, 100h, 80h,40h, 20h, 10h, 8h, 4h, 2h, 1hTHE ELEMENTARY FUNCTIONS;; circular-implementation of the circular routine, a subset of the CORDIC devicescircularproc uses bx cx dx di si, x:word, y:word, z:wordlocalsmallx:qword, smally:qword, smallz:qword, i:byte,shifter:worddi, word ptr smallxsi, word ptr xcx, 4repleamovmovmovswdi, word ptr smallysi, word ptr ycx, 4repleamovmovmovswdi, word ptr smallzsi, word ptr zcx, 4repleamovmovmovswsubmovmovmovtwist:submovmovax, axbyte ptr i, albx, axcx, ax;i=Oax, axal, iword ptr shifter, axmovmovmovmovmovsi,ax,bx,cx,dx,cmpword ptr shifter, 0load_smallxje;load input x, y, and zwordwordwordwordwordptrptrptrptrptrx[si][si][2][si][4][si][6];multiply by 2ˆ-i243NUMERICAL METHODSshiftx:sarrcrrcrrcrdecjnzload_smallx:movmovmovmovsubmovmovmovmovmovmovmovcmpjeshifty:sarrcrrcrrcrdecjnzload_smally:movmovmovmovget_atan:submovshlshl244;note the arithmetic shift;for sign extension;negative powers of two;require right shiftsdx, 1cx, 1bx, 1ax, 1word ptr shiftershiftxwordwordwordwordptrptrptrptrsmallx, axsmallx [2], bxsmallx [4], cxsmallx [6], dxax, axal, iword ptr shifter, axsi, word ptryax, word ptr [si]bx, word ptr [si][2]cx, word ptr [si][4]dx, word ptr [si][6];return x to smallx;get yword ptr shifter, 0load_smally;multiply by 2^-idx, 1cx, 1bx, 1ax, 1word ptr shiftershifty;note the arithmetic shift;for sign extension;take to a negative powerwordwordwordword;return to smallybx,bl,bx,bx,ptrptrptrptrbxi11smally, axsmally[2], bxsmally[4], cxsmally[6], dx;have to point into a dword;tableTHE ELEMENTARY FUNCTIONSleamovsi, word ptr atan_arrayax, word ptr [si] [bx]movdx, word ptr [si] [bx][2]movmovsubmovmovword ptrword ptrax, axword ptrword ptrmovmovsi, word ptr zax, word ptr [si][6]orjnsax, axpositivesmallz, axsmallz [2], d x;use the negative power;of two as a pointer;to get proper atan;z=atan[i]smallz [4], axsmallz [6], axtest_Z:negative:movmovmovmovax,bx,cx,dx,wordwordwordwordptrptrptrptrsmallysmally[2]smally[4]smally[6]movaddadcadcadcdi, wordword ptrword ptrword ptrword ptrptr x[di], ax[di][2], bx[di][41, cx[di][6], dxmovmovmovmovax,bx,cx,dx,ptrptrptrptrmovsubsbbsbbsbbdi, wordword ptrword ptrword ptrword ptrmovmovax, word ptr small Zbx, word ptr small z[2]wordwordwordword;the sign of z determines;whether the arguments;are summed or subtracted;smally is added xsmallxsmallx[2]smallx[4]smallx[6]ptr y[di], ax[di][2], bx[di][4], cx[di][6], dx;smallx is added toy245NUMERICAL METHODSmovmovcx, word ptr smallz [4]dx, word ptr smallz [6]movaddadcadcadcdi, wordword ptrword ptrword ptrword ptrjmpfor_nextpositive:movmovmovmovmovsub246ax, wordbx, wordcx, worddx, worddi, wordword ptrptr z[di], ax[di][2], bx[di][4], cx[di][6], dxptr smallyptr smally[2]ptr smally[4]ptr smally[6]ptr x[di], axsbbsbbsbbword ptr [di][2], bxword ptr [di][4], cxword ptr [di][6], dxmovmovmovmovmovaddadcadcadcax, word ptr smallxbx, word ptr smallx[2]cx, word ptr smallx[4]dx, word ptr smallx[6]di, word ptryword ptr [di], axword ptr [di] [2], bxword ptr [di][4], cxword ptr [di][6], dxmovmovmovmovmovax,bx,cx,dx,di,subsbbsbbword ptr [di], axword ptr [di][2], bxword ptr [di][4], cxwordwordwordwordwordptrptrptrptrptrsmallzsmallz[2]smallz[4]smallz[6]z;and smallz is added to z;z was positive, so;smally is subtracted;from y;smallx is added to y;and smallz is subtracted;fromzTHE ELEMENTARY FUNCTIONSsbbfor_next:inccmpjajmpword ptr [di][6], dxbyte ptr ibyte ptr i, 32circular_exittwist;bump exponent on each passcircular_exit:retcircularendpPolynomial EvaluationsOne of the most popular and most accurate ways to develop the transcendentalsis evaluation of a power series.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
