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