In 1956, engineers at Convair were asked to develop a digital computer to replace the B-58 Hustler bomber’s existing analog navigation system [1]. This was quite a challenge: the first commercial silicon bipolar junction transistors were only announced in 1954, so fully transistorized computers were still very new [2]. In addition, navigation requires the calculation of trigonometric functions in real-time, something existing algorithms running on the transistors of the day were much too slow to accomplish.

Existing analog navigation systems made use of electromechanical resolvers to calculate rotations of vectors and determine their magnitude and angle. There was no digital equivalent to a resolver so “computation of navigation problems required either series expansions, approximations or table lookup of individual trigonometric functions” [2]. The COrdinate Rotation DIgital Computer (CORDIC) was developed to replace the resolvers [3].

Jack Volder’s initial inspiration for the CORDIC algorithm came from “the massaging of the basic angle addition equations to obtain the following interesting equations” [1]. If then:

He realized that these operations could be repeated for increasing powers of which would allow the rotation of vector through any angle. It may not be immediately apparent from looking at these equations, but this is an implementation of trigonometric functions using only addition and subtraction, bit shifts, a small lookup table, and a single multiplication. In the sections that follow, we’ll attempt to unravel this.

Coordinate rotation

Let’s set the stage: we start with the Cartesian coordinates and which we want to rotate by some angle to determine a new set of coordinates and . These coordinates can be interpreted as the components of a vector with magnitude and angle which is rotated by to obtain .

Coordinate rotation

Figure 1. Graphical representation of coordinate rotation

Mathematically speaking, this can be represented as

which is pretty straightforward; however, we want to eliminate the trigonometric functions and develop an algorithm a digital computer can solve quickly.

We’ll begin by considering and . These equations involve the addition of angles, for which we have identities:

Making use of these identities, we obtain

into which we can substitude our original equations for and ,

and pull a from each equation to obtain

Now we let and to get

and make use of another identity,

Since , this simplifies our equations nicely:

And just like that, the trigonometric functions are gone! Well, not quite , so unless we already know the tangent of our rotation angle , we haven’t really solved anything. What’s worse, it also requires a square root, division, and a fair bit of multiplication which is no good at all.

A table for two

Suppose, for a moment, that we had a small table containing values of and . If these values were chosen carefully, any rotation angle could be determined by iteratively applying rotations from the table until the desired rotation angle was achieved. Instead of calculating and directly, we’ll calculate and , increasing until and where is the number of iterations. We’ll call our starting coordinates and .

We’ll use a very simple algorithm: we’ll set to our desired rotation angle, and then add or subtract each angle from the table to such that . If we’ll subtract the table angle, and if we’ll add the table angle. For mathematical convenience, we’ll define the direction of rotation as which can only be . Our equations are thus rewritten in an iterative form:

The first simplification we’ll make is to address the ugly inverse square root term, which we will call :

Since and are both linear combinations of and , it’s fairly easy to see what will happen as increases:

The coefficient in front of th term, , can clearly be expressed as the product of a sequence:

This means we can either add a new column to our table with each value, or if is fixed, we can simply multiply the final results and by . Let’s put it aside for now, and simply ignore :

Now we’re down to two muliplications per iteration, plus two final (or initial) multiplications to account for . This is significantly better than before, but we can do better still.

10 types of people

In the 1959 paper, Volder states that “these algorithms are suitable only for use with a binary code [which] possibly accounts for their late appearance as a numerical computing technique” [3]. So far what we’ve accomplished isn’t really that optimized for binary (as opposed to decimal) computation. We can change that by selecting a clever value for :

And just like that, we’ve arrived at the CORDIC. For a binary computer, multiplying any number by is simply dividing it by a multiple of 2 which is equivalent to bit shifting right by . For example, suppose and :

Our equations are now:

At each iteration, we either add or subtract a bit shifted version of the previous iteration’s values and add or subtract a value from a table. No multiplication, no division, no square roots, no trigonometric functions.

Vectoring

Figure 2. pseudorotations converge to a scaled version of our desired rotation

Each iteration is not a true rotation of the original vector due to the coefficient we’ve ignored. Instead we’re adding a scaled, orthogonal version of the previous iteration to itself to calculate the current iteration. We’ll refer to this as pseudorotation. It’s apparent in figure 2 that these pseudorotations do indeed converge on our desired rotation angle and that their magnitudes converge to some multiple of the original vector.

Figure 2 shows six iterations that seem to approach our desired rotation pretty nicely. Is six iterations, in general, enough?

A table for 10

Let’s try to determine how many iterations we’ll need and thus how large our table is. Assuming our numbers are all represented as -bit two’s complement, there’s going to be a fixed limit for how many iterations we can do. The range of numbers that can be represented is to . The absolute value of the largest possible number we can represent is thus or . As increases, we shift right by more and more bits. We’re therefore going to hit at some point. For an -bit number, once we’ve shifted right times, we’re done:

Note that this is arithmetic shifting which maintains the sign bit and rounds towards negative infinity. The result is clear, though. When (an 8-bit signed number), nothing changes past iteration . This sets an upper bound on our table size. What does such a table look like?

That’s only 8 numbers to store since there’s no need to store . For an -bit algorithm, we need to store a table of -bit numbers, or bits total. We can’t forget , though, which is another number we’ll need to store.

Let’s go back to our definition of and substitute our chosen value for :

It turns out this product converges pretty quickly:

By dropping the coefficients, our final results and are larger than and . In some circumstances this isn’t an issue. It’s more likely, though, that we’ll need to multiply and by . However, this is just a single multiplication per component.

Putting it all together, we have the following equations:

Turning it around

Now that we’ve developed the rotation algorithm, going in the other direction is straightforward: Given some and , suppose we want the corresponding vector ’s magnitude and angle . This is referred to as translation or vectoring. The way to approach this is to use the rotation algorithm, but instead of iterating until , we rotate the vector such that . Once complete, and .

Vectoring

Figure 3. Graphical representation of vectoring

Our iterative equations are the same,

but the algorithm is modified slightly. We let since we don’t have any angle information to start with, and instead of choosing based on , we’ll choose it based on : if we’ll rotate clockwise by setting , and if we’ll rotate counterclockwise by setting .

The limits

Starting our table at places a limit on our possible rotation angles. The maximum angle we can rotate through is given by

For a 16-bit number, . This means that our existing algorithm can’t actually rotate beyond . The simplest way to solve this is to perform an initial coarse rotation of the input vector:

This is essentially applying a rotation matrix to the input vector to prerotate it by . The prerotation is also subtracted or added to our initial angle to bring it into the required range. And it’s pretty much that simple.

References

[1] J.E. Volder, “The Birth of CORDIC,” J. VLSI Signal Processing, vol. 25, n. 2, pp. 101–105, June 2000. [Online].

[3] Texas Instruments, “Semiconductor Timeline - 1954: First commercial silicon transistor.” [Online].

[2] J.E. Volder, “The CORDIC Computing Technique,” Proc. Western Joint Computer Conf., pp. 257–261, 1959. [Online].