C - Taylor Series for Sine / Cosine

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

#define TSLIM 20
#define EMAX 1.0e-8

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950
#endif

/** sin with taylor series expansion to n = TSLIM
 *  (no function reliance, quadrants handled)
 */
double sinenfq (const double deg)
{
    double fp = deg - (int)deg,     /* save fractional part of deg */
        qdeg = (int)deg % 360,      /* get equivalent 0-359 deg angle */
        rad, sine_deg;              /* radians, sine_deg */
    int pos_quad = 1,               /* positive quadrant flag 1,2  */
        sign = -1;                  /* taylor series term sign */

    qdeg += fp;                     /* add fractional part back to angle */

    /* get equivalent 0-90 degree angle, set pos_quad flag */
    if (90 < qdeg && qdeg <= 180)           /* in 2nd quadrant */
        qdeg = 180 - qdeg;
    else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
        qdeg = qdeg - 180;
        pos_quad = 0;
    }
    else if (270 < qdeg && qdeg <= 360) {   /* in 4th quadrant */
        qdeg = 360 - qdeg;
        pos_quad = 0;
    }

    rad = qdeg * M_PI / 180.0;      /* convert to radians */
    sine_deg = rad;                 /* save copy for computation */

    /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
    for (int n = 3; n < TSLIM; n += 2, sign *= -1) {
        double p = rad;
        uint64_t f = n;

        for (int i = 1; i < n; i++) {   /* pow & nfact */
            p *= rad;
            f *= i;
        }

        sine_deg += sign * p / f;       /* Taylor-series term */
    }

    return pos_quad ? sine_deg : -sine_deg;
}

/** cos with taylor series expansion to n = TSLIM
 *  (no function reliance, quadrants handled)
 */
double cosenfq (const double deg)
{
    double fp = deg - (int)deg,     /* save fractional part of deg */
        qdeg = (int)deg % 360,      /* get equivalent 0-359 deg angle */
        rad, cose_deg = 1.0;        /* radians, cose_deg */
    int pos_quad = 1,               /* positive quadrant flag 1,4  */
        sign = -1;                  /* taylor series term sign */

    qdeg += fp;                     /* add fractional part back to angle */

    /* get equivalent 0-90 degree angle, set pos_quad flag */
    if (90 < qdeg && qdeg <= 180) {         /* in 2nd quadrant */
        qdeg = 180 - qdeg;
        pos_quad = 0;
    }
    else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
        qdeg = qdeg - 180;
        pos_quad = 0;
    }
    else if (270 < qdeg && qdeg <= 360)     /* in 4th quadrant */
        qdeg = 360 - qdeg;

    rad = qdeg * M_PI / 180.0;      /* convert to radians */

    /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
    for (int n = 2; n < TSLIM; n += 2, sign *= -1) {
        double p = rad;
        uint64_t f = n;

        for (int i = 1; i < n; i++) {   /* pow & nfact */
            p *= rad;
            f *= i;
        }

        cose_deg += sign * p / f;       /* Taylor-series term */
    }

    return pos_quad ? cose_deg : -cose_deg;
}

int main (int argc, char **argv) {

    double a = (argc > 1) ? strtod (argv[1], NULL) : 187.2,
        coslibc = cos (a * M_PI / 180.0),
        sinlibc = sin (a * M_PI / 180.0),
        sints = sinenfq (a),
        costs = cosenfq (a),
        err = fabs (sinlibc - sints),
        cerr = fabs (coslibc - costs);

    if (err > EMAX)
        fprintf (stderr, "sine error exceeds limit of %e\n"
        "%.2f    %11.8f    %11.8f    %e\n",
        EMAX, a, sinlibc, sints, err);

    if (cerr > EMAX)
        fprintf (stderr, "cose error exceeds limit of %e\n"
        "%.2f    %11.8f    %11.8f    %e\n",
        EMAX, a, coslibc, costs, cerr);

    printf ("angle        libc_sin         sinets       libc_cos         cosets"
            "         sineerr        coseerr\n");
    printf ("%6.2f    %11.8f    %11.8f    %11.8f    %11.8f    %e    %e\n",
            a, sinlibc, sints, coslibc, costs, err, cerr);

    return 0;
}