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 |