OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.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 | 15 |
16 #include "libm.h" | 16 #include "libm.h" |
17 | 17 |
18 static const float | 18 static const float pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */ |
19 pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */ | 19 pio2_lo = 7.5497894159e-08, /* 0x33a22168 */ |
20 pio2_lo = 7.5497894159e-08, /* 0x33a22168 */ | 20 pS0 = 1.6666586697e-01, pS1 = -4.2743422091e-02, pS2 = -8.6563630030e-03, |
21 pS0 = 1.6666586697e-01, | 21 qS1 = -7.0662963390e-01; |
22 pS1 = -4.2743422091e-02, | |
23 pS2 = -8.6563630030e-03, | |
24 qS1 = -7.0662963390e-01; | |
25 | 22 |
26 static float R(float z) | 23 static float R(float z) { |
27 { | 24 float_t p, q; |
28 » float_t p, q; | 25 p = z * (pS0 + z * (pS1 + z * pS2)); |
29 » p = z*(pS0+z*(pS1+z*pS2)); | 26 q = 1.0f + z * qS1; |
30 » q = 1.0f+z*qS1; | 27 return p / q; |
31 » return p/q; | |
32 } | 28 } |
33 | 29 |
34 float acosf(float x) | 30 float acosf(float x) { |
35 { | 31 float z, w, s, c, df; |
36 » float z,w,s,c,df; | 32 uint32_t hx, ix; |
37 » uint32_t hx,ix; | |
38 | 33 |
39 » GET_FLOAT_WORD(hx, x); | 34 GET_FLOAT_WORD(hx, x); |
40 » ix = hx & 0x7fffffff; | 35 ix = hx & 0x7fffffff; |
41 » /* |x| >= 1 or nan */ | 36 /* |x| >= 1 or nan */ |
42 » if (ix >= 0x3f800000) { | 37 if (ix >= 0x3f800000) { |
43 » » if (ix == 0x3f800000) { | 38 if (ix == 0x3f800000) { |
44 » » » if (hx >> 31) | 39 if (hx >> 31) |
45 » » » » return 2*pio2_hi + 0x1p-120f; | 40 return 2 * pio2_hi + 0x1p-120f; |
46 » » » return 0; | 41 return 0; |
47 » » } | 42 } |
48 » » return 0/(x-x); | 43 return 0 / (x - x); |
49 » } | 44 } |
50 » /* |x| < 0.5 */ | 45 /* |x| < 0.5 */ |
51 » if (ix < 0x3f000000) { | 46 if (ix < 0x3f000000) { |
52 » » if (ix <= 0x32800000) /* |x| < 2**-26 */ | 47 if (ix <= 0x32800000) /* |x| < 2**-26 */ |
53 » » » return pio2_hi + 0x1p-120f; | 48 return pio2_hi + 0x1p-120f; |
54 » » return pio2_hi - (x - (pio2_lo-x*R(x*x))); | 49 return pio2_hi - (x - (pio2_lo - x * R(x * x))); |
55 » } | 50 } |
56 » /* x < -0.5 */ | 51 /* x < -0.5 */ |
57 » if (hx >> 31) { | 52 if (hx >> 31) { |
58 » » z = (1+x)*0.5f; | 53 z = (1 + x) * 0.5f; |
59 » » s = sqrtf(z); | 54 s = sqrtf(z); |
60 » » w = R(z)*s-pio2_lo; | 55 w = R(z) * s - pio2_lo; |
61 » » return 2*(pio2_hi - (s+w)); | 56 return 2 * (pio2_hi - (s + w)); |
62 » } | 57 } |
63 » /* x > 0.5 */ | 58 /* x > 0.5 */ |
64 » z = (1-x)*0.5f; | 59 z = (1 - x) * 0.5f; |
65 » s = sqrtf(z); | 60 s = sqrtf(z); |
66 » GET_FLOAT_WORD(hx,s); | 61 GET_FLOAT_WORD(hx, s); |
67 » SET_FLOAT_WORD(df,hx&0xfffff000); | 62 SET_FLOAT_WORD(df, hx & 0xfffff000); |
68 » c = (z-df*df)/(s+df); | 63 c = (z - df * df) / (s + df); |
69 » w = R(z)*s+c; | 64 w = R(z) * s + c; |
70 » return 2*(df+w); | 65 return 2 * (df + w); |
71 } | 66 } |
OLD | NEW |