/* sincos.c
 *
 * From Eyal Lebedinsky, eyal@ise.canberra.edu.au
 * 8 Jan 94.
 *
 * fractions kept such that 1.0 is 16*1024, so the range is [-2.0,2.0)
 * angles kept such that -180 degrees is 0x8000, so the range is [-180,180)
 *
 * You must call set_tring() first. I actually dump the tables and then
 * use them as initialized static arrays, no need to initialize. Floating
 * point is used ONLY during the initialization phase, <math.h> only used
 * then.
 *
 * to compile, try:
 * gcc -o sincos sincos.c -lm
*/

#include <stdio.h>
#include <math.h>

#define NEAR
#define FAR
#define FASTCALL
#define ANGLE		short

#define FSCALE		14		/* fraction bits in sine/cos etc. */
#define FONE		(1<<FSCALE)
#define	FnDIV(n,x,y)	((int)(((long)(x) << (n)) / (y)))
#define	fdiv(x,y)	FnDIV (FSCALE, (x), (y))
#define D90		0x4000
#define D180		0x8000
#define DEG(x)		((ANGLE)((x)*(long)D90/90))

/* The following belongs in your application header to provide access to these
 * functions.
*/

#define	SIN(a)		my_sin (a)
#define	COS(a)		my_sin ((ANGLE)((a)+D90))
#define	ASIN(i)		my_asin(i)			/* [-90...+90] */
#define	ACOS(i)		((ANGLE)(D90-my_asin(i)))	/* [0...180] */
#define	ATAN(x,y)	my_atan((x),(y))

extern short	FAR FASTCALL my_sin (ANGLE angle);
extern ANGLE	FAR FASTCALL my_asin (int sin);
extern ANGLE	FAR FASTCALL my_atan (int y, int x);


/* This is part of my ifuncs.c module.
*/

/* angles are stored such that -pi..0..pi radians is 0x8000..0..0x7fff
*/

#define	P3_4	((1024/4)*3)
 
static short	NEAR sin_tab[1024+1] = {0};		/* 0..90 */
static ANGLE	NEAR asin_tab0[1024+1] = {0};		/* 0......1/4 */
static ANGLE	NEAR asin_tab1[P3_4+1] = {0};		/* 0......3/4 */
static ANGLE	NEAR asin_tab2[P3_4+1] = {0};		/* 3/4....15/16 */
static ANGLE	NEAR asin_tab3[1024+1] = {0};		/* 15/16..1 */
static ANGLE	NEAR atan_tab1[P3_4+1] = {0};		/* 0......3/4 */
static ANGLE	NEAR atan_tab2[P3_4+1] = {0};		/* 3/4....15/16 */
static ANGLE	NEAR atan_tab3[1024+1] = {0};		/* 15/16..1 */

static short FAR
my_round (double f)
{
	return ((short)((f<0) ? f-0.5 : f+0.5));
}

/* Fill in the trig tables.
*/
static void FAR
set_trig (void)
{
	double	s, c, t, d, sd, cd, pi;
	int	i;

	pi = atan(1.0)*4;
	d = pi/2048;			/* 0...pi */
	sd = sin (d);
	cd = cos (d);
	s = 0.0;
	c = 1.0*FONE;

	for (i = 0; i <= 1024; ++i) {
		sin_tab[i] = my_round (s);
		t = c*sd+s*cd;
		c = c*cd-s*sd;
		s = t;
	}

	for (i = 0; i <= 1024; ++i)
		asin_tab0[i] = (ANGLE)my_round (asin (i/16.0/1024.0)
				*0x08000/pi);
	for (i = 0; i <= P3_4; ++i)
		asin_tab1[i] = (ANGLE)my_round (asin (i*3.0/4.0/P3_4)
				*0x08000/pi);
	for (i = 0; i <= P3_4; ++i)
		asin_tab2[i] = (ANGLE)my_round (asin (
					3.0/4.0+i*(15.0/16.0-3.0/4.0)/P3_4)
					*0x08000/pi);
	for (i = 0; i <= 1024; ++i)
		asin_tab3[i] = (ANGLE)my_round (asin (
					15.0/16.0+i*(1.0-15.0/16.0)/1024)
					*0x08000/pi);

	for (i = 0; i <= P3_4; ++i)
		atan_tab1[i] = (ANGLE)my_round (atan (i*3.0/4.0/P3_4)
				*0x08000/pi);
	for (i = 0; i <= P3_4; ++i)
		atan_tab2[i] = (ANGLE)my_round (atan (
					3.0/4.0+i*(15.0/16.0-3.0/4.0)/P3_4)
					*0x08000/pi);
	for (i = 0; i <= 1024; ++i)
		atan_tab3[i] = (ANGLE)my_round (atan (
					15.0/16.0+i*(1.0-15.0/16.0)/1024)
					*0x08000/pi);
}

#undef P3_4

extern short FAR FASTCALL
my_sin (register ANGLE d)
{
	register int	ind, f, l, h;

	if (d < 0) {
		d = -d;
		if (d < 0)		/* -180 degrees! */
			return (0);
		ind = 1;
	} else
		ind = 0;
	if (d > D90)
		d = D180 - d;

	f = d&15;			/* d%16 */
	d >>= 4;			/* d/16 */
	l = sin_tab[d];
	h = sin_tab[d+1] - l;
	if (h < 0)
		l -= (8-h*f)>>4;	/* (8-h*f)/16; */
	else
		l += (8+h*f)>>4;	/* (8+h*f)/16; */

	return ((short)(ind ? -l : l));
}

extern ANGLE FAR FASTCALL
my_asin (int d)
{
	register int	ind, f;
	register ANGLE	l, h;

	if (d < 0) {
		d = -d;
		ind = 1;
	} else
		ind = 0;

	if (d < FONE/16) {		/* 0..1/16 */
		l = asin_tab0[d];
	} else if (d < (FONE/4)*3) {	/* 1/16..3/4 */
		f = d&15;
		d = d>>4;
		l = asin_tab1[d];
		h = asin_tab1[d+1] - l;
		if (h < 0)
			l -= (8-h*f)>>4;
		else
			l += (8+h*f)>>4;
	} else if (d < (FONE/16)*15) {	/* 3/4..15/16 */
		d -= (FONE/4)*3;	/* 0..1/4 */
		f = d&3;
		d = d>>2;
		l = asin_tab2[d];
		h = asin_tab2[d+1] - l;
		if (h < 0)
			l -= (2-h*f)/4;
		else
			l += (2+h*f)/4;
	} else {			/* 15/16..1 */
		if (d > FONE)
			d = FONE;
		d -= (FONE/16)*15;	/* 0..1/16 */
		l = asin_tab3[d];
	}
	return (ind ? -l : l);
}

extern ANGLE FAR FASTCALL
my_atan (int y, int x)
{
	int	neg, f, i, d;
	ANGLE	l, h;

	if (y < 0) {
		y = -y;
		neg = 1;
	} else
		neg = 0;

	if (x < 0) {
		x = -x;
		neg += 2;
	}

	if (y > x) {
		d = y;
		y = x;
		x = d;
		neg += 4;
	}

	if (x == 0)
		return (0);

	d = fdiv (y, x);		/* temp */

	if (d < (FONE/4)*3) {		/* 0..3/4 */
		f = d&15;
		i = d>>4;
		l = atan_tab1[i];
		h = atan_tab1[i+1] - l;
		if (h < 0)
			l -= (8-h*f)>>4;
		else
			l += (8+h*f)>>4;
	} else if (d < (FONE/16)*15) {	/* 3/4..15/16 */
		d -= (FONE/4)*3;
		f = d&3;
		i = d>>2;
		l = atan_tab2[i];
		h = atan_tab2[i+1] - l;
		if (h < 0)
			l -= (2-h*f)>>2;
		else
			l += (2+h*f)>>2;
	} else {			/* 15/16..1 */
		if (d > FONE)
			d = FONE;
		d -= (FONE/16)*15;
		l = atan_tab3[d];
	}

	switch (neg) {
	case 0:				/* ++ */
		return (l);
	case 1:				/* -+ */
		return (-l);
	case 2:				/* +- */
		return (D180-l);
	case 3:				/* -- */
		return (D180+l);
	case 4:				/* ++ rev */
		return (D90-l);
	case 5:				/* -+ rev */
		return (-D90+l);
	case 6:				/* +- rev */
		return (D90+l);
	case 7:				/* -- rev */
		return (-D90-l);
	}
	return (0);			/* never reached */
}

main ()
{
	ANGLE	a;
	int	c, s;

	set_trig ();

	a = DEG (30);
	s = SIN (a);
	c = COS (a);

#define FFONE	((float)FONE)	/* just for nice printing */
#define FD90	(90.0/FONE)

	printf ("sin(30)=%g cos(30)=%g\n", s/FFONE, c/FFONE);
	printf ("asin(s)=%g acos(c)=%g\n", ASIN(s)*FD90, ACOS(c)*FD90);
	printf ("atan(s/c)=%g\n", ATAN(s, c)*FD90);
	printf ("atan(c/s)=%g\n", ATAN(c, s)*FD90);
	printf ("atan(1/1)=%g\n", ATAN(1, 1)*FD90);

	exit (0);
}