ken (04/20/83)
The CORDIC algorithms are superb for implementing trigonometric,
inverse trigonometric, and hyperbolic functions.
rotate() rotates a rectangular vector, but can be used to convert from
polar to rectangular coordinates. It is given rectangular coordinates
and an angle, and it returns rectangular coordinates.
polarize() converts from rectangular to polar coordinates. It is given
rectangular coordinates (x, y) and returns polar coordinates (r, theta).
polarize() is more of an arctan2(x, y) than an arctan(x/y), but gives
an angle around the whole circle, rather than just the right-half
plane.
BRADS is my name for binary radians, the natural measure of angle on a
16-bit machine: 32768 is a half circle, and 65536 is a whole circle.
You get your choice of degrees, radians, or brads, or whatever other
measure you want. The difference is in the arctangent table.
The scheme is as follows: Rotate a vector through successively smaller
angles until it meets a certain condition. For polarize(), it is until
the Y-coordinate is zero. For rotate(), it is until the rotation angle
is zero. The papers prove convergence. Approximately one bit of
accuracy is got for each iteration through the loop. If you aren't
interested in the magnitude of the vector, you can eliminate the final
multiplication by COSCALE to save time.
These were implemented on the PDP-11, and hence assume that an int is
16 bits.
References:
Chen, T. C.
"Automatic Computation of Exponentials, Logarithms, Ratios, and Square Roots",
IBM J. Res. Develop., July 1972, pp. 380-388
Despain, A. M.
"Fourier Transform Computers Using CORDIC Iterations",
IEEE Trans. Comp., vol. C-23, no. 10, Oct. 1974, pp. 993-1001
Turkowski, K. E.
"Anti-Aliasing through the Use of Coordinate Transformations",
ACM Transactions on Graphics,
Vol. 1, No. 3, July 1982, pp. 215-234
Volder, J. E.
"The CORDIC Trigonometric Computing Technique",
IRE Transactions on Electronic Computers,
vol. EC-8, no. 3, Sept. 1959, pp. 330-334
Walther, J. S.
"A Unified Algorithm for Elementary Functions",
Proc. AFIPS 1971 Spring Joint Computer Conference,
pp. 379-385
/********************************************************/
cat << 'EOF' > polarize.c
# define COSCALE 042606 /* 0.271572 */
# ifdef DEGREES
# define QUARTER (90 << 6) /* Before shifting left by two */
static int arctantab[16] = {
16239,11520,6801,3593,1824,916,458,229,115,57,29,14,7,4,2,1
};
# else
# ifdef RADIANS
# define QUARTER ((int)(3.141592654 / 2.0 * 4096.0)) /* Before shifting */
static int arctantab[16] = {
18140,12868,7596,4014,2037,1023,512,256,128,64,32,16,8,4,2,1
};
# else
# define BRADS BRADS
# define QUARTER 16384 /* Before shifting by 1 */
static int arctantab[16] = {
23096,16384,9672,5110,2594,1302,652,326,163,81,41,20,10,5,3,1
};
# endif
# endif
struct rectangular {
int x;
int y;
};
struct polar {
int magnitude;
int angle;
};
/* This program changes from rectangular to polar coordinates by rotating
* the vector several times until it lies within epsilon of 0 degrees. A
* fixed magnitude correction is done at the end, as is a rounding of the
* accumulated angle.
*
* To rotate a vector through an angle of theta, we calculate:
*
* x' = x cos(theta) - y sin(theta)
* y' = x sin(theta) + y cos(theta)
*
* The rotate() routine performs multiple rotations of the form:
*
* x[i+1] = cos(theta[i]) { x[i] - y[i] tan(theta[i]) }
* y[i+1] = cos(theta[i]) { y[i] + x[i] tan(theta[i]) }
*
* with the constraint that tan(theta[i]) = pow(2, -i), which can be implemented
* by shifting. We always shift by either a positive or negative angle, so the
* convergence has the ringing property. Since the cosine is always positive
* for positive and negative angles between -90 and 90 degrees, a predictable
* magnitude scaling occurs at each step, and can be compensated for instead
* at the end of the iterations by a composite scaling of the product of all
* the cos(theta[i])'s.
*
* Several modifications to the technically elegant program were made to
* accomodate other desirable properties.
*
* We have accommodated any of the various ways of representing an angle:
* 0 to 360 degrees, -180 to 180 degrees, or even possibly -360 to 0 degrees
* (correspondingly with radians). After the angle has been brought to the
* -90 to 90 degree range, the range is expanded to have the maximum accuracy
* in the forthcoming calculations. For degrees, this means that the input
* angle has 6 bits of fractional accuracy, and is increased to 8 bits for
* calculation purposes. For radians, we have 12 bits of fractional accuracy
* on input, and 14 bits during calculations. For brads (the natural way of
* representing angles on a binary computer), we have 16 bits to represent the
* full 360 degrees on input, and 16 bits to represent +/- 90 degrees during
* calculations.
*
* The other modification made was to accommodate any radius we might encounter
* in terms of sixteenths of pixels. There should be enough range to accommodate
* any vectors that are legal, even in 4xPAL.
*
* An empirical analysis shows that the error is about 1/2 LSB until one feeds
* in a vector with a magnitude of 6464 (404 pixels in terms of sixteenths);
* after that, the least significant bit sometimes flips up or down. This gives
* an accuracy of approximately 0.015%.
*/
struct polar
polarize(vector)
struct rectangular vector;
{
register int i, theta;
long x0, y0, x1;
struct polar tempolar;
struct doubleword {
int hiword;
int loword;
};
x0.loword = 0;
x0.hiword = vector.x;
y0.loword = 0;
y0.hiword = vector.y;
/* Get angle between -90 and 90 degrees */
if (x0.hiword < 0) {
x0.hiword = -x0.hiword;
if ((y0.hiword = -y0.hiword) <= 0)
tempolar.angle = 2 * QUARTER;
else
tempolar.angle = - 2 * QUARTER;
}
else
tempolar.angle = 0;
theta = 0;
/* Prescale vectors to allow for full dynamic range in sixteenths */
x0 >>= 1;
y0 >>= 1;
for (i = 1; i > -15; --i) {
if (y0 < 0) {
x1 = x0 - (y0 << i);
y0 = y0 + (x0 << i);
x0 = x1;
theta -= arctantab[1 - i];
}
else {
x1 = x0 + (y0 << i);
y0 = y0 - (x0 << i);
x0 = x1;
theta += arctantab[1 - i];
}
}
x0 += 0100000L; /* Round x0 */
/* Scale the rotated vectors to keep the magnitude the same,
* compensating for the prescaling. */
tempolar.magnitude = frmul(x0.hiword, COSCALE);
# ifdef BRADS
tempolar.angle += (theta + 1) >> 1;
# else
tempolar.angle += (theta + 2) >> 2;
# endif
return(tempolar);
}
'EOF'
cat << 'EOF' > rotate.c
# define COSCALE 042606 /* 0.271572 */
# ifdef DEGREES
# define QUARTER (90 << 6) /* Before shifting left by two */
static int arctantab[16] = {
16239,11520,6801,3593,1824,916,458,229,115,57,29,14,7,4,2,1
};
# else
# ifdef RADIANS
# define QUARTER ((int)(3.141592654 / 2.0 * 4096.0)) /* Before shifting */
static int arctantab[16] = {
18140,12868,7596,4014,2037,1023,512,256,128,64,32,16,8,4,2,1
};
# else
# define BRADS BRADS
# define QUARTER 16384 /* Before shifting by 1 */
static int arctantab[16] = {
23096,16384,9672,5110,2594,1302,652,326,163,81,41,20,10,5,3,1
};
# endif
# endif
struct coordinate {
int x;
int y;
};
/* To rotate a vector through an angle of theta, we calculate:
*
* x' = x cos(theta) - y sin(theta)
* y' = x sin(theta) + y cos(theta)
*
* The rotate() routine performs multiple rotations of the form:
*
* x[i+1] = cos(theta[i]) { x[i] - y[i] tan(theta[i]) }
* y[i+1] = cos(theta[i]) { y[i] + x[i] tan(theta[i]) }
*
* with the constraint that tan(theta[i]) = pow(2, -i), which can be implemented
* by shifting. We always shift by either a positive or negative angle, so the
* convergence has the ringing property. Since the cosine is always positive
* for positive and negative angles between -90 and 90 degrees, a predictable
* magnitude scaling occurs at each step, and can be compensated for instead
* at the end of the iterations by a composite scaling of the product of all
* the cos(theta[i])'s.
*
* Several modifications to the technically elegant program were made to
* accomodate other desirable properties.
*
* We have accommodated any of the various ways of representing an angle:
* 0 to 360 degrees, -180 to 180 degrees, or even possibly -360 to 0 degrees
* (correspondingly with radians). After the angle has been brought to the
* -90 to 90 degree range, the range is expanded to have the maximum accuracy
* in the forthcoming calculations. For degrees, this means that the input
* angle has 6 bits of fractional accuracy, and is increased to 8 bits for
* calculation purposes. For radians, we have 12 bits of fractional accuracy
* on input, and 14 bits during calculations. For brads (the natural way of
* representing angles on a binary computer), we have 16 bits to represent the
* full 360 degrees on input, and 16 bits to represent +/- 90 degrees during
* calculations.
*
* The other modification made was to accommodate any radius we might encounter
* in terms of sixteenths of pixels. There should be enough range to accommodate
* any vectors that are legal, even in 4xPAL.
*
* An empirical analysis shows that the error is about 1/2 LSB until one feeds
* in a vector with a magnitude of 6464 (404 pixels in terms of sixteenths);
* after that, the least significant bit sometimes flips up or down. This gives
* an accuracy of approximately 0.015%.
*/
struct coordinate
rotate(vector, theta)
struct coordinate vector;
register int theta;
{
register int i;
long x0, y0, x1;
struct doubleword {
int hiword;
int loword;
};
x0.loword = 0;
x0.hiword = vector.x;
y0.loword = 0;
y0.hiword = vector.y;
/* Get angle between -90 and 90 degrees */
# ifdef BRADS
if (theta < -QUARTER) {
# else
while (theta < -QUARTER) {
# endif
x0.hiword = -x0.hiword;
y0.hiword = -y0.hiword;
theta += 2 * QUARTER;
}
# ifdef BRADS
if (theta > QUARTER) {
# else
while (theta > QUARTER) {
# endif
x0.hiword = -x0.hiword;
y0.hiword = -y0.hiword;
theta -= 2 * QUARTER;
}
# ifdef BRADS
/* Expand accuracy of angle specification. */
theta <<= 1;
# else
theta <<= 2;
# endif
/* Prescale vectors to allow for full dynamic range in sixteenths */
x0 >>= 1;
y0 >>= 1;
# ifdef DEBUG
printf("starting: x = %d, y = %d, theta = %d\n", x0.hiword,
y0.hiword, theta);
# endif
for (i = 1; i > -15; --i) {
if (theta < 0) {
x1 = x0 + (y0 << i);
y0 = y0 - (x0 << i);
x0 = x1;
theta += arctantab[1 - i];
}
else {
x1 = x0 - (y0 << i);
y0 = y0 + (x0 << i);
x0 = x1;
theta -= arctantab[1 - i];
}
# ifdef DEBUG
printf("x = %d, y = %d, theta = %d\n", x0.hiword, y0.hiword, theta);
# endif
}
/* Round x0 and y0 */
x0 += 0100000L;
y0 += 0100000L;
/* Scale the rotated vectors to keep the magnitude the same,
* compensating for the prescaling. */
vector.x = frmul(x0.hiword, COSCALE);
vector.y = frmul(y0.hiword, COSCALE);
return(vector);
}
'EOF'