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 >
Wrap
Text File
|
1992-06-18
|
18KB
|
533 lines
/* CRDTOOLS.C : Tools to exercise CORDIC calculations for sine/cosine.
* Make in Borland C++'s internal environment.
* (c) 1991 Michael Bertrand.
*/
#include <stdio.h> /* printf, scanf */
#include <math.h> /* sqrt, atan, sin, cos, fabs */
#include <conio.h> /* clrscr, gotoxy, getch, cputs, etc. */
#include <dos.h> /* struct time, gettime */
typedef unsigned int WORD;
typedef struct time TIME;
typedef struct
{
int x;
int y;
} POINT;
typedef struct
{
double u;
double v;
} DOUB_PT;
char Menu(void);
void SingleAngleDegree(void);
void SingleAngleCordic(void);
void ListSeries(void);
void CalcStatistics(void);
void TimeTest(void);
int HundSecs(TIME *tp, TIME *sp);
void CalcHexPtsCLib(POINT center, POINT vertex);
void CalcHexPtsCORDIC(POINT center, POINT vertex);
void SinCosSetup(void);
void SinCos(WORD theta, int *sin, int *cos);
#define TRUE 1
#define ESC 0x1B
/* Quadrants for angles. */
#define QUAD1 1
#define QUAD2 2
#define QUAD3 3
#define QUAD4 4
/* NBITS is number of bits used for CORDIC units. */
#define NBITS 14
/* NUM_PTS is number of vertices in polygon. */
#define NUM_PTS 6
/* Macros to convert angles to different units. */
#define DEG_TO_CORDIC(x) ((WORD) ((double)x/90.0*CordicBase + 0.5))
#define DEG_TO_RADIAN(x) (x/180.0*M_PI)
#define CORDIC_TO_RADIAN(x) ((double)x/CordicBase*M_PI/2)
/* NUM_TIMES is number of times to iterate for timing test. */
#define NUM_TIMES 1000
int ArcTan[NBITS]; /* angular constants : arctans */
int xInit; /* initial x projection */
WORD CordicBase; /* base for CORDIC units */
WORD HalfBase; /* CordicBase / 2 */
WORD Quad2Boundary; /* 2 * CordicBase */
WORD Quad3Boundary; /* 3 * CordicBase */
POINT HexPts[NUM_PTS+1]; /* to hold calculated polygon points */
void main(void)
{
SinCosSetup();
do
{
switch(Menu())
{
case '1' : case 'd': case 'D':
SingleAngleDegree();
break;
case '2' : case 'c': case 'C':
SingleAngleCordic();
break;
case '3' : case 'l': case 'L':
ListSeries();
break;
case '4' : case 's': case 'S':
CalcStatistics();
break;
case '5' : case 't': case 'T':
TimeTest();
break;
case '6' : case 'q': case 'Q':
clrscr();
return;
} /* switch */
} while (TRUE);
}
char Menu(void)
/*
USE : Display menu and get user choice.
RET : User choice.
*/
{
struct text_info ti; /* to get default text color */
gettextinfo(&ti);
clrscr();
textattr(LIGHTGRAY);
gotoxy(26, 1); cputs("╔═════ CORDIC Tests ═════╗");
gotoxy(26, 2); cputs("║ ║");
gotoxy(26, 3); cputs("║ 1) One Angle (Degree) ║");
gotoxy(26, 4); cputs("║ ║");
gotoxy(26, 5); cputs("║ 2) One Angle (CORDIC) ║");
gotoxy(26, 6); cputs("║ ║");
gotoxy(26, 7); cputs("║ 3) List Series ║");
gotoxy(26, 8); cputs("║ ║");
gotoxy(26, 9); cputs("║ 4) Statistics ║");
gotoxy(26,10); cputs("║ ║");
gotoxy(26,11); cputs("║ 5) Time Test ║");
gotoxy(26,12); cputs("║ ║");
gotoxy(26,13); cputs("║ 6) Quit ║");
gotoxy(26,14); cputs("║ ║");
gotoxy(26,15); cputs("║ Enter choice : ║");
gotoxy(26,16); cputs("╚════════════════════════╝");
textattr(RED);
gotoxy(43, 3); putch('D');
gotoxy(43, 5); putch('C');
gotoxy(32, 7); putch('L');
gotoxy(32, 9); putch('S');
gotoxy(32,11); putch('T');
gotoxy(32,13); putch('Q');
gotoxy(46, 15); /* put cursor inside box */
textattr(ti.attribute); /* restore default text color */
return(getch()); /* get char from user, return to main() */
}
void SingleAngleDegree(void)
/*
USE : Calculate sin/cos for single angle, in degrees.
NOTE: Also compares with actual values from C library functions.
*/
{
WORD theta; /* user-entered angle, in degrees */
int sinTheta; /* calculated sin in CORDIC units */
int cosTheta; /* calculated cos in CORDIC units */
double dSin; /* error in CORDIC calculation */
double dCos; /* error in CORDIC calculation */
gotoxy(25, 20);
printf("Enter angle, in degrees : ");
scanf("%u", &theta);
SinCos(DEG_TO_CORDIC(theta), &sinTheta, &cosTheta);
dSin = sin(DEG_TO_RADIAN(theta)) - (double) sinTheta / CordicBase;
dCos = cos(DEG_TO_RADIAN(theta)) - (double) cosTheta / CordicBase;
gotoxy(18, 22);
printf("%u deg : dSin = %8.5f dCos = %8.5f", theta, dSin, dCos);
gotoxy(24, 25);
printf("Press a key to continue ......");
getch();
}
void SingleAngleCordic(void)
/*
USE : Calculate sin/cos for single angle, in CORDIC units.
NOTE: Also compares with actual values from C library functions.
*/
{
WORD theta; /* user-entered angle, in CORDIC units */
int sinTheta; /* calculated sin in CORDIC units */
int cosTheta; /* calculated cos in CORDIC units */
double dSin; /* error in CORDIC calculation */
double dCos; /* error in CORDIC calculation */
gotoxy(22, 20);
printf("Enter angle, in CORDIC units : ");
scanf("%u", &theta);
SinCos(theta, &sinTheta, &cosTheta);
dSin = sin(CORDIC_TO_RADIAN(theta)) - (double) sinTheta / CordicBase;
dCos = cos(CORDIC_TO_RADIAN(theta)) - (double) cosTheta / CordicBase;
gotoxy(18, 22);
printf("%u CORDIC : dSin = %8.5f dCos = %8.5f", theta, dSin, dCos);
gotoxy(24, 25);
printf("Press a key to continue ......");
getch();
}
void ListSeries(void)
/*
USE : Calculate sin/cos for series of values (0 deg - 89 deg).
NOTE: Also compares with actual values from C library functions.
*/
{
WORD theta; /* increment from 0 to 89 (degrees) */
int sinTheta; /* calculated sin in CORDIC units */
int cosTheta; /* calculated cos in CORDIC units */
double dSin; /* error in CORDIC calculation */
double dCos; /* error in CORDIC calculation */
clrscr();
gotoxy(14,1); printf("╔═══════════════════════════════════════════════╗");
gotoxy(14,2); printf("║ Press <SPACE> to continue, <ESC> to quit. ║");
gotoxy(14,3); printf("╚═══════════════════════════════════════════════╝");
gotoxy(1, 5);
for (theta = 0; theta < 90; theta++)
{
SinCos(DEG_TO_CORDIC(theta), &sinTheta, &cosTheta);
dSin = sin(DEG_TO_RADIAN(theta)) - (double) sinTheta / CordicBase;
dCos = cos(DEG_TO_RADIAN(theta)) - (double) cosTheta / CordicBase;
printf("%3u deg : dSin = %8.5f dCos = %8.5f\n", theta, dSin, dCos);
if (getch() == ESC) return;
}
}
void CalcStatistics(void)
/*
USE : Calculate average error and worst error for all angles.
NOTE: Get worstError = 0.00064, avgError = 0.00011 for NBITS = 14
This is 13+ bits of precision for the average.
*/
{
WORD theta; /* to index thru all angles, in CORDIC units */
double thetaRad; /* theta in radians */
int sinTheta; /* calculated sin in CORDIC units */
int cosTheta; /* calculated cos in CORDIC units */
double dSin; /* error in CORDIC calculation (absVal) */
double dCos; /* error in CORDIC calculation (absVal) */
double accumError; /* accumulated error between CORDIC/actual sin/cos */
double avgError; /* average error per sin/cos */
double worstError; /* worst error for all sin/cos */
gotoxy(25,18); printf("╔══════════════════════════╗");
gotoxy(25,19); printf("║ Press <ESC> to quit. ║");
gotoxy(25,20); printf("╚══════════════════════════╝");
gotoxy(1, 21);
worstError = 0.0;
accumError = 0.0;
/* Update worstError and accumError as progress thru entire series. */
for (theta = 0; theta < CordicBase; theta++)
{
SinCos(theta, &sinTheta, &cosTheta);
thetaRad = CORDIC_TO_RADIAN(theta);
dSin = fabs(sin(thetaRad) - (double) sinTheta / CordicBase);
dCos = fabs(cos(thetaRad) - (double) cosTheta / CordicBase);
if (dSin > worstError) worstError = dSin;
if (dCos > worstError) worstError = dCos;
accumError += (dSin + dCos);
/* Put '.' and check for ESC every 256th time. */
if ((theta & 0x00FF) == 0)
{
putch('.');
if (kbhit() && (getch() == ESC)) return;
}
} /* for (theta ... */
/* Calculate avgError and display statistics. */
avgError = 0.5 * accumError / CordicBase;
gotoxy(15, 23);
printf("Worst error = %8.5f Average error = %8.5f", worstError, avgError);
gotoxy(24, 25);
printf("Press a key to continue ......");
getch();
}
void TimeTest(void)
/*
USE : Timing test for CORDIC vs C Library sine/cosine.
NOTE: Call CalcHexPtsCORDIC() and CalcHexPtsCLib() each NUM_TIMES times.
These routines use the CORDIC/long integer and C Library/floating
point methods respectively for rotating to calculate hexagon pts.
*/
{
int numTimes; /* counter to iterate for timing */
TIME t0, t1; /* to use with gettime() */
POINT center; /* center of hexagon */
POINT vertex; /* vertex of hexagon */
center.x = 320; center.y = 240;
vertex.x = 400; vertex.y = 305;
gotoxy(25,18); printf("╔══════════════════════════╗");
gotoxy(25,19); printf("║ Timing ............. ║");
gotoxy(25,20); printf("╚══════════════════════════╝");
gotoxy(1, 22);
gettime(&t0);
for (numTimes = 0; numTimes < NUM_TIMES; numTimes++)
CalcHexPtsCORDIC(center, vertex);
gettime(&t1);
printf("%d iterations, CORDIC sin/cos = %d hundSecs\n", NUM_TIMES, HundSecs(&t0, &t1));
gettime(&t0);
for (numTimes = 0; numTimes < NUM_TIMES; numTimes++)
CalcHexPtsCLib(center, vertex);
gettime(&t1);
printf("%d iterations, C's sin/cos = %d hundSecs\n", NUM_TIMES, HundSecs(&t0, &t1));
gotoxy(24, 25);
printf("Press a key to continue ......");
getch();
}
int HundSecs(TIME *tp, TIME *sp)
/*
USE : Calculate hundredths of a second between two times.
IN : tp = ptr to 1st time
sp = ptr to 2nd time
RET : hundredths of a second elapsed
*/
{
long t, s; /* to convert from structures to hundredth of second */
t = 100L * (3600L*tp->ti_hour + 60*tp->ti_min + tp->ti_sec) + tp->ti_hund;
s = 100L * (3600L*sp->ti_hour + 60*sp->ti_min + sp->ti_sec) + sp->ti_hund;
return ((int) (s - t));
}
void CalcHexPtsCLib(POINT center, POINT vertex)
/*
USE: Calculate array of hexagon vertices using C Lib calls.
IN: center = center of hexagon.
vertex = one of the hexagon vertices
NOTE: Loads global array HexPoints[] with other 5 vertices of the
regular hexagon. Uses C's floating point routines for trig
and floating point calculations.
*/
{
double sinTheta; /* sine of central angle */
double cosTheta; /* cosine of central angle */
DOUB_PT centerD; /* center of rotation, as doubles */
DOUB_PT rotPt; /* rotated point, as doubles */
DOUB_PT delta; /* rotPt - center (a radial spoke) */
int corner; /* index for vertices of polygon */
/* 60 deg = 1.047197551 radians : use constants for accurate timing. */
sinTheta = sin(1.047197551);
/* sqrt(1 - sin * sin) is faster than cos */
cosTheta = sqrt(1 - sinTheta * sinTheta);
/* Convert center to doubles and store. */
centerD.u = (double) center.x;
centerD.v = (double) center.y;
/* Set initial and final point to incoming vertex; convert vertex
to doubles and store in rotPt. */
rotPt.u = (double) (HexPts[0].x = HexPts[NUM_PTS].x = vertex.x);
rotPt.v = (double) (HexPts[0].y = HexPts[NUM_PTS].y = vertex.y);
/* Go clockwise around circle calculating hexagon points. */
for (corner = 1; corner < NUM_PTS; corner++)
{
/* Calculate the radial spoke : rotPt - center. */
delta.u = rotPt.u - centerD.u;
delta.v = rotPt.v - centerD.v;
/* Generate new rotPt by rotating around center. */
rotPt.u = (delta.u * cosTheta - delta.v * sinTheta) + centerD.u;
rotPt.v = (delta.u * sinTheta + delta.v * cosTheta) + centerD.v;
/* Round rotPt and store in array. */
HexPts[corner].x = (int) (rotPt.u + 0.5);
HexPts[corner].y = (int) (rotPt.v + 0.5);
}
}
void CalcHexPtsCORDIC(POINT center, POINT vertex)
/*
USE: Calculate array of hexagon vertices using CORDIC calculations.
IN: center = center of hexagon.
vertex = one of the hexagon vertices
NOTE: Loads global array HexPoints[] with other 5 vertices of the
regular hexagon. Uses CORDIC routines for trig and long integer
calculations.
*/
{
int sinTheta; /* sine of central angle */
int cosTheta; /* cosine of central angle */
int corner; /* index for vertices of polygon */
POINT delta; /* vertex - center (a radial spoke) */
/* 60 deg = 10923 CORDIC units. */
SinCos(10923, &sinTheta, &cosTheta);
/* Set initial and final point to incoming vertex. */
HexPts[0].x = HexPts[NUM_PTS].x = vertex.x;
HexPts[0].y = HexPts[NUM_PTS].y = vertex.y;
/* Go clockwise around circle calculating hexagon points. */
for (corner = 1; corner < NUM_PTS; corner++)
{
/* Calculate the radial spoke : vertex - center. */
delta.x = vertex.x - center.x;
delta.y = vertex.y - center.y;
/* Generate new vertex by rotating around center. */
vertex.x = (int) (((long) delta.x * cosTheta - (long) delta.y * sinTheta + HalfBase) >> NBITS) + center.x;
vertex.y = (int) (((long) delta.x * sinTheta + (long) delta.y * cosTheta + HalfBase) >> NBITS) + center.y;
/* Store new vertex in array. */
HexPts[corner].x = vertex.x;
HexPts[corner].y = vertex.y;
}
}
void SinCosSetup(void)
/*
USE : Load globals used by SinCos().
OUT : Loads globals used in SinCos() :
CordicBase = base for CORDIC units
HalfBase = Cordicbase / 2
Quad2Boundary = 2 * CordicBase
Quad3Boundary = 3 * CordicBase
ArcTan[] = the arctans of 1/(2^i)
xInit = initial value for x projection, sufficiently
less than CordicBase to exactly compensate
for the expansion in each rotation).
NOTE: Must be called once only to initialize before calling SinCos().
*/
{
int i; /* to index ArcTan[] */
double f; /* to calculate initial x projection */
long powr; /* for powers of 2 : up to 2^(2*(NBITS-1)) */
CordicBase = 1 << NBITS;
HalfBase = CordicBase >> 1;
Quad2Boundary = CordicBase << 1;
Quad3Boundary = CordicBase + Quad2Boundary;
/* ArcTan's are the arctans of diminishingly small angles. */
powr = 1;
for (i = 0; i < NBITS; i++)
{
ArcTan[i] = (int)(atan (1.0/powr) / (M_PI / 2) * CordicBase + 0.5);
powr <<= 1;
}
/* xInit is initial value of x projection to compensate for expansions.
f = 1/sqrt(2/1 * 5/4 * 17/16 * 65/64 * ...). Normalize as an NBITS
binary fraction (multiply by CordicBase) and store in xInit.
Get f = 0.607253 and xInit = 9949 = 0x26DD for NBITS = 14. */
f = 1.0;
powr = 1;
for (i = 0; i < NBITS; i++)
{
f = (f * (powr + 1)) / powr;
powr <<= 2;
}
f = 1.0/sqrt(f);
xInit = (int) (CordicBase * f + 0.5);
}
void SinCos(WORD theta, int *sin, int *cos)
/*
USE : Calculate sin and cos with integer CORDIC routine.
IN : theta = angle whose sin and cos we want (in CORDIC angle units)
OUT : sin = ptr to calculated sin (in CORDIC fixed point units)
cos = ptr to calculated cos (in CORDIC fixed point units)
NOTE: The incoming angle theta is in CORDIC angle units, which subdivide
the unit circle into 64K parts, with 0 deg = 0, 90 deg = 16384
(16384 = CordicBase), 180 deg = 32768, 270 deg = 49152, etc. The
calculated sine and cosine are in CORDIC fixed point units : the
calculated value must be considered a fraction of 16384 (CordicBase).
*/
{
int quadrant; /* quadrant of incoming angle */
int z; /* incoming angle translated to 1st quadrant */
int i; /* to index rotations : one per bit */
int x, y; /* projections onto axes */
int x1, y1; /* projections of rotated vector */
/* Determine quadrant of incoming angle, translate to 1st quadrant.
Note use of previously calculated values CordicBase, Quad2Boundary,
Quad3Boundary for speed. */
if (theta < CordicBase)
{
quadrant = QUAD1;
z = (int) theta;
}
else if (theta < Quad2Boundary)
{
quadrant = QUAD2;
z = (int) (Quad2Boundary - theta);
}
else if (theta < Quad3Boundary)
{
quadrant = QUAD3;
z = (int) (theta - Quad2Boundary);
}
else
{
quadrant = QUAD4;
z = - ((int) theta);
}
/* Initialize projections. */
x = xInit;
y = 0;
/* Negate z, so same rotations taking angle z to 0 will take
(x, y) = (xInit, 0) to (*cos, *sin). */
z = -z;
/* Rotate NBITS times. */
for (i = 0; i < NBITS; i++)
{
if (z < 0)
{
/* Counter-clockwise rotation : move z closer to 0 deg. */
z += ArcTan[i];
y1 = y + (x >> i);
x1 = x - (y >> i);
}
else
{
/* Clockwise rotation : move z closer to 0 deg. */
z -= ArcTan[i];
y1 = y - (x >> i);
x1 = x + (y >> i);
}
/* Put new projections into (x,y) for next go. */
x = x1;
y = y1;
} /* for i */
/* Attach signs depending on quadrant. */
*cos = (quadrant == QUAD1 || quadrant == QUAD4) ? x : -x;
*sin = (quadrant == QUAD1 || quadrant == QUAD2) ? y : -y;
}