Octave (or Matlab) & C

A bit of a complicated build process, but kind of a novel approach and the results were encouraging.

The approach is to generate approximating quadratic polynomials for each degree. So degree = [0, 1), degree = [1, 2), ..., degree = [359, 360) will each have a different polynomial.

Octave - building part

Octave is publicly available - Google download octave.

This determines the best fit quadratic polynomial for every degree.

Save as build-fast-trig.m:

format long;
for d = 0:359
    x = (d-1):0.1:(d+1);
    y = sin(x / 360 * 2 * pi);
    polyfit(x, y, 2)
endfor

C - building part

This converts doubles in text format to native binary format on your system.

Save as build-fast-trig.c:

#include <stdio.h>

int main()
{
    double d[3];

    while (scanf("%lf %lf %lf", d, d + 1, d + 2) == 3)
        fwrite(d, sizeof(double), 3, stdout);

    return 0;
}

Compile:

gcc -o build-fast-trig build-fast-trig.c

Generating the coefficients file

Run:

octave build-fast-trig.m | grep '^ ' | ./build-fast-trig > qcoeffs.dat

Now we have qcoeffs.dat as the data file to use for the actual program.

C - fast-trig part

Save as fast-trig.c:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define INPUT    "qcoeffs.dat"

#define DEGREES    360

typedef struct {double a, b, c;} QCOEFFS;

double normalize(double d)
{
    if (d < 0.0)
        d += ceil(d / -(double)DEGREES) * (double)DEGREES;

    if (d >= (double)DEGREES)
        d -= floor(d / (double)DEGREES) * (double)DEGREES;

    return d;
}

int main()
{
    FILE *f;
    time_t tm;
    double d;
    QCOEFFS qc[DEGREES];

    if (!(f = fopen(INPUT, "rb")) || fread(qc, sizeof(QCOEFFS), DEGREES, f) < DEGREES)
    {
        fprintf(stderr, "Problem with %s - aborting.", INPUT);
        return EXIT_FAILURE;
    }
    fclose(f);

    tm = -clock();

    while (scanf("%lf", &d) > 0)
    {
        int i;
        double e, f;

        /* sin */
        d = normalize(d);
        i = (int)d;
        e = (qc[i].a * d + qc[i].b) * d + qc[i].c;

        /* cos */
        d = normalize((double)DEGREES / 4.0 - d);
        i = (int)d;
        f = (qc[i].a * d + qc[i].b) * d + qc[i].c;

        /* tan = sin / cos */

        /* output - format closest to specs, AFAICT */
        if (d != 0.0 && d != 180.0)
            printf("%.6e %.6e %.6e\n", e, f, e / f);
        else
            printf("%.6e %.6e n\n", e, f);
    }

    tm += clock();

    fprintf(stderr, "time: %.3fs\n", (double)tm/(double)CLOCKS_PER_SEC);    

    return EXIT_SUCCESS;
}

Compile:

gcc -o fast-trig fast-trig.c -lm

Run:

./fast-trig < trig.in > trig.out

It will read from trig.in, save to trig.out and print to console the elapsed time with millisecond precision.

Depending on the testing methods used it may fail on certain input, e.g.:

$ ./fast-trig 
0
-6.194924e-19 1.000000e+00 -6.194924e-19

The correct output should be 0.000000e+00 1.000000e+00 0.000000e+00. If the results are validated using strings, the input will fail, if they are validated using absolute error, e.g. fabs(actual - result) < 1e-06, the input will pass.

The maximum absolute error for sin and cos was ≤ 3e-07. For tan, because the result isn't limited to ± 1 and you can divide a relatively large number by a relatively small number, the absolute error could be larger. From -1 ≤ tan(x) ≤ +1, the maximum absolute error was ≤ 4e-07. For tan(x) > 1 and tan(x) < -1, the maximum relative error, e.g. fabs((actual - result) / actual) was usually < 1e-06 until you get in the area of (90 ± 5) or (270 ± 5) degrees, then the error gets worse.

In testing, the average time per single input was (1.053 ± 0.007) µs, which on my machine was about 0.070 µs faster than native sin and cos, tan being defined the same way.