sin_approx.cpp

/**
 * Code generated from https://thomasahle.com/fast-polynomials/
 * License: 0BSD. Use freely; no attribution required.
 *
 * P(x) = -189823930421/200000000000000*x^6 + 51096785741/5000000000000*x^5 -
 *   38058038929/20000000000000*x^4 - 82835861483/500000000000*x^3 -
 *   122999528541/500000000000000*x^2 + 10000204137/10000000000*x
 * Field: Q, evaluated in double precision.
 * This paper: 5 multiplications (Horner: 6), 5 preprocessed constants.
 * Reference: T. D. Ahle and J. B. T. Knudsen, "Fast Evaluation of Polynomials with
 *   Rational Preprocessing", arXiv:2609.06022, 2026. https://arxiv.org/abs/2609.06022
 * Compile: cc -O3 -march=native <this file>  (MSVC: cl /O2 <this file>)
 *
 * The chain is exact over Q; here the constants are rounded to the nearest double
 * (shortest decimal that round-trips to the correctly rounded double).
 * Optionally compile with -ffp-contract=fast (/fp:contract): it lets the compiler fuse
 * multiply-adds (FMA), which is faster on most cores and changes the rounding slightly.
 */

/* preprocessed constants (the paper's alpha_i) */
static constexpr double P_alpha[5] = {
    336.61213901928437,   // alpha0 = 12129187133850356881807611/36033124560476649237241
    -322.8400306051717,   // alpha1 = -11632935035904245021900840/36033124560476649237241
    -11.767195817234478,  // alpha2 = -2233695360061/189823930421
    112633.11340391554,   // alpha3 = 731206324838775917501415024480502203917129982677457/6491930327954127581370475409866412135506460405
    45559.0589700695,     // alpha4 = 295766236640844717748711468969424757151464475459137/6491930327954127581370475409866412135506460405
};

constexpr double eval_P(double x) {
    /* P_5 base */
    double y = x * x;                                    // y = x * x
    double z = (y + P_alpha[0]) * (x + y + P_alpha[1]);  // z = (y + alpha0) * (x + y + alpha1)
    double t = (x + P_alpha[2]) * (z + P_alpha[3]);      // t = (x + alpha2) * (z + alpha3)
    /* even lift P_6 = x·P_5 + α_0 */
    double u = (t + P_alpha[4]) * x;                     // u = (t + alpha4) * x
    double P = u;                                        // P = u
    return P * -0.000949119652098215;  // leading coefficient
}