tinytrig.c

CD_高级PALM编程

/********************************************************************
* FILE:				TinyTrig.c
*
* DESCRIPTION:	Various modest-resolution trignometric functions.
*
* VERSION:		1.0
**********************************************************************/

#undef COUNTRY	// Because of OS 3.1 & OS 3.5 includes

#include <PalmOS.h>
#include <PalmCompatibility.h>
#include 	"TinyTrig.h"	

/* 
* Lookup tables for radians, sine, cosine, and tangent values.
* The values are stored as 32-bit integers, with a scaling of 1.oe+5, to
* reduce memory requirements and speed up computations.
* There is basically one entry per degree, between 0 and 90 degrees, but 
* because the radians values are more sensitive at the beginning and end,
* there are extra values in all of the tables. These tables are used strictly
* for computing the inverse functions asin, acos, and atan. There are
* approximation equations for computing sin, cos, and tan values.
*/

// Radians table

	UInt32 radians [ TRIG_TAB_SIZE ] = {
	
		0,            349, 	  698,	    1047,      1396, 			// 0    - .8 in 20ths
		1745,     2094,     2443,      2793,      3142,			// 1    - 1.8 in 20ths
		3491,     5236,     6981,      8727,  					// 3    - 5
		10472,   12217,   13963,   15708,    17453,  		// 6    - 10
		19199,   20944,   22689,   24435,    26180,  		// 11   - 15
		27925,   29671,   31416,   33161,    34907,  		// 16   - 20
		36652,   38397,   40143,   41888,    43633,  		// 21   - 25
		45379,   47124,   48869,   50615,    52360,  		// 26   - 30
		54105,   55850,   57596,   59341,    61086,			// 31   - 35
		62832,   64577,   66322,   68068,    69813,  		// 36   - 40
		71558,   73304,   75049,   76794,    78540,  		// 41   - 45
		80285,   82030,   83776,   85521,    87266,  		// 46   - 50
		89012,   90757,   92502,   94248,    95993, 		// 51   - 55
		97738,   99484,   101229, 102974, 104720, 		// 56   - 60
		106465, 108210, 109956, 111701, 113446, 		// 61   - 65
		115192, 116937, 118682, 120428, 122173, 		// 66   - 70
		123918, 125664, 127409, 129154, 130900, 		// 71   - 75
		132645, 134390, 136136, 137881, 139626, 		// 76   - 80
		141372, 143117, 144862, 146608, 148353, 		// 81   - 85
		150098, 151844, 153589, 							// 86   - 88
		153938, 154287, 154636, 154985, 155334,		// 88.2 - 89 in 20ths
		155683, 156032, 156381, 156730, 157079		// 89.2 - 90 in 20th	
	};

// sin table

	UInt32 sines [ TRIG_TAB_SIZE ] = {
	
	        0,          349,      698,      1047,     1396, 				// 0    - .8 in 20ths
	    	1745,    2094,    2443,    2792,     3141,				// 1    - 1.8 in 20ths
	    	3490,    5234,    6976,     8716, 						// 2    - 5
	   	10453,  12187,  13917,  15643,  17365, 			// 6    - 10
	   	19081,  20791,  22495,  24192,  25882, 			// 11   - 15
	   	27564,  29237,  30902,  32557,  34202, 			// 16   - 20
	   	35837,  37461,  39073,  40674,  42262, 			// 21   - 25
	   	43837,  45399,  46947,  48481,  50000, 			// 26   - 30
	   	51504,  52992,  54464,  55919,  57358, 			// 31   - 35
		58778,  60181,  61566,  62932,  64279, 			// 36   - 40
		65606,  66913,  68200,  69466,  70711, 			// 41   - 45
		71934,  73135,  74314,  75471,  76604, 			// 46   - 50			
		77715,  78801,  79864,  80902,  81915, 			// 51   - 55
		82904,  83867,  84805,  85717,  86602, 			// 56   - 60
		87462,  88295,  89101,  89879,  90631, 			// 61   - 65
		91355,  92050,  92718,  93358,  93969, 			// 66   - 70
		94552,  95106,  95630,  96126,  96593, 			// 71   - 75
		97030,  97437,  97815,  98163,  98481, 			// 76   - 80
		98769,  99027,  99255,  99452,  99619,  			// 81   - 85
		99756,  99863,  99939,  							// 86   - 88
		99951,  99961,  99970,  99978,  99985, 			// 88.2 - 89 in 20ths
		99990,  99995,  99998,  99999, 100000			// 89.2 - 90 in 20ths
	};

// cos table

	UInt32 cosines [ TRIG_TAB_SIZE ] = {
	
		100000,  99999,  99998,  99995,  99990, 			// 0    - .8 in 20ths
		99985,  99978,  99970,  99961,  99951,			// 1    - 1.8 in 20ths
		99939,  99863,  99756,  99619, 					// 3    - 5
		99452,  99255,  99027,  98769,  98481, 			// 6    - 10
		98163,  97815,  97437,  97030,  96593, 			// 11   - 15
		96126,  95631,  95106,  94552,  93969, 			// 16   - 20
		93358,  92718,  92051,  91355,  90631, 			// 21   - 25
		89879,  89101,  88295,  87462,  86603, 			// 26   - 30
		85717,  84805,  83867,  82904,  81915, 			// 31   - 35
		80902,  79864,  78801,  77715,  76604, 			// 36   - 40
		75471,  74315,  73135,  71934,  70711, 			// 41   - 45
		69466,  68200,  66913,  65606,  64279, 			// 46   - 50			
		62932,  61566,  60182,  58779,  57358, 			// 51   - 55
		55919,  54464,  52992,  51504,  50000, 			// 56   - 60
		48481,  46947,  45399,  43837,  42262, 			// 61   - 65
		40674,  39073,  37461,  35837,  34202, 			// 66   - 70
		32557,  30902,  29237,  27564,  25882, 			// 71   - 75
		24192,  22495,  20791,  19081,  17365, 			// 76   - 80
		15644,  13917,  12187,  10453,   8716,  			// 81   - 85
		 6976,   5234,   3490,  								// 86   - 88
		 3141,   2792,   2443,   2094,   1745, 				// 88.2 - 89 in 20ths
		 1396,   1047,    698,    349,      0					// 89.2 - 90 in 20ths
	};

// tan table

	UInt32  tangents [ TRIG_TAB_SIZE ] = {
	
		0,               349,          698,          1047,           1396, 		// 0    - .8 in 20ths
		1746,	    2095,        2444,         2793,          3143,		// 1    - 1.8 in 20ths
		3492,	    5241,        6993,         8749,  					// 3    - 5
		10510,     12278,      14054,      15838,        17633,  		// 6    - 10
		19438,     21256,      23087,      24933,        26795,  		// 11   - 15
		28675,     30573,      32492,      34433,        36397,  		// 16   - 20
		38386,     40403,      42447,      44523,        46631,  		// 21   - 25
		48773,     50952,      53171,      55431,        57735,  		// 26   - 30
		60086,     62487,      64941,      67451,        70021,		// 31   - 35
		72654,     75355,      78128,      80978,        83910,  		// 36   - 40
		86929,     90040,      93251,      96569,        100000,  	// 41   - 45
		103553,   107237,    111061,    115037,     119175,  	// 46   - 50
		123490,   127994,    132704,    137638,     142815, 	// 51   - 55
	 	148256,   153986,    160033,    166428,     173205, 	// 56   - 60
	 	180404,   188072,    196261,    205030,     214450, 	// 61   - 65
	 	224603,   235585,    247508,    260508,     274747, 	// 66   - 70
	 	290420,   307767,    327084,    348740,     373203, 	// 71   - 75
	 	401076,   433145,    470460,    514452,      567124, 	// 76   - 80
	 	631370,   711531,    814426,    951424,      1142990, 	// 81   - 85
		1430040, 1908060,  2863530, 						// 86   - 88
		3181900, 3579910,  4091510,  4773610,   5728550,	// 88.2 - 89 in 20ths
		7160780, 9547760, 14320900, 28637600, 55200000	// 89.2 - 90 in 20th	
	};


/*****************************************************
 * FUNCTION:    	TTLookup
 *
 * DESCRIPTION: 	Find the two consecutive values in a  table
 *					that bound x, and then interpolate another
 *					value using a second table. 
 *
 * RETURNED:    	The result.
*****************************************************/
static double TTLookup // ( out ) the interpolated result
(
	double 	x, 	// ( in ) the value to look up
	UInt32 	lkupTab	[ TRIG_TAB_SIZE ],	// ( in ) lookup table to use
	UInt32		rsltTab	[ TRIG_TAB_SIZE ] 	// ( in ) the result table to use
)
{
	UInt8 		tabloc = 0, increment = 1;
	double	result = -1.0, lkup1, lkup2, rslt1, rslt2, term1, term2, term3;
	UInt32		lval;

// Convert the incoming double to a UInt32 for table lookup
	
	lval = ( UInt32 ) ( SCALE_FACTOR * x ); 
	
// Find the trig table entries that bound the value being searched for

	tabloc = SearchTable ( lkupTab, lval );
	
// Get the interpolation variables in doubles to make the math more accurate.
				
	rslt1 = (( double ) rsltTab [ tabloc ] ) / SCALE_FACTOR;
	rslt2 = (( double ) rsltTab [ tabloc + increment ] ) / SCALE_FACTOR;
			
	lkup1 = (( double ) lkupTab   [ tabloc ] ) / SCALE_FACTOR;
	lkup2 = (( double ) lkupTab   [ tabloc + increment ] ) / SCALE_FACTOR;

// Do the interpolation and get the answer in radians. 
// Use the interpolation formula:
//
//		result = result  +  abs ( lkup1 - Ival )  *  abs ( rslt1 - rslt2 )
//				  --------------------------------------------
//							            abs ( lkup1 - lkup2 )
				
	term1 = _abs ( lkup1 - ( double ) x );
	term2 = _abs ( rslt1 - rslt2 );
	term3 = _abs ( lkup1 - lkup2 );		
	result = ( term1 * term2 ) / term3;
	
// Add or subtract the adjustment, depending on the lookup
// table's direction--ascending or descending.

	if ( rsltTab [ 0 ] < rsltTab [ 1 ] )
		result = rslt1 + result;
	else
		result = rslt1 - result;
	
	return ( result );
}


/*****************************************************
 * FUNCTION:    	SearchTable
 *
 * DESCRIPTION: 	Binary search a table of size TRIG_TAB_SIZE for
 *					value. Return a pointer to the lower of two 
 *					table positions that bound it. 
 *
 * RETURNED:    	An index into the table of the lower boundary position.
*****************************************************/
static UInt8 SearchTable	// ( out ) the found location index
(
	UInt32		table [ TRIG_TAB_SIZE ], 	// ( in ) the table to search
	UInt32		value	// ( in ) the value to look for
)
{
	UInt8 		mid, lower = 0, upper = TRIG_TAB_SIZE - 1;
	Int8		increment = 1;
	UInt32 	val1, val2;
	Boolean	found = false;

// Determine if ascending or descending table
	
	if ( table [ 0 ] > table [ 1 ] ) 
		increment = -1;

// Binary search thru table for a pair of bounding values
	
	while (! found )
	{

// Get next search location
	
		mid = ( UInt8 ) (( lower + upper ) / 2 );	
		val1 = table [ mid ];
		val2 = table [ mid + increment ];

// Check for bounds
		
		if (( value >= val1 ) && ( value <= val2 )) 
			found = true;

// Not found, move search location
		
		else
		{
			if ( increment == 1 )	// ascending table
			{
				if ( value < val1 ) 
					upper = mid;		
				else
					lower = mid;
			}
			else // descending table
			{					
				if ( value < val1 ) 
					lower = mid;		
				else 
					upper = mid;
			}
		}		
	}
	
// Return the proper pointer to the lower of the two values

	if ( increment == 1 ) 
		return  mid;
	else
		return (( UInt8 ) ( mid - 1 ));
}


/*****************************************************
 * FUNCTION:    	ScaleNearZero
 *
 * DESCRIPTION: 	Scale a radians number that is outside the 
 *					range - scale <--> scale to something within that range.
 *					Assume scale is either PI or 2PI
 *
 * RETURNED:    	The scaled result.
*****************************************************/
static double		ScaleNearZero // ( out ) the scaled result
(
	double	x,		// ( in ) the number to scale
	double	scale	// ( in ) the range to scale it to
)
{
	double 	adjCntD, result;
	UInt32		adjCnt;

// Get the number of scale increments in the input value.

	adjCntD 	= scale;
	adjCntD 	= x / adjCntD;
	adjCntD 	= _abs ( adjCntD );
	adjCnt 	= ( UInt32 ) adjCntD;
	
// Adjust increment count for scale factor == PI

	if ( scale == PI )
		adjCnt = adjCnt * 2;
	
// Adjust the input value accordingly.

	if ( x > 0.0 ) 
		result 	= x - ( scale * ( double ) adjCnt );
	else
		result 	= x + ( scale * ( double ) adjCnt );
		
	return ( result );
}

/***********************************************************************
 * FUNCTION:    	_abs
 *
 * DESCRIPTION: 	Calculate the absolute value of a double number.
 *
 * RETURNED:    	The result.
 ***********************************************************************/
double _abs		// ( out ) the result
( 
	double	x	// ( in ) The number to get the absolute value of.
)
{
	if ( x < 0.0 )
		return ( - x );
	else
		return ( x );
}


/***********************************************************************
 * FUNCTION:    	_sqrt
 *
 * DESCRIPTION: 	Calculate the square root of a double number using
 *					successive approximations. The loop stops when the 
 *					difference between two successive iterations is smaller
 *					than the Tiny Trig accuracy threshold.
 *
 * RETURNED:    	The square root or zero if the number is smaller
 *					than the minimum value handled by the library.
 ***********************************************************************/
double _sqrt		// ( out ) the square root
( 
	double	x	// ( in ) The number to get the square root of
)
{
	double 	result = 0.0, prevResult, diff = 0.0;

// Eliminate negatives and small values outside the library range.       
	
	if ( x >= TT_MIN_VALUE )
	{	
		result 			= x;	

// Calculate until two successive guesses are less than
// the library's minimum value. # iterations ranges from about
// 10 (for very small and very large #s) to about 5 for numbers
// close to zero.

		do
		{
			prevResult 	= result;
			result 			= ( result + ( x / result )) / 2.0;
			diff				= _abs ( prevResult - result );
		}
		while ( diff > TT_ACCURACY );
	}	
	return result;
}


/*****************************************************
 * FUNCTION:    	_sin
 *
 * DESCRIPTION: 	Calculate the sine of a double radians number
 *					using  the power series sin x = x - 
 *					( x**3 / 3! ) + ( x**5 / 5! ) - ( x**7 / 7! ) + ...
 *
 * RETURNED:    	The result.
*****************************************************/
double _sin		// ( out ) the result
(
	double x		// ( in ) the number to get the sine of
)
{