| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_atan.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 SunPro, a Sun Microsystems, Inc. business. | 6 * Developed at SunPro, 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 11 matching lines...) Expand all Loading... |
| 22 * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) | 22 * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) |
| 23 * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) | 23 * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) |
| 24 * | 24 * |
| 25 * Constants: | 25 * Constants: |
| 26 * The hexadecimal values are the intended ones for the following | 26 * The hexadecimal values are the intended ones for the following |
| 27 * constants. The decimal values may be used, provided that the | 27 * constants. The decimal values may be used, provided that the |
| 28 * compiler will convert from decimal to binary accurately enough | 28 * compiler will convert from decimal to binary accurately enough |
| 29 * to produce the hexadecimal values shown. | 29 * to produce the hexadecimal values shown. |
| 30 */ | 30 */ |
| 31 | 31 |
| 32 | |
| 33 #include "libm.h" | 32 #include "libm.h" |
| 34 | 33 |
| 35 static const double atanhi[] = { | 34 static const double atanhi[] = { |
| 36 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ | 35 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ |
| 37 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ | 36 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ |
| 38 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ | 37 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ |
| 39 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ | 38 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ |
| 40 }; | 39 }; |
| 41 | 40 |
| 42 static const double atanlo[] = { | 41 static const double atanlo[] = { |
| 43 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ | 42 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ |
| 44 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ | 43 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ |
| 45 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ | 44 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ |
| 46 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ | 45 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ |
| 47 }; | 46 }; |
| 48 | 47 |
| 49 static const double aT[] = { | 48 static const double aT[] = { |
| 50 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ | 49 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ |
| 51 -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ | 50 -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ |
| 52 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ | 51 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ |
| 53 -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ | 52 -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ |
| 54 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ | 53 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ |
| 55 -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ | 54 -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ |
| 56 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ | 55 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ |
| 57 -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ | 56 -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ |
| 58 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ | 57 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ |
| 59 -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ | 58 -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ |
| 60 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ | 59 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ |
| 61 }; | 60 }; |
| 62 | 61 |
| 63 double atan(double x) | 62 double atan(double x) { |
| 64 { | 63 double_t w, s1, s2, z; |
| 65 » double_t w,s1,s2,z; | 64 uint32_t ix, sign; |
| 66 » uint32_t ix,sign; | 65 int id; |
| 67 » int id; | |
| 68 | 66 |
| 69 » GET_HIGH_WORD(ix, x); | 67 GET_HIGH_WORD(ix, x); |
| 70 » sign = ix >> 31; | 68 sign = ix >> 31; |
| 71 » ix &= 0x7fffffff; | 69 ix &= 0x7fffffff; |
| 72 » if (ix >= 0x44100000) { /* if |x| >= 2^66 */ | 70 if (ix >= 0x44100000) { /* if |x| >= 2^66 */ |
| 73 » » if (isnan(x)) | 71 if (isnan(x)) |
| 74 » » » return x; | 72 return x; |
| 75 » » z = atanhi[3] + 0x1p-120f; | 73 z = atanhi[3] + 0x1p-120f; |
| 76 » » return sign ? -z : z; | 74 return sign ? -z : z; |
| 77 » } | 75 } |
| 78 » if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ | 76 if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ |
| 79 » » if (ix < 0x3e400000) { /* |x| < 2^-27 */ | 77 if (ix < 0x3e400000) { /* |x| < 2^-27 */ |
| 80 » » » if (ix < 0x00100000) | 78 if (ix < 0x00100000) |
| 81 » » » » /* raise underflow for subnormal x */ | 79 /* raise underflow for subnormal x */ |
| 82 » » » » FORCE_EVAL((float)x); | 80 FORCE_EVAL((float)x); |
| 83 » » » return x; | 81 return x; |
| 84 » » } | 82 } |
| 85 » » id = -1; | 83 id = -1; |
| 86 » } else { | 84 } else { |
| 87 » » x = fabs(x); | 85 x = fabs(x); |
| 88 » » if (ix < 0x3ff30000) { /* |x| < 1.1875 */ | 86 if (ix < 0x3ff30000) { /* |x| < 1.1875 */ |
| 89 » » » if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */ | 87 if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */ |
| 90 » » » » id = 0; | 88 id = 0; |
| 91 » » » » x = (2.0*x-1.0)/(2.0+x); | 89 x = (2.0 * x - 1.0) / (2.0 + x); |
| 92 » » » } else { /* 11/16 <= |x| < 19/16 */ | 90 } else { /* 11/16 <= |x| < 19/16 */ |
| 93 » » » » id = 1; | 91 id = 1; |
| 94 » » » » x = (x-1.0)/(x+1.0); | 92 x = (x - 1.0) / (x + 1.0); |
| 95 » » » } | 93 } |
| 96 » » } else { | 94 } else { |
| 97 » » » if (ix < 0x40038000) { /* |x| < 2.4375 */ | 95 if (ix < 0x40038000) { /* |x| < 2.4375 */ |
| 98 » » » » id = 2; | 96 id = 2; |
| 99 » » » » x = (x-1.5)/(1.0+1.5*x); | 97 x = (x - 1.5) / (1.0 + 1.5 * x); |
| 100 » » » } else { /* 2.4375 <= |x| < 2^66 */ | 98 } else { /* 2.4375 <= |x| < 2^66 */ |
| 101 » » » » id = 3; | 99 id = 3; |
| 102 » » » » x = -1.0/x; | 100 x = -1.0 / x; |
| 103 » » » } | 101 } |
| 104 » » } | 102 } |
| 105 » } | 103 } |
| 106 » /* end of argument reduction */ | 104 /* end of argument reduction */ |
| 107 » z = x*x; | 105 z = x * x; |
| 108 » w = z*z; | 106 w = z * z; |
| 109 » /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ | 107 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ |
| 110 » s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); | 108 s1 = z * (aT[0] + |
| 111 » s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); | 109 w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); |
| 112 » if (id < 0) | 110 s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); |
| 113 » » return x - x*(s1+s2); | 111 if (id < 0) |
| 114 » z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x); | 112 return x - x * (s1 + s2); |
| 115 » return sign ? -z : z; | 113 z = atanhi[id] - (x * (s1 + s2) - atanlo[id] - x); |
| 114 return sign ? -z : z; |
| 116 } | 115 } |
| OLD | NEW |