CORDIC algorithms: How do calculators calculate?

motivated by a comment by Doug K.
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 = 10^{2}
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)10^{2}+(2)10^{1}+(3)10^{0}+(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)2^{2}+(1)2^{1}+(0)2^{0}+(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:
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 + tan^{2} θ) so that:
x' = (x  y tan θ) / sqrt(1 + tan^{2} θ)
y' = (y + x tanθ) / sqrt(1 + tan^{2} θ)
Now we restrict the successive angles of rotation to those for which tan(angle) = 2^{n}.
Then we'd get:
x_{n+1} = (x_{n}  y_{n} 2^{n}) / sqrt(1 + 2^{2n})
x_{n+1} = (y_{n} + x_{n} 2^{n}) / sqrt(1 + 2^{2n})
We begin with a couple of numbers (in binary, of course):
x_{0}, y_{0}.
Then we calculate a second pair via:
x_{1} = (x_{0}  y_{0}/2^{1}) / sqrt(1+2^{1})
y_{1} = (y_{0} + x_{0}/2^{1}) / sqrt(1+2^{1})
Then another pair:
x_{2} = (x_{1}  y_{1}/2^{2}) / sqrt(1+2^{2})
y_{2} = (y_{1} + x_{1}/2^{2}) / 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 (x_{0}, y_{0}).
If we start with the pair (1,0) and continue the ritual we started above, we'd get (as n ∞):
x_{n} 0.5753
y_{n} 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.
 
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 (x_{0}, y_{0}) at π/6 and rotate clockwise this time ... until y = 0.
To make sure we don't go past the xaxis (with our minirotations), 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 xaxis), we change direction again.
>And the next angle of rotation is even smaller, right?
Exactly! In fact, we'd be using the iterative process:
x_{n+1} = (x_{n}  d_{n} y_{n} 2^{n}) / sqrt(1 + 2^{2n})
x_{n+1} = (y_{n} + d_{n} x_{n} 2^{n}) / sqrt(1 + 2^{2n})
where d_{n} changes sign if y changes sign (so we rotate in the opposite direction when y crosses the xaxis).
Note that x_{n+1}^{2} + y_{n+1}^{2} = x_{n}^{2} + y_{n}^{2}
That way the computer will oscillate about the xaxis, closing in on the exact xvalue. Example, starting at an angle of π/6
Here's a closeup of the "ending", illustrating the convergence to the final xvalue:

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.
Note that we can write the magic formulas for CORDIC rotation as:
x_{n+1} = (x_{n}  d_{n} y_{n} 2^{n}) / sqrt(1 + 2^{2n}) = K_{n} (x_{n}  d_{n} y_{n} 2^{n})
x_{n+1} = (y_{n} + d_{n} x_{n} 2^{n}) / sqrt(1 + 2^{2n}) = K_{n} (y_{n} + d_{n} x_{n} 2^{n})
where K_{n} = 1/sqrt(1 + 2^{2n}) and d_{n} determines the direction of rotation. (d_{n} = 1 for counterclockwise, d_{n} = 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 d_{n}.
In general we use the above equations and stop at a given angle ... changing the sign of d_{n} = ±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/2^{0} = π/4) before taking smaller steps.
Note that the first arctan value in Table 1 generates the angle π/4. We should then start with K_{0} = 1 / sqrt(1 + 2^{0}) = 1/sqrt(2).
Indeed, the equations for CORDIC rotation insist that you use K_{0} with x_{0}, y_{0}
and K_{1} with x_{1}, y_{1}, etc..
Successive iterations multiply the Ks:
K_{0}K_{1}K_{2} ... K_{n} = 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:
[ K_{0}K_{1}K_{2} ... K_{n}] 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 Kfactor 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 x_{0} = 1 and y_{0} = 0 and a specified angle θ, the final values are: x_{n} = cos θ, y_{n} = 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 handheld 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:
x_{n+1} = K_{n} (x_{n} + d_{n} y_{n} 2^{n})
x_{n+1} = K_{n} (y_{n} + d_{n} x_{n} 2^{n})
where the incremental "angles" are K_{n} = arctanh(2^{n}) and d_{n} determines the direction of "rotation".
As before, values of K_{n} would be stored and x = cosh β, y = sinh β and, starting at (1,0), x_{n}^{2}  y_{n}^{2} = 1.
Note, however, that K_{1}+K_{2}+K_{3}+... 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[ (z1 ) / (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.
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 B58 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 hardwired 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.
