ActiveState Code

Recipe 576792: Polar-to-rectangular conversions using CORDIC

Demonstrate CORDIC algorithm for computing coordinate conversions using only adds and shifts (plus one multiply in the outer loop).

Python 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
def to_polar(x, y):
    'Rectangular to polar conversion using ints scaled by 100000. Angle in degrees.'
    theta = 0
    for i, adj in enumerate((4500000, 2656505, 1403624, 712502, 357633, 178991, 89517, 44761)):
        sign = 1 if y < 0 else -1
        x, y, theta = x - sign*(y >> i) , y + sign*(x >> i), theta - sign*adj
    return theta, x * 60726 // 100000

def to_rect(r, theta):
    'Polar to rectangular conversion using ints scaled by 100000. Angle in degrees.'
    x, y = 60726 * r // 100000, 0
    for i, adj in enumerate((4500000, 2656505, 1403624, 712502, 357633, 178991, 89517, 44761)):
        sign = 1 if theta > 0 else -1
        x, y, theta = x - sign*(y >> i) , y + sign*(x >> i), theta - sign*adj
    return x, y

if __name__ == '__main__':
    print(to_rect(471700, 5799460))     # r=4.71700  theta=57.99460
    print(to_polar(250000, 400000))     # x=2.50000  y=4.00000

Discussion

Impressive algorithmic technology from yesteryear. Computes cosine, sine, tangent, hypoteneuse, and arctangent using only simple scaled integer arithmetic (a few constants, bit shifts, sign flips, and integer adds in the inner loop and one scaled-multiply in the outer loop).

Here's how to extract the components from the coordination transformations:

cos, sin = to_rect(100000, angle)
tan = sin * 100000 // cos
hypot, arctan = to_polar(x, y)

It's easy to make new constants for other scaling factors, for higher levels of precision, or for using radians instead of degrees:

>>> from operator import mul
>>> from math import degrees, atan

>>> scale = 10**5       # can by any large number including powers of two
>>> n = 8               # more iterations give more accuracy
>>> round(scale * reduce(mul, [1/(1+4**-i)**0.5 for i in range(n)]))
60726
>>> [round(scale * degrees(atan(2**-i))) for i in range(n)]
[4500000, 2656505, 1403624, 712502, 357633, 178991, 89517, 44761]

Sign in to comment