OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */ |
2 /* | 2 /* |
3 * ==================================================== | 3 * ==================================================== |
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 * | 5 * |
6 * Developed at SunPro, a Sun Microsystems, Inc. business. | 6 * Developed at SunPro, a Sun Microsystems, Inc. business. |
7 * Permission to use, copy, modify, and distribute this | 7 * Permission to use, copy, modify, and distribute this |
8 * software is freely granted, provided that this notice | 8 * software is freely granted, provided that this notice |
9 * is preserved. | 9 * is preserved. |
10 * ==================================================== | 10 * ==================================================== |
11 */ | 11 */ |
12 | 12 |
13 #define _GNU_SOURCE | 13 #define _GNU_SOURCE |
14 #include "libm.h" | 14 #include "libm.h" |
15 | 15 |
16 void sincos(double x, double *sin, double *cos) | 16 void sincos(double x, double* sin, double* cos) { |
17 { | 17 double y[2], s, c; |
18 » double y[2], s, c; | 18 uint32_t ix; |
19 » uint32_t ix; | 19 unsigned n; |
20 » unsigned n; | |
21 | 20 |
22 » GET_HIGH_WORD(ix, x); | 21 GET_HIGH_WORD(ix, x); |
23 » ix &= 0x7fffffff; | 22 ix &= 0x7fffffff; |
24 | 23 |
25 » /* |x| ~< pi/4 */ | 24 /* |x| ~< pi/4 */ |
26 » if (ix <= 0x3fe921fb) { | 25 if (ix <= 0x3fe921fb) { |
27 » » /* if |x| < 2**-27 * sqrt(2) */ | 26 /* if |x| < 2**-27 * sqrt(2) */ |
28 » » if (ix < 0x3e46a09e) { | 27 if (ix < 0x3e46a09e) { |
29 » » » /* raise inexact if x!=0 and underflow if subnormal */ | 28 /* raise inexact if x!=0 and underflow if subnormal */ |
30 » » » FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f); | 29 FORCE_EVAL(ix < 0x00100000 ? x / 0x1p120f : x + 0x1p120f); |
31 » » » *sin = x; | 30 *sin = x; |
32 » » » *cos = 1.0; | 31 *cos = 1.0; |
33 » » » return; | 32 return; |
34 » » } | 33 } |
35 » » *sin = __sin(x, 0.0, 0); | 34 *sin = __sin(x, 0.0, 0); |
36 » » *cos = __cos(x, 0.0); | 35 *cos = __cos(x, 0.0); |
37 » » return; | 36 return; |
38 » } | 37 } |
39 | 38 |
40 » /* sincos(Inf or NaN) is NaN */ | 39 /* sincos(Inf or NaN) is NaN */ |
41 » if (ix >= 0x7ff00000) { | 40 if (ix >= 0x7ff00000) { |
42 » » *sin = *cos = x - x; | 41 *sin = *cos = x - x; |
43 » » return; | 42 return; |
44 » } | 43 } |
45 | 44 |
46 » /* argument reduction needed */ | 45 /* argument reduction needed */ |
47 » n = __rem_pio2(x, y); | 46 n = __rem_pio2(x, y); |
48 » s = __sin(y[0], y[1], 1); | 47 s = __sin(y[0], y[1], 1); |
49 » c = __cos(y[0], y[1]); | 48 c = __cos(y[0], y[1]); |
50 » switch (n&3) { | 49 switch (n & 3) { |
51 » case 0: | 50 case 0: |
52 » » *sin = s; | 51 *sin = s; |
53 » » *cos = c; | 52 *cos = c; |
54 » » break; | 53 break; |
55 » case 1: | 54 case 1: |
56 » » *sin = c; | 55 *sin = c; |
57 » » *cos = -s; | 56 *cos = -s; |
58 » » break; | 57 break; |
59 » case 2: | 58 case 2: |
60 » » *sin = -s; | 59 *sin = -s; |
61 » » *cos = -c; | 60 *cos = -c; |
62 » » break; | 61 break; |
63 » case 3: | 62 case 3: |
64 » default: | 63 default: |
65 » » *sin = -c; | 64 *sin = -c; |
66 » » *cos = s; | 65 *cos = s; |
67 » » break; | 66 break; |
68 » } | 67 } |
69 } | 68 } |
OLD | NEW |