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