# Understanding the CORDIC algorithm with Python

Started by May 6, 2009
```Hello, I've been trying to figure out the way the CORDIC algorithm
works and I am not sure if this is the right place to post this, so
feel free to point me in the right direction to get help.

I usually code a quick prototype in python for verification and
understanding of an algorithm I need to implement. For those of you
who don't know Python but know C, you should be able to follow it
since python has a similar syntax to C. I've been struggling with this
somehow I messing something up.

import math

# I know CORDIC is only valid for inputs between
# -pi/2 and pi/2; I am not exactly sure what I need
# to do add to make it where any input is acceptable
# I believe is keep on adding/subtracting pi, but I don't
# get why this works?
def cordic_trig(beta,N=40):
# in hardware, put this in a table.
def K_vals(n):
K = []
acc = 1.0
for i in range(0, n):
acc = acc * (1.0/math.sqrt(1 + 2.0**(-2*i)))
K.append(acc)
return K

#K = K_vals(N)
K = 0.6072529350088812561694
# emulation for hardware lookup table
atans = [math.atan(2.0**(-i)) for i in range(0,N)]
#print K
#print atans
x = 1
y = 0

for i in range(0,N):
d = 1.0
if beta < 0:
d = -1.0

x = (x - (d*(2.0**(-i))*y))
y = ((d*(2.0**(-i))*x) + y)
# in hardware put the atan values in a table
beta = beta - (d*math.atan(2**(-i)))
return (K*x, K*y)

if __name__ == '__main__':
beta = math.pi/6.0
print "Actual cos(%f) = %f" % (beta, math.cos(beta))
print "Actual sin(%f) = %f" % (beta, math.sin(beta))
cos_val, sin_val = cordic_trig(beta)
print "CORDIC cos(%f) = %f" % (beta, cos_val)
print "CORDIC sin(%f) = %f" % (beta, sin_val)
```
```On Wed, 06 May 2009 14:33:14 -0700, Chico wrote:

> Hello, I've been trying to figure out the way the CORDIC algorithm
> works and I am not sure if this is the right place to post this, so
> feel free to point me in the right direction to get help.
>
> I usually code a quick prototype in python for verification and
> understanding of an algorithm I need to implement. For those of you
> who don't know Python but know C, you should be able to follow it
> since python has a similar syntax to C. I've been struggling with this
> somehow I messing something up.

>         x = (x - (d*(2.0**(-i))*y))
>         y = ((d*(2.0**(-i))*x) + y)

You're calculating y[i+1] from x[i+1], when it should be x[i].

Rewriting the above as:

(x,y) = (x - (d*(2.0**(-i))*y), (d*(2.0**(-i))*x) + y)

to update x and y concurrently gives the correct result.

```
```On May 6, 7:41=A0pm, Nobody <nob...@nowhere.com> wrote:
> On Wed, 06 May 2009 14:33:14 -0700, Chico wrote:
> > Hello, I've been trying to figure out the way the CORDIC algorithm
> > works and I am not sure if this is the right place to post this, so
> > feel free to point me in the right direction to get help.
>
> > I usually code a quick prototype in python for verification and
> > understanding of an algorithm I need to implement. For those of you
> > who don't know Python but know C, you should be able to follow it
> > since python has a similar syntax to C. I've been struggling with this
> > CORDIC algorithm, seems straight forward from the wikipedia page, but
> > somehow I messing something up.
> > =A0 =A0 =A0 =A0 x =3D (x - (d*(2.0**(-i))*y))
> > =A0 =A0 =A0 =A0 y =3D ((d*(2.0**(-i))*x) + y)
>
> You're calculating y[i+1] from x[i+1], when it should be x[i].
>
> Rewriting the above as:
>
> =A0 =A0 =A0 =A0 =A0(x,y) =3D (x - (d*(2.0**(-i))*y), (d*(2.0**(-i))*x) + =
y)
>
> to update x and y concurrently gives the correct result.

Thank you so much for the help (whoever you are), those are little
subtle bugs that can drive you crazy. Now I know why our Prof. in
Functional programming was making a big deal out of assignment
statements! :)

As for the domain between -pi/2 and pi/2 do I just keep adding/
subtracting until I get it within that domain? Thanks.
```
```Chico wrote:
> On May 6, 7:41 pm, Nobody <nob...@nowhere.com> wrote:
>> On Wed, 06 May 2009 14:33:14 -0700, Chico wrote:
>>> Hello, I've been trying to figure out the way the CORDIC algorithm
>>> works and I am not sure if this is the right place to post this, so
>>> feel free to point me in the right direction to get help.
>>> I usually code a quick prototype in python for verification and
>>> understanding of an algorithm I need to implement. For those of you
>>> who don't know Python but know C, you should be able to follow it
>>> since python has a similar syntax to C. I've been struggling with this
>>> somehow I messing something up.
>>>         x = (x - (d*(2.0**(-i))*y))
>>>         y = ((d*(2.0**(-i))*x) + y)
>> You're calculating y[i+1] from x[i+1], when it should be x[i].
>>
>> Rewriting the above as:
>>
>>          (x,y) = (x - (d*(2.0**(-i))*y), (d*(2.0**(-i))*x) + y)
>>
>> to update x and y concurrently gives the correct result.
>
> Thank you so much for the help (whoever you are), those are little
> subtle bugs that can drive you crazy. Now I know why our Prof. in
> Functional programming was making a big deal out of assignment
> statements! :)
>
> As for the domain between -pi/2 and pi/2 do I just keep adding/
> subtracting until I get it within that domain? Thanks.

Usually with embedded trig algorithms you just implement a single
quadrant (0 to pi/2), but you can use a half-circle if you want.  And
you might need to negate the angle, or the result - I'll let you figure
out the details yourself.

This is the first time I've heard someone describe Python syntax as
being "similar to C" - it's about as different as you could get while
still remaining an imperative programming language.  But a well-written
Python program should be fairly easily understood by C programmers.

When posting sample code in a newsgroup (or other public forum), it's a
good idea to cut out unnecessarily or old code - keep it to a minimum.
Including the nested function K_vals(), for example, serves only to confuse.
```
```On Wed, 06 May 2009 19:45:04 -0700, Chico wrote:

> As for the domain between -pi/2 and pi/2 do I just keep adding/
> subtracting until I get it within that domain? Thanks.

For any trigonometric function f, start by using the 2*pi periodicity,
f(x+2*n*pi) = f(x) to coerce the value into the range 0..2*pi, e.g.
x = x % 2*pi or x = x - 2*pi * floor(x / 2*pi) (but check how % handles
negative values).

Then, for values outside of the range 0..pi/2, use the identities:

sin(x) = sin(pi - x)    for pi/2 <= x <= pi
sin(x) = -sin(x - pi)   for pi <= x <= 3*pi/2
sin(x) = -sin(2*pi - x) for 3*pi/2 <= x <= 2*pi

cos(x) = -cos(pi - x)   for pi/2 <= x <= pi
cos(x) = -cos(x - pi)   for pi <= x <= 3*pi/2
cos(x) = cos(2*pi - x)  for 3*pi/2 <= x <= 2*pi

tan(x) = -tan(pi - x)   for pi/2 <= x <= pi
tan(x) = tan(x - pi)    for pi <= x <= 3*pi/2
tan(x) = -tan(2*pi - x) for 3*pi/2 <= x <= 2*pi

This ensures that the argument is always in the range 0..pi/2.

IOW:

0..pi/2         n/a             n/a
pi/2..pi        x = pi - x      cos = -cos, tan = -tan
pi..3*pi/2      x = x - pi      sin = -sin, cos = -cos
3*pi/2..2*pi    x = 2*pi - x    sin = -sin, tan = -tan

```
```In article <pan.2009.05.07.12.31.53.141000@nowhere.com>,
nobody@nowhere.com says...
> On Wed, 06 May 2009 19:45:04 -0700, Chico wrote:
>
> > As for the domain between -pi/2 and pi/2 do I just keep adding/
> > subtracting until I get it within that domain? Thanks.
>
> For any trigonometric function f, start by using the 2*pi periodicity,
> f(x+2*n*pi) = f(x) to coerce the value into the range 0..2*pi, e.g.
> x = x % 2*pi or x = x - 2*pi * floor(x / 2*pi) (but check how % handles
> negative values).
>
> Then, for values outside of the range 0..pi/2, use the identities:
>
>         sin(x) = sin(pi - x)    for pi/2 <= x <= pi
>         sin(x) = -sin(x - pi)   for pi <= x <= 3*pi/2
>         sin(x) = -sin(2*pi - x) for 3*pi/2 <= x <= 2*pi
>
>         cos(x) = -cos(pi - x)   for pi/2 <= x <= pi
>         cos(x) = -cos(x - pi)   for pi <= x <= 3*pi/2
>         cos(x) = cos(2*pi - x)  for 3*pi/2 <= x <= 2*pi
>
>         tan(x) = -tan(pi - x)   for pi/2 <= x <= pi
>         tan(x) = tan(x - pi)    for pi <= x <= 3*pi/2
>         tan(x) = -tan(2*pi - x) for 3*pi/2 <= x <= 2*pi
>
> This ensures that the argument is always in the range 0..pi/2.
>
> IOW:
>
>         0..pi/2         n/a             n/a
>         pi/2..pi        x = pi - x      cos = -cos, tan = -tan
>         pi..3*pi/2      x = x - pi      sin = -sin, cos = -cos
>         3*pi/2..2*pi    x = 2*pi - x    sin = -sin, tan = -tan

Also known as the ASTC rule

All, Sin, Tan, Cos

--
Paul Carpenter          | paul@pcserviceselectronics.co.uk
<http://www.pcserviceselectronics.co.uk/>    PC Services
<http://www.pcserviceselectronics.co.uk/fonts/> Timing Diagram Font
<http://www.gnuh8.org.uk/>  GNU H8 - compiler & Renesas H8/H8S/H8 Tiny
<http://www.badweb.org.uk/> For those web sites you hate
```
```On May 7, 1:31=A0pm, Nobody <nob...@nowhere.com> wrote:
> On Wed, 06 May 2009 19:45:04 -0700, Chico wrote:
> > As for the domain between -pi/2 and pi/2 do I just keep adding/
> > subtracting until I get it within that domain? Thanks.
>
> For any trigonometric function f, start by using the 2*pi periodicity,
> f(x+2*n*pi) =3D f(x) to coerce the value into the range 0..2*pi, e.g.
> x =3D x % 2*pi or x =3D x - 2*pi * floor(x / 2*pi) (but check how % handl=
es
> negative values).
>
> Then, for values outside of the range 0..pi/2, use the identities:
>
> =A0 =A0 =A0 =A0 sin(x) =3D sin(pi - x) =A0 =A0for pi/2 <=3D x <=3D pi
> =A0 =A0 =A0 =A0 sin(x) =3D -sin(x - pi) =A0 for pi <=3D x <=3D 3*pi/2
> =A0 =A0 =A0 =A0 sin(x) =3D -sin(2*pi - x) for 3*pi/2 <=3D x <=3D 2*pi
>
> =A0 =A0 =A0 =A0 cos(x) =3D -cos(pi - x) =A0 for pi/2 <=3D x <=3D pi
> =A0 =A0 =A0 =A0 cos(x) =3D -cos(x - pi) =A0 for pi <=3D x <=3D 3*pi/2
> =A0 =A0 =A0 =A0 cos(x) =3D cos(2*pi - x) =A0for 3*pi/2 <=3D x <=3D 2*pi
>
> =A0 =A0 =A0 =A0 tan(x) =3D -tan(pi - x) =A0 for pi/2 <=3D x <=3D pi
> =A0 =A0 =A0 =A0 tan(x) =3D tan(x - pi) =A0 =A0for pi <=3D x <=3D 3*pi/2
> =A0 =A0 =A0 =A0 tan(x) =3D -tan(2*pi - x) for 3*pi/2 <=3D x <=3D 2*pi
>
> This ensures that the argument is always in the range 0..pi/2.
>
> IOW:
>
> =A0 =A0 =A0 =A0 Quadrant =A0 =A0 =A0 =A0pre =A0 =A0 =A0 =A0 =A0 =A0 post
> =A0 =A0 =A0 =A0 0..pi/2 =A0 =A0 =A0 =A0 n/a =A0 =A0 =A0 =A0 =A0 =A0 n/a
> =A0 =A0 =A0 =A0 pi/2..pi =A0 =A0 =A0 =A0x =3D pi - x =A0 =A0 =A0cos =3D -=
cos, tan =3D -tan
> =A0 =A0 =A0 =A0 pi..3*pi/2 =A0 =A0 =A0x =3D x - pi =A0 =A0 =A0sin =3D -si=
n, cos =3D -cos
> =A0 =A0 =A0 =A0 3*pi/2..2*pi =A0 =A0x =3D 2*pi - x =A0 =A0sin =3D -sin, t=
an =3D -tan

At school we learned that the trig functions that are positive in each

0..pi/2         sin, cos, tan  All
pi/2..pi        sin            Stations
pi..3*pi/2      tan            To
3*pi/2..2*pi    cos            Crewe

Might mean something if you are English ;->
```
```Chico wrote:
> As for the domain between -pi/2 and pi/2 do I just keep adding/
> subtracting until I get it within that domain? Thanks.

If you represent angle such that 0 to pi maps to 0 to 2^n then getting
within the domain only requires a bit-mask (no iteration).

Bob
```
```s0lstice wrote:
> Nobody <nob...@nowhere.com> wrote:
>
... snip ...
>
>>     0..pi/2         n/a             n/a
>>     pi/2..pi        x = pi - x      cos = -cos, tan = -tan
>>     pi..3*pi/2      x = x - pi      sin = -sin, cos = -cos
>>     3*pi/2..2*pi    x = 2*pi - x    sin = -sin, tan = -tan
>
> At school we learned that the trig functions that are positive
>
> 0..pi/2         sin, cos, tan  All
> pi/2..pi        sin            Stations
> pi..3*pi/2      tan            To
> 3*pi/2..2*pi    cos            Crewe

Much simpler, and not language dependent, is just think of the
signs applied to the x and y axis when you draw a rt triangle in
it.  x is positive with displacement to right, y with displacement
up.  Now all you need to remember is what sin, cos, tan mean, i.e.
the ratio of which sides.

--
[mail]: Chuck F (cbfalconer at maineline dot net)
[page]: <http://cbfalconer.home.att.net>

```
```On Thu, 07 May 2009 16:51:10 -0400, CBFalconer <cbfalconer@yahoo.com>
wrote:

><snip>
>Much simpler, and not language dependent, is just think of the
>signs applied to the x and y axis when you draw a rt triangle in
>it.  x is positive with displacement to right, y with displacement
>up.  Now all you need to remember is what sin, cos, tan mean, i.e.
>the ratio of which sides.

Your basic "unit circle" definitions.  Sometimes, I wonder if it is
even taught, anymore.  Few younger folks in today's college calculus
classes seem to even know about it -- and my latest experience with
them was this last Winter term.  Yet it's what I always go back to if
I want to gain some insight on a new question (or one I have long
since forgotten about.)  You can't go wrong with that and Pythagorean.
You also get all your half-angle and double-angle equivalences, etc.
It's all there to be had with merely a gander.

Also, there are symmetries to exploit beyond the 180 or 90 degree ones
that I've seen mentioned in this thread.  At a minimum, octants (or 45