Xref: utzoo misc.wanted:4170 comp.lang.c:16889
Path: utzoo!attcan!uunet!mcvax!hp4nl!philmds!leo
From: leo@philmds.UUCP (Leo de Wit)
Newsgroups: misc.wanted,comp.lang.c
Subject: Re: Wanted: fast, low-precision trig functions for C
Message-ID: <974@philmds.UUCP>
Date: 12 Mar 89 15:08:34 GMT
References: <1256@dukeac.UUCP> <1989Feb28.170925.17264@utzoo.uucp> <1274@dukeac.UUCP> <973@philmds.UUCP>
Reply-To: leo@philmds.UUCP (Leo de Wit)
Organization: Philips I&E DTS Eindhoven
Lines: 108

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.