Xref: unixg.ubc.ca alt.fractals:2839 bit.listserv.frac-l:1252 Newsgroups: alt.fractals,bit.listserv.frac-l Path: unixg.ubc.ca!ubc-cs!destroyer!caen!uakari.primate.wisc.edu!usenet.coe.montana.edu!news.u.washington.edu!milton.u.washington.edu!jesjones From: jesjones@milton.u.washington.edu (Jesse Jones) Subject: Fractal Paper #1: Cordic Math Message-ID: <1992Sep17.224950.11642@u.washington.edu> Summary: Describes a fast method of computing elementary functions. Keywords: Fractal Cordic Integer Math Elementary Functions Sender: news@u.washington.edu (USENET News System) Organization: University of Washington, Seattle Date: Thu, 17 Sep 1992 22:49:50 GMT Lines: 345 I've had around fifteen requests for my papers on implementing a fractal program so I'm just going to repost them. I'll start off with two new papers describing cordic math... Cordic Math 1.0d1 by Jesse Jones INTRODUCTION This file discusses a fast algorithm for calculating the elementary functions using cordic math. The advantage of this method over more traditional approaches like power series [1] include: 1) The same algorithm can be used to calculate multiplication, division, sin, cos, tan, arctan, sinh, cosh, tanh, arctanh, ln, exp, and square-root. 2) The algorithm uses only shifts and adds so it will usually be faster then an equivalent power series. 3) Sin/Cos and Sinh/Cosh are calculated simultaneously (this is especially helpful when dealing with complex trig functions). DESCRIPTION The algorithm is based on coordinate rotation, i.e. an initial vector P0 is rotated until some condition is met. The vector P0 is determined by the coordinate system, m, and two initial points, x0 and y0. The coodinate systems are: circular (m = 1), linear (m = 0), and hyperbolic (m = -1). P0 is then given by: Length = Sqrt(x*x + y*y) Eq. 1 Angle = ArcTan(y/x) m = 0 or 1 ArcTanh(y/x) m = -1 The rotated vector, P1, is determined by the points x1 and y1: x1 = x0 + m*y0*delta0 Eq. 2 y1 = y0 - x0*delta0 where delta is an arbitrary power of two. This has the effect of both rotating and scaling the original vector, P0: Length1 = Length0*Scale0 Eq. 3 Angle1 = Angle0 - Alpha0 when Scale and Alpha are defined thusly: Scale = Sqrt(1 + m*delta*delta) Eq. 4 ArcTan(delta) m = 1 Alpha = delta m = 0 ArcTanh(delta) m = -1 The total change in Angle is given by: z1 = z0 + Alpha0 Eq. 5 The system of equations to be iterated is then: x1 = x0 + m*y0*delta0 Eq. 6 y1 = y0 - x0*delta0 z1 = z0 + Alpha0 After n+1 iterations we wind up with these values for x, y, and z: xn = K*[x0*Cos(A*Sqrt(m)) + y0*Sin(A*Sqrt(m))*Sqrt(m)] Eq.7 yn = K*[y0*Cos(A*Sqrt(m)) - x0*Sin(A*Sqrt(m))/Sqrt(m)] zn = z0 + A where: A = Alpha0 + Alpha1 + Alpha2 + ... + AlphaN Eq. 8 K = Scale0*Scale1*Scale2*...*ScaleN (note that Cos(i*z) = Cosh(z) and Sin(i*z) = i*Sinh(z)). In the actual implementation of the algorithm the magnitude of each delta is a fixed power of two. The sign is chosen to force either Angle or z to zero. Alpha can be stored in a small look-up table. CONVERGENCE For the system in Eq. 6 to converge to the desired result the delta values must be chosen appropiately. To speed the algorithm delta should be a power of two so that a shift can be used instead of a multiplication, e.g. delta = 2^(-k) where k is given below: m k Table 1 -1 1, 2, 3, 4, 4, 5, ... 0 1, 2, 3, 4, 5, 6, ... 1 0, 1, 2, 3, 4, 5, ... Note that for m = -1 the following shifts are repeated: {4, 13, 40, 121, ...., i, 3i+1}. If these shifts are not repeated sinh and cosh will have (at best) only four significant bits for small arguments. The number of significant bits will be roughly equal to the number of iterations. The sole exception is for functions returning large results. For example if Eq. 6 is iterated 30 times the cordic algorithm returns only 19 significant bits for Tan(1.545) = 38.75. For smaller results the algorithm returns 28 or more significant bits. The system in Eq. 6 will converge to different values depending on x0, y0, z0, and m. The values are as follows: Circular (m = 1) Table 2 z -> 0 x' = K*(x*Cos(z) - y*Sin(z)) y' = K*(y*Cos(z) + x*Sin(z)) z' = 0 A -> 0 x' = K*Sqrt(x*x + y*y) y' = 0 z' = z + ArcTan(y/x) Linear (m = 0) z -> 0 x' = x y' = y + x*z z' = 0 A -> 0 x' = x y' = 0 z' = z + y/x Hyperbolic (m = -1) z -> 0 x' = K*(x*Cosh(z) + y*Sinh(z)) y' = K*(y*Cosh(z) + x*Sinh(z)) z' = 0 A -> 0 x' = K*Sqrt(x*x - y*y) y' = 0 z' = z + ArcTanh(y/x) where K is as in Eq.8 and x' = xn, x = x0, etc. So, to calculate sin/cos we simply let: x0 = 1/K (which will be constant for a given number of iterations), y0 = 0, and z = the angle to be evaluated. Xn will be Cos(z) and yn will be Sin(z). The identities below can be used to generate additional functions: Tan(z) = Sin(z)/Cos(z) Table 3 Tanh(z) = Sinh(z)/Cosh(z) Exp(z) = Sinh(z) + Cosh(z) Ln(w) = 2*ArcTanh(y/x) where x = w+1 and y = w-1 Sqrt(w) = Sqrt(x*x - y*y) where x = w+1/4 and y = w-1/4 The domain of all the functions but ArcTan is limited to small values so the input arguments must be reduced. For example sin and cos converge only for ABS(z) <= 1.74. To reduce z we can divide by PI/2 to obtain an integer quotient, Q, and a remainder, D. The remainder is then fed into the algorithm. The reductions work as follows: Sin(D) if Q mod 4 = 0 Table 4 Sin(Q*PI/2 + D) = Cos(D) if Q mod 4 = 1 -Sin(D) if Q mod 4 = 2 -Cos(D) if Q mod 4 = 3 Cos(D) if Q mod 4 = 0 Cos(Q*PI/2 + D) = -Sin(D) if Q mod 4 = 1 -Cos(D) if Q mod 4 = 2 Sin(D) if Q mod 4 = 3 Tan(Q*PI/2 + D) = Sin(Q*PI/2 + D)/Cos(Q*PI/2 + D) Sinh(Q*Ln(2)+D) = 2^(Q-1)*[Cosh(D) + Sinh(D) - (Cosh(D) - Sinh(D))/2^(2*Q)] Cosh(Q*Ln(2)+D) = 2^(Q-1)*[Cosh(D) + Sinh(D) + (Cosh(D) - Sinh(D))/2^(2*Q)] Tanh(Q*Ln(2)+D) = Sinh(Q*Ln(2)+D)/Cosh(Q*Ln(2)+D) ArcTanh(1 - M/2^E) = ArcTanh(T) + Ln(2)*E/2 where T = (2 - M - M/2^E)/(2 + M - M/2^E) Exp(Q*Ln(2)+D) = 2^Q*(Cosh(D) + Sinh(D)) Ln(M*2^E) = Ln(M) + E*Ln(2) where 0.5 <= M < 1.0 Sqrt(M*2^E) = 2^(E/2)*Sqrt(M) if E mod 2 = 0 0.5 <= M < 1.0 2^[(E+1)/2]*Sqrt(M/2) if E mod 2 = 1 0.25 <= M < 0.4 ALGORITHM Putting all the pieces together we wind up with an algorithm like this: CONST MaxBits = 32; (* Maximum bits of accuracy desired. *) TYPE Reply = RECORD x, y, z: REAL END; VAR Alpha: ARRAY [-1..1, 0..MaxBits] OF REAL; (* look up table *) Scale: ARRAY [-1..1] OF REAL; PROCEDURE InitLookUp; VAR pow2: REAL; m, k: INTEGER; next: INTEGER; BEGIN pow2 := 1.0; FOR k := 0 TO MaxBits DO Alpha[1, k] := ArcTan(pow2); Alpha[0, k] := pow2; IF k = 0 THEN Alpha[-1, k] := 0.0; (* Can't take ArcTanh of 1.0 *) ELSE Alpha[-1, k] := ArcTanh(pow2); END; pow2 := pow2/2.0; END; FOR m := -1 TO 1 DO k := ORD(m <> 1); next := 4; Scale[m] := 1.0; REPEAT Scale[m] := Scale[m]*Sqrt(1.0 + Scalb(-2*k, FLOAT(m))); IF (m = -1) AND (k = next) THEN next := 3*k + 1; ELSE INC(k); END; UNTIL k > NoBits; Scale[m] := 1.0/Scale[m]; END; END InitLookUp; PROCEDURE DoCordic (x, y, z: REAL; m: INTEGER; reduceAngle: BOOLEAN): Reply; CONST NoBits = 16; (* want 16 bits of accuracy *) VAR xn, yn, zn: REAL; start, k : INTEGER; nextK : INTEGER; positive : BOOLEAN; reply : Reply; BEGIN start := ORD(m <> 1); (* start = 0 if m=1, else start = 1 *) k := start; nextK := 4; REPEAT IF reduceAngle THEN positive := y > 0.0; ELSE positive := z < 0.0; END; IF positive THEN xn := x + Scalb(-k, FLOAT(m)*y); (* Scalb(k, x) = x*2^k *) yn := y - Scalb(-k, x); zn := z + Alpha[m, k] ELSE xn := x - Scalb(-k, FLOAT(m)*y); (* note that m = 0, -1, or 1! *) yn := y + Scalb(-k, x); zn := z - Alpha[m, k] END; x := xn; y := yn; z := zn; IF (m = -1) AND (k = nextK) THEN nextK := 3*k + 1; ELSE INC(k); END; UNTIL k > NoBits; reply.x := xn; reply.y := yn; reply.z := zn; RETURN reply; END DoCordic; (* x = multiple*interval + delta *) PROCEDURE Reduce (x, interval: REAL; VAR multiple: INTEGER; VAR delta: REAL); BEGIN x := ABS(x); IF x >= interval THEN multiple := TRUNC(x/interval); delta := x - FLOAT(multiple)*interval; ELSE multiple := 0; delta := x; END; END Reduce; PROCEDURE Sin (x: REAL): REAL; VAR answer : Reply; quadrant: INTEGER; angle : REAL; result : REAL; BEGIN Reduce(x, PI/2.0, quadrant, angle); answer := DoCordic(Scale[1], 0.0, angle, 1, FALSE); CASE quadrant MOD 4 OF 0: result := answer.y; | (* result = Sin(angle) *) 1: result := answer.x; | (* result = Cos(angle) *) 2: result := -answer.y; | (* result = -Sin(angle) *) 3: result := -answer.x; | (* result = -Cos(angle) *) END; IF x < 0.0 THEN result := -result END; RETURN result; END Sin; REFERENCES I should mention here that I have not read several of these references. Others I have only skimmed. The only one I have closely studied is [4]. [1] Cody and Waite, "Software Manual for Elementary Functions". Discusses the classical methods for calculating elementary functions (i.e. power series). Plauger gives this book an excellent recommendation. [2] Knuth, "SemiNumerical Algorithms". Presents algorithms for multiple-word floating point aritmetic and for an integer power function. [3] Plauger, "Standard C Library". Includes an implementation of the standard C math functions using algorithms from [1]. [4] Swartzlander, Earl Jr editor, "Computer Arithmetic" pg. 272, "A Unified Algorithm for Elementary Functions". This is a good explanation of the Cordic algorithm. The algorithm outlined here is based on this article. [5] "Computer Approximations". Discusses properties of power series using computers and their applications to elementary functions. [6] IEEE Transactions on Computers, Jan 1991. "Expanding the Range of Convergence of the Cordic Algorithm." Outlines an efficient method for extending the range of the Cordic routines so that arguments don't have to be reduced. jesjones@milton.u.washington.edu