| 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 |