tcamp@dukeac.UUCP (Ted A. Campbell) (02/28/89)
I am writing a graphics-based program that has to calculate
several thousand map coordinates. The slowest part of this
program by far is its trigonometric functions--it takes
3 minutes and 45 seconds to calculate 3000 sets ofd coordinates
on an XT clone at 10 mhz (no coprocessor), which is the lowest
resolution map I can make.
I wondered if there were any fast, low-precision approaches
to the standard trig functions, since the compilers I use
(Microsoft QuickC and AT&T Unix C) boast of their high
level of precision?
Any suggestions will be appreciated. Thanks.
Ted A. Campbell
Duke Divinity School
{ethos, ecsgate}!dukeac!numen!tcampbill@twwells.uucp (T. William Wells) (02/28/89)
In article <1256@dukeac.UUCP> tcamp@dukeac.UUCP (Ted A. Campbell) writes:
: I wondered if there were any fast, low-precision approaches
: to the standard trig functions, since the compilers I use
: (Microsoft QuickC and AT&T Unix C) boast of their high
: level of precision?
Try table lookup and linear interpolation.
---
Bill
{ uunet | novavax } !twwells!billhenry@utzoo.uucp (Henry Spencer) (03/01/89)
In article <1256@dukeac.UUCP> tcamp@dukeac.UUCP (Ted A. Campbell) writes: >...I wondered if there were any fast, low-precision approaches >to the standard trig functions... Consider table lookup, followed (if necessary) by linear interpolation between adjacent table values (most functions look fairly linear on a small scale). -- The Earth is our mother; | Henry Spencer at U of Toronto Zoology our nine months are up. | uunet!attcan!utzoo!henry henry@zoo.toronto.edu
tcamp@dukeac.UUCP (Ted A. Campbell) (03/09/89)
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 );
}
}leo@philmds.UUCP (Leo de Wit) (03/11/89)
In article <1274@dukeac.UUCP> tcamp@dukeac.UUCP (Ted A. Campbell) writes: | | |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: [] |(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. Simple. Given x values x0, x1 (x1 >= x0) and their corresponding y values f(x0) = y0, f(x1) = y1, linear interpolation between these two points just says: replace the curve f(x) by the straight line between (x0,y0) and (x1,y1) for x values in the interval [x0,x1]. l(x) - y0 y1 - y0 --------- = ---------- x - x0 x1 - x0 or l(x) = y0 + (x - x0) * (y1 - y0) / (x1 - x0) where l(x) is a linear approximation for f(x) on the interval [x0,x1]. Alternatively one can use linear approximation using the derivative (?) of the functions at hand, since they are easy available: for sin(x) it is cos(x) and for cos(x) it is -sin(x). This is what I used in the example below. This is also called a Taylor series: f(x) = f(x0) + (x - x0) * f'(x0) + .... |(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. I hope you don't mind me critisizing (part of) this code, any criticism on my code would also be appreciated. [] |vpt_init() | { | int i; | for ( i = 0; i < 91; i++ ) | { | sin_table[ i ] = sin( (double) DEG_RAD * i ); | } | } Since we're talking approximations here, there is no need to calculate all sinuses for your table using sin(). The example I'll give shows how it can be done. |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; This can take a loooooong time if work (i if you want) is a large negative number. work == 3600000 for instance would take 10000 iterations. I would suggest in this case work = 360 - (-work % 360). Note also that for initially negative numbers work = work % 360 is not needed (because work will already be >= 0 and < 360). Another point is that C float to int conversion does truncation instead of rounding (and if I understand the dpANS correctly, is now no longer implementation defined for negative numbers). So you should preferably use work = i + 0.5 for positive i and work = i - 0.5 for negative i. [] | else if ( ( work >= 90 ) && ( work < 180 )) | { | sign = 1; | target = 90 - ( work % 90 ); | } Since work is >= 90 and < 180, work % 90 will amount to work - 90, which is perhaps a preferable operation (subtraction is a more basic operation for many processors), especially since this leads to target = 180 - work; There are many other examples where this could be used in the code. | return sign * sin_table[ target ]; I'd suggest you use a boolean for sign, and avoid doing a multiplication: return (sign) ? -sin_table[target] : sin_table[target]; And below is how I would do it (perhaps). sintab[] and costab[] I would probably generate by a C program, thus avoiding having to do the initialization. In this case however, init_tab() initiates the sin and cos tables, using the gonio rules: sin(a+b) = sin(a) * cos(b) + cos(a) * sin(b); cos(a+b) = cos(a) * cos(b) - sin(a) * sin(b); Although these may lead to cumulative errors, these errors are quite small in respect to any approximations. The dsin() and dcos() macros are the simple, table-driven macros, the macros dsin2() and dcos2() are a first order approximation (the error in them is typically a factor 100 smaller). Perhaps dsin2() and dcos2() should better be functions, thus avoiding to recalculate x % 360. Be warned! dsin(), dcos, dsin2(), and dcos2() all use their argument several times, so avoid using side effects or (complex) calculations, or implement them as functions if you must. --------- c u t h e r e ---------- #include <math.h> /* One degree in radians */ #define ONE_DEG (M_PI / 180) /* sin, cos for pos. x and neg. x (x in degrees), to be used by the macros * that follow. */ #define dsinp(x) sintab[(int)((x) + 0.5) % 180] #define dsinn(x) -sintab[(int)(-(x) + 0.5) % 180] #define dcosp(x) costab[(int)((x) + 0.5) % 180] #define dcosn(x) costab[(int)(-(x) + 0.5) % 180] /* Zero order approximation: simple table lookup (x in degrees) */ #define dsin(x) (((x) >= 0) ? dsinp(x) : dsinn(x)) #define dcos(x) (((x) >= 0) ? dcosp(x) : dcosn(x)) /* First order approximation (x in degrees) */ #define dsin2(x) (((x) >= 0) \ ? (dsinp(x) + ONE_DEG * ((x) - (int)((x) + 0.5)) * dcosp(x)) \ : (dsinn(x) + ONE_DEG * ((x) - (int)((x) - 0.5)) * dcosn(x))) #define dcos2(x) (((x) >= 0) \ ? (dcosp(x) + ONE_DEG * ((x) - (int)((x) + 0.5)) * dsinp(x)) \ : (dcosn(x) - ONE_DEG * ((x) - (int)((x) - 0.5)) * dsinn(x))) static double sintab[180], costab[180]; static void init_tab(), demo(); main() { int i; init_tab(); demo(); } /* Shows abberation for dsin(), dcos(), dsin2() and dcos2() in the interval * [-56,-55] */ static void demo() { double d; for (d = -55.; d >= -56.; d = d - 0.1) { printf("%g: %g %g\n", d,sin(d * ONE_DEG) - dsin(d), cos(d * ONE_DEG) - dcos(d)); printf("%g: %g %g\n", d,sin(d * ONE_DEG) - dsin2(d), cos(d * ONE_DEG) - dcos2(d)); } } /* This builds sintab,costab with errors of magnitude e-16, which can * be ignored without problem (the approximations introduce much * larger errors). */ static void init_tab() { int i; double sin1 = sin(ONE_DEG), cos1 = cos(ONE_DEG); sintab[0] = 0.; costab[0] = 1.; for (i = 1; i <= 45; i++) { sintab[i] = sintab[180 - i] = costab[90 - i] = cos1 * sintab[i - 1] + sin1 * costab[i - 1]; costab[90 + i] = -sintab[i]; costab[i] = sintab[90 + i] = sintab[90 - i] = cos1 * costab[i - 1] - sin1 * sintab[i - 1]; costab[180 - i] = -costab[i]; } } ------- e n d s h e r e ---------- Leo. B.T.W. News had objections against some of the newsgroups in the list, so I had to strip it a bit to get it out. Well ... (Newsgroups was: misc.wanted,comp.lang.c,triangle.general,micro.general,triangle.graphics,micro.ibm)
leo@philmds.UUCP (Leo de Wit) (03/12/89)
Sorry for the followup, but my original posting contained an error with
respect to the sin / cos periods. The following one should be ok.
And below is how I would do it (perhaps). sintab[] and costab[] I would
probably generate by a C program, thus avoiding having to do the
initialization. In this case however, init_tab() initiates the sin and
cos tables, using the gonio rules:
sin(a+b) = sin(a) * cos(b) + cos(a) * sin(b);
cos(a+b) = cos(a) * cos(b) - sin(a) * sin(b);
Although these may lead to cumulative errors, these errors are quite small
in respect to the approximations used.
The dsin() and dcos() macros are the simple, table-driven macros, the
macros dsin2() and dcos2() are a first order approximation (the error in
them is typically a factor 100 smaller).
Perhaps dsin2() and dcos2() should better be functions, thus avoiding to
recalculate x % 360.
Be warned! dsin(), dcos, dsin2(), and dcos2() all use their argument
several times, so avoid using side effects or (complex) calculations,
or implement them as functions if you must.
--------- c u t h e r e ----------
#include <math.h>
/* One degree in radians */
#define ONE_DEG (M_PI / 180)
/* sin, cos for pos. x and neg. x (x in degrees), to be used by the macros
* that follow.
*/
#define dsinp(x) sintab[(int)((x) + 0.5) % 360]
#define dsinn(x) -sintab[(int)(-(x) + 0.5) % 360]
#define dcosp(x) costab[(int)((x) + 0.5) % 360]
#define dcosn(x) costab[(int)(-(x) + 0.5) % 360]
/* Zero order approximation: simple table lookup (x in degrees) */
#define dsin(x) (((x) >= 0) ? dsinp(x) : dsinn(x))
#define dcos(x) (((x) >= 0) ? dcosp(x) : dcosn(x))
/* First order approximation (x in degrees) */
#define dsin2(x) (((x) >= 0) \
? (dsinp(x) + ONE_DEG * ((x) - (int)((x) + 0.5)) * dcosp(x)) \
: (dsinn(x) + ONE_DEG * ((x) - (int)((x) - 0.5)) * dcosn(x)))
#define dcos2(x) (((x) >= 0) \
? (dcosp(x) + ONE_DEG * ((x) - (int)((x) + 0.5)) * dsinp(x)) \
: (dcosn(x) - ONE_DEG * ((x) - (int)((x) - 0.5)) * dsinn(x)))
static double sintab[360], costab[360];
static void init_tab(), demo();
main()
{
init_tab();
demo();
}
/* Shows abberation for dsin(), dcos(), dsin2() and dcos2() in the interval
* [-56,-55]
*/
static void demo()
{
double d;
for (d = -55.; d >= -56.; d = d - 0.1) {
printf("%g: %g %g\n", d,sin(d * ONE_DEG) - dsin(d),
cos(d * ONE_DEG) - dcos(d));
printf("%g: %g %g\n", d,sin(d * ONE_DEG) - dsin2(d),
cos(d * ONE_DEG) - dcos2(d));
}
}
/* This builds sintab,costab with errors of size 10^-16, which can
* be ignored without problem (the approximations introduce much
* larger errors).
*/
static void init_tab()
{
int i;
double sin1 = sin(ONE_DEG), cos1 = cos(ONE_DEG);
sintab[0] = sintab[180] = costab[90] = costab[270] = 0.;
costab[0] = sintab[90] = 1.;
costab[180] = sintab[270] = -1.;
for (i = 1; i <= 45; i++) {
sintab[i] = sintab[180 - i]
= costab[90 - i]
= costab[270 + i]
= cos1 * sintab[i - 1] + sin1 * costab[i - 1];
sintab[180 + i] = sintab[360 - i]
= costab[90 + i]
= costab[270 - i]
= -sintab[i];
costab[i] = costab[360 - i]
= sintab[90 + i]
= sintab[90 - i]
= cos1 * costab[i - 1] - sin1 * sintab[i - 1];
costab[180 + i] = costab[180 - i]
= sintab[270 - i]
= sintab[270 + i]
= -costab[i];
}
}
------- e n d s h e r e ----------
Leo.