Computer, Digitally Rotate the Coordinates!

Computer, digitally rotate the coordinates! I can almost hear Captain Jean-Luc Picard delivering this command from the captain’s chair aboard the Enterprise. If this reference doesn’t ring a bell, then congratulations, you’re young! There’s a generational gap between us. But while our pop culture touchstones may differ, we’re united by something far older than either of us: the timeless beauty of mathematics.

The concept of numbers laid the foundation for algebra. Performing computations with Roman numerals was notoriously difficult, but the decimal number system, already in use in India between the 5th and 7th centuries CE, made calculations far more efficient. This system eventually spread to the Islamic world, and by the 9th century, Al-Khwarizmi formalized algebra, introducing systematic methods for solving equations. The word algorithm itself originates from the Latinized mispronunciation of his name by European mathematicians.

Algebra enabled the structured study of functions, a concept further advanced by René Descartes and Leonhard Euler in the 17th and 18th centuries. The idea of power series expansions laid the groundwork for calculus, which was formally developed by Isaac Newton and Gottfried Wilhelm Leibniz in the late 17th century. However, over 300 years earlier, in the 14th century, the Indian mathematician Madhava of Sangamagrama had already derived power series expansions for sine and cosine functions. The ingenious method Madhava used to derive these expansions deserves an article of its own.

The cumulative effect of these mathematical advancements over the centuries has led to exponential growth in the field, particularly in recent history. In fact, most of us can only grasp the mathematics developed until the last century — unless, of course, we are trained, practicing mathematicians.

With the invention of computers in the 20th century, one important problem was efficiently computing the values of fundamental functions such as the exponential, logarithm, sine, cosine, tangent, and arctangent. While Taylor series expansions provide a way to calculate these functions, they require multiple high-precision multiplications and divisions, which were computationally expensive on early digital computers. Additionally, different functions require separately coded power series expansions.

To address this problem, the CORDIC (COordinate Rotation DIgital Computer) algorithm was developed around 1960 by Jack Volder. It provides an iterative method to compute these functions efficiently using only shifts and additions, as it operates with powers of 2 as factors and this makes it particularly well-suited for digital hardware.

As the name suggests, what the algorithm fundamentally does is performing rotations. In 2D, rotations are represented by multiplication with exp(jθ). Given a complex number (x+jy), multiplying by exp(jθ) and restricting θ to only inverse tangents of negative powers of 2 results in the following set of equations:

The original point (x,y) is on some circle whose radius is given by x²+y². We can show that the recursive equaions x[n+1] = x[n]-y[n]/2^n and y[n+1]=y[n]+x[n]/2^n are also on a scaled circle (see image below). The scaling at every iteration is (1+1/2^(2n)). If we adjust for this scaling at every iteration, then (x[n],y[n]) at every iteration n will lie on the same circle as the original point.

This factor (1+1/2^(2n)) is called the CORDIC gain and the cumulative gain across all iterations can be compensated for at the end by just one multiplication. The factor Kn (see image below) to compensate for the CORDIC gain An converges to 0.6073 as n increases.

Any θ can be written as ∑d[n]θ[n] where d[n] is {-1,1} and θ[n] are a finite set of angles that form a geometrically decreasing sequence. This is similar to how any real number can be approximated using a binary expansion. In this case, we are summing angles instead of power of 2. We then have the following set of equations for the CORDIC iteration.

While the formualtion itself is generic and should work for any angle, in practice, the input angle is restricted to the first quadrant and the result for angles in other quadrants can be easily derived, either via rotations or via reflections.

To rotate any (x,y) by an angle θ, we initialize θ[0]=θ and iterate until θ[n] approaches 0. Each iteration rotates (x[n],y[n]) by exp(jθ[n]). At every iteration, d[n] is chosen to be -1 of θ[n] < 0 or 1 otherwise. The cosine and sine of the angle θ can be computed by rotating (1,0) by θ.

If we instead choose d[n] to get y[n] to approach 0, we are effectively rotating the (x,y) vector to the x-axis. The accumulated angle θ[n] should then give us the angle of (x,y) which is the arctangent(y/x). Once the vector (x,y) is rotated to x-axis, the value x[n] has the magnitude of the vector (x,y). This mode thus converts from the cartesian representation (x,y) to the polar representation.

For arcsine, we need to find the angle at which (1,0) after rotation would give us the y-component that equals the given y. We then choose d[n] = 1 if y[n] < y, -1 otherwise. But since x[n] and y[n] are scaled at every iteration due to the CORDIC gain, after n iterations, we get θ[n]=arcsin(y/An).

The unit circle x²+y² = 1 is parameterized by x = cos(θ) and y = sin(θ). Given our familiarity with the Euler equation exp(jθ) = cos(θ) + j sin(θ), we have cos(θ)=(exp(jθ)+exp(-jθ))/2 and sin(θ)=(exp(jθ)-exp(-jθ))/2j. The equation of the hyperbola is x²-y² = 1. What could be the parameterization for the hyperbola? If we try x(t) = (exp(t)+exp(-t))/2 and y(t) = (exp(t)-exp(-t))/2, we see that these satisfy x(t)²-y(t)² = 1.

We define cosh(t) = (exp(t)+exp(-t))/2, sinh(t) = (exp(t)-exp(-t))/2 and tanh(t) = sinh(t)/cosh(t). This allows us to replace the trignometric functions in the CORDIC iteration with the hyperbolic equivalents. Note that the negative sign inside the matrix is removed in the hyperbolic case (see image above). The hyperbolic transformation can be interpreted as a diagonal scaling operation in the Hadamard transform basis, as shown by the equation below. If the Hadamard transform is used to rotate the co-ordinate axis, we then have u = (x+y)/2 and v = (x-y)/2 and the hyperbola equation then becomes uv = 1. This form of hyperbola is parameterized by u=exp(t) and v=exp(-t).

Why are we doing this again with hyperbola equation? Because we can then use the CORDIC iterations to compute cosh(t) and sinh(t) using the rotation mode, which then allows us to get exp(t) as cosh(t) + sinh(t). As with the circle case, we see that the points x[n], y[n] continue to lie on a hyperbola.

If we instead use the arctangent mode, we can get the hyperbolic arctangent and with y=α-1 and x=α+1, that gives us the logarithm, as shown by the equations below.

Just like the arctangent also providing the magnitude of the vector (x,y), the hyperbolic arctangent computation too provides the magnitude of the semi-transverse axis of the hyperbola. This allows us to compute square root, as shown by the equations below.

Instead of using invese tangent of 1/2^n or inverse hyperbolic tangent of 1/2^n to update θ[n] or z[n], if we simply use 1/2^n, then we have the linear mode and rotation in this case becomes the bitwise shift and add multiplication. If we operate in arctangent mode (also called the vectoring mode), we get division via bitwise shift and subtract.

The following table lists all the CORDIC modes and the functions. For an excellent summary of CORDIC, take a look at Ray Andraka, “A survey of CORDIC algorithms for FPGA based computers.”

While the hyperbolic mode indirectly computes the natual logarithm, we can come up with a much more straightforward iteration by realizing that any number x can be written as a product of smaller numbers: x = ∏a[n]. If we choose a[n] = 1+d[n]/2^n, we can find a sequence of d[n] = {-1,0,1} such that yx = 1. We therefore have ln(y) = -ln(x) = -ln(∏a[n])= ∑ln(a[n]). This means, if we find the d[n] such that x = ∏(1+d[n]/2^n) results in xy becoming 1, we have the ln(y) from this product.

Similarly, if can write exp(x) as exp(∑d[n]/2^n) and this is same as ∏exp(d[n]/2^n). For small values of exponents, exp(d[n]/2^n) ≈ 1+d[n]/2^n and this gives us exp(x) ≈ ∏(1+d[n]/2^n). This mean x = ln(∏(1+d[n]/2^n)) = ∑ln(1+d[n]/2^n).

The following image shows these equations more clearly. If we choose d[n] such that we get x ≈ ∑ln(1+d[n]/2^n), this means ∏(1+d[n]/2^n) ≈ exp(x). And if we choose d[n] to get x ≈ ∏(1+d[n]/2^n), then ∑ln(1+d[n]/2^n) ≈ ln(x).

We can now define the following iterations of the BKM algorithm (named after its inventors Bajard, Kla and Muller):

If we drive En to 1, meaning we select d[n] such that ∏(1+d[n]/2^n) ≈ 1/E0, then with L0=0, we get Ln = ln(E0). This is called the L-mode. At every iteration n, if En > 1, we choose d[n] = -1 to decrease E[n+1], otherwise choose d[n] = 1.

If we drive Ln to 0, meaning we select d[n] such that ∑ln(1+d[n]/2^n) ≈ L0, then with E0=1, we get En = exp(L0). This is called the E-mode. At every iteration n, we choose d[n]={-1,0,1} such that |Ln-ln(1+d[n]/2^n)| is minimized.

With exp(x) = exp(∑d[n]/2^n), there is nothing to stop us from having a complex exponent and we would get exp(x+jy) = exp(∑(dr[n]+jdi[n])/2^n). This approximation exp(x+jy) ≈ ∏(1+(dr[n]+jdi[n])/2^n) is valid too. This means x+jy ≈ ∑ln(1+(dr[n]+jdi[n])/2^n). But what does logarithm of a complex number mean? We have r exp(jθ) as our complex number, then its logarithm is ln(r)+j(θ+2πn) for any integer n. Note that the same complex number has many θ+2πn values the logarithm maps it to. So the complex logarithm is multivalued function that is a one to many map.

We can extend the BKM algorithm for complex numbers by allowing d[n]=dr[n]+jdi[n] with both dr[n] and di[n] taking values {-1,0,1}. Possible values for d[n] then are {−1−j,−1,−1+j,−j,0,j,1−j,1,1+j}.

In E-mode, we can have L0 as complex and drive Ln to 0 to get En = E0·exp(L0) = E0·exp(L0r)·exp(jL0i). This allows us to do both exponential scaling and rotation at the same time!

Similarly, in L-mode, we can have E0 as complex and drive En to 0 to get Ln=ln(E0) = ln(|E0|) + j∡E0. Thus we can compute the logarithm and arctangent at the same time!

And that brings us to the end. Hopefully, this article provides you enough insights to enable you, dear reader, to confidently make your computer digitally rotate the coordinates!