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