CORDIC algorithms: How do calculators calculate?
motivated by a comment by Doug K.


Binary arithmetic

If you had a number, like 123.456 and you wanted to multiply by 100 you'd just ...

>Move the decimal point to the right!
Exactly. Move it 2 places to the right, since 100 = 102
And to multiply by 1/10000 = 10-4, you'd move the decimal point 4 places to the left.
These have to be the world's easiest multiplications, eh? That's because we work with numbers to the base 10.
That is: 123.456 = (1)102+(2)101+(3)100+(4)10-1+(5)10-2+(6)10-3.

Aah, but if you were a computer, working with the binary number system (base = 2), your numbers would look like:
110.101 = (1)22+(1)21+(0)20+(1)2-1+(0)2-2+(1)2-3.

Now what's the worlds simplest multiplication? It'd be ...

>Multiply by 2.
Well, multiplications by powers of 2.

Okay, now suppose a computer wanted to do some elaborate calculations.
If it could arrange that all multiplications were by powers of 2, it'd be FAST!
In fact, in a binary computer, moving the decimal point is just a matter of "shifting" the number.

>Huh?
Imagine the number 110.101 stored in a 16 bit computer like so:
        000000110.1010000
To divide by 2, the computer would just shift the decimal to the left one place (or shift the number one place to the right), getting:
        000000011.0101000
Neat, eh?

>Yeah, but how does it calculate?
Let's do an example:


CORDIC example

Consider the equations for rotation:
      x' = (x cos θ - y sin θ)
      y' = (y cos θ + x sinθ)

These can be rewritten as::
      x' = cos θ (x - y tan θ)
      y' = cos θ (y + x tanθ)

Note that cos θ = 1/sqrt(1 + tan2 θ) so that:
      x' = (x - y tan θ) / sqrt(1 + tan2 θ)
      y' = (y + x tanθ) / sqrt(1 + tan2 θ)

Now we restrict the successive angles of rotation to those for which tan(angle) = 2-n.
Then we'd get:
      xn+1 = (xn - yn 2-n) / sqrt(1 + 2-2n)
      xn+1 = (yn + xn 2-n) / sqrt(1 + 2-2n)

We begin with a couple of numbers (in binary, of course):
      x0, y0.
Then we calculate a second pair via:
      x1 = (x0 - y0/21) / sqrt(1+2-1)
      y1 = (y0 + x0/21) / sqrt(1+2-1)
Then another pair:
      x2 = (x1 - y1/22) / sqrt(1+2-2)
      y2 = (y1 + x1/22) / sqrt(1+2-2)
Continuing ...

>Wait! Where are you going?
Patience. However, notice that there's just some adding, subtracting and no "real" multiplication ... just shifting.
Note, too, that the procedure is simple: it's iterative, repeating the same ritual over and over again.
Note that we're actually applying successive rotations to the initial point (x0, y0).
If we start with the pair (1,0) and continue the ritual we started above, we'd get (as n ∞):
      xn 0.5753
      yn 0.8180

In other words, we've rotated the point (1, 0) by an angle ... so it's become (0.5753, 0.8180)
Normally, to rotate by an angle θ, we'd use the magic formulas
      x' = x cos θ - y sin θ
      y' = y cos θ + x sin θ

That'd involve calculating functions like cos θ and sin θ and that ain't easy.
For our computer, it just did some adding, subtracting and ...

>And shifting. Yes, I see, but what's the angle of rotation and do you get the same result and what's CORDIC?    
Patience. We're getting there.


CORDIC rotation

In the above example, we were rotating by smaller and smaller angles.

If we add all the successive rotations, starting at (1, 0), we'd get the total angle of rotation. That is:
      Total Angle = arctan(2-1) + arctan(2-2) + arctan(2-3) + arctan(2-4) + ... = 0.9579.

>Is that in radians?
Is there any other measure?
In degrees, it'd be (180/ π)(0.9579) or about 55 degrees.

>And did you know that beforehand? Suppose I gave you the angle, like maybe 30 degrees.
I assume you mean π/6.
Okay, let's start with (x0, y0) at π/6 and rotate clockwise this time ... until y = 0.
To make sure we don't go past the x-axis (with our mini-rotations), we'd keep an eye on the y value and rotate counterclockwise if y < 0.
Of course, we're making microscopic rotations so we'd rotate in that opposite direction by an angle whose tangent is some very small number: 2-n.
If y goes positive (meaning we've gone past the x-axis), we change direction again.

>And the next angle of rotation is even smaller, right?
Exactly! In fact, we'd be using the iterative process:
      xn+1 = (xn - dn yn 2-n) / sqrt(1 + 2-2n)
      xn+1 = (yn + dn xn 2-n) / sqrt(1 + 2-2n)
where dn changes sign if y changes sign (so we rotate in the opposite direction when y crosses the x-axis).
Note that xn+12 + yn+12 = xn2 + yn2

That way the computer will oscillate about the x-axis, closing in on the exact x-value. Example, starting at an angle of π/6
Here's a closeup of the "ending", illustrating the convergence to the final x-value:
Successive iterates
starting at (0.866, 0.500):

(0.866, 0.500)
(0.998, 0.060)
(0.983, -0.184)
(0.998, -0.061)
(1.000, 0.002)
(1.000, -0.029)
(1.000, -0.014)
(1.000, -0.006)
(1.000, -0.002)
(1.000, 0.000)


The total angle of rotation is (surprise!) π/6 and the distance of iterates from the origin is always (surprise!) 1 and it only took a few iterations.
>And this with just adding, subtracting and shifting, eh?
Yes.
>Is that's it?
Hardly.


Observations

Note that we can write the magic formulas for CORDIC rotation as:
      xn+1 = (xn - dn yn 2-n) / sqrt(1 + 2-2n) = Kn (xn - dn yn 2-n)
      xn+1 = (yn + dn xn 2-n) / sqrt(1 + 2-2n) = Kn (yn + dn xn 2-n)
where Kn = 1/sqrt(1 + 2-2n) and dn determines the direction of rotation. (dn = 1 for counterclockwise, dn = -1 for clockwise.)

In our first example, starting at (1,0), we moved counterclockwise until the sequence of iterates converged.
In our second example, starting at π/6, we moved clockwise and discovered that we could stop by introducing dn.
In general we use the above equations and stop at a given angle ... changing the sign of dn = ±1 if we pass the given angle.

Normally, one would start at (1,0) and rotate to a given angle.
Alas, the largest angle is given by the first example, namely (about) 55 degrees.
To accommodate angles larger than this (for example, π/3), we can start with a giant step (by angles of arctan(1+1/20 = π/4) before taking smaller steps.
Note that the first arctan value in Table 1 generates the angle π/4. We should then start with K0 = 1 / sqrt(1 + 2-0) = 1/sqrt(2).
Indeed, the equations for CORDIC rotation insist that you use K0 with x0, y0 and K1 with x1, y1, etc..
Successive iterations multiply the Ks:
      K0K1K2 ... Kn = 1/ sqrt[2(1+2-1)(1+2-2)(1+2-3)...(1+2-n)].

In fact, regardless of the angle involved in the rotation, we have the same set of Ks.
In fact:
      [ K0K1K2 ... Kn] 0.6072529350088812561694... as n
Since this don't hardly change with the angle, the computer could simply store this number and regurgitate when it's needed.
In fact, the computer could iterate without the Ks and apply the K-factor at the end of the iterations.

Further, note that we must keep checking the angle to see if we've rotated enough or have gone beyond the desired angle of rotation      
(in which case we change the direction of rotation).
That involves keeping track of
      Total Angle = arctan(2-0) + arctan(2-1) + arctan(2-2) + arctan(2-3) + ...

That suggests that we'd have to evaluate a whole bunch of arctan functions.
In fact, the computer would simply store a bunch of arctan values.
(Since it'd rarely take more than a few dozen iterations for 10 decimal place accuracy, only a few dozen values need be stored.)

> I wish you had said that earlier. I figured you were confused and ...
Listen! I'm just learning this stuff.

stored
arctan values
0.785398...
0.463648...
0.244979...
0.124355...
0.062419...
0.031240...
0.015624...
0.007812...
0.003906...
0.001953...
0.000977...
0.000488...
0.000244...
0.000122...
0.000061...
0.000031...
0.000015...
0.000008...
0.000004...
...
Table 1
Finally we note that, starting at x0 = 1 and y0 = 0 and a specified angle θ, the final values are: xn = cos θ, yn = sin θ.

>Aha! So it can evaluate sines and cosines!
Yes ... and more. Inverse trig function, logs and exponentials, etc. etc. The guys you'd normally find on a hand-held calculator.
You just have to be very clever ... to writing iterative equations where the iterates converge to the desired value and ...

>And you forget about multiplying, right?
Right. For example, if you want arcsin θ, you can use the above algorithm, rotating until y = θ.
Then y = θ = sin β hence the (eventual) angle of rotation is β = arcsin θ.

Further, one can generate the hyperbolic function, tanh(), using the magic equations:
      xn+1 = Kn (xn + dn yn 2-n)
      xn+1 = Kn (yn + dn xn 2-n)
where the incremental "angles" are Kn = arctanh(2-n) and dn determines the direction of "rotation".
As before, values of Kn would be stored and x = cosh β, y = sinh β and, starting at (1,0), xn2 - yn2 = 1.
Note, however, that K1+K2+K3+... is bounded above, so there's a limit to the size of the "angle".
Note, too, that tanh(∞) = 1 so arctanh (2-0) = arctanh (1) = ∞

We have:
      tanh θ = (eθ - e) / (eθ + e)
      cosh θ = (eθ + e) / 2
      sinh θ = (eθ - e) / 2
Hence cosh θ + sinh θ = eθ, so exponentials can be calculated.

Further, we can show that log(z) = 2 arctanh[ (z-1 ) / (z+1) ] for any (respectable) z.

The spreadsheet illustrates all this. That is, it'll calculate sinh θ, cosh θ, tanh θ, and eθ and log θ (to the base e).

>Using just adding and shifting?
No, of course not. The spreadsheet just illustrates the convergence of the iterative schemes and plots pretty graphs.

Note that the hyperbolic functions are closely related to the circular trig functions ... hence the similarity in CORDIC machinations.

>So what's the CORDIC story?
Good question.


Who is CORDIC?

Apparently, Jack E. Volder saw the following equation in the 1946 edition of the CRC Handbook of Chemistry and Physics:

In 1959 Volder described the CORDIC method. (COordinate Rotation DIgital Computer)
It was developed at Convair to replace the analog ritual in the B-58 bomber's navigation computer.
The idea was to do digital calculation without having any multiplication facility.
It was originally implemented for computers using binary numbers (as we noted above), but has since been used in decimal systems, including pocket calculators.
It eventually was hard-wired into many computer Central Processing Units (CPUs) including the Intel 486 (the first x86 chip that used more than a million transistors).

There's a spreadsheet to play with. Click on the picture to download.