fast_filter_trig.c

typedef float mp_float_t;

typedef struct { mp_float_t c, s; } sincos_result;

// Compute sin and cos of (x+1)*pi/4, -1 < x < 1 to about 4.5 places
/* polynomials from
```py
import numpy as np
from numpy.polynomial.polynomial import Polynomial

x = np.linspace(0, 1, 10000)   

y = np.cos(x * np.pi / 2)
cos_fit = Polynomial.fit(x, y, deg=5)
cos_fit = cos_fit - cos_fit(1)
cos_fit = cos_fit / cos_fit(0)
print(cos_fit.coef)
```
*/
void fast_offset_scaled_sincos( mp_float_t x, sincos_result *result) {
    mp_float_t x2 = x*x, x3 = x2*x, x4 = x2*x2, x5 = x2*x3;
    mp_float_t c0 = 0.70708592,
        c1x = -0.55535724 * x,
        c2x2 = -0.21798592 * x2,
        c3x3 = 0.05707685 * x3,
        c4x4 = 0.0109 * x4,
        c5x5 = -0.00171961 * x5;

    mp_float_t evens = c4x4 + c2x2 + c0, odds = c5x5 + c3x3 + c1x;
    result->c = evens + odds;
    result->s = evens - odds;
}

#if defined(STANDALONE)
#include <math.h>
#include <stdio.h>
int main() {
    for(int i=0; i<=1000; i++) {
        mp_float_t arg = (i / 1000. - .5) * 2;
        mp_float_t theta = (arg + 1) * M_PI / 4;

        sincos_result result;
        fast_offset_scaled_sincos(arg, &result);
        mp_float_t p = cos(theta);
        printf("%f %f %f %f %f\n", arg, theta, result.c, p, result.c - p);
    }
}
#endif