```In article <2479@orca.TEK.COM> brucec@orca.UUCP (Bruce Cohen) writes: >In article <6019@iuvax.UUCP> viking@iuvax writes: >> >>The major enhancements that Jim added to the program were to recode it in >>C and to use "a CORDIC rotator instead of trig functions". I want to talk >>with Jim about his program....where is he?! >> >>Alternately, does anyone know anything about "CORDIC rotators"? I'd like to >>see the source to his modification of Michiel's program. The algorithm used >>was based on one developed by LucasFilm Ltd. and published in the Sept. '84 >>Scientific American. The idea is really simple. A rotation can be expressed as: [ x' ] [ cos t -sin y ] [ x ] [ y' ] = [ sin t cos t ] [ y ] Factor out cos t, giving: [ x' ] [ 1 -tan y ] [ x ] [ y' ] = cos t [ tan t 1 ] [ y ] Now let t be expressed as a signed sum of -1 -i t = Sum { d tan 2 }, where d = {+1, -1} i i i This yields (please forgive the cruddy typesetting; Mathtype doesn't work on UNIX): -i [ x' ] [ 1 -d 2 ] [ x ] [ ] [ i ] [ ] [ ] = Sum { C } Sum { [ -i ] [ ] } [ y' ] i i i [ d 2 1 ] [ y ] i The first sum is a constant. The second sum is a series of shifts and adds (or subtracts). This is illustrated by the following working code, which performs all of the following operations: rotate a vector polar to rectangular conversion rectangular to polar conversion simultaneous calculation of sine and cosine. atan2() Note that the normalization and denormalization is done mainly to increase accuracy; if you're more concerned about speed, you'd probably hack this out entirely or just shift by a fixed amount. You can do similar sorts of things with hyperbolic functions, thereby being able to do exponential and logarithms. I have never written such code, although I would like to have it, especially to calculate gamma correction tables. ------------------------------------------------------------ # define RADIANS 1 # define COSCALE 0x22c2dd1c /* 0.271572 */ static long arctantab = { # ifdef DEGREES /* MS 10 integral bits */ # define QUARTER (90 << 22) 266065460, 188743680, 111421900, 58872272, 29884485, 15000234, 7507429, 3754631, 1877430, 938729, 469366, 234683, 117342, 58671, 29335, 14668, 7334, 3667, 1833, 917, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1, 0, 0, # else # ifdef RADIANS /* MS 4 integral bits */ # define QUARTER ((int)(3.141592654 / 2.0 * (1 << 28))) 297197971, 210828714, 124459457, 65760959, 33381290, 16755422, 8385879, 4193963, 2097109, 1048571, 524287, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 0, 0, # else # define BRADS 1 # define QUARTER (1 << 30) 756808418, 536870912, 316933406, 167458907, 85004756, 42667331, 21354465, 10679838, 5340245, 2670163, 1335087, 667544, 333772, 166886, 83443, 41722, 20861, 10430, 5215, 2608, 1304, 652, 326, 163, 81, 41, 20, 10, 5, 3, 1, 1, # endif # endif }; /* 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. */ static pseudorotate(px, py, theta) long *px, *py; register long theta; { register int i; register long x, y, xtemp; register long *arctanptr; x = *px; y = *py; /* Get angle between -90 and 90 degrees */ while (theta < -QUARTER) { x = -x; y = -y; theta += 2 * QUARTER; } while (theta > QUARTER) { x = -x; y = -y; theta -= 2 * QUARTER; } /* Initial pseudorotation, with left shift */ arctanptr = arctantab; if (theta < 0) { xtemp = x + (y << 1); y = y - (x << 1); x = xtemp; theta += *arctanptr++; } else { xtemp = x - (y << 1); y = y + (x << 1); x = xtemp; theta -= *arctanptr++; } /* Subsequent pseudorotations, with right shifts */ for (i = 0; i <= 28; i++) { if (theta < 0) { xtemp = x + (y >> i); y = y - (x >> i); x = xtemp; theta += *arctanptr++; } else { xtemp = x - (y >> i); y = y + (x >> i); x = xtemp; theta -= *arctanptr++; } } *px = x; *py = y; } static pseudopolarize(argx, argy) long *argx, *argy; { register long theta; register long yi, i; register long x, y; register long *arctanptr; x = *argx; y = *argy; /* Get the vector into the right half plane */ theta = 0; if (x < 0) { x = -x; y = -y; theta = 2 * QUARTER; } if (y > 0) theta = - theta; arctanptr = arctantab; if (y < 0) { /* Rotate positive */ yi = y + (x << 1); x = x - (y << 1); y = yi; theta -= *arctanptr++; /* Subtract angle */ } else { /* Rotate negative */ yi = y - (x << 1); x = x + (y << 1); y = yi; theta += *arctanptr++; /* Add angle */ } for (i = 0; i <= 28; i++) { if (y < 0) { /* Rotate positive */ yi = y + (x >> i); x = x - (y >> i); y = yi; theta -= *arctanptr++; } else { /* Rotate negative */ yi = y - (x >> i); x = x + (y >> i); y = yi; theta += *arctanptr++; } } *argx = x; *argy = theta; } /* Fxprenorm() block normalizes the arguments to a magnitude suitable for * CORDIC pseudorotations. The returned value is the block exponent. */ static int fxprenorm(argx, argy) long *argx, *argy; { register long x, y; register int shiftexp; int signx, signy; shiftexp = 0; /* Block normalization exponent */ signx = signy = 1; if ((x = *argx) < 0) { x = -x; signx = -signx; } if ((y = *argy) < 0) { y = -y; signy = -signy; } /* Prenormalize vector for maximum precision */ if (x < y) { /* |y| > |x| */ while (y < (1 << 27)) { x <<= 1; y <<= 1; shiftexp--; } while (y > (1 << 28)) { x >>= 1; y >>= 1; shiftexp++; } } else { /* |x| > |y| */ while (x < (1 << 27)) { x <<= 1; y <<= 1; shiftexp--; } while (x > (1 << 28)) { x >>= 1; y >>= 1; shiftexp++; } } *argx = (signx < 0) ? -x : x; *argy = (signy < 0) ? -y : y; return(shiftexp); } /* Return a unit vector corresponding to theta. * sin and cos are fixed-point fractions. */ fxunitvec(cos, sin, theta) long *cos, *sin, theta; { *cos = COSCALE >> 1; /* Save a place for the integer bit */ *sin = 0; pseudorotate(cos, sin, theta); /* Saturate results to fractions less than 1, shift out integer bit */ if (*cos >= (1 << 30)) *cos = 0x7FFFFFFF; /* Just shy of 1 */ else if (*cos <= -(1 << 30)) *cos = 0x80000001; /* Just shy of -1*/ else *cos <<= 1; if (*sin >= (1 << 30)) *sin = 0x7FFFFFFF; /* Just shy of 1 */ else if (*sin <= -(1 << 30)) *sin = 0x80000001; /* Just shy of -1*/ else *sin <<= 1; } /* Fxrotate() rotates the vector (argx, argy) by the angle theta. */ fxrotate(argx, argy, theta) long *argx, *argy; long theta; { register long x, y; int shiftexp; long frmul(); if (((*argx || *argy) == 0) || (theta == 0)) return; shiftexp = fxprenorm(argx, argy); /* Prenormalize vector */ pseudorotate(argx, argy, theta); /* Perform CORDIC pseudorotation */ x = frmul(*argx, COSCALE); /* Compensate for CORDIC enlargement */ y = frmul(*argy, COSCALE); if (shiftexp < 0) { /* Denormalize vector */ *argx = x >> -shiftexp; *argy = y >> -shiftexp; } else { *argx = x << shiftexp; *argy = y << shiftexp; } } fxatan2(x, y) long x, y; { if ((x || y) == 0) return(0); fxprenorm(&x, &y); /* Prenormalize vector for maximum precision */ pseudopolarize(&x, &y); /* Convert to polar coordinates */ return(y); } fxpolarize(argx, argy) long *argx, *argy; { int shiftexp; long frmul(); if ((*argx || *argy) == 0) { *argx = 0; /* Radius */ *argy = 0; /* Angle */ return; } /* Prenormalize vector for maximum precision */ shiftexp = fxprenorm(argx, argy); /* Perform CORDIC conversion to polar coordinates */ pseudopolarize(argx, argy); /* Scale radius to undo pseudorotation enlargement factor */ *argx = frmul(*argx, COSCALE); /* Denormalize radius */ *argx = (shiftexp < 0) ? (*argx >> -shiftexp) : (*argx << shiftexp); } # ifdef vax long frmul(a, b) /* Just for testing */ long a, b; { /* This routine should be written in assembler to calculate the * high part of the product, i.e. the product when the operands * are considered fractions. */ return((a >> 16) * (b >> 15)); } # endif ##### This is the end of my posting ##### -- Ken Turkowski @ Apple Computer, Inc., Cupertino, CA UUCP: {mtxinu,sun,decwrl,nsc,voder}!apple!turk CSNET: turk@Apple.CSNET ARPA: turk%Apple@csnet-relay.ARPA ``` Discuss this article in the forums Date this article was posted to GameDev.net: 7/16/1999 (Note that this date does not necessarily correspond to the date the article was written) See Also: Miscellaneous © 1999-2011 Gamedev.net. All rights reserved. Terms of Use Privacy Policy Comments? Questions? Feedback? Click here!