4th-order polynomial cosine approximation

isin.c

#include <stdint.h>
#include <stdbool.h>
#include <math.h>

// 4th-order polynomial cosine approximation
// Derived from
// http://web.archive.org/web/20210224233553/https://www.coranac.com/2009/07/sines/
// Changes:
// - constants split out into variables
// - 32-bits-per-turn angles as input
// - adjustable fixed-point output
// - first-class cosine
int32_t icos(uint32_t angle) {
    // approximation only covers "middle" two quadrants, from -90deg to +90deg
    bool flip = angle >= 0x40000000 && angle < 0xC0000000;
    if (flip)
        angle += 0x80000000;
    int32_t y = *(int32_t*)&angle;

    const uint8_t fractional_digits = 24;

    // divide by a quarter turn, then square
    // same as square, then divide by 1/16th turn
    int32_t z2 = ((int64_t)y*(int64_t)y)>>(28+32-fractional_digits);
    
    // these constants are not actually 64-bit, but it saves us from casting
    const int64_t one = 1<<fractional_digits;
    // naive coefficients
    //const int64_t pi_4 = (one>>2)*M_PI;
    //const int64_t b = 2*one - pi_4;
    //const int64_t c = one - pi_4;
    // coefficients optimized for average error
    const int64_t three_over_pi = one*3/M_PI;
    const int64_t c = 5*(one-three_over_pi);
    const int64_t b = c + one;
    
    int32_t C_4;
    C_4 = (z2*c)>>fractional_digits;
    C_4 = b - C_4;
    C_4 = (z2*(int64_t)C_4)>>fractional_digits;
    C_4 = one - C_4;
    if (flip)
        return -1*C_4;
    else
        return C_4;
}

int32_t isin(uint32_t angle) {
    return icos(angle-0x40000000);
}

#include <stdio.h>
int main(int argc, char** argv) {
    for (int i = 0; i < 256; i++) {
        int fixedpoint_result = icos(i<<24);
        int floatingpoint_result = cos(i*M_PI/128.0)*0x1000000;
        printf("% 3d 0x%08X 0x%08X %d\n", i, fixedpoint_result, floatingpoint_result, fixedpoint_result-floatingpoint_result);
    }
    for (int i = 0; i < 256; i++) {
        int fixedpoint_result = isin(i<<24);
        int floatingpoint_result = sin(i*M_PI/128.0)*0x1000000;
        printf("% 3d 0x%08X 0x%08X %d\n", i, fixedpoint_result, floatingpoint_result, fixedpoint_result-floatingpoint_result);
    }
    return 0;
}