Exploring the CORDIC algorithm with the WP-34S

[Thomas Klemm]

Introduction
In a recent thread somebody mentioned John Stephen Walther which lead me to his paper: A Unified Algorithm for Elementary Functions. This made me look again at the CORDIC algorithm.

I like to hit a key repeatedly and look what happens:

Quote: Try this during a boring lecture: Set your calculator to radians mode and then repeatedly press the cos button until you obtain the fixed point.

Structure and Interpretation of Computer Programs
Footnote 57

Rotation
Since CORDIC is all about rotation I wanted to use complex multiplication. The WP-34S provides a command that simulates this with the stack:

Code:

[cmplx][times]

Let's assume we want to calculate arg(4, 3), that is the angle of the complex number z = 4 + 3i.
How would we fill up the stack?
For reasons that become apparent later I want to keep the imaginary part of z in register X. Thus X and Y are swapped and the usual rotation becomes counter-clockwise.

We use for instance w = 1 + 0.1i.
This will rotate z by tan^-1(0.1) = 5.7105931375°.

T: 0.1
Z: 1
Y: 4
X: 3

We run ©× multiple times until the register X turns negative:

\(\begin{matrix}
Y: & 4 & 4.3 & 4.56 & 4.777 & 4.9484 & 5.07203 & 5.146176 & 5.1696017 \\
X: & 3 & 2.6 & 2.17 & 1.714 & 1.2363 & 0.74146 & 0.234257 & -0.2803606
\end{matrix}\)

But now we've gone one step too far. Therefore let's get the last value back:

Code:

[cmplx]x[<->] L

Thus in total we executed the rotation 6 times which is about 34.263558825°.

Let's write a little program for that. We use register A as counter:

Code:

[cmplx][times]
INC A
x>0?
BACK 003
[cmplx]x[<->] L
DEC A

Now you can probably see why I wanted to keep the imaginary part of z in register X: we have to check whether it is still positive.

Iteration
What's the next step of the algorithm? We have to reduce the angle of rotation and use w = 1 + 0.01i instead.
Thus we're going to shift the imaginary value of w one digit to the right using SDR:

Code:

R[^]
SDR 001
R[v]

This is the small program that we have so far:

Code:

0001 [cmplx][times]
0002 INC A
0003 x>0?
0004 BACK 003
0005 [cmplx]x[<->] L
0006 DEC A
0007 R[^]
0008 SDR 001
0009 R[v]
0010 END

Initialize the stack and registers with:

A: 0
T: 1
Z: 1
Y: 4
X: 3

And then run the program multiple times to see what happens:

\(\begin{matrix}
& X & Y & A \\
0: & 3.0000 & 4.0000 & 0 \\
1: & 0.2343 & 5.1462 & 6 \\
2: & 0.0283 & 5.1525 & 10 \\
3: & 0.0025 & 5.1525 & 15 \\
4: & 0.0005 & 5.1525 & 19
\end{matrix}\)

From the counter A we can conclude that the total angle of rotation was:
6tan^-1(0.1) + 4tan^-1(0.01) + 5tan^-1(0.001) + 4tan^-1(0.0001) = 36.8647107295°

Therefore we should probably reset the counter A to 0. But before we have to multiply it by the corresponding angle and add it to the total of angles.

Initialization of Constants
This little program will fill the registers with the corresponding angles:

Code:

0001 # 013
0002 SDR 003
0003 STO J
0004 # 001
0005 ATAN
0006 STO[->]J
0007 RCL L
0008 SDR 001
0009 ISG J
0010 BACK 005
0011 END

Prototype
Since we don't want to loose the values of the stack we have to use SSIZE8 for further calculations.
Please keep in mind that now both registers A and B are part of the stack.
We keep using J as the loop control variable.

Code:

0001 [cmplx][times]
0002 INC A
0003 x>0?
0004 BACK 003
0005 [cmplx]x[<->] L
0006 DEC A
0007 R[^]
0008 SDR 001
0009 R[v]
0010 SSIZE8
0011 RCL A
0012 STO- B
0013 RCL[times][->]J
0014 STO+ C
0015 R[v]
0016 SSIZE4
0017 ISG J
0018 BACK 017
0019 RCL B
0020 END

Initialize the stack and registers with:

J: 0.013
B: 0
A: 0
T: 1
Z: 1
Y: 4
X: 3

And then run the program to get: 36.8698976458
Compare this to tan^-1(0.75).

Finalisation
We're not finished yet: it's a hassle that we have to initialize the stack and registers for each run.
Thus let's do that beforehand. Add a nice label as well:

Code:

0001 LBL'ARG'
0002 # 013
0003 SDR 003
0004 STO J
0005 CLx
0006 STO A
0007 STO B
0008 R[v]
0009 # 001
0010 ENTER[^]
0011 [cmplx]x[<->] Z
0012 [cmplx][times]
0013 INC A
0014 x>0?
0015 BACK 003
0016 [cmplx]x[<->] L
0017 DEC A
0018 R[^]
0019 SDR 001
0020 R[v]
0021 SSIZE8
0022 RCL A
0023 STO- B
0024 RCL[times][->]J
0025 STO+ C
0026 R[v]
0027 SSIZE4
0028 ISG J
0029 BACK 017
0030 RCL B
0031 END

This final version of the program can be used to calculate arg(4, 3):

4 ENTER
3
XEQ'ARG'
36.8698976458

Addendum
On request of Tugdual I attached a Python program with examples from the aforementioned paper.