November 1992/The CORDIC Method for Faster sin and cos Calculations/Listing 1
/* CORDIC.C : Demonstrates the integer-based CORDIC
* system for calculating sines and cosines. The
* vertices of a regular hexagon are calculated using
* CORDIC trig and the hexagon itself is rotated.
* Make in Borland C++'s internal environment.
* (c) 1991 Michael Bertrand.
*/
#include <stdio.h> /* printf */
#include <math.h> /* sqrt, atan */
#include <conio.h> /* getch */
#include <graphics.h> /* BGI */
typedef struct
{
int x;
int y;
} POINT;
typedef unsigned int WORD;
void CalcHexPtsCORDIC(POINT center, POINT vertex);
void DrawHexagon(POINT center, POINT vertex);
void DrawCross(POINT pt, int colr);
void SinCosSetup(void);
void SinCos(WORD theta, int *sin, int *cos);
#define ESC 0x1B
/* Quadrants for angles. */
#define QUAD1 1
#define QUAD2 2
#define QUAD3 3
#define QUAD4 4
/* NBITS is number of bits for CORDIC units. */
#define NBITS 14
/* NUM_PTS is number of vertices in polygon. */
#define NUM_PTS 6
int ArcTan[NBITS]; /* angles for rotation */
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]; /* calculated poly points */
void main(void)
{
int driver; /* for initgraph() */
int mode; /* for initgraph() */
WORD theta; /* CORDIC angle */
int sine; /* sine of CORDIC angle */
int cosine; /* cosine of CORDIC angle */
POINT center; /* center of hexagon */
POINT vertex; /* hexagon's original base vertex */
POINT vertex1; /* hexagon's changing base vertex */
POINT del; /* vertex - center (radial spoke) */
int radius; /* radius of circumscribing circle */
driver = VGA; /* for initgraph() */
mode = VGAHI; /* mode 0x12 : 640x480 16 color */
if (registerbgidriver(EGAVGA_driver) < 0)
{
printf("couldn't find VGA driver"); return;
}
initgraph(&driver, &mode, NULL);
printf("Press ENTER to rotate, ESC to quit.");
center.x = 320; center.y = 240;
vertex.x = 470; vertex.y = 240;
radius = vertex.x - center.x;
/* Calculate the radial spoke : vertex - center. */
del.x = vertex.x - center.x;
del.y = vertex.y - center.y;
setwritemode(XOR_PUT);
/* Draw circumscribing circle. */
setcolor(RED);
circle(center.x, center.y, radius);
/* Draw small cross at center. */
DrawCross(center, YELLOW);
setcolor(WHITE);
/* Setup CORDIC system, initialize theta = 0. */
SinCosSetup();
theta = 0;
/* Draw initial hexagon. */
DrawHexagon(center, vertex1 = vertex);
/* Rotate hexagon. vertex is fixed; vertex1 rotates
* clockwise around vertex in increments of 650
* CORDIC units (3.57 deg). CORDIC sines/cosines are
* used to find vertex1. DrawHexagon() also uses
* CORDIC sines/cosines to calculate the remaining
* vertices for each hexagon so they can be drawn.
*/
while (getch() ! = ESC)
{
/* Erase last hexagon. */
DrawHexagon(center, vertex1);
/* Inc theta by 650 CORDIC units (3.57 deg). */
theta += 650;
SinCos(theta, &sine, &cosine);
/* Calc new vertex by rotating around center. */
vertex1.x = (int)
(((long) del.x * cosine - (long) del.y * sine +
HalfBase) >> NBIIS) + center.x;
vertexl.y = (int)
(((long) del.x * sine + (long) del.y * cosine +
HalfBase) >> NBITS) + center.y;
/* Draw new hexagon. */
DrawHexagon(center, vertex1);
}
closegraph();
}
void CalcHexPtsCORDIC(POINT center, POINT vertex)
/*
USE: Calc array of hex vertices using CORDIC calcs.
IN: center = center of hexagon.
vertex = One of the hexagon vertices
NOTE: Loads global array HexPoints[] with other 5
vertices of regular hexagon. Uses CORDIC routines
for trig and long integer calculations.
*/
{
int sine; /* sine of central angle */
int cosine; /* cosine of central angle */
int corner; /* index for vertices of polygon */
POINT del; /* vertex - center (radial spoke) */
/* 60 deg = 10923 CORDIC units. */
SinCos(10923, &sine, &cosine);
/* 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 to calc hex points. */
for (corner = 1; corner < NUM_PTS; corner++)
{
/* Calculate the radial spoke : vertex - center. */
del.x = vertex.x - center.x;
del.y = vertex.y - center.y;
/* calc new vertex by rotating around center. */
vertex.x = (int)
(((long) del.x * cosine - (long) del.y * sine +
HalfBase) >> NBITS) + center.x;
vertex.y = (int)
(((long) del.x * sine + (long) del.y * cosine +
HalfBase) >> NBITS) + center.y;
/* Store new vertex in array. */
HexPts[corner].x = vertex.x;
HexPts[corner].y = vertex.y;
}
}
void DrawHexagon(POINT center, POINT vertex)
/*
USE: Draw Hexagon given center and one vertex.
IN: center = center of hexagon.
vertex = one of the hexagon vertices
NOTE: Call CalcHexPtsCORDIC() to load global array
HexPts[] with hexagon vertices.
*/
{
CalcHexPtsCORDIC(center, vertex);
drawpoly(NUM_PTS+1, (int far *)HexPts);
}
void DrawCross(POINT pt, int colr)
/*
USE: Draw cross on screen at pt with given color.
*/
{
int oldColor;
setwritemode(COPY_PUT);
oldColor = getcolor();
setcolor(colr);
moveto(pt.x - 2, pt.y); lineto(pt.x + 2, pt.y);
moveto(pt.x, pt.y - 2); lineto(pt.x, pt.y + 2);
setcolor(oldColor);
setwritemode(XOR_PUT);
}
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
NOTE: Must be called once only to initialize before
calling SinCos(). xInit is sufficiently less than
CordicBase to exactly compensate for the expansion
in each rotation.
*/
{
int i; /* to index ArcTan[] */
double f; /* to calc initial x projection */
long powr; /* powers of 2 up to 2^(2*(NBITS-1)) */
CordicBase = 1 << NBITS;
HalfBase = CordicBase >> 1;
Quad2Boundary = CordicBase << 1;
Quad3Boundary = CordicBase + Quad2Boundary;
/* ArcTan's are 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 * ...
* 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 : Calc sin and cos with integer CORDIC routine.
IN : theta = incoming angle (in CORDIC angle units)
OUT : sin = ptr to sin (in CORDIC fixed point units)
cos = ptr to cos (in CORDIC fixed point units)
NOTE: The incoming angle theta is in CORDIC angle
units, which subdivide the circle into 64K parts,
with 0 deg = 0, 90 deg = 16384 (CordicBase), 180 deg
= 32768, 270 deg = 49152, etc. The calculated sine
and cosine are in CORDIC fixed point units : an int
considered as a fraction of 16384 (CordicBase).
*/
{
int quadrant; /* quadrant of incoming angle */
int z; /* incoming angle moved to 1st quad */
int i; /* to index rotations : one per bit */
int x, y; /* projections onto axes */
int xl, yl; /* projections of rotated vector */
/* Determine quadrant of incoming angle, translate to
* 1st quadrant. Note use of previously calculated
* values CordicBase, etc. 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. */
z += ArcTan[i];
y1 = y + (x >> i);
x1 = x - (y >> i);
}
else
{
/* Clockwise rotation. */
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;
}
/* End of File */
