home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_10_11 / crdtools.c < prev    next >
Text File  |  1992-06-18  |  18KB  |  533 lines
  1. /*  CRDTOOLS.C : Tools to exercise CORDIC calculations for sine/cosine.
  2.  *   Make in Borland C++'s internal environment.
  3.  *   (c) 1991 Michael Bertrand.
  4.  */
  5.  
  6. #include <stdio.h>  /* printf, scanf */
  7. #include <math.h>   /* sqrt, atan, sin, cos, fabs */
  8. #include <conio.h>  /* clrscr, gotoxy, getch, cputs, etc. */
  9. #include <dos.h>    /* struct time, gettime */
  10.  
  11. typedef unsigned int WORD;
  12. typedef struct time TIME;
  13.  
  14. typedef struct
  15. {
  16.   int x;
  17.   int y;
  18. } POINT;
  19.  
  20. typedef struct
  21. {
  22.   double u;
  23.   double v;
  24. } DOUB_PT;
  25.  
  26. char Menu(void);
  27. void SingleAngleDegree(void);
  28. void SingleAngleCordic(void);
  29. void ListSeries(void);
  30. void CalcStatistics(void);
  31. void TimeTest(void);
  32. int  HundSecs(TIME *tp, TIME *sp);
  33. void CalcHexPtsCLib(POINT center, POINT vertex);
  34. void CalcHexPtsCORDIC(POINT center, POINT vertex);
  35. void SinCosSetup(void);
  36. void SinCos(WORD theta, int *sin, int *cos);
  37.  
  38. #define TRUE  1
  39. #define ESC   0x1B
  40.  
  41. /* Quadrants for angles. */
  42. #define QUAD1 1
  43. #define QUAD2 2
  44. #define QUAD3 3
  45. #define QUAD4 4
  46.  
  47. /* NBITS is number of bits used for CORDIC units. */
  48. #define NBITS 14
  49.  
  50. /* NUM_PTS is number of vertices in polygon. */
  51. #define NUM_PTS 6
  52.  
  53. /* Macros to convert angles to different units. */
  54. #define DEG_TO_CORDIC(x) ((WORD) ((double)x/90.0*CordicBase + 0.5))
  55. #define DEG_TO_RADIAN(x) (x/180.0*M_PI)
  56. #define CORDIC_TO_RADIAN(x) ((double)x/CordicBase*M_PI/2)
  57.  
  58. /* NUM_TIMES is number of times to iterate for timing test. */
  59. #define NUM_TIMES 1000
  60.  
  61. int  ArcTan[NBITS];      /* angular constants : arctans */
  62. int  xInit;              /* initial x projection */
  63. WORD CordicBase;         /* base for CORDIC units */
  64. WORD HalfBase;           /* CordicBase / 2 */
  65. WORD Quad2Boundary;      /* 2 * CordicBase */
  66. WORD Quad3Boundary;      /* 3 * CordicBase */
  67. POINT HexPts[NUM_PTS+1]; /* to hold calculated polygon points */
  68.  
  69. void main(void)
  70. {
  71.   SinCosSetup();
  72.  
  73.   do
  74.     {
  75.     switch(Menu())
  76.       {
  77.       case '1' : case 'd':  case 'D':
  78.         SingleAngleDegree();
  79.         break;
  80.       case '2' : case 'c':  case 'C':
  81.         SingleAngleCordic();
  82.         break;
  83.       case '3' : case 'l':  case 'L':
  84.         ListSeries();
  85.         break;
  86.       case '4' : case 's':  case 'S':
  87.         CalcStatistics();
  88.         break;
  89.       case '5' : case 't':  case 'T':
  90.         TimeTest();
  91.         break;
  92.       case '6' : case 'q':  case 'Q':
  93.         clrscr();
  94.         return;
  95.       }  /* switch */
  96.     } while (TRUE);
  97. }
  98.  
  99. char Menu(void)
  100. /*
  101.    USE  : Display menu and get user choice.
  102.    RET  : User choice.
  103. */
  104. {
  105.   struct text_info ti;    /* to get default text color */
  106.  
  107.   gettextinfo(&ti);
  108.   clrscr();
  109.   textattr(LIGHTGRAY);
  110.   gotoxy(26, 1); cputs("╔═════ CORDIC Tests ═════╗");
  111.   gotoxy(26, 2); cputs("║                        ║");
  112.   gotoxy(26, 3); cputs("║  1) One Angle (Degree) ║");
  113.   gotoxy(26, 4); cputs("║                        ║");
  114.   gotoxy(26, 5); cputs("║  2) One Angle (CORDIC) ║");
  115.   gotoxy(26, 6); cputs("║                        ║");
  116.   gotoxy(26, 7); cputs("║  3) List Series        ║");
  117.   gotoxy(26, 8); cputs("║                        ║");
  118.   gotoxy(26, 9); cputs("║  4) Statistics         ║");
  119.   gotoxy(26,10); cputs("║                        ║");
  120.   gotoxy(26,11); cputs("║  5) Time Test          ║");
  121.   gotoxy(26,12); cputs("║                        ║");
  122.   gotoxy(26,13); cputs("║  6) Quit               ║");
  123.   gotoxy(26,14); cputs("║                        ║");
  124.   gotoxy(26,15); cputs("║  Enter choice :        ║");
  125.   gotoxy(26,16); cputs("╚════════════════════════╝");
  126.  
  127.   textattr(RED);
  128.   gotoxy(43, 3); putch('D');
  129.   gotoxy(43, 5); putch('C');
  130.   gotoxy(32, 7); putch('L');
  131.   gotoxy(32, 9); putch('S');
  132.   gotoxy(32,11); putch('T');
  133.   gotoxy(32,13); putch('Q');
  134.  
  135.   gotoxy(46, 15);           /* put cursor inside box */
  136.   textattr(ti.attribute);   /* restore default text color */
  137.   return(getch());          /* get char from user, return to main() */
  138. }
  139.  
  140. void SingleAngleDegree(void)
  141. /*
  142.   USE : Calculate sin/cos for single angle, in degrees.
  143.   NOTE: Also compares with actual values from C library functions.
  144. */
  145. {
  146.   WORD   theta;      /* user-entered angle, in degrees */
  147.   int    sinTheta;   /* calculated sin in CORDIC units */
  148.   int    cosTheta;   /* calculated cos in CORDIC units */
  149.   double dSin;       /* error in CORDIC calculation */
  150.   double dCos;       /* error in CORDIC calculation */
  151.  
  152.   gotoxy(25, 20);
  153.   printf("Enter angle, in degrees : ");
  154.   scanf("%u", &theta);
  155.  
  156.   SinCos(DEG_TO_CORDIC(theta), &sinTheta, &cosTheta);
  157.   dSin = sin(DEG_TO_RADIAN(theta)) - (double) sinTheta / CordicBase;
  158.   dCos = cos(DEG_TO_RADIAN(theta)) - (double) cosTheta / CordicBase;
  159.   gotoxy(18, 22);
  160.   printf("%u deg : dSin = %8.5f   dCos = %8.5f", theta, dSin, dCos);
  161.   gotoxy(24, 25);
  162.   printf("Press a key to continue ......");
  163.   getch();
  164. }
  165.  
  166. void SingleAngleCordic(void)
  167. /*
  168.   USE : Calculate sin/cos for single angle, in CORDIC units.
  169.   NOTE: Also compares with actual values from C library functions.
  170. */
  171. {
  172.   WORD   theta;      /* user-entered angle, in CORDIC units */
  173.   int    sinTheta;   /* calculated sin in CORDIC units */
  174.   int    cosTheta;   /* calculated cos in CORDIC units */
  175.   double dSin;       /* error in CORDIC calculation */
  176.   double dCos;       /* error in CORDIC calculation */
  177.  
  178.   gotoxy(22, 20);
  179.   printf("Enter angle, in CORDIC units : ");
  180.   scanf("%u", &theta);
  181.  
  182.   SinCos(theta, &sinTheta, &cosTheta);
  183.   dSin = sin(CORDIC_TO_RADIAN(theta)) - (double) sinTheta / CordicBase;
  184.   dCos = cos(CORDIC_TO_RADIAN(theta)) - (double) cosTheta / CordicBase;
  185.   gotoxy(18, 22);
  186.   printf("%u CORDIC : dSin = %8.5f   dCos = %8.5f", theta, dSin, dCos);
  187.   gotoxy(24, 25);
  188.   printf("Press a key to continue ......");
  189.   getch();
  190. }
  191.  
  192. void ListSeries(void)
  193. /*
  194.   USE : Calculate sin/cos for series of values (0 deg - 89 deg).
  195.   NOTE: Also compares with actual values from C library functions.
  196. */
  197. {
  198.   WORD   theta;      /* increment from 0 to 89 (degrees) */
  199.   int    sinTheta;   /* calculated sin in CORDIC units */
  200.   int    cosTheta;   /* calculated cos in CORDIC units */
  201.   double dSin;       /* error in CORDIC calculation */
  202.   double dCos;       /* error in CORDIC calculation */
  203.  
  204.   clrscr();
  205.   gotoxy(14,1); printf("╔═══════════════════════════════════════════════╗");
  206.   gotoxy(14,2); printf("║   Press <SPACE> to continue, <ESC> to quit.   ║");
  207.   gotoxy(14,3); printf("╚═══════════════════════════════════════════════╝");
  208.   gotoxy(1, 5);
  209.  
  210.   for (theta = 0; theta < 90; theta++)
  211.     {
  212.     SinCos(DEG_TO_CORDIC(theta), &sinTheta, &cosTheta);
  213.     dSin = sin(DEG_TO_RADIAN(theta)) - (double) sinTheta / CordicBase;
  214.     dCos = cos(DEG_TO_RADIAN(theta)) - (double) cosTheta / CordicBase;
  215.     printf("%3u deg : dSin = %8.5f   dCos = %8.5f\n", theta, dSin, dCos);
  216.     if (getch() == ESC) return;
  217.     }
  218. }
  219.  
  220. void CalcStatistics(void)
  221. /*
  222.   USE : Calculate average error and worst error for all angles.
  223.   NOTE: Get worstError = 0.00064, avgError = 0.00011 for NBITS = 14
  224.         This is 13+ bits of precision for the average.
  225. */
  226. {
  227.   WORD   theta;      /* to index thru all angles, in CORDIC units */
  228.   double thetaRad;   /* theta in radians */
  229.   int    sinTheta;   /* calculated sin in CORDIC units */
  230.   int    cosTheta;   /* calculated cos in CORDIC units */
  231.   double dSin;       /* error in CORDIC calculation (absVal) */
  232.   double dCos;       /* error in CORDIC calculation (absVal) */
  233.   double accumError; /* accumulated error between CORDIC/actual sin/cos */
  234.   double avgError;   /* average error per sin/cos */
  235.   double worstError; /* worst error for all sin/cos */
  236.  
  237.   gotoxy(25,18); printf("╔══════════════════════════╗");
  238.   gotoxy(25,19); printf("║   Press <ESC> to quit.   ║");
  239.   gotoxy(25,20); printf("╚══════════════════════════╝");
  240.   gotoxy(1, 21);
  241.  
  242.   worstError = 0.0;
  243.   accumError = 0.0;
  244.  
  245.   /* Update worstError and accumError as progress thru entire series. */
  246.   for (theta = 0; theta < CordicBase; theta++)
  247.     {
  248.     SinCos(theta, &sinTheta, &cosTheta);
  249.     thetaRad = CORDIC_TO_RADIAN(theta);
  250.     dSin = fabs(sin(thetaRad) - (double) sinTheta / CordicBase);
  251.     dCos = fabs(cos(thetaRad) - (double) cosTheta / CordicBase);
  252.     if (dSin > worstError) worstError = dSin;
  253.     if (dCos > worstError) worstError = dCos;
  254.     accumError += (dSin + dCos);
  255.     /* Put '.' and check for ESC every 256th time. */
  256.     if ((theta & 0x00FF) == 0)
  257.       {
  258.       putch('.');
  259.       if (kbhit() && (getch() == ESC)) return;
  260.       }
  261.     }  /* for (theta ... */
  262.  
  263.   /* Calculate avgError and display statistics. */
  264.   avgError = 0.5 * accumError / CordicBase;
  265.   gotoxy(15, 23);
  266.   printf("Worst error = %8.5f  Average error = %8.5f", worstError, avgError);
  267.   gotoxy(24, 25);
  268.   printf("Press a key to continue ......");
  269.   getch();
  270. }  
  271.  
  272. void TimeTest(void)
  273. /*
  274.   USE : Timing test for CORDIC vs C Library sine/cosine.
  275.   NOTE: Call CalcHexPtsCORDIC() and CalcHexPtsCLib() each NUM_TIMES times.
  276.         These routines use the CORDIC/long integer and C Library/floating
  277.         point methods respectively for rotating to calculate hexagon pts.
  278. */
  279. {
  280.   int   numTimes;    /* counter to iterate for timing */
  281.   TIME  t0, t1;      /* to use with gettime() */
  282.   POINT center;      /* center of hexagon */
  283.   POINT vertex;      /* vertex of hexagon */
  284.  
  285.   center.x = 320;  center.y = 240;
  286.   vertex.x = 400;  vertex.y = 305;
  287.  
  288.   gotoxy(25,18); printf("╔══════════════════════════╗");
  289.   gotoxy(25,19); printf("║   Timing .............   ║");
  290.   gotoxy(25,20); printf("╚══════════════════════════╝");
  291.   gotoxy(1, 22);
  292.  
  293.   gettime(&t0);
  294.   for (numTimes = 0; numTimes < NUM_TIMES; numTimes++)
  295.     CalcHexPtsCORDIC(center, vertex);
  296.   gettime(&t1);
  297.   printf("%d iterations, CORDIC sin/cos = %d hundSecs\n", NUM_TIMES, HundSecs(&t0, &t1));
  298.  
  299.   gettime(&t0);
  300.   for (numTimes = 0; numTimes < NUM_TIMES; numTimes++)
  301.     CalcHexPtsCLib(center, vertex);
  302.   gettime(&t1);
  303.   printf("%d iterations, C's sin/cos = %d hundSecs\n", NUM_TIMES, HundSecs(&t0, &t1));
  304.   gotoxy(24, 25);
  305.   printf("Press a key to continue ......");
  306.   getch();
  307. }
  308.  
  309. int HundSecs(TIME *tp, TIME *sp)
  310. /*
  311.   USE : Calculate hundredths of a second between two times.
  312.   IN  : tp = ptr to 1st time
  313.         sp = ptr to 2nd time
  314.   RET : hundredths of a second elapsed
  315. */
  316. {
  317.   long t, s;  /* to convert from structures to hundredth of second */
  318.  
  319.   t = 100L * (3600L*tp->ti_hour + 60*tp->ti_min + tp->ti_sec) + tp->ti_hund;
  320.   s = 100L * (3600L*sp->ti_hour + 60*sp->ti_min + sp->ti_sec) + sp->ti_hund;
  321.  
  322.   return ((int) (s - t));
  323. }
  324.  
  325. void CalcHexPtsCLib(POINT center, POINT vertex)
  326. /*
  327.   USE:  Calculate array of hexagon vertices using C Lib calls.
  328.   IN:   center = center of hexagon.
  329.         vertex = one of the hexagon vertices
  330.   NOTE: Loads global array HexPoints[] with other 5 vertices of the
  331.         regular hexagon. Uses C's floating point routines for trig
  332.         and floating point calculations.
  333. */
  334. {
  335.   double  sinTheta;   /* sine of central angle */
  336.   double  cosTheta;   /* cosine of central angle */
  337.   DOUB_PT centerD;    /* center of rotation, as doubles */
  338.   DOUB_PT rotPt;      /* rotated point, as doubles */
  339.   DOUB_PT delta;      /* rotPt - center (a radial spoke) */
  340.   int     corner;     /* index for vertices of polygon */
  341.  
  342.   /* 60 deg = 1.047197551 radians : use constants for accurate timing. */
  343.   sinTheta = sin(1.047197551);
  344.   /* sqrt(1 - sin * sin) is faster than cos */
  345.   cosTheta = sqrt(1 - sinTheta * sinTheta);
  346.  
  347.   /* Convert center to doubles and store. */
  348.   centerD.u = (double) center.x;
  349.   centerD.v = (double) center.y;
  350.  
  351.   /* Set initial and final point to incoming vertex; convert vertex
  352.      to doubles and store in rotPt. */
  353.   rotPt.u = (double) (HexPts[0].x = HexPts[NUM_PTS].x = vertex.x);
  354.   rotPt.v = (double) (HexPts[0].y = HexPts[NUM_PTS].y = vertex.y);
  355.  
  356.   /* Go clockwise around circle calculating hexagon points. */
  357.   for (corner = 1; corner < NUM_PTS; corner++)
  358.     {
  359.     /* Calculate the radial spoke : rotPt - center. */
  360.     delta.u = rotPt.u - centerD.u;
  361.     delta.v = rotPt.v - centerD.v;
  362.     /* Generate new rotPt by rotating around center. */
  363.     rotPt.u = (delta.u * cosTheta - delta.v * sinTheta) + centerD.u;
  364.     rotPt.v = (delta.u * sinTheta + delta.v * cosTheta) + centerD.v;
  365.     /* Round rotPt and store in array. */
  366.     HexPts[corner].x = (int) (rotPt.u + 0.5);
  367.     HexPts[corner].y = (int) (rotPt.v + 0.5);
  368.     }
  369. }
  370.  
  371. void CalcHexPtsCORDIC(POINT center, POINT vertex)
  372. /*
  373.   USE:  Calculate array of hexagon vertices using CORDIC calculations.
  374.   IN:   center = center of hexagon.
  375.         vertex = one of the hexagon vertices
  376.   NOTE: Loads global array HexPoints[] with other 5 vertices of the
  377.         regular hexagon. Uses CORDIC routines for trig and long integer
  378.         calculations.
  379. */
  380. {
  381.   int   sinTheta;     /* sine of central angle */
  382.   int   cosTheta;     /* cosine of central angle */
  383.   int   corner;       /* index for vertices of polygon */
  384.   POINT delta;        /* vertex - center (a radial spoke) */
  385.  
  386.   /* 60 deg = 10923 CORDIC units. */
  387.   SinCos(10923, &sinTheta, &cosTheta);
  388.  
  389.   /* Set initial and final point to incoming vertex. */
  390.   HexPts[0].x = HexPts[NUM_PTS].x = vertex.x;
  391.   HexPts[0].y = HexPts[NUM_PTS].y = vertex.y;
  392.  
  393.   /* Go clockwise around circle calculating hexagon points. */
  394.   for (corner = 1; corner < NUM_PTS; corner++)
  395.     {
  396.     /* Calculate the radial spoke : vertex - center. */
  397.     delta.x = vertex.x - center.x;
  398.     delta.y = vertex.y - center.y;
  399.     /* Generate new vertex by rotating around center. */
  400.     vertex.x = (int) (((long) delta.x * cosTheta - (long) delta.y * sinTheta + HalfBase) >> NBITS) + center.x;
  401.     vertex.y = (int) (((long) delta.x * sinTheta + (long) delta.y * cosTheta + HalfBase) >> NBITS) + center.y;
  402.     /* Store new vertex in array. */
  403.     HexPts[corner].x = vertex.x;
  404.     HexPts[corner].y = vertex.y;
  405.     }
  406. }
  407.  
  408. void SinCosSetup(void)
  409. /*
  410.   USE : Load globals used by SinCos().
  411.   OUT : Loads globals used in SinCos() :
  412.         CordicBase    = base for CORDIC units
  413.         HalfBase      = Cordicbase / 2
  414.         Quad2Boundary = 2 * CordicBase
  415.         Quad3Boundary = 3 * CordicBase
  416.         ArcTan[]      = the arctans of 1/(2^i)
  417.         xInit         = initial value for x projection, sufficiently
  418.                         less than CordicBase to exactly compensate
  419.                         for the expansion in each rotation).
  420.   NOTE: Must be called once only to initialize before calling SinCos().
  421. */
  422. {
  423.   int i;        /* to index ArcTan[] */
  424.   double f;     /* to calculate initial x projection */
  425.   long powr;    /* for powers of 2 : up to 2^(2*(NBITS-1)) */
  426.  
  427.   CordicBase = 1 << NBITS;
  428.   HalfBase = CordicBase >> 1;
  429.   Quad2Boundary = CordicBase << 1;
  430.   Quad3Boundary = CordicBase + Quad2Boundary;
  431.  
  432.   /* ArcTan's are the arctans of diminishingly small angles. */
  433.   powr = 1;
  434.   for (i = 0; i < NBITS; i++)
  435.     {
  436.     ArcTan[i] = (int)(atan (1.0/powr) / (M_PI / 2) * CordicBase + 0.5);
  437.     powr <<= 1;
  438.     }
  439.  
  440.   /* xInit is initial value of x projection to compensate for expansions.
  441.      f = 1/sqrt(2/1 * 5/4 * 17/16 * 65/64 * ...). Normalize as an NBITS
  442.      binary fraction (multiply by CordicBase) and store in xInit.
  443.      Get f = 0.607253 and xInit = 9949 = 0x26DD for NBITS = 14. */
  444.   f = 1.0;
  445.   powr = 1;
  446.   for (i = 0; i < NBITS; i++)
  447.     {
  448.     f = (f * (powr + 1)) / powr;
  449.     powr <<= 2;
  450.     }
  451.   f = 1.0/sqrt(f);
  452.   xInit = (int) (CordicBase * f + 0.5);
  453. }
  454.  
  455. void SinCos(WORD theta, int *sin, int *cos)
  456. /*
  457.   USE : Calculate sin and cos with integer CORDIC routine.
  458.   IN  : theta = angle whose sin and cos we want (in CORDIC angle units)
  459.   OUT : sin = ptr to calculated sin (in CORDIC fixed point units)
  460.         cos = ptr to calculated cos (in CORDIC fixed point units)
  461.   NOTE: The incoming angle theta is in CORDIC angle units, which subdivide
  462.         the unit circle into 64K parts, with 0 deg = 0, 90 deg = 16384
  463.         (16384 = CordicBase), 180 deg = 32768, 270 deg = 49152, etc. The
  464.         calculated sine and cosine are in CORDIC fixed point units : the
  465.         calculated value must be considered a fraction of 16384 (CordicBase).
  466. */
  467. {
  468.   int quadrant;    /* quadrant of incoming angle */
  469.   int z;           /* incoming angle translated to 1st quadrant */
  470.   int i;           /* to index rotations : one per bit */
  471.   int x, y;        /* projections onto axes */
  472.   int x1, y1;      /* projections of rotated vector */
  473.  
  474.   /* Determine quadrant of incoming angle, translate to 1st quadrant.
  475.      Note use of previously calculated values CordicBase, Quad2Boundary,
  476.      Quad3Boundary for speed. */
  477.   if (theta < CordicBase)
  478.     {
  479.     quadrant = QUAD1;
  480.     z = (int) theta;
  481.     }
  482.   else if (theta < Quad2Boundary)
  483.     {
  484.     quadrant = QUAD2;
  485.     z = (int) (Quad2Boundary - theta);
  486.     }
  487.   else if (theta < Quad3Boundary)
  488.     {
  489.     quadrant = QUAD3;
  490.     z = (int) (theta - Quad2Boundary);
  491.     }
  492.   else
  493.     {
  494.     quadrant = QUAD4;
  495.     z = - ((int) theta);
  496.     }
  497.  
  498.   /* Initialize projections. */
  499.   x = xInit;
  500.   y = 0;
  501.  
  502.   /* Negate z, so same rotations taking angle z to 0 will take
  503.      (x, y) = (xInit, 0) to (*cos, *sin). */
  504.   z = -z;
  505.  
  506.   /* Rotate NBITS times. */
  507.   for (i = 0; i < NBITS; i++)
  508.     {
  509.     if (z < 0)
  510.       {
  511.       /* Counter-clockwise rotation : move z closer to 0 deg. */
  512.       z += ArcTan[i];
  513.       y1 = y + (x >> i);
  514.       x1 = x - (y >> i);
  515.       }
  516.     else
  517.       {
  518.       /* Clockwise rotation : move z closer to 0 deg. */
  519.       z -= ArcTan[i];
  520.       y1 = y - (x >> i);
  521.       x1 = x + (y >> i);
  522.       }
  523.  
  524.     /* Put new projections into (x,y) for next go. */
  525.     x = x1;
  526.     y = y1;
  527.     }  /* for i */
  528.  
  529.   /* Attach signs depending on quadrant. */
  530.   *cos = (quadrant == QUAD1 || quadrant == QUAD4) ? x : -x;
  531.   *sin = (quadrant == QUAD1 || quadrant == QUAD2) ? y : -y;
  532. }
  533.