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 */
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;
32.     }
33.     else if (270 < qdeg && qdeg <= 360) {   /* in 4th quadrant */
34.         qdeg = 360 - qdeg;
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) {
44.         uint64_t f = n;
45.
46.         for (int i = 1; i < n; i++) {   /* pow & nfact */
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 */
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;
74.     }
75.     else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
76.         qdeg = qdeg - 180;
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) {
87.         uint64_t f = n;
88.
89.         for (int i = 1; i < n; i++) {   /* pow & nfact */
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