OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */ |
2 /* | 2 /* |
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 */ | 4 */ |
5 /* | 5 /* |
6 * ==================================================== | 6 * ==================================================== |
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 * | 8 * |
9 * Developed at SunPro, a Sun Microsystems, Inc. business. | 9 * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 * Permission to use, copy, modify, and distribute this | 10 * Permission to use, copy, modify, and distribute this |
11 * software is freely granted, provided that this notice | 11 * software is freely granted, provided that this notice |
12 * is preserved. | 12 * is preserved. |
13 * ==================================================== | 13 * ==================================================== |
14 */ | 14 */ |
15 #include "libm.h" | 15 #include "libm.h" |
16 | 16 |
17 static const double | 17 static const double pio2 = 1.570796326794896558e+00; |
18 pio2 = 1.570796326794896558e+00; | |
19 | 18 |
20 static const float | 19 static const float |
21 /* coefficients for R(x^2) */ | 20 /* coefficients for R(x^2) */ |
22 pS0 = 1.6666586697e-01, | 21 pS0 = 1.6666586697e-01, |
23 pS1 = -4.2743422091e-02, | 22 pS1 = -4.2743422091e-02, pS2 = -8.6563630030e-03, qS1 = -7.0662963390e-01; |
24 pS2 = -8.6563630030e-03, | |
25 qS1 = -7.0662963390e-01; | |
26 | 23 |
27 static float R(float z) | 24 static float R(float z) { |
28 { | 25 float_t p, q; |
29 » float_t p, q; | 26 p = z * (pS0 + z * (pS1 + z * pS2)); |
30 » p = z*(pS0+z*(pS1+z*pS2)); | 27 q = 1.0f + z * qS1; |
31 » q = 1.0f+z*qS1; | 28 return p / q; |
32 » return p/q; | |
33 } | 29 } |
34 | 30 |
35 float asinf(float x) | 31 float asinf(float x) { |
36 { | 32 double s; |
37 » double s; | 33 float z; |
38 » float z; | 34 uint32_t hx, ix; |
39 » uint32_t hx,ix; | |
40 | 35 |
41 » GET_FLOAT_WORD(hx, x); | 36 GET_FLOAT_WORD(hx, x); |
42 » ix = hx & 0x7fffffff; | 37 ix = hx & 0x7fffffff; |
43 » if (ix >= 0x3f800000) { /* |x| >= 1 */ | 38 if (ix >= 0x3f800000) { /* |x| >= 1 */ |
44 » » if (ix == 0x3f800000) /* |x| == 1 */ | 39 if (ix == 0x3f800000) /* |x| == 1 */ |
45 » » » return x*pio2 + 0x1p-120f; /* asin(+-1) = +-pi/2 with i
nexact */ | 40 return x * pio2 + 0x1p-120f; /* asin(+-1) = +-pi/2 with inexact */ |
46 » » return 0/(x-x); /* asin(|x|>1) is NaN */ | 41 return 0 / (x - x); /* asin(|x|>1) is NaN */ |
47 » } | 42 } |
48 » if (ix < 0x3f000000) { /* |x| < 0.5 */ | 43 if (ix < 0x3f000000) { /* |x| < 0.5 */ |
49 » » /* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */ | 44 /* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */ |
50 » » if (ix < 0x39800000 && ix >= 0x00800000) | 45 if (ix < 0x39800000 && ix >= 0x00800000) |
51 » » » return x; | 46 return x; |
52 » » return x + x*R(x*x); | 47 return x + x * R(x * x); |
53 » } | 48 } |
54 » /* 1 > |x| >= 0.5 */ | 49 /* 1 > |x| >= 0.5 */ |
55 » z = (1 - fabsf(x))*0.5f; | 50 z = (1 - fabsf(x)) * 0.5f; |
56 » s = sqrt(z); | 51 s = sqrt(z); |
57 » x = pio2 - 2*(s+s*R(z)); | 52 x = pio2 - 2 * (s + s * R(z)); |
58 » if (hx >> 31) | 53 if (hx >> 31) |
59 » » return -x; | 54 return -x; |
60 » return x; | 55 return x; |
61 } | 56 } |
OLD | NEW |