OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_acos.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_acos.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 SunSoft, a Sun Microsystems, Inc. business. | 6 * Developed at SunSoft, 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 * ==================================================== |
(...skipping 17 matching lines...) Expand all Loading... |
28 * | 28 * |
29 * Special cases: | 29 * Special cases: |
30 * if x is NaN, return x itself; | 30 * if x is NaN, return x itself; |
31 * if |x|>1, return NaN with invalid signal. | 31 * if |x|>1, return NaN with invalid signal. |
32 * | 32 * |
33 * Function needed: sqrt | 33 * Function needed: sqrt |
34 */ | 34 */ |
35 | 35 |
36 #include "libm.h" | 36 #include "libm.h" |
37 | 37 |
38 static const double | 38 static const double pio2_hi = |
39 pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ | 39 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ |
40 pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ | 40 pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ |
41 pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ | 41 pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ |
42 pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ | 42 pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ |
43 pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ | 43 pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ |
44 pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ | 44 pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ |
45 pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ | 45 pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ |
46 pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ | 46 pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ |
47 qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ | 47 qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ |
48 qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ | 48 qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ |
49 qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ | 49 qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ |
50 qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ | 50 qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ |
51 | 51 |
52 static double R(double z) | 52 static double R(double z) { |
53 { | 53 double_t p, q; |
54 » double_t p, q; | 54 p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); |
55 » p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); | 55 q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); |
56 » q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4))); | 56 return p / q; |
57 » return p/q; | |
58 } | 57 } |
59 | 58 |
60 double acos(double x) | 59 double acos(double x) { |
61 { | 60 double z, w, s, c, df; |
62 » double z,w,s,c,df; | 61 uint32_t hx, ix; |
63 » uint32_t hx,ix; | |
64 | 62 |
65 » GET_HIGH_WORD(hx, x); | 63 GET_HIGH_WORD(hx, x); |
66 » ix = hx & 0x7fffffff; | 64 ix = hx & 0x7fffffff; |
67 » /* |x| >= 1 or nan */ | 65 /* |x| >= 1 or nan */ |
68 » if (ix >= 0x3ff00000) { | 66 if (ix >= 0x3ff00000) { |
69 » » uint32_t lx; | 67 uint32_t lx; |
70 | 68 |
71 » » GET_LOW_WORD(lx,x); | 69 GET_LOW_WORD(lx, x); |
72 » » if ((ix-0x3ff00000 | lx) == 0) { | 70 if ((ix - 0x3ff00000 | lx) == 0) { |
73 » » » /* acos(1)=0, acos(-1)=pi */ | 71 /* acos(1)=0, acos(-1)=pi */ |
74 » » » if (hx >> 31) | 72 if (hx >> 31) |
75 » » » » return 2*pio2_hi + 0x1p-120f; | 73 return 2 * pio2_hi + 0x1p-120f; |
76 » » » return 0; | 74 return 0; |
77 » » } | 75 } |
78 » » return 0/(x-x); | 76 return 0 / (x - x); |
79 » } | 77 } |
80 » /* |x| < 0.5 */ | 78 /* |x| < 0.5 */ |
81 » if (ix < 0x3fe00000) { | 79 if (ix < 0x3fe00000) { |
82 » » if (ix <= 0x3c600000) /* |x| < 2**-57 */ | 80 if (ix <= 0x3c600000) /* |x| < 2**-57 */ |
83 » » » return pio2_hi + 0x1p-120f; | 81 return pio2_hi + 0x1p-120f; |
84 » » return pio2_hi - (x - (pio2_lo-x*R(x*x))); | 82 return pio2_hi - (x - (pio2_lo - x * R(x * x))); |
85 » } | 83 } |
86 » /* x < -0.5 */ | 84 /* x < -0.5 */ |
87 » if (hx >> 31) { | 85 if (hx >> 31) { |
88 » » z = (1.0+x)*0.5; | 86 z = (1.0 + x) * 0.5; |
89 » » s = sqrt(z); | 87 s = sqrt(z); |
90 » » w = R(z)*s-pio2_lo; | 88 w = R(z) * s - pio2_lo; |
91 » » return 2*(pio2_hi - (s+w)); | 89 return 2 * (pio2_hi - (s + w)); |
92 » } | 90 } |
93 » /* x > 0.5 */ | 91 /* x > 0.5 */ |
94 » z = (1.0-x)*0.5; | 92 z = (1.0 - x) * 0.5; |
95 » s = sqrt(z); | 93 s = sqrt(z); |
96 » df = s; | 94 df = s; |
97 » SET_LOW_WORD(df,0); | 95 SET_LOW_WORD(df, 0); |
98 » c = (z-df*df)/(s+df); | 96 c = (z - df * df) / (s + df); |
99 » w = R(z)*s+c; | 97 w = R(z) * s + c; |
100 » return 2*(df+w); | 98 return 2 * (df + w); |
101 } | 99 } |
OLD | NEW |