drankinatty

C - Taylor Series for Sine / Cosine

Oct 25th, 2020 (edited)
1,492
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <stdint.h>
  4. #include <math.h>
  5.  
  6. #define TSLIM 20
  7. #define EMAX 1.0e-8
  8.  
  9. #ifndef M_PI
  10. #define M_PI 3.14159265358979323846264338327950
  11. #endif
  12.  
  13. /** sin with taylor series expansion to n = TSLIM
  14.  *  (no function reliance, quadrants handled)
  15.  */
  16. double sinenfq (const double deg)
  17. {
  18.     double fp = deg - (int)deg,     /* save fractional part of deg */
  19.         qdeg = (int)deg % 360,      /* get equivalent 0-359 deg angle */
  20.         rad, sine_deg;              /* radians, sine_deg */
  21.     int pos_quad = 1,               /* positive quadrant flag 1,2  */
  22.         sign = -1;                  /* taylor series term sign */
  23.  
  24.     qdeg += fp;                     /* add fractional part back to angle */
  25.  
  26.     /* get equivalent 0-90 degree angle, set pos_quad flag */
  27.     if (90 < qdeg && qdeg <= 180)           /* in 2nd quadrant */
  28.         qdeg = 180 - qdeg;
  29.     else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
  30.         qdeg = qdeg - 180;
  31.         pos_quad = 0;
  32.     }
  33.     else if (270 < qdeg && qdeg <= 360) {   /* in 4th quadrant */
  34.         qdeg = 360 - qdeg;
  35.         pos_quad = 0;
  36.     }
  37.  
  38.     rad = qdeg * M_PI / 180.0;      /* convert to radians */
  39.     sine_deg = rad;                 /* save copy for computation */
  40.  
  41.     /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
  42.     for (int n = 3; n < TSLIM; n += 2, sign *= -1) {
  43.         double p = rad;
  44.         uint64_t f = n;
  45.  
  46.         for (int i = 1; i < n; i++) {   /* pow & nfact */
  47.             p *= rad;
  48.             f *= i;
  49.         }
  50.  
  51.         sine_deg += sign * p / f;       /* Taylor-series term */
  52.     }
  53.  
  54.     return pos_quad ? sine_deg : -sine_deg;
  55. }
  56.  
  57. /** cos with taylor series expansion to n = TSLIM
  58.  *  (no function reliance, quadrants handled)
  59.  */
  60. double cosenfq (const double deg)
  61. {
  62.     double fp = deg - (int)deg,     /* save fractional part of deg */
  63.         qdeg = (int)deg % 360,      /* get equivalent 0-359 deg angle */
  64.         rad, cose_deg = 1.0;        /* radians, cose_deg */
  65.     int pos_quad = 1,               /* positive quadrant flag 1,4  */
  66.         sign = -1;                  /* taylor series term sign */
  67.  
  68.     qdeg += fp;                     /* add fractional part back to angle */
  69.  
  70.     /* get equivalent 0-90 degree angle, set pos_quad flag */
  71.     if (90 < qdeg && qdeg <= 180) {         /* in 2nd quadrant */
  72.         qdeg = 180 - qdeg;
  73.         pos_quad = 0;
  74.     }
  75.     else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
  76.         qdeg = qdeg - 180;
  77.         pos_quad = 0;
  78.     }
  79.     else if (270 < qdeg && qdeg <= 360)     /* in 4th quadrant */
  80.         qdeg = 360 - qdeg;
  81.  
  82.     rad = qdeg * M_PI / 180.0;      /* convert to radians */
  83.  
  84.     /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
  85.     for (int n = 2; n < TSLIM; n += 2, sign *= -1) {
  86.         double p = rad;
  87.         uint64_t f = n;
  88.  
  89.         for (int i = 1; i < n; i++) {   /* pow & nfact */
  90.             p *= rad;
  91.             f *= i;
  92.         }
  93.  
  94.         cose_deg += sign * p / f;       /* Taylor-series term */
  95.     }
  96.  
  97.     return pos_quad ? cose_deg : -cose_deg;
  98. }
  99.  
  100. int main (int argc, char **argv) {
  101.  
  102.     double a = (argc > 1) ? strtod (argv[1], NULL) : 187.2,
  103.         coslibc = cos (a * M_PI / 180.0),
  104.         sinlibc = sin (a * M_PI / 180.0),
  105.         sints = sinenfq (a),
  106.         costs = cosenfq (a),
  107.         err = fabs (sinlibc - sints),
  108.         cerr = fabs (coslibc - costs);
  109.  
  110.     if (err > EMAX)
  111.         fprintf (stderr, "sine error exceeds limit of %e\n"
  112.         "%.2f    %11.8f    %11.8f    %e\n",
  113.         EMAX, a, sinlibc, sints, err);
  114.  
  115.     if (cerr > EMAX)
  116.         fprintf (stderr, "cose error exceeds limit of %e\n"
  117.         "%.2f    %11.8f    %11.8f    %e\n",
  118.         EMAX, a, coslibc, costs, cerr);
  119.  
  120.     printf ("angle        libc_sin         sinets       libc_cos         cosets"
  121.             "         sineerr        coseerr\n");
  122.     printf ("%6.2f    %11.8f    %11.8f    %11.8f    %11.8f    %e    %e\n",
  123.             a, sinlibc, sints, coslibc, costs, err, cerr);
  124.  
  125.     return 0;
  126. }
  127.  
RAW Paste Data