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