Mike teaches Mathematics and programming at Madison Area Technical College, Madison, WI 53704.
CORDIC stands for Coordinate Rotation Digital Computer, an early device implementing fast integer sine and cosine calculations (Volder, 1959). Whenever trigonometry functions must be evaluated repeatedly, as in computer graphics, integer methods, such as CORDIC, should be considered. While integer methods are less accurate than the C Library functions sin and cos, the improved speed makes the tradeoff quite acceptable in some applications.
CORDIC Units
The key to the CORDIC method is the representation of both angles and trigonometric ratios as integers. In this implementation a 16-bit unsigned integer represents the angles around the circle as shown in Table 1.
With CORDIC angle units, or CAUs, the circle divides into 64K parts instead of 360 parts using degrees. Each degree measures about 182 CAU.
Sine and cosine values are represented as signed integers, with an implicit denominator of 16384 (CordicBase in the accompanying program). Calculated sines and cosines lie in the range -16384 to +16384, corresponding to a trigonometric ratio between -1 and + 1. Table 2 contains sample correspondents in this fixed point scheme.
Suppose your application receives a value of 100 and must multiply it by sin(54°) to produce the nearest integer (a realistic example from computer graphics where input and output are pixel locations). The standard approach, using the C Standard Library sin call, amounts to:
100 x sin(54°) = 100 x 0.8090 = 80.90 -> 81Both the floating-point multiplication and the sin are expensive.
The CORDIC version of this calculation substitutes a fast sin and long integer multiplication (where 54° = 9830 CAU and sin(9830 CAU) = 13255 CORDIC fixed-point units):
100 x sin(9830) = (100 x 13255 + 8192) / 16384 = 81.40 -> 81The 8192 is added for integer rounding (8192/16384 = 0.5). The division by 16384 is done with an inexpensive right shift. CORDIC's speed is due to the fast sine calculation and the complete avoidance of floating-point calculations.
CORDIC Special Angles
The CORDIC algorithm depends on representing a given angle by a set of special angles, {arctan(2-i)} = {arctan(1), arctan(1/2), arctan(1/4), ...}:
arctan( 1) = 45.00° = 8192 CAU
arctan( 1/2) = 26.57° = 4836 CAU
arctan( 1/4) = 14.04° = 2555 CAU
arctan( 1/8) = 7.13° = 1297 CAU
arctan(1/16) = 3.58° = 651 CAU
.
.
.
Using these special angles, 54 is represented to a finer and finer precision as follows:
54° = 45.00 = 45.00 54° = 45.00 + 26.57 = 71.57 54° = 45.00 + 26.57 - 14.04 = 57.53 54° = 45.00 + 26.57 - 14.04 - 7.13 = 50.40 54° = 45.00 + 26.57 - 14.04 - 7.13 + 3.58 = 53.98This approximation has a physical interpretation. Think in terms of a vector 16384 units long and emanating from the origin in a standard x-y coordinate system. Starting at 0° along the positive x axis, the vector rotates through each of the special angles one step at a time. Rotation at each step may be clockwise or counter-clockwise, whichever is needed to bring the vector closer to 54°. The special angles represent rotations by smaller and smaller amounts, with positive signifying counter-clockwise rotation and negative signifying clockwise rotation.
Starting at 0° along the positive x axis, +45.00° is a counter-clockwise rotation into the middle of the first quadrant. The next step again rotates counter-clockwise by another 26.57°, but this results in 71.57°, which overshoots the mark. The third rotation is therefore clockwise, signified by the minus 14.04° bringing the result back to 57.53°. This continues as many times as there are bits in 16384 — 14 times.
The x and y coordinates of the rotating vector are the cosine and sine of the vector's angle at each step, assuming that the vector's length, 16384, doesn't change during rotation. As the vector's angle approaches 54° more closely, its x and y coordinates provide better approximations of cos(54°) and sin(54°).
CORDIC Equations
Figure 1 illustrates the geometry of a counter-clockwise rotation at the ith step. Rotate vector p through an angle of arctan(2-i) to new vector p1 such that the indicated angle at p is a right angle. The task is to express the new (cosine, sine) approximations (x1, y1) in terms of the old ones (x, y). The two shaded triangles are similar because they are right triangles with acute angles that are equal. In the right triangle connecting p and p1 to the origin, the two legs R* and R are in the proportion 1 : 2i by the definition of arctan(2-i)) :
2<SUP>-i</SUP> = tan(arctan(2<SUP>-i</SUP>)) = (R*)/RR and R* are the hypotenuses of the two similar triangles, so these triangles are in the proportion 1 : 2-i. This implies that the shorter horizontal leg of the small triangle is y/(2i) and the longer, vertical leg is x/(2i), leading to the counter-clockwise equations:
ccw rotation x1 = x - y/(2<SUP>i</SUP>)
y1 = y + x/(2<SUP>i</SUP>)
Rotating p1 clockwise leads to the clockwise equations, identical except for a sign reversal:
cw rotation x1 = x + y/(2<SUP>i</SUP>)
y1 = y - x/(2<SUP>i</SUP>)
These are fast operations involving integer addition, subtraction and shifting.
Expansion Factor
A problem arises with p1, which is further from the origin than p was at each step, so the vector does not simply rotate around a circle of radius 16384, but also expands. The expansion can be exactly measured by applying the Pythagorean theorem to the right triangle containing p and p1:
p1<SUP>2</SUP> = R<SUP>2</SUP> + (R/2<SUP>i</SUP>)<SUP>2</SUP> = R<SUP>2</SUP> + R<SUP>2</SUP>/2<SUP>2i</SUP> = R<SUP>2</SUP> * (1 + 1 / 2<SUP>2i</SUP>) = R<SUP>2</SUP> * (2<SUP>2i</SUP> + 1) / 2<SUP>2i</SUP>or:
The first rotation (45°) expands the vector by a factor of
the second rotation (26.57°) expands it further by a factor of
the third rotation (14.04°) expands by
and so on. Each of the 14 rotations entails such an expansion, although the later ones are negligible.
The net effect is an expansion by a factor of 1.414 x 1.118 x 1.031 x ... = 1.646760. The original vector will expand to 1.646760 x 16384 = 26981 in the course of rotating. To offset this expansion, the original vector is contracted before starting the process to bring its final length, after the rotation/expansions, back to 16384. Instead of initializing the vector to 16384, it is initialized to xInit = 16384/1.646760 = 9949, which expands to 16384 after 14 steps. At the last step the vector's x and y coordinates will be the cosine and sine in CORDIC fixed-point units (based on 16384).
Implementation Notes
The central routine is SinCos (See Listing 1) , which calculates the sine and cosine of an incoming angle. Both incoming angle and calculated sine and cosine are assumed to be in CORDIC units. SinCosSetup initializes needed global variables, including the special arctan angles and xInit, the initial contracted vector length. SinCosSetup must be called once only for initialization before calling SinCos. The CORDIC algorithm in SinCos works on first quadrant angles only (0 — 90, or 0 — 16383 CAU). SinCos translates angles from other quadrants into the first quadrant before applying the algorithm. The calculated sine and cosine will be correct, except possibly for the sign, which is adjusted before returning from the routine.
A hexagon rotator is included to demonstrate the CORDIC system (see Figure 2) . A center point and initial vertex remain fixed while another point, vertex1, rotates clockwise around the original vertex in increments of 650 CAU (3.57). For each vertex1, the five associated regular hexagon vertices are calculated and the hexagon is drawn using Borland C++'s line drawing commands. The routine calculating the vertices, CalcHexPts, calculates the sine and cosine of 60° (10923 CAU) only once and calculates each vertex from the previous one. This practice cannot be used for continual rotating because of accumulating errors from integer rounding as well as inexact sines and cosines.
Timing tests show the CORDIC method to be over 20 times faster than the standard method when calculating the vertices of a regular hexagon, a characteristic computer graphics task. The savings are due to the elimination of floating-point calculations as well as fast sine and cosine evaluation.
An exhaustive test of all 16384 CAUs shows that the worst error in a sine or cosine is 0.00064 and the average error is 0.00011. This is over 13 bits of accuracy on average, or better than one part in 8000, quite adequate for many screen-related computer graphics applications.
The original papers on this topic read very well, and are accessible through an excellent reprint by IEEE (Computer Arithmetic). Also helpful are articles cited below in Graphics Gems. Both volumes are essential for computer graphics workers.
References
Linhardt, R. J., and H. S. Miller 1969. "Digit-by-Digit Transcendental Function Computation". RCA Rev. 30:209-247. (Reprinted in E. Swartzlander (ed.) 1990. Computer Arithmetic. IEEE Computer Society Press, 233-271.)
Ritter, Jack. "Fast 2D-3D Rotation". In A. S. Glassner 1990. Graphics Gems. Academic Press, 440-441.
Turkowski, Ken. "Fixed-Point Trigonometry with CORDIC Iterations". In A. S. Glassner 1990. Graphics Gems. Academic Press, 494-497.
Volder, Jack E. 1959. "The CORDIC Trigonometric Computing Technique". IRE Trans. Electron. Comput. EC-8:330-334. (Reprinted in E. Swartzlander (ed.) 1990. Computer Arithmetic. IEEE Computer Society Press, 226-230.)
