1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
5 | * |
6 | * This program is free software; you can redistribute it and/or modify |
7 | * it under the terms of the GNU Lesser General Public License as published by |
8 | * the Free Software Foundation; either version 2.1 of the License, or |
9 | * (at your option) any later version. |
10 | * |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | * GNU Lesser General Public License for more details. |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License |
17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
18 | */ |
19 | /********************************************************************/ |
20 | /* */ |
21 | /* MODULE_NAME: dosincos.c */ |
22 | /* */ |
23 | /* */ |
24 | /* FUNCTIONS: dubsin */ |
25 | /* dubcos */ |
26 | /* docos */ |
27 | /* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h */ |
28 | /* sincos.tbl */ |
29 | /* */ |
30 | /* Routines compute sin() and cos() as Double-Length numbers */ |
31 | /********************************************************************/ |
32 | |
33 | |
34 | |
35 | #include "endian.h" |
36 | #include "mydefs.h" |
37 | #include <dla.h> |
38 | #include "dosincos.h" |
39 | #include <math_private.h> |
40 | |
41 | #ifndef SECTION |
42 | # define SECTION |
43 | #endif |
44 | |
45 | extern const union |
46 | { |
47 | int4 i[880]; |
48 | double x[440]; |
49 | } __sincostab attribute_hidden; |
50 | |
51 | /***********************************************************************/ |
52 | /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */ |
53 | /* as Double-Length number and store it at array v .It computes it by */ |
54 | /* arithmetic action on Double-Length numbers */ |
55 | /*(x+dx) between 0 and PI/4 */ |
56 | /***********************************************************************/ |
57 | |
58 | void |
59 | SECTION |
60 | __dubsin (double x, double dx, double v[]) |
61 | { |
62 | double r, s, c, cc, d, dd, d2, dd2, e, ee, |
63 | sn, ssn, cs, ccs, ds, dss, dc, dcc; |
64 | #ifndef DLA_FMS |
65 | double p, hx, tx, hy, ty, q; |
66 | #endif |
67 | mynumber u; |
68 | int4 k; |
69 | |
70 | u.x = x + big.x; |
71 | k = u.i[LOW_HALF] << 2; |
72 | x = x - (u.x - big.x); |
73 | d = x + dx; |
74 | dd = (x - d) + dx; |
75 | /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ |
76 | MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc); |
77 | sn = __sincostab.x[k]; /* */ |
78 | ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */ |
79 | cs = __sincostab.x[k + 2]; /* */ |
80 | ccs = __sincostab.x[k + 3]; /* */ |
81 | /* Taylor series for sin ds=sin(t) */ |
82 | MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
83 | ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s); |
84 | MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
85 | ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s); |
86 | MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
87 | MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
88 | ADD2 (ds, dss, d, dd, ds, dss, r, s); |
89 | |
90 | /* Taylor series for cos dc=cos(t) */ |
91 | MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
92 | ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s); |
93 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
94 | ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s); |
95 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
96 | ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s); |
97 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
98 | |
99 | MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc); |
100 | MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
101 | SUB2 (e, ee, dc, dcc, e, ee, r, s); |
102 | ADD2 (e, ee, sn, ssn, e, ee, r, s); /* e+ee=sin(x+dx) */ |
103 | |
104 | v[0] = e; |
105 | v[1] = ee; |
106 | } |
107 | /**********************************************************************/ |
108 | /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ |
109 | /* as Double-Length number and store it in array v .It computes it by */ |
110 | /* arithmetic action on Double-Length numbers */ |
111 | /*(x+dx) between 0 and PI/4 */ |
112 | /**********************************************************************/ |
113 | |
114 | void |
115 | SECTION |
116 | __dubcos (double x, double dx, double v[]) |
117 | { |
118 | double r, s, c, cc, d, dd, d2, dd2, e, ee, |
119 | sn, ssn, cs, ccs, ds, dss, dc, dcc; |
120 | #ifndef DLA_FMS |
121 | double p, hx, tx, hy, ty, q; |
122 | #endif |
123 | mynumber u; |
124 | int4 k; |
125 | u.x = x + big.x; |
126 | k = u.i[LOW_HALF] << 2; |
127 | x = x - (u.x - big.x); |
128 | d = x + dx; |
129 | dd = (x - d) + dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */ |
130 | MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc); |
131 | sn = __sincostab.x[k]; /* */ |
132 | ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */ |
133 | cs = __sincostab.x[k + 2]; /* */ |
134 | ccs = __sincostab.x[k + 3]; /* */ |
135 | MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
136 | ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s); |
137 | MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
138 | ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s); |
139 | MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
140 | MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
141 | ADD2 (ds, dss, d, dd, ds, dss, r, s); |
142 | |
143 | MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
144 | ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s); |
145 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
146 | ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s); |
147 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
148 | ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s); |
149 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
150 | |
151 | MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc); |
152 | MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
153 | |
154 | MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
155 | ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s); |
156 | MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
157 | ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s); |
158 | MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
159 | MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc); |
160 | ADD2 (ds, dss, d, dd, ds, dss, r, s); |
161 | MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
162 | ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s); |
163 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
164 | ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s); |
165 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
166 | ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s); |
167 | MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
168 | MUL2 (sn, ssn, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc); |
169 | MUL2 (dc, dcc, cs, ccs, dc, dcc, p, hx, tx, hy, ty, q, c, cc); |
170 | ADD2 (e, ee, dc, dcc, e, ee, r, s); |
171 | SUB2 (cs, ccs, e, ee, e, ee, r, s); |
172 | |
173 | v[0] = e; |
174 | v[1] = ee; |
175 | } |
176 | /**********************************************************************/ |
177 | /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ |
178 | /* as Double-Length number and store it in array v */ |
179 | /**********************************************************************/ |
180 | void |
181 | SECTION |
182 | __docos (double x, double dx, double v[]) |
183 | { |
184 | double y, yy, p, w[2]; |
185 | if (x > 0) |
186 | { |
187 | y = x; yy = dx; |
188 | } |
189 | else |
190 | { |
191 | y = -x; yy = -dx; |
192 | } |
193 | if (y < 0.5 * hp0.x) /* y< PI/4 */ |
194 | { |
195 | __dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1]; |
196 | } |
197 | else if (y < 1.5 * hp0.x) /* y< 3/4 * PI */ |
198 | { |
199 | p = hp0.x - y; /* p = PI/2 - y */ |
200 | yy = hp1.x - yy; |
201 | y = p + yy; |
202 | yy = (p - y) + yy; |
203 | if (y > 0) |
204 | { |
205 | __dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1]; |
206 | } |
207 | /* cos(x) = sin ( 90 - x ) */ |
208 | else |
209 | { |
210 | __dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1]; |
211 | } |
212 | } |
213 | else /* y>= 3/4 * PI */ |
214 | { |
215 | p = 2.0 * hp0.x - y; /* p = PI- y */ |
216 | yy = 2.0 * hp1.x - yy; |
217 | y = p + yy; |
218 | yy = (p - y) + yy; |
219 | __dubcos (y, yy, w); |
220 | v[0] = -w[0]; |
221 | v[1] = -w[1]; |
222 | } |
223 | } |
224 | |