Xref: utzoo misc.wanted:4143 comp.lang.c:16815
Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!ames!xanth!mcnc!ecsvax!dukeac!tcamp
From: tcamp@dukeac.UUCP (Ted A. Campbell)
Newsgroups: misc.wanted,comp.lang.c,triangle.general,micro.general,triangle.graphics,micro.ibm
Subject: Re: Wanted: fast, low-precision trig functions for C
Summary: Responses (source code included)
Message-ID: <1274@dukeac.UUCP>
Date: 8 Mar 89 17:51:33 GMT
References: <1256@dukeac.UUCP> <1989Feb28.170925.17264@utzoo.uucp>
Organization: Academic Computing, Duke University, Durham, NC
Lines: 220



Thanks to all who responded to my query about high-speed, low
precision trig functions for use in graphics applications.  
I received numerous replies, along the following lines: 

(a) Most suggested the development of fast look-up tables, 
perhaps with some variation of interpolation between near values
to increase precision.  

(b) Some suggested more elaborate numerical algorithms which, 
alas, I don't understand.  

(c) Some suggest the purchase of a co-processor, which perhaps 
someday I shall do.  

The attached code is based on a suggestion for a lookup table and 
sample code sent to me from Fraser Orr at the Department of Computer
Science, University of Glasgow.  (If he is listening--THANKS--I doubt
that I could find a path to his machine.)  The code below allows 
me to use the functions vpt_sin(), vpt_cos(), and vpt_tan() at either
low precision (high speed) or high precision (low speed).  Note:

(a) The global variable vpt_level selects precision level:  if set to 1,
high-speed look-up tables are used; if set to any other values, the 
routines will convert to radians and utilize the normal sin(), cos(), 
and tan() functions.  

(b) The function vpt_init() must be called before any of these routines 
are used at level 1.  

(c) I have not utilized interpolation here.  If those who understand 
mathematics can give me a model for an interpolation function, 
we could incorporate it as a second level between 1 and the slower 
routines.  

(d) These routines have been tested with AT&T Unix C compiler (7300) 
and with the Microsoft QuickC (tm) compiler.  

Thanks again to all who responded.  

Ted A. Campbell
Duke Divinity School 
{ethos, ecsgate}!dukeac!numen!tcamp

/****************************************************************

	VPT.C		Variable Precision Trigonometric Functions

****************************************************************/

#include "math.h"

extern double vpt_sin(), vpt_cos(), vpt_tan();

#define	DEG_RAD	1.745329252e-2
#define	RAD_DEG	5.729577951e1

double 	sin_table[ 91 ];
int	vpt_level = 1;

/* #define TEST */

#ifdef	TEST

#include "stdio.h"

char test_buf[ 64 ];

main()
	{
	int cont;
	static double x;
	vpt_init();
	cont = 1;
	while ( cont == 1 )
		{
		printf( "Enter precision level (1 or 2 ): " );
		fflush( stdout );
		gets( test_buf );
		if ( test_buf[ 0 ] == 0 )
			{
			return 0;
			}
		sscanf( test_buf, "%d", &vpt_level );
		printf( "Enter a number: " );
		fflush( stdout );
		gets( test_buf );
		if ( test_buf[ 0 ] == 0 )
			{
			return 0;
			}
		sscanf( test_buf, "%lf", &x );
		printf( "sin = %f, cos = %f, tan = %f \n", 
			vpt_sin( x ),
			vpt_cos( x ), 
			vpt_tan( x ) );
		}
	}

#endif

vpt_init()
	{	
	int i;
	for ( i = 0; i < 91; i++ )
		{
		sin_table[ i ] = sin( (double) DEG_RAD * i );
		}
	}

double 
vpt_sin( i )
	double i;
	{
	int sign, target, work;
	switch( vpt_level )
		{
		case 1:
			work = i;
			while ( work < 0 )
				{
				work += 360;
				}
			work = work % 360;
			if ( ( work >= 0 ) && ( work < 90 ))
				{
				sign = 1;
				target = work;
				}
			else if ( ( work >= 90 ) && ( work < 180 ))
				{
				sign = 1;
				target = 90 - ( work % 90 );
				}
			else if ( ( work >= 180 ) && ( work < 270 ))
				{
				sign = -1;
				target = work % 90;
				}
			else if ( ( work >= 270 ) && ( work < 360 ))
				{
				sign = -1;
				target = 90 - ( work % 90 );
				}
			else
				{
				return 0;
				}
			return sign * sin_table[ target ];
		default:
			return sin( DEG_RAD * i );
		}
	}

double
vpt_cos( i )
	double i;
	{
	int sign, target, work;
	switch( vpt_level )
		{
		case 1:
			work = i;
			while ( work < 0 )
				{
				work += 360;
				}
			work = work % 360;
			if ( ( work >= 0 ) && ( work < 90 ))
				{
				sign   = 1;
				target = 90 - work;
				}
			else if ( ( work >= 90 ) && ( work < 180 ))
				{
				sign   = -1;
				target = work % 90;
				}
			else if ( ( work >= 180 ) && ( work < 270 ))
				{
				sign   = -1;
				target = 90 - ( work % 90 );
				}
			else if ( ( work >= 270 ) && ( work < 360 ))
				{
				sign   = 1;
				target = work % 90;
				}
			else
				{
				return 0;
				}
			return sign * sin_table[ target ];
		default: 
			return cos( DEG_RAD * i );
		}
	}

double
vpt_tan( i )
	double i;
	{
	double c;
	switch( vpt_level )
		{
		case 1:
			c = vpt_cos( i );
			if ( c == 0 )
				{
				return 0;
				}
			else
				{
				return vpt_sin( i ) / vpt_cos( i );
				}
		default: 
			return tan( DEG_RAD * i );
		}
	}