]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
b168057a | 4 | * Copyright (C) 2001-2015 Free Software Foundation, Inc. |
f7eac6eb | 5 | * |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
f7eac6eb | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
f7eac6eb | 15 | * |
e4d82761 | 16 | * You should have received a copy of the GNU Lesser General Public License |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f7eac6eb | 18 | */ |
e4d82761 UD |
19 | /****************************************************************************/ |
20 | /* */ | |
21 | /* MODULE_NAME:usncs.c */ | |
22 | /* */ | |
23 | /* FUNCTIONS: usin */ | |
24 | /* ucos */ | |
25 | /* slow */ | |
26 | /* slow1 */ | |
27 | /* slow2 */ | |
28 | /* sloww */ | |
29 | /* sloww1 */ | |
30 | /* sloww2 */ | |
31 | /* bsloww */ | |
32 | /* bsloww1 */ | |
33 | /* bsloww2 */ | |
34 | /* cslow2 */ | |
35 | /* csloww */ | |
36 | /* csloww1 */ | |
37 | /* csloww2 */ | |
38 | /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ | |
39 | /* branred.c sincos32.c dosincos.c mpa.c */ | |
40 | /* sincos.tbl */ | |
41 | /* */ | |
42 | /* An ultimate sin and routine. Given an IEEE double machine number x */ | |
43 | /* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */ | |
44 | /* Assumption: Machine arithmetic operations are performed in */ | |
45 | /* round to nearest mode of IEEE 754 standard. */ | |
46 | /* */ | |
47 | /****************************************************************************/ | |
f7eac6eb | 48 | |
f7eac6eb | 49 | |
0c59a196 | 50 | #include <errno.h> |
e4d82761 UD |
51 | #include "endian.h" |
52 | #include "mydefs.h" | |
53 | #include "usncs.h" | |
54 | #include "MathLib.h" | |
0e9be4db | 55 | #include <math.h> |
1ed0291c | 56 | #include <math_private.h> |
804360ed | 57 | #include <fenv.h> |
e4d82761 | 58 | |
4aafb73c | 59 | /* Helper macros to compute sin of the input values. */ |
196f7f5d | 60 | #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx)) |
4aafb73c | 61 | |
196f7f5d | 62 | #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1) |
4aafb73c SP |
63 | |
64 | /* The computed polynomial is a variation of the Taylor series expansion for | |
65 | sin(a): | |
66 | ||
67 | a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2 | |
68 | ||
69 | The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so | |
70 | on. The result is returned to LHS and correction in COR. */ | |
8d561986 | 71 | #define TAYLOR_SIN(xx, a, da, cor) \ |
4aafb73c SP |
72 | ({ \ |
73 | double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \ | |
74 | double res = (a) + t; \ | |
75 | (cor) = ((a) - res) + t; \ | |
76 | res; \ | |
77 | }) | |
78 | ||
79 | /* This is again a variation of the Taylor series expansion with the term | |
80 | x^3/3! expanded into the following for better accuracy: | |
81 | ||
82 | bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3 | |
83 | ||
84 | The correction term is dx and bb + aa = -1/3! | |
85 | */ | |
86 | #define TAYLOR_SLOW(x0, dx, cor) \ | |
87 | ({ \ | |
88 | static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ \ | |
89 | double xx = (x0) * (x0); \ | |
90 | double x1 = ((x0) + th2_36) - th2_36; \ | |
196f7f5d | 91 | double y = aa * x1 * x1 * x1; \ |
4aafb73c SP |
92 | double r = (x0) + y; \ |
93 | double x2 = ((x0) - x1) + (dx); \ | |
196f7f5d SP |
94 | double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2) \ |
95 | * (x0) + aa * x2 * x2 * x2 + (dx)); \ | |
4aafb73c SP |
96 | t = (((x0) - r) + y) + t; \ |
97 | double res = r + t; \ | |
98 | (cor) = (r - res) + t; \ | |
99 | res; \ | |
100 | }) | |
101 | ||
b348e1e3 SP |
102 | #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \ |
103 | ({ \ | |
104 | int4 k = u.i[LOW_HALF] << 2; \ | |
105 | sn = __sincostab.x[k]; \ | |
106 | ssn = __sincostab.x[k + 1]; \ | |
107 | cs = __sincostab.x[k + 2]; \ | |
108 | ccs = __sincostab.x[k + 3]; \ | |
109 | }) | |
110 | ||
31d3cc00 UD |
111 | #ifndef SECTION |
112 | # define SECTION | |
113 | #endif | |
114 | ||
af968f62 UD |
115 | extern const union |
116 | { | |
117 | int4 i[880]; | |
118 | double x[440]; | |
119 | } __sincostab attribute_hidden; | |
120 | ||
e4d82761 | 121 | static const double |
6dbe713d SP |
122 | sn3 = -1.66666666666664880952546298448555E-01, |
123 | sn5 = 8.33333214285722277379541354343671E-03, | |
124 | cs2 = 4.99999999999999999999950396842453E-01, | |
125 | cs4 = -4.16666666666664434524222570944589E-02, | |
126 | cs6 = 1.38888874007937613028114285595617E-03; | |
127 | ||
7a74607f SP |
128 | static const double t22 = 0x1.8p22; |
129 | ||
6dbe713d SP |
130 | void __dubsin (double x, double dx, double w[]); |
131 | void __docos (double x, double dx, double w[]); | |
09544cbc SP |
132 | double __mpsin (double x, double dx, bool reduce_range); |
133 | double __mpcos (double x, double dx, bool reduce_range); | |
6dbe713d SP |
134 | static double slow (double x); |
135 | static double slow1 (double x); | |
136 | static double slow2 (double x); | |
137 | static double sloww (double x, double dx, double orig); | |
84ba214c | 138 | static double sloww1 (double x, double dx, double orig, int m); |
6dbe713d SP |
139 | static double sloww2 (double x, double dx, double orig, int n); |
140 | static double bsloww (double x, double dx, double orig, int n); | |
141 | static double bsloww1 (double x, double dx, double orig, int n); | |
142 | static double bsloww2 (double x, double dx, double orig, int n); | |
143 | int __branred (double x, double *a, double *aa); | |
144 | static double cslow2 (double x); | |
145 | static double csloww (double x, double dx, double orig); | |
84ba214c | 146 | static double csloww1 (double x, double dx, double orig, int m); |
6dbe713d SP |
147 | static double csloww2 (double x, double dx, double orig, int n); |
148 | ||
392dd2de SP |
149 | /* Given a number partitioned into U and X such that U is an index into the |
150 | sin/cos table, this macro computes the cosine of the number by combining | |
151 | the sin and cos of X (as computed by a variation of the Taylor series) with | |
152 | the values looked up from the sin/cos table to get the result in RES and a | |
153 | correction value in COR. */ | |
154 | static double | |
155 | do_cos (mynumber u, double x, double *corp) | |
156 | { | |
157 | double xx, s, sn, ssn, c, cs, ccs, res, cor; | |
158 | xx = x * x; | |
159 | s = x + x * xx * (sn3 + xx * sn5); | |
160 | c = xx * (cs2 + xx * (cs4 + xx * cs6)); | |
161 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
162 | cor = (ccs - s * ssn - cs * c) - sn * s; | |
163 | res = cs + cor; | |
164 | cor = (cs - res) + cor; | |
165 | *corp = cor; | |
166 | return res; | |
167 | } | |
168 | ||
169 | /* A more precise variant of DO_COS where the number is partitioned into U, X | |
170 | and DX. EPS is the adjustment to the correction COR. */ | |
171 | static double | |
172 | do_cos_slow (mynumber u, double x, double dx, double eps, double *corp) | |
173 | { | |
174 | double xx, y, x1, x2, e1, e2, res, cor; | |
175 | double s, sn, ssn, c, cs, ccs; | |
176 | xx = x * x; | |
177 | s = x * xx * (sn3 + xx * sn5); | |
178 | c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); | |
179 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
180 | x1 = (x + t22) - t22; | |
181 | x2 = (x - x1) + dx; | |
182 | e1 = (sn + t22) - t22; | |
183 | e2 = (sn - e1) + ssn; | |
184 | cor = (ccs - cs * c - e1 * x2 - e2 * x) - sn * s; | |
185 | y = cs - e1 * x1; | |
186 | cor = cor + ((cs - y) - e1 * x1); | |
187 | res = y + cor; | |
188 | cor = (y - res) + cor; | |
189 | if (cor > 0) | |
190 | cor = 1.0005 * cor + eps; | |
191 | else | |
192 | cor = 1.0005 * cor - eps; | |
193 | *corp = cor; | |
194 | return res; | |
195 | } | |
196 | ||
197 | /* Given a number partitioned into U and X and DX such that U is an index into | |
198 | the sin/cos table, this macro computes the sine of the number by combining | |
199 | the sin and cos of X (as computed by a variation of the Taylor series) with | |
200 | the values looked up from the sin/cos table to get the result in RES and a | |
201 | correction value in COR. */ | |
202 | static double | |
203 | do_sin (mynumber u, double x, double dx, double *corp) | |
204 | { | |
205 | double xx, s, sn, ssn, c, cs, ccs, cor, res; | |
206 | xx = x * x; | |
207 | s = x + (dx + x * xx * (sn3 + xx * sn5)); | |
208 | c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); | |
209 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
210 | cor = (ssn + s * ccs - sn * c) + cs * s; | |
211 | res = sn + cor; | |
212 | cor = (sn - res) + cor; | |
213 | *corp = cor; | |
214 | return res; | |
215 | } | |
216 | ||
217 | /* A more precise variant of res = do_sin where the number is partitioned into U, X | |
218 | and DX. EPS is the adjustment to the correction COR. */ | |
219 | static double | |
220 | do_sin_slow (mynumber u, double x, double dx, double eps, double *corp) | |
221 | { | |
222 | double xx, y, x1, x2, c1, c2, res, cor; | |
223 | double s, sn, ssn, c, cs, ccs; | |
224 | xx = x * x; | |
225 | s = x * xx * (sn3 + xx * sn5); | |
226 | c = xx * (cs2 + xx * (cs4 + xx * cs6)); | |
227 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
228 | x1 = (x + t22) - t22; | |
229 | x2 = (x - x1) + dx; | |
230 | c1 = (cs + t22) - t22; | |
231 | c2 = (cs - c1) + ccs; | |
232 | cor = (ssn + s * ccs + cs * s + c2 * x + c1 * x2 - sn * x * dx) - sn * c; | |
233 | y = sn + c1 * x1; | |
234 | cor = cor + ((sn - y) + c1 * x1); | |
235 | res = y + cor; | |
236 | cor = (y - res) + cor; | |
237 | if (cor > 0) | |
238 | cor = 1.0005 * cor + eps; | |
239 | else | |
240 | cor = 1.0005 * cor - eps; | |
241 | *corp = cor; | |
242 | return res; | |
243 | } | |
244 | ||
6cce25f8 SP |
245 | /* Reduce range of X and compute sin of a + da. K is the amount by which to |
246 | rotate the quadrants. This allows us to use the same routine to compute cos | |
247 | by simply rotating the quadrants by 1. */ | |
248 | static inline double | |
249 | __always_inline | |
975195e4 | 250 | reduce_and_compute (double x, unsigned int k) |
6cce25f8 | 251 | { |
975195e4 | 252 | double retval = 0, a, da; |
6cce25f8 SP |
253 | unsigned int n = __branred (x, &a, &da); |
254 | k = (n + k) % 4; | |
255 | switch (k) | |
256 | { | |
257 | case 0: | |
258 | if (a * a < 0.01588) | |
259 | retval = bsloww (a, da, x, n); | |
260 | else | |
261 | retval = bsloww1 (a, da, x, n); | |
262 | break; | |
263 | case 2: | |
264 | if (a * a < 0.01588) | |
265 | retval = bsloww (-a, -da, x, n); | |
266 | else | |
267 | retval = bsloww1 (-a, -da, x, n); | |
268 | break; | |
269 | ||
270 | case 1: | |
271 | case 3: | |
272 | retval = bsloww2 (a, da, x, n); | |
273 | break; | |
274 | } | |
275 | return retval; | |
276 | } | |
277 | ||
e4d82761 UD |
278 | /*******************************************************************/ |
279 | /* An ultimate sin routine. Given an IEEE double machine number x */ | |
280 | /* it computes the correctly rounded (to nearest) value of sin(x) */ | |
281 | /*******************************************************************/ | |
31d3cc00 UD |
282 | double |
283 | SECTION | |
6dbe713d SP |
284 | __sin (double x) |
285 | { | |
286 | double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1, | |
287 | xn2; | |
288 | mynumber u, v; | |
289 | int4 k, m, n; | |
290 | double retval = 0; | |
291 | ||
292 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); | |
293 | ||
294 | u.x = x; | |
295 | m = u.i[HIGH_HALF]; | |
296 | k = 0x7fffffff & m; /* no sign */ | |
297 | if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ | |
5eea0404 | 298 | retval = x; |
e4d82761 | 299 | /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ |
6dbe713d SP |
300 | else if (k < 0x3fd00000) |
301 | { | |
302 | xx = x * x; | |
4aafb73c SP |
303 | /* Taylor series. */ |
304 | t = POLYNOMIAL (xx) * (xx * x); | |
6dbe713d SP |
305 | res = x + t; |
306 | cor = (x - res) + t; | |
307 | retval = (res == res + 1.07 * cor) ? res : slow (x); | |
6dbe713d | 308 | } /* else if (k < 0x3fd00000) */ |
e4d82761 | 309 | /*---------------------------- 0.25<|x|< 0.855469---------------------- */ |
6dbe713d SP |
310 | else if (k < 0x3feb6000) |
311 | { | |
196f7f5d SP |
312 | u.x = (m > 0) ? big + x : big - x; |
313 | y = (m > 0) ? x - (u.x - big) : x + (u.x - big); | |
6dbe713d SP |
314 | xx = y * y; |
315 | s = y + y * xx * (sn3 + xx * sn5); | |
316 | c = xx * (cs2 + xx * (cs4 + xx * cs6)); | |
b348e1e3 SP |
317 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); |
318 | if (m <= 0) | |
319 | { | |
320 | sn = -sn; | |
321 | ssn = -ssn; | |
322 | } | |
6dbe713d SP |
323 | cor = (ssn + s * ccs - sn * c) + cs * s; |
324 | res = sn + cor; | |
325 | cor = (sn - res) + cor; | |
326 | retval = (res == res + 1.096 * cor) ? res : slow1 (x); | |
6dbe713d | 327 | } /* else if (k < 0x3feb6000) */ |
e4d82761 UD |
328 | |
329 | /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ | |
6dbe713d SP |
330 | else if (k < 0x400368fd) |
331 | { | |
332 | ||
196f7f5d | 333 | y = (m > 0) ? hp0 - x : hp0 + x; |
6dbe713d SP |
334 | if (y >= 0) |
335 | { | |
196f7f5d SP |
336 | u.x = big + y; |
337 | y = (y - (u.x - big)) + hp1; | |
6dbe713d SP |
338 | } |
339 | else | |
340 | { | |
196f7f5d SP |
341 | u.x = big - y; |
342 | y = (-hp1) - (y + (u.x - big)); | |
6dbe713d | 343 | } |
392dd2de | 344 | res = do_cos (u, y, &cor); |
6dbe713d | 345 | retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x); |
6dbe713d | 346 | } /* else if (k < 0x400368fd) */ |
e4d82761 UD |
347 | |
348 | /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ | |
6dbe713d SP |
349 | else if (k < 0x419921FB) |
350 | { | |
196f7f5d SP |
351 | t = (x * hpinv + toint); |
352 | xn = t - toint; | |
6dbe713d | 353 | v.x = t; |
196f7f5d | 354 | y = (x - xn * mp1) - xn * mp2; |
6dbe713d | 355 | n = v.i[LOW_HALF] & 3; |
196f7f5d | 356 | da = xn * mp3; |
6dbe713d SP |
357 | a = y - da; |
358 | da = (y - a) - da; | |
0e9be4db | 359 | eps = fabs (x) * 1.2e-30; |
6dbe713d SP |
360 | |
361 | switch (n) | |
362 | { /* quarter of unit circle */ | |
363 | case 0: | |
364 | case 2: | |
365 | xx = a * a; | |
366 | if (n) | |
367 | { | |
368 | a = -a; | |
369 | da = -da; | |
370 | } | |
371 | if (xx < 0.01588) | |
372 | { | |
4aafb73c | 373 | /* Taylor series. */ |
8d561986 | 374 | res = TAYLOR_SIN (xx, a, da, cor); |
6dbe713d SP |
375 | cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; |
376 | retval = (res == res + cor) ? res : sloww (a, da, x); | |
f7eac6eb | 377 | } |
6dbe713d SP |
378 | else |
379 | { | |
380 | if (a > 0) | |
84ba214c | 381 | m = 1; |
e4d82761 | 382 | else |
6dbe713d SP |
383 | { |
384 | m = 0; | |
84ba214c | 385 | a = -a; |
5ff8d60e | 386 | da = -da; |
6dbe713d | 387 | } |
84ba214c SP |
388 | u.x = big + a; |
389 | y = a - (u.x - big); | |
392dd2de | 390 | res = do_sin (u, y, da, &cor); |
6dbe713d SP |
391 | cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; |
392 | retval = ((res == res + cor) ? ((m) ? res : -res) | |
84ba214c | 393 | : sloww1 (a, da, x, m)); |
e4d82761 | 394 | } |
6dbe713d SP |
395 | break; |
396 | ||
397 | case 1: | |
398 | case 3: | |
399 | if (a < 0) | |
400 | { | |
401 | a = -a; | |
402 | da = -da; | |
403 | } | |
196f7f5d SP |
404 | u.x = big + a; |
405 | y = a - (u.x - big) + da; | |
392dd2de | 406 | res = do_cos (u, y, &cor); |
6dbe713d SP |
407 | cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; |
408 | retval = ((res == res + cor) ? ((n & 2) ? -res : res) | |
409 | : sloww2 (a, da, x, n)); | |
6dbe713d SP |
410 | break; |
411 | } | |
6dbe713d | 412 | } /* else if (k < 0x419921FB ) */ |
e4d82761 UD |
413 | |
414 | /*---------------------105414350 <|x|< 281474976710656 --------------------*/ | |
6dbe713d SP |
415 | else if (k < 0x42F00000) |
416 | { | |
196f7f5d SP |
417 | t = (x * hpinv + toint); |
418 | xn = t - toint; | |
6dbe713d SP |
419 | v.x = t; |
420 | xn1 = (xn + 8.0e22) - 8.0e22; | |
421 | xn2 = xn - xn1; | |
196f7f5d | 422 | y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); |
6dbe713d | 423 | n = v.i[LOW_HALF] & 3; |
196f7f5d | 424 | da = xn1 * pp3; |
6dbe713d SP |
425 | t = y - da; |
426 | da = (y - t) - da; | |
196f7f5d | 427 | da = (da - xn2 * pp3) - xn * pp4; |
6dbe713d SP |
428 | a = t + da; |
429 | da = (t - a) + da; | |
430 | eps = 1.0e-24; | |
431 | ||
432 | switch (n) | |
433 | { | |
434 | case 0: | |
435 | case 2: | |
436 | xx = a * a; | |
437 | if (n) | |
438 | { | |
439 | a = -a; | |
440 | da = -da; | |
441 | } | |
442 | if (xx < 0.01588) | |
443 | { | |
4aafb73c | 444 | /* Taylor series. */ |
8d561986 | 445 | res = TAYLOR_SIN (xx, a, da, cor); |
6dbe713d SP |
446 | cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; |
447 | retval = (res == res + cor) ? res : bsloww (a, da, x, n); | |
e4d82761 | 448 | } |
6dbe713d SP |
449 | else |
450 | { | |
1cadc858 | 451 | double t; |
6dbe713d SP |
452 | if (a > 0) |
453 | { | |
454 | m = 1; | |
1cadc858 | 455 | t = a; |
6dbe713d SP |
456 | db = da; |
457 | } | |
458 | else | |
459 | { | |
460 | m = 0; | |
1cadc858 | 461 | t = -a; |
6dbe713d SP |
462 | db = -da; |
463 | } | |
1cadc858 SP |
464 | u.x = big + t; |
465 | y = t - (u.x - big); | |
392dd2de | 466 | res = do_sin (u, y, db, &cor); |
6dbe713d SP |
467 | cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; |
468 | retval = ((res == res + cor) ? ((m) ? res : -res) | |
469 | : bsloww1 (a, da, x, n)); | |
6dbe713d SP |
470 | } |
471 | break; | |
472 | ||
473 | case 1: | |
474 | case 3: | |
475 | if (a < 0) | |
476 | { | |
477 | a = -a; | |
478 | da = -da; | |
479 | } | |
196f7f5d SP |
480 | u.x = big + a; |
481 | y = a - (u.x - big) + da; | |
392dd2de | 482 | res = do_cos (u, y, &cor); |
6dbe713d SP |
483 | cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; |
484 | retval = ((res == res + cor) ? ((n & 2) ? -res : res) | |
485 | : bsloww2 (a, da, x, n)); | |
6dbe713d SP |
486 | break; |
487 | } | |
488 | } /* else if (k < 0x42F00000 ) */ | |
e4d82761 UD |
489 | |
490 | /* -----------------281474976710656 <|x| <2^1024----------------------------*/ | |
6dbe713d | 491 | else if (k < 0x7ff00000) |
975195e4 | 492 | retval = reduce_and_compute (x, 0); |
6dbe713d SP |
493 | |
494 | /*--------------------- |x| > 2^1024 ----------------------------------*/ | |
495 | else | |
496 | { | |
497 | if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) | |
498 | __set_errno (EDOM); | |
499 | retval = x / x; | |
6dbe713d | 500 | } |
804360ed | 501 | |
6dbe713d | 502 | return retval; |
e4d82761 UD |
503 | } |
504 | ||
505 | ||
506 | /*******************************************************************/ | |
507 | /* An ultimate cos routine. Given an IEEE double machine number x */ | |
508 | /* it computes the correctly rounded (to nearest) value of cos(x) */ | |
509 | /*******************************************************************/ | |
510 | ||
31d3cc00 UD |
511 | double |
512 | SECTION | |
6dbe713d | 513 | __cos (double x) |
e4d82761 | 514 | { |
392dd2de | 515 | double y, xx, res, t, cor, xn, a, da, db, eps, xn1, |
6dbe713d SP |
516 | xn2; |
517 | mynumber u, v; | |
518 | int4 k, m, n; | |
e4d82761 | 519 | |
804360ed JM |
520 | double retval = 0; |
521 | ||
eb92c487 | 522 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); |
804360ed | 523 | |
e4d82761 UD |
524 | u.x = x; |
525 | m = u.i[HIGH_HALF]; | |
6dbe713d | 526 | k = 0x7fffffff & m; |
e4d82761 | 527 | |
5eea0404 | 528 | /* |x|<2^-27 => cos(x)=1 */ |
6dbe713d | 529 | if (k < 0x3e400000) |
5eea0404 | 530 | retval = 1.0; |
6dbe713d SP |
531 | |
532 | else if (k < 0x3feb6000) | |
533 | { /* 2^-27 < |x| < 0.855469 */ | |
0e9be4db | 534 | y = fabs (x); |
196f7f5d SP |
535 | u.x = big + y; |
536 | y = y - (u.x - big); | |
392dd2de | 537 | res = do_cos (u, y, &cor); |
6dbe713d | 538 | retval = (res == res + 1.020 * cor) ? res : cslow2 (x); |
6dbe713d SP |
539 | } /* else if (k < 0x3feb6000) */ |
540 | ||
541 | else if (k < 0x400368fd) | |
542 | { /* 0.855469 <|x|<2.426265 */ ; | |
0e9be4db | 543 | y = hp0 - fabs (x); |
196f7f5d SP |
544 | a = y + hp1; |
545 | da = (y - a) + hp1; | |
6dbe713d SP |
546 | xx = a * a; |
547 | if (xx < 0.01588) | |
548 | { | |
8d561986 | 549 | res = TAYLOR_SIN (xx, a, da, cor); |
6dbe713d SP |
550 | cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31; |
551 | retval = (res == res + cor) ? res : csloww (a, da, x); | |
6dbe713d SP |
552 | } |
553 | else | |
554 | { | |
555 | if (a > 0) | |
556 | { | |
557 | m = 1; | |
6dbe713d SP |
558 | } |
559 | else | |
560 | { | |
561 | m = 0; | |
84ba214c | 562 | a = -a; |
5ff8d60e | 563 | da = -da; |
6dbe713d | 564 | } |
84ba214c SP |
565 | u.x = big + a; |
566 | y = a - (u.x - big); | |
392dd2de | 567 | res = do_sin (u, y, da, &cor); |
6dbe713d SP |
568 | cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31; |
569 | retval = ((res == res + cor) ? ((m) ? res : -res) | |
84ba214c | 570 | : csloww1 (a, da, x, m)); |
6dbe713d | 571 | } |
e4d82761 | 572 | |
6dbe713d | 573 | } /* else if (k < 0x400368fd) */ |
e4d82761 | 574 | |
e4d82761 | 575 | |
6dbe713d SP |
576 | else if (k < 0x419921FB) |
577 | { /* 2.426265<|x|< 105414350 */ | |
196f7f5d SP |
578 | t = (x * hpinv + toint); |
579 | xn = t - toint; | |
6dbe713d | 580 | v.x = t; |
196f7f5d | 581 | y = (x - xn * mp1) - xn * mp2; |
6dbe713d | 582 | n = v.i[LOW_HALF] & 3; |
196f7f5d | 583 | da = xn * mp3; |
6dbe713d SP |
584 | a = y - da; |
585 | da = (y - a) - da; | |
0e9be4db | 586 | eps = fabs (x) * 1.2e-30; |
6dbe713d SP |
587 | |
588 | switch (n) | |
589 | { | |
590 | case 1: | |
591 | case 3: | |
592 | xx = a * a; | |
593 | if (n == 1) | |
594 | { | |
595 | a = -a; | |
596 | da = -da; | |
597 | } | |
598 | if (xx < 0.01588) | |
599 | { | |
8d561986 | 600 | res = TAYLOR_SIN (xx, a, da, cor); |
6dbe713d SP |
601 | cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; |
602 | retval = (res == res + cor) ? res : csloww (a, da, x); | |
6dbe713d SP |
603 | } |
604 | else | |
605 | { | |
606 | if (a > 0) | |
607 | { | |
608 | m = 1; | |
6dbe713d SP |
609 | } |
610 | else | |
611 | { | |
612 | m = 0; | |
84ba214c | 613 | a = -a; |
5ff8d60e | 614 | da = -da; |
6dbe713d | 615 | } |
84ba214c SP |
616 | u.x = big + a; |
617 | y = a - (u.x - big); | |
392dd2de | 618 | res = do_sin (u, y, da, &cor); |
6dbe713d SP |
619 | cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; |
620 | retval = ((res == res + cor) ? ((m) ? res : -res) | |
84ba214c | 621 | : csloww1 (a, da, x, m)); |
6dbe713d SP |
622 | } |
623 | break; | |
624 | ||
625 | case 0: | |
626 | case 2: | |
627 | if (a < 0) | |
628 | { | |
629 | a = -a; | |
630 | da = -da; | |
631 | } | |
196f7f5d SP |
632 | u.x = big + a; |
633 | y = a - (u.x - big) + da; | |
392dd2de | 634 | res = do_cos (u, y, &cor); |
6dbe713d SP |
635 | cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; |
636 | retval = ((res == res + cor) ? ((n) ? -res : res) | |
637 | : csloww2 (a, da, x, n)); | |
6dbe713d SP |
638 | break; |
639 | } | |
640 | } /* else if (k < 0x419921FB ) */ | |
e4d82761 | 641 | |
6dbe713d SP |
642 | else if (k < 0x42F00000) |
643 | { | |
196f7f5d SP |
644 | t = (x * hpinv + toint); |
645 | xn = t - toint; | |
6dbe713d SP |
646 | v.x = t; |
647 | xn1 = (xn + 8.0e22) - 8.0e22; | |
648 | xn2 = xn - xn1; | |
196f7f5d | 649 | y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); |
6dbe713d | 650 | n = v.i[LOW_HALF] & 3; |
196f7f5d | 651 | da = xn1 * pp3; |
6dbe713d SP |
652 | t = y - da; |
653 | da = (y - t) - da; | |
196f7f5d | 654 | da = (da - xn2 * pp3) - xn * pp4; |
6dbe713d SP |
655 | a = t + da; |
656 | da = (t - a) + da; | |
657 | eps = 1.0e-24; | |
658 | ||
659 | switch (n) | |
660 | { | |
661 | case 1: | |
662 | case 3: | |
663 | xx = a * a; | |
664 | if (n == 1) | |
665 | { | |
666 | a = -a; | |
667 | da = -da; | |
668 | } | |
669 | if (xx < 0.01588) | |
670 | { | |
8d561986 | 671 | res = TAYLOR_SIN (xx, a, da, cor); |
6dbe713d SP |
672 | cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; |
673 | retval = (res == res + cor) ? res : bsloww (a, da, x, n); | |
6dbe713d SP |
674 | } |
675 | else | |
676 | { | |
1cadc858 | 677 | double t; |
6dbe713d SP |
678 | if (a > 0) |
679 | { | |
680 | m = 1; | |
1cadc858 | 681 | t = a; |
6dbe713d SP |
682 | db = da; |
683 | } | |
684 | else | |
685 | { | |
686 | m = 0; | |
1cadc858 | 687 | t = -a; |
6dbe713d SP |
688 | db = -da; |
689 | } | |
1cadc858 SP |
690 | u.x = big + t; |
691 | y = t - (u.x - big); | |
392dd2de | 692 | res = do_sin (u, y, db, &cor); |
6dbe713d SP |
693 | cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; |
694 | retval = ((res == res + cor) ? ((m) ? res : -res) | |
695 | : bsloww1 (a, da, x, n)); | |
6dbe713d SP |
696 | } |
697 | break; | |
698 | ||
699 | case 0: | |
700 | case 2: | |
701 | if (a < 0) | |
702 | { | |
703 | a = -a; | |
704 | da = -da; | |
705 | } | |
196f7f5d SP |
706 | u.x = big + a; |
707 | y = a - (u.x - big) + da; | |
392dd2de | 708 | res = do_cos (u, y, &cor); |
6dbe713d SP |
709 | cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; |
710 | retval = ((res == res + cor) ? ((n) ? -res : res) | |
711 | : bsloww2 (a, da, x, n)); | |
6dbe713d SP |
712 | break; |
713 | } | |
714 | } /* else if (k < 0x42F00000 ) */ | |
715 | ||
6cce25f8 | 716 | /* 281474976710656 <|x| <2^1024 */ |
6dbe713d | 717 | else if (k < 0x7ff00000) |
975195e4 | 718 | retval = reduce_and_compute (x, 1); |
e4d82761 | 719 | |
6dbe713d SP |
720 | else |
721 | { | |
722 | if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) | |
723 | __set_errno (EDOM); | |
724 | retval = x / x; /* |x| > 2^1024 */ | |
e4d82761 UD |
725 | } |
726 | ||
804360ed | 727 | return retval; |
e4d82761 UD |
728 | } |
729 | ||
730 | /************************************************************************/ | |
731 | /* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */ | |
732 | /* precision and if still doesn't accurate enough by mpsin or dubsin */ | |
733 | /************************************************************************/ | |
734 | ||
31d3cc00 UD |
735 | static double |
736 | SECTION | |
6dbe713d SP |
737 | slow (double x) |
738 | { | |
4aafb73c SP |
739 | double res, cor, w[2]; |
740 | res = TAYLOR_SLOW (x, 0, cor); | |
6dbe713d SP |
741 | if (res == res + 1.0007 * cor) |
742 | return res; | |
743 | else | |
744 | { | |
0e9be4db | 745 | __dubsin (fabs (x), 0, w); |
6dbe713d SP |
746 | if (w[0] == w[0] + 1.000000001 * w[1]) |
747 | return (x > 0) ? w[0] : -w[0]; | |
748 | else | |
09544cbc | 749 | return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); |
6dbe713d | 750 | } |
e4d82761 | 751 | } |
6dbe713d | 752 | |
e4d82761 | 753 | /*******************************************************************************/ |
c5d5d574 | 754 | /* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */ |
e4d82761 UD |
755 | /* and if result still doesn't accurate enough by mpsin or dubsin */ |
756 | /*******************************************************************************/ | |
757 | ||
31d3cc00 UD |
758 | static double |
759 | SECTION | |
6dbe713d SP |
760 | slow1 (double x) |
761 | { | |
e4d82761 | 762 | mynumber u; |
392dd2de | 763 | double w[2], y, cor, res; |
0e9be4db | 764 | y = fabs (x); |
196f7f5d SP |
765 | u.x = big + y; |
766 | y = y - (u.x - big); | |
392dd2de SP |
767 | res = do_sin_slow (u, y, 0, 0, &cor); |
768 | if (res == res + cor) | |
6dbe713d SP |
769 | return (x > 0) ? res : -res; |
770 | else | |
771 | { | |
0e9be4db | 772 | __dubsin (fabs (x), 0, w); |
6dbe713d SP |
773 | if (w[0] == w[0] + 1.000000005 * w[1]) |
774 | return (x > 0) ? w[0] : -w[0]; | |
775 | else | |
09544cbc | 776 | return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); |
6dbe713d | 777 | } |
e4d82761 | 778 | } |
6dbe713d | 779 | |
e4d82761 | 780 | /**************************************************************************/ |
af968f62 | 781 | /* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */ |
e4d82761 UD |
782 | /* and if result still doesn't accurate enough by mpsin or dubsin */ |
783 | /**************************************************************************/ | |
31d3cc00 UD |
784 | static double |
785 | SECTION | |
6dbe713d SP |
786 | slow2 (double x) |
787 | { | |
e4d82761 | 788 | mynumber u; |
392dd2de | 789 | double w[2], y, y1, y2, cor, res, del; |
7a74607f | 790 | |
0e9be4db | 791 | y = fabs (x); |
196f7f5d | 792 | y = hp0 - y; |
6dbe713d SP |
793 | if (y >= 0) |
794 | { | |
196f7f5d SP |
795 | u.x = big + y; |
796 | y = y - (u.x - big); | |
797 | del = hp1; | |
6dbe713d SP |
798 | } |
799 | else | |
800 | { | |
196f7f5d SP |
801 | u.x = big - y; |
802 | y = -(y + (u.x - big)); | |
803 | del = -hp1; | |
6dbe713d | 804 | } |
392dd2de SP |
805 | res = do_cos_slow (u, y, del, 0, &cor); |
806 | if (res == res + cor) | |
6dbe713d SP |
807 | return (x > 0) ? res : -res; |
808 | else | |
809 | { | |
0e9be4db | 810 | y = fabs (x) - hp0; |
196f7f5d SP |
811 | y1 = y - hp1; |
812 | y2 = (y - y1) - hp1; | |
6dbe713d SP |
813 | __docos (y1, y2, w); |
814 | if (w[0] == w[0] + 1.000000005 * w[1]) | |
815 | return (x > 0) ? w[0] : -w[0]; | |
816 | else | |
09544cbc | 817 | return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); |
6dbe713d | 818 | } |
e4d82761 | 819 | } |
6dbe713d | 820 | |
e4d82761 UD |
821 | /***************************************************************************/ |
822 | /* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/ | |
823 | /* to use Taylor series around zero and (x+dx) */ | |
824 | /* in first or third quarter of unit circle.Routine receive also */ | |
825 | /* (right argument) the original value of x for computing error of */ | |
826 | /* result.And if result not accurate enough routine calls mpsin1 or dubsin */ | |
827 | /***************************************************************************/ | |
828 | ||
31d3cc00 UD |
829 | static double |
830 | SECTION | |
6dbe713d SP |
831 | sloww (double x, double dx, double orig) |
832 | { | |
4aafb73c | 833 | double y, t, res, cor, w[2], a, da, xn; |
8d561986 | 834 | mynumber v; |
e4d82761 | 835 | int4 n; |
4aafb73c | 836 | res = TAYLOR_SLOW (x, dx, cor); |
8d561986 | 837 | if (cor > 0) |
0e9be4db | 838 | cor = 1.0005 * cor + fabs (orig) * 3.1e-30; |
8d561986 | 839 | else |
0e9be4db | 840 | cor = 1.0005 * cor - fabs (orig) * 3.1e-30; |
8d561986 | 841 | |
6dbe713d SP |
842 | if (res == res + cor) |
843 | return res; | |
844 | else | |
845 | { | |
846 | (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); | |
847 | if (w[1] > 0) | |
0e9be4db | 848 | cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; |
6dbe713d | 849 | else |
0e9be4db | 850 | cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; |
6dbe713d SP |
851 | |
852 | if (w[0] == w[0] + cor) | |
853 | return (x > 0) ? w[0] : -w[0]; | |
854 | else | |
855 | { | |
196f7f5d SP |
856 | t = (orig * hpinv + toint); |
857 | xn = t - toint; | |
6dbe713d | 858 | v.x = t; |
196f7f5d | 859 | y = (orig - xn * mp1) - xn * mp2; |
6dbe713d | 860 | n = v.i[LOW_HALF] & 3; |
196f7f5d | 861 | da = xn * pp3; |
6dbe713d SP |
862 | t = y - da; |
863 | da = (y - t) - da; | |
196f7f5d | 864 | y = xn * pp4; |
6dbe713d SP |
865 | a = t - y; |
866 | da = ((t - a) - y) + da; | |
867 | if (n & 2) | |
868 | { | |
869 | a = -a; | |
870 | da = -da; | |
871 | } | |
872 | (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); | |
873 | if (w[1] > 0) | |
0e9be4db | 874 | cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; |
6dbe713d | 875 | else |
0e9be4db | 876 | cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; |
6dbe713d SP |
877 | |
878 | if (w[0] == w[0] + cor) | |
879 | return (a > 0) ? w[0] : -w[0]; | |
880 | else | |
09544cbc | 881 | return __mpsin (orig, 0, true); |
6dbe713d | 882 | } |
e4d82761 | 883 | } |
e4d82761 | 884 | } |
6dbe713d | 885 | |
e4d82761 UD |
886 | /***************************************************************************/ |
887 | /* Routine compute sin(x+dx) (Double-Length number) where x in first or */ | |
888 | /* third quarter of unit circle.Routine receive also (right argument) the */ | |
889 | /* original value of x for computing error of result.And if result not */ | |
890 | /* accurate enough routine calls mpsin1 or dubsin */ | |
891 | /***************************************************************************/ | |
892 | ||
31d3cc00 UD |
893 | static double |
894 | SECTION | |
84ba214c | 895 | sloww1 (double x, double dx, double orig, int m) |
6dbe713d | 896 | { |
e4d82761 | 897 | mynumber u; |
392dd2de | 898 | double w[2], y, cor, res; |
6dbe713d | 899 | |
84ba214c SP |
900 | u.x = big + x; |
901 | y = x - (u.x - big); | |
0e9be4db | 902 | res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
6dbe713d SP |
903 | |
904 | if (res == res + cor) | |
84ba214c | 905 | return (m > 0) ? res : -res; |
6dbe713d SP |
906 | else |
907 | { | |
84ba214c | 908 | __dubsin (x, dx, w); |
6dbe713d SP |
909 | |
910 | if (w[1] > 0) | |
0e9be4db | 911 | cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
6dbe713d | 912 | else |
0e9be4db | 913 | cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
6dbe713d SP |
914 | |
915 | if (w[0] == w[0] + cor) | |
84ba214c | 916 | return (m > 0) ? w[0] : -w[0]; |
6dbe713d | 917 | else |
09544cbc | 918 | return __mpsin (orig, 0, true); |
6dbe713d | 919 | } |
e4d82761 | 920 | } |
6dbe713d | 921 | |
e4d82761 UD |
922 | /***************************************************************************/ |
923 | /* Routine compute sin(x+dx) (Double-Length number) where x in second or */ | |
924 | /* fourth quarter of unit circle.Routine receive also the original value */ | |
925 | /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ | |
926 | /* accurate enough routine calls mpsin1 or dubsin */ | |
927 | /***************************************************************************/ | |
928 | ||
31d3cc00 UD |
929 | static double |
930 | SECTION | |
6dbe713d SP |
931 | sloww2 (double x, double dx, double orig, int n) |
932 | { | |
e4d82761 | 933 | mynumber u; |
392dd2de | 934 | double w[2], y, cor, res; |
6dbe713d | 935 | |
5ff8d60e SP |
936 | u.x = big + x; |
937 | y = x - (u.x - big); | |
0e9be4db | 938 | res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
6dbe713d SP |
939 | |
940 | if (res == res + cor) | |
941 | return (n & 2) ? -res : res; | |
942 | else | |
943 | { | |
5ff8d60e | 944 | __docos (x, dx, w); |
6dbe713d SP |
945 | |
946 | if (w[1] > 0) | |
0e9be4db | 947 | cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
6dbe713d | 948 | else |
0e9be4db | 949 | cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
6dbe713d SP |
950 | |
951 | if (w[0] == w[0] + cor) | |
952 | return (n & 2) ? -w[0] : w[0]; | |
953 | else | |
09544cbc | 954 | return __mpsin (orig, 0, true); |
6dbe713d | 955 | } |
e4d82761 | 956 | } |
6dbe713d | 957 | |
e4d82761 UD |
958 | /***************************************************************************/ |
959 | /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ | |
960 | /* is small enough to use Taylor series around zero and (x+dx) */ | |
961 | /* in first or third quarter of unit circle.Routine receive also */ | |
962 | /* (right argument) the original value of x for computing error of */ | |
963 | /* result.And if result not accurate enough routine calls other routines */ | |
964 | /***************************************************************************/ | |
965 | ||
31d3cc00 UD |
966 | static double |
967 | SECTION | |
6dbe713d SP |
968 | bsloww (double x, double dx, double orig, int n) |
969 | { | |
4aafb73c SP |
970 | double res, cor, w[2]; |
971 | ||
972 | res = TAYLOR_SLOW (x, dx, cor); | |
6dbe713d SP |
973 | cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; |
974 | if (res == res + cor) | |
975 | return res; | |
976 | else | |
977 | { | |
978 | (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); | |
979 | if (w[1] > 0) | |
980 | cor = 1.000000001 * w[1] + 1.1e-24; | |
981 | else | |
982 | cor = 1.000000001 * w[1] - 1.1e-24; | |
983 | if (w[0] == w[0] + cor) | |
984 | return (x > 0) ? w[0] : -w[0]; | |
985 | else | |
09544cbc | 986 | return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
6dbe713d | 987 | } |
e4d82761 UD |
988 | } |
989 | ||
990 | /***************************************************************************/ | |
991 | /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ | |
992 | /* in first or third quarter of unit circle.Routine receive also */ | |
993 | /* (right argument) the original value of x for computing error of result.*/ | |
994 | /* And if result not accurate enough routine calls other routines */ | |
995 | /***************************************************************************/ | |
996 | ||
31d3cc00 UD |
997 | static double |
998 | SECTION | |
6dbe713d SP |
999 | bsloww1 (double x, double dx, double orig, int n) |
1000 | { | |
1001 | mynumber u; | |
392dd2de | 1002 | double w[2], y, cor, res; |
6dbe713d | 1003 | |
0e9be4db | 1004 | y = fabs (x); |
196f7f5d SP |
1005 | u.x = big + y; |
1006 | y = y - (u.x - big); | |
6dbe713d | 1007 | dx = (x > 0) ? dx : -dx; |
392dd2de | 1008 | res = do_sin_slow (u, y, dx, 1.1e-24, &cor); |
6dbe713d SP |
1009 | if (res == res + cor) |
1010 | return (x > 0) ? res : -res; | |
1011 | else | |
1012 | { | |
0e9be4db | 1013 | __dubsin (fabs (x), dx, w); |
6dbe713d SP |
1014 | |
1015 | if (w[1] > 0) | |
1016 | cor = 1.000000005 * w[1] + 1.1e-24; | |
1017 | else | |
1018 | cor = 1.000000005 * w[1] - 1.1e-24; | |
1019 | ||
1020 | if (w[0] == w[0] + cor) | |
1021 | return (x > 0) ? w[0] : -w[0]; | |
1022 | else | |
09544cbc | 1023 | return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
6dbe713d | 1024 | } |
e4d82761 UD |
1025 | } |
1026 | ||
1027 | /***************************************************************************/ | |
1028 | /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ | |
1029 | /* in second or fourth quarter of unit circle.Routine receive also the */ | |
1030 | /* original value and quarter(n= 1or 3)of x for computing error of result. */ | |
1031 | /* And if result not accurate enough routine calls other routines */ | |
1032 | /***************************************************************************/ | |
1033 | ||
31d3cc00 UD |
1034 | static double |
1035 | SECTION | |
6dbe713d SP |
1036 | bsloww2 (double x, double dx, double orig, int n) |
1037 | { | |
1038 | mynumber u; | |
392dd2de | 1039 | double w[2], y, cor, res; |
6dbe713d | 1040 | |
0e9be4db | 1041 | y = fabs (x); |
196f7f5d SP |
1042 | u.x = big + y; |
1043 | y = y - (u.x - big); | |
6dbe713d | 1044 | dx = (x > 0) ? dx : -dx; |
392dd2de | 1045 | res = do_cos_slow (u, y, dx, 1.1e-24, &cor); |
6dbe713d SP |
1046 | if (res == res + cor) |
1047 | return (n & 2) ? -res : res; | |
1048 | else | |
1049 | { | |
0e9be4db | 1050 | __docos (fabs (x), dx, w); |
6dbe713d SP |
1051 | |
1052 | if (w[1] > 0) | |
1053 | cor = 1.000000005 * w[1] + 1.1e-24; | |
1054 | else | |
1055 | cor = 1.000000005 * w[1] - 1.1e-24; | |
1056 | ||
1057 | if (w[0] == w[0] + cor) | |
1058 | return (n & 2) ? -w[0] : w[0]; | |
1059 | else | |
09544cbc | 1060 | return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); |
6dbe713d | 1061 | } |
e4d82761 UD |
1062 | } |
1063 | ||
1064 | /************************************************************************/ | |
1065 | /* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */ | |
1066 | /* precision and if still doesn't accurate enough by mpcos or docos */ | |
1067 | /************************************************************************/ | |
1068 | ||
31d3cc00 UD |
1069 | static double |
1070 | SECTION | |
6dbe713d SP |
1071 | cslow2 (double x) |
1072 | { | |
e4d82761 | 1073 | mynumber u; |
392dd2de | 1074 | double w[2], y, cor, res; |
6dbe713d | 1075 | |
0e9be4db | 1076 | y = fabs (x); |
196f7f5d SP |
1077 | u.x = big + y; |
1078 | y = y - (u.x - big); | |
392dd2de SP |
1079 | res = do_cos_slow (u, y, 0, 0, &cor); |
1080 | if (res == res + cor) | |
e4d82761 | 1081 | return res; |
6dbe713d SP |
1082 | else |
1083 | { | |
0e9be4db | 1084 | y = fabs (x); |
6dbe713d SP |
1085 | __docos (y, 0, w); |
1086 | if (w[0] == w[0] + 1.000000005 * w[1]) | |
1087 | return w[0]; | |
1088 | else | |
09544cbc | 1089 | return __mpcos (x, 0, false); |
6dbe713d | 1090 | } |
e4d82761 UD |
1091 | } |
1092 | ||
1093 | /***************************************************************************/ | |
1094 | /* Routine compute cos(x+dx) (Double-Length number) where x is small enough*/ | |
1095 | /* to use Taylor series around zero and (x+dx) .Routine receive also */ | |
1096 | /* (right argument) the original value of x for computing error of */ | |
1097 | /* result.And if result not accurate enough routine calls other routines */ | |
1098 | /***************************************************************************/ | |
1099 | ||
1100 | ||
31d3cc00 UD |
1101 | static double |
1102 | SECTION | |
6dbe713d SP |
1103 | csloww (double x, double dx, double orig) |
1104 | { | |
4aafb73c | 1105 | double y, t, res, cor, w[2], a, da, xn; |
8d561986 | 1106 | mynumber v; |
e4d82761 | 1107 | int4 n; |
6dbe713d | 1108 | |
6dbe713d | 1109 | /* Taylor series */ |
6c9642ed | 1110 | res = TAYLOR_SLOW (x, dx, cor); |
6dbe713d SP |
1111 | |
1112 | if (cor > 0) | |
0e9be4db | 1113 | cor = 1.0005 * cor + fabs (orig) * 3.1e-30; |
6dbe713d | 1114 | else |
0e9be4db | 1115 | cor = 1.0005 * cor - fabs (orig) * 3.1e-30; |
6dbe713d SP |
1116 | |
1117 | if (res == res + cor) | |
1118 | return res; | |
1119 | else | |
1120 | { | |
1121 | (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); | |
1122 | ||
1123 | if (w[1] > 0) | |
0e9be4db | 1124 | cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; |
6dbe713d | 1125 | else |
0e9be4db | 1126 | cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; |
6dbe713d SP |
1127 | |
1128 | if (w[0] == w[0] + cor) | |
1129 | return (x > 0) ? w[0] : -w[0]; | |
1130 | else | |
1131 | { | |
196f7f5d SP |
1132 | t = (orig * hpinv + toint); |
1133 | xn = t - toint; | |
6dbe713d | 1134 | v.x = t; |
196f7f5d | 1135 | y = (orig - xn * mp1) - xn * mp2; |
6dbe713d | 1136 | n = v.i[LOW_HALF] & 3; |
196f7f5d | 1137 | da = xn * pp3; |
6dbe713d SP |
1138 | t = y - da; |
1139 | da = (y - t) - da; | |
196f7f5d | 1140 | y = xn * pp4; |
6dbe713d SP |
1141 | a = t - y; |
1142 | da = ((t - a) - y) + da; | |
1143 | if (n == 1) | |
1144 | { | |
1145 | a = -a; | |
1146 | da = -da; | |
1147 | } | |
1148 | (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); | |
1149 | ||
1150 | if (w[1] > 0) | |
0e9be4db | 1151 | cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; |
6dbe713d | 1152 | else |
0e9be4db | 1153 | cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; |
6dbe713d SP |
1154 | |
1155 | if (w[0] == w[0] + cor) | |
1156 | return (a > 0) ? w[0] : -w[0]; | |
1157 | else | |
09544cbc | 1158 | return __mpcos (orig, 0, true); |
6dbe713d | 1159 | } |
e4d82761 | 1160 | } |
e4d82761 UD |
1161 | } |
1162 | ||
1163 | /***************************************************************************/ | |
1164 | /* Routine compute sin(x+dx) (Double-Length number) where x in first or */ | |
1165 | /* third quarter of unit circle.Routine receive also (right argument) the */ | |
1166 | /* original value of x for computing error of result.And if result not */ | |
1167 | /* accurate enough routine calls other routines */ | |
1168 | /***************************************************************************/ | |
1169 | ||
31d3cc00 UD |
1170 | static double |
1171 | SECTION | |
84ba214c | 1172 | csloww1 (double x, double dx, double orig, int m) |
6dbe713d | 1173 | { |
e4d82761 | 1174 | mynumber u; |
392dd2de | 1175 | double w[2], y, cor, res; |
6dbe713d | 1176 | |
84ba214c SP |
1177 | u.x = big + x; |
1178 | y = x - (u.x - big); | |
0e9be4db | 1179 | res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
6dbe713d SP |
1180 | |
1181 | if (res == res + cor) | |
84ba214c | 1182 | return (m > 0) ? res : -res; |
6dbe713d SP |
1183 | else |
1184 | { | |
84ba214c | 1185 | __dubsin (x, dx, w); |
6dbe713d | 1186 | if (w[1] > 0) |
0e9be4db | 1187 | cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
6dbe713d | 1188 | else |
0e9be4db | 1189 | cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
6dbe713d | 1190 | if (w[0] == w[0] + cor) |
84ba214c | 1191 | return (m > 0) ? w[0] : -w[0]; |
6dbe713d | 1192 | else |
09544cbc | 1193 | return __mpcos (orig, 0, true); |
6dbe713d | 1194 | } |
e4d82761 UD |
1195 | } |
1196 | ||
1197 | ||
1198 | /***************************************************************************/ | |
1199 | /* Routine compute sin(x+dx) (Double-Length number) where x in second or */ | |
1200 | /* fourth quarter of unit circle.Routine receive also the original value */ | |
1201 | /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ | |
1202 | /* accurate enough routine calls other routines */ | |
1203 | /***************************************************************************/ | |
1204 | ||
31d3cc00 UD |
1205 | static double |
1206 | SECTION | |
6dbe713d SP |
1207 | csloww2 (double x, double dx, double orig, int n) |
1208 | { | |
e4d82761 | 1209 | mynumber u; |
392dd2de | 1210 | double w[2], y, cor, res; |
6dbe713d | 1211 | |
5ff8d60e SP |
1212 | u.x = big + x; |
1213 | y = x - (u.x - big); | |
0e9be4db | 1214 | res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); |
6dbe713d SP |
1215 | |
1216 | if (res == res + cor) | |
1217 | return (n) ? -res : res; | |
1218 | else | |
1219 | { | |
5ff8d60e | 1220 | __docos (x, dx, w); |
6dbe713d | 1221 | if (w[1] > 0) |
0e9be4db | 1222 | cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); |
6dbe713d | 1223 | else |
0e9be4db | 1224 | cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); |
6dbe713d SP |
1225 | if (w[0] == w[0] + cor) |
1226 | return (n) ? -w[0] : w[0]; | |
1227 | else | |
09544cbc | 1228 | return __mpcos (orig, 0, true); |
6dbe713d | 1229 | } |
e4d82761 UD |
1230 | } |
1231 | ||
af968f62 | 1232 | #ifndef __cos |
ca58f1db | 1233 | weak_alias (__cos, cos) |
af968f62 UD |
1234 | # ifdef NO_LONG_DOUBLE |
1235 | strong_alias (__cos, __cosl) | |
1236 | weak_alias (__cos, cosl) | |
1237 | # endif | |
1238 | #endif | |
1239 | #ifndef __sin | |
ca58f1db | 1240 | weak_alias (__sin, sin) |
af968f62 | 1241 | # ifdef NO_LONG_DOUBLE |
ca58f1db UD |
1242 | strong_alias (__sin, __sinl) |
1243 | weak_alias (__sin, sinl) | |
af968f62 | 1244 | # endif |
cccda09f | 1245 | #endif |