#include "math.h"

__attribute__((always_inline))
static float fminf(float a, float b) {
    return a < b ? a : b;
}

/* Approach 1: polynomial; input in radians. */
static float sin_poly_radians(float x)
{
    x *= 0.159155f;          /* approximately 1 / (2*pi) */
    x -= __builtin_floorf(x);          /* reduce to [0, 1) */
    float xx = x * x;
    float y = -6.87897f;
    y = y * xx + 33.7755f;
    y = y * xx - 72.5257f;
    y = y * xx + 80.5874f;
    y = y * xx - 41.2408f;
    y = y * xx + 6.28077f;
    return x * y;
}

// Approach 2: rational approximation; input in turns (1.0 = 2*pi radians)
// Calculate sin(tau x) (where x is in the range [0, 1]) using Bhaskara I's approximation
// https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximation_formula
// Accurate to about 0.2%
static float sin_rational_tau(float x)
{
    float frac = x - __builtin_floorf(x);
    float sign = 0.5f - frac;
    float y = 2.0f * fminf(frac, 1.0f - frac);
    float q = y * (1.0f - y);
    float approximation = (16.0f * q) / (5.0f - 4.0f * q);
    return __builtin_copysignf(approximation, sign);
}

/* Approach 3: folded polynomial; input in turns (1.0 = 2*pi radians). */
static float sin_folded_poly_tau(float x)
{
    union { uint32_t i; float f; } coeffs[6] = {
        { 0x4261c73fu }, { 0x4113e0d9u }, { 0xc22aaa4bu },
        { 0x3da6b273u }, { 0x40c90140u }, { 0x36d0957eu }
    };

    float frac = x - __builtin_floorf(x);
    float x1 = __builtin_fabsf(frac - 0.5f);
    float t = fminf(0.5f - x1, x1); /* fold to [0, 0.25] */

    float y = 0.0f;
    for (int i = 0; i < 6; ++i) {
        y = y * t + coeffs[i].f;
    }

    return __builtin_copysignf(y, 0.5f - frac);
}

float sintau(float x) {
    return sin_rational_tau(x);
}