Wanted: fast, low-precision trig functions for C

[Leo de Wit]

In article <1274 at dukeac.UUCP> tcamp at 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)