home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsi / fixedtrig.c < prev    next >
Text File  |  1992-09-15  |  2KB  |  74 lines
  1. /* 
  2. Fixed-Point Trigonometry with CORDIC Iterations
  3. by Ken Turkowski
  4. from "Graphics Gems", Academic Press, 1990
  5.  
  6. provided by user:
  7.     frmul(a,b)=(a*b)>>31, high part of 64-bit product
  8. */
  9.  
  10. #define COSCALE 0x22c2dd1c /* 0.271572 */
  11. #define QUARTER ((int)(3.141592654 / 2.0 * (1 << 28)))
  12. static long arctantab[32] = {  /* MS 4 integral bits for radians */
  13.     297197971, 210828714, 124459457, 65760959, 33381290, 16755422,
  14.     8385879, 4193963, 2097109, 1048571, 524287, 262144, 131072,
  15.     65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64,
  16.     32, 16, 8, 4, 2, 1, 0, 0,
  17. };
  18.  
  19. CordicRotate(px, py, theta)
  20. long *px, *py;
  21. register long theta;    /* Assume that abs(theta) <= pi */
  22. {
  23.     register int i;
  24.     register long x = *px, y = *py, xtemp;
  25.     register long *arctanptr = arctantab;
  26.  
  27.     /* The -1 may need to be pulled out and done as a left shift */
  28.     for (i = -1; i <= 28; i++) {
  29.         if (theta < 0) {
  30.             xtemp = x + (y >> i);
  31.             y     = y - (x >> i);
  32.             x = xtemp;
  33.             theta += *arctanptr++;
  34.         } else {
  35.             xtemp = x - (y >> i);    
  36.             y     = y + (x >> i);
  37.             x = xtemp;    
  38.             theta -= *arctanptr++;
  39.         }
  40.     }
  41.  
  42.     *px = frmul(x, COSCALE); /* Compensate for CORDIC enlargement */
  43.     *py = frmul(y, COSCALE); /* frmul(a,b)=(a*b)>>31, high part  */
  44.                  /* of 64-bit product */
  45. }
  46.  
  47.  
  48.  
  49.  
  50. CordicPolarize(argx, argy)
  51. long *argx, *argy;    /* We assume these are already in the */
  52.                     /*  right half plane */
  53. {
  54.     register long theta = 0, yi, i;
  55.     register long x = *argx, y = *argy;
  56.     register long *arctanptr = arctantab;
  57.     for (i = -1; i <= 28; i++) {
  58.         if (y < 0) {        /* Rotate positive */
  59.             yi = y + (x >> i);
  60.             x  = x - (y >> i);
  61.             y  = yi;
  62.             theta -= *arctanptr++;
  63.         } else {         /* Rotate negative */
  64.             yi = y - (x >> i);
  65.             x  = x + (y >> i);
  66.             y  = yi;
  67.             theta += *arctanptr++;
  68.         }
  69.     }
  70.  
  71.     *argx = frmul(x, COSCALE);
  72.     *argy = theta;
  73. }
  74.