brucec@orca.TEK.COM (Bruce Cohen) (02/11/88)
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. CORDIC (stands for COordinate Rotation DIgital Computer) algorithms are a class of algorithms for computing approximations to transcendental functions quickly. They typically converge towards the solution at a rate of one bit of precision per iteration of the loop. The trig functions, and maybe some of the other transcendentals, on HP calculators are implemented with CORDICS, for instance. The technique has been around since at least the '50s (1959 was the earliest citation I found in about two minutes of search, it was claimed to be the original paper), so I doubt that that part of the algorithm is what LucasFilm claims as theirs. The basis of the algorithm is a fast loop which rotates a vector using only shifts and adds. There's a nice description of how it works in an article by Ken Turkowski in the July, 1982 issue of ACM Transctions on Graphics. He even shows the kernel of the algorithm in C. I don't have time right now to write a dissertation on the algorithm. If you have access to Ken's article, read it. If Ken is out there listening as he used to be, he might be willing to write something up now, or let me paraphrase some of his article in a day or two, when I have time. Umm ... now that I think about it, someone (possibly Ken) sent out C sources for a set of CORDIC functions a couple of years ago. Anyone have those squirreled away? --------------------------------------------------------------------------- "The authorities are excellent at amassing facts, though they do not always use them to advantage." Sherlock Holmes, "The Naval Treaty" My opinions are my own; no-one else seems to want them. Bruce Cohen bang-syntax: {the real world}...!tektronix!ruby!brucec at-syntax: brucec@ruby.TEK.COM overland: Tektronix Inc., M/S 61-028, P.O. Box 1000, Wilsonville, OR 97070
jimc@iscuva.ISCS.COM (Jim Cathey) (02/12/88)
In article <2479@orca.TEK.COM> brucec@orca.UUCP (Bruce Cohen) writes: >In article <6019@iuvax.UUCP> viking@iuvax writes: >>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. >CORDIC (stands for COordinate Rotation DIgital Computer) algorithms are a >... The technique has been around >since at least the '50s (1959 was the earliest citation I found in about >two minutes of search, it was claimed to be the original paper), so I >doubt that that part of the algorithm is what LucasFilm claims as theirs. The LucasFilm reference is for the part of the Fractals program that generates the mountainous landscapes, not the CORDIC routine. The SA article was the 'original' source for that information. My biggest improvement to de Panne's program was to use CORDIC rather than the Rect->Polar+theta->Rect method he used to flip the database up to a nice angle for viewing. CORDIC has always been on my list of 'elegant' algorithms. +----------------+ ! II CCCCCC ! Jim Cathey ! II SSSSCC ! ISC Systems Corp. ! II CC ! TAF-C8; Spokane, WA 99220 ! IISSSS CC ! UUCP: uunet!iscuva!jimc ! II CCCCCC ! (509) 927-5757 +----------------+ "With excitement like this, who is needing enemas?"
olson@endor.harvard.edu (Eric K. Olson) (02/12/88)
I believe Ken Turkowski now works for the Advanced Technology Group at Apple Computer. I know he has an account on husc6!sri-unix!garth!apple, a/k/a apple.CSNET, a/k/a 192.5.58.12, but I was unable to finger him there. Best bet is to call Apple's main number at 408-996-1010. -Eric Eric K. Olson olson@endor.harvard.edu harvard!endor!olson D0760 (Name) (ArpaNet) (UseNet) (AppleLink)
root@cca.ucsf.edu (Computer Center) (02/12/88)
In article <2479@orca.TEK.COM>, brucec@orca.TEK.COM (Bruce Cohen) writes: ... > Umm ... now that I think about it, someone (possibly Ken) sent out C sources > for a set of CORDIC functions a couple of years ago. Anyone have those > squirreled away? > There was a posting in volume 1 of comp.sources.misc 8706/5 (fifth posting in June 1987 if I understand the notation). Thos Sumner (thos@cca.ucsf.edu) BITNET: thos@ucsfcca (The I.G.) (...ucbvax!ucsfcgl!cca.ucsf!thos) OS|2 -- an Operating System for puppets. #include <disclaimer.std>
turk@apple.UUCP (Ken "Turk" Turkowski) (02/18/88)
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[32] = { # 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