OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_atanl.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_atanl.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 * ==================================================== |
11 */ | 11 */ |
12 /* | 12 /* |
13 * See comments in atan.c. | 13 * See comments in atan.c. |
14 * Converted to long double by David Schultz <das@FreeBSD.ORG>. | 14 * Converted to long double by David Schultz <das@FreeBSD.ORG>. |
15 */ | 15 */ |
16 | 16 |
17 #include "libm.h" | 17 #include "libm.h" |
18 | 18 |
19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
20 long double atanl(long double x) | 20 long double atanl(long double x) { |
21 { | 21 return atan(x); |
22 » return atan(x); | |
23 } | 22 } |
24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | 23 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
25 | 24 |
26 #if LDBL_MANT_DIG == 64 | 25 #if LDBL_MANT_DIG == 64 |
27 #define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | (u.i.m>>55 & 0xff)) | 26 #define EXPMAN(u) ((u.i.se & 0x7fff) << 8 | (u.i.m >> 55 & 0xff)) |
28 | 27 |
29 static const long double atanhi[] = { | 28 static const long double atanhi[] = { |
30 » 4.63647609000806116202e-01L, | 29 4.63647609000806116202e-01L, 7.85398163397448309628e-01L, |
31 » 7.85398163397448309628e-01L, | 30 9.82793723247329067960e-01L, 1.57079632679489661926e+00L, |
32 » 9.82793723247329067960e-01L, | |
33 » 1.57079632679489661926e+00L, | |
34 }; | 31 }; |
35 | 32 |
36 static const long double atanlo[] = { | 33 static const long double atanlo[] = { |
37 » 1.18469937025062860669e-20L, | 34 1.18469937025062860669e-20L, -1.25413940316708300586e-20L, |
38 » -1.25413940316708300586e-20L, | 35 2.55232234165405176172e-20L, -2.50827880633416601173e-20L, |
39 » 2.55232234165405176172e-20L, | |
40 » -2.50827880633416601173e-20L, | |
41 }; | 36 }; |
42 | 37 |
43 static const long double aT[] = { | 38 static const long double aT[] = { |
44 » 3.33333333333333333017e-01L, | 39 3.33333333333333333017e-01L, -1.99999999999999632011e-01L, |
45 » -1.99999999999999632011e-01L, | 40 1.42857142857046531280e-01L, -1.11111111100562372733e-01L, |
46 » 1.42857142857046531280e-01L, | 41 9.09090902935647302252e-02L, -7.69230552476207730353e-02L, |
47 » -1.11111111100562372733e-01L, | 42 6.66661718042406260546e-02L, -5.88158892835030888692e-02L, |
48 » 9.09090902935647302252e-02L, | 43 5.25499891539726639379e-02L, -4.70119845393155721494e-02L, |
49 » -7.69230552476207730353e-02L, | 44 4.03539201366454414072e-02L, -2.91303858419364158725e-02L, |
50 » 6.66661718042406260546e-02L, | 45 1.24822046299269234080e-02L, |
51 » -5.88158892835030888692e-02L, | |
52 » 5.25499891539726639379e-02L, | |
53 » -4.70119845393155721494e-02L, | |
54 » 4.03539201366454414072e-02L, | |
55 » -2.91303858419364158725e-02L, | |
56 » 1.24822046299269234080e-02L, | |
57 }; | 46 }; |
58 | 47 |
59 static long double T_even(long double x) | 48 static long double T_even(long double x) { |
60 { | 49 return aT[0] + |
61 » return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + | 50 x * (aT[2] + |
62 » » x * (aT[8] + x * (aT[10] + x * aT[12]))))); | 51 x * (aT[4] + |
| 52 x * (aT[6] + x * (aT[8] + x * (aT[10] + x * aT[12]))))); |
63 } | 53 } |
64 | 54 |
65 static long double T_odd(long double x) | 55 static long double T_odd(long double x) { |
66 { | 56 return aT[1] + |
67 » return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + | 57 x * (aT[3] + x * (aT[5] + x * (aT[7] + x * (aT[9] + x * aT[11])))); |
68 » » x * (aT[9] + x * aT[11])))); | |
69 } | 58 } |
70 #elif LDBL_MANT_DIG == 113 | 59 #elif LDBL_MANT_DIG == 113 |
71 #define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | u.i.top>>8) | 60 #define EXPMAN(u) ((u.i.se & 0x7fff) << 8 | u.i.top >> 8) |
72 | 61 |
73 const long double atanhi[] = { | 62 const long double atanhi[] = { |
74 » 4.63647609000806116214256231461214397e-01L, | 63 4.63647609000806116214256231461214397e-01L, |
75 » 7.85398163397448309615660845819875699e-01L, | 64 7.85398163397448309615660845819875699e-01L, |
76 » 9.82793723247329067985710611014666038e-01L, | 65 9.82793723247329067985710611014666038e-01L, |
77 » 1.57079632679489661923132169163975140e+00L, | 66 1.57079632679489661923132169163975140e+00L, |
78 }; | 67 }; |
79 | 68 |
80 const long double atanlo[] = { | 69 const long double atanlo[] = { |
81 » 4.89509642257333492668618435220297706e-36L, | 70 4.89509642257333492668618435220297706e-36L, |
82 » 2.16795253253094525619926100651083806e-35L, | 71 2.16795253253094525619926100651083806e-35L, |
83 » -2.31288434538183565909319952098066272e-35L, | 72 -2.31288434538183565909319952098066272e-35L, |
84 » 4.33590506506189051239852201302167613e-35L, | 73 4.33590506506189051239852201302167613e-35L, |
85 }; | 74 }; |
86 | 75 |
87 const long double aT[] = { | 76 const long double aT[] = { |
88 » 3.33333333333333333333333333333333125e-01L, | 77 3.33333333333333333333333333333333125e-01L, |
89 » -1.99999999999999999999999999999180430e-01L, | 78 -1.99999999999999999999999999999180430e-01L, |
90 » 1.42857142857142857142857142125269827e-01L, | 79 1.42857142857142857142857142125269827e-01L, |
91 » -1.11111111111111111111110834490810169e-01L, | 80 -1.11111111111111111111110834490810169e-01L, |
92 » 9.09090909090909090908522355708623681e-02L, | 81 9.09090909090909090908522355708623681e-02L, |
93 » -7.69230769230769230696553844935357021e-02L, | 82 -7.69230769230769230696553844935357021e-02L, |
94 » 6.66666666666666660390096773046256096e-02L, | 83 6.66666666666666660390096773046256096e-02L, |
95 » -5.88235294117646671706582985209643694e-02L, | 84 -5.88235294117646671706582985209643694e-02L, |
96 » 5.26315789473666478515847092020327506e-02L, | 85 5.26315789473666478515847092020327506e-02L, |
97 » -4.76190476189855517021024424991436144e-02L, | 86 -4.76190476189855517021024424991436144e-02L, |
98 » 4.34782608678695085948531993458097026e-02L, | 87 4.34782608678695085948531993458097026e-02L, |
99 » -3.99999999632663469330634215991142368e-02L, | 88 -3.99999999632663469330634215991142368e-02L, |
100 » 3.70370363987423702891250829918659723e-02L, | 89 3.70370363987423702891250829918659723e-02L, |
101 » -3.44827496515048090726669907612335954e-02L, | 90 -3.44827496515048090726669907612335954e-02L, |
102 » 3.22579620681420149871973710852268528e-02L, | 91 3.22579620681420149871973710852268528e-02L, |
103 » -3.03020767654269261041647570626778067e-02L, | 92 -3.03020767654269261041647570626778067e-02L, |
104 » 2.85641979882534783223403715930946138e-02L, | 93 2.85641979882534783223403715930946138e-02L, |
105 » -2.69824879726738568189929461383741323e-02L, | 94 -2.69824879726738568189929461383741323e-02L, |
106 » 2.54194698498808542954187110873675769e-02L, | 95 2.54194698498808542954187110873675769e-02L, |
107 » -2.35083879708189059926183138130183215e-02L, | 96 -2.35083879708189059926183138130183215e-02L, |
108 » 2.04832358998165364349957325067131428e-02L, | 97 2.04832358998165364349957325067131428e-02L, |
109 » -1.54489555488544397858507248612362957e-02L, | 98 -1.54489555488544397858507248612362957e-02L, |
110 » 8.64492360989278761493037861575248038e-03L, | 99 8.64492360989278761493037861575248038e-03L, |
111 » -2.58521121597609872727919154569765469e-03L, | 100 -2.58521121597609872727919154569765469e-03L, |
112 }; | 101 }; |
113 | 102 |
114 static long double T_even(long double x) | 103 static long double T_even(long double x) { |
115 { | 104 return ( |
116 » return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * (aT[8] + | 105 aT[0] + |
117 » » x * (aT[10] + x * (aT[12] + x * (aT[14] + x * (aT[16] + | 106 x * (aT[2] + |
118 » » x * (aT[18] + x * (aT[20] + x * aT[22]))))))))))); | 107 x * (aT[4] + |
| 108 x * (aT[6] + |
| 109 x * (aT[8] + |
| 110 x * (aT[10] + |
| 111 x * (aT[12] + |
| 112 x * (aT[14] + |
| 113 x * (aT[16] + |
| 114 x * (aT[18] + |
| 115 x * (aT[20] + |
| 116 x * aT[22]))))))))))); |
119 } | 117 } |
120 | 118 |
121 static long double T_odd(long double x) | 119 static long double T_odd(long double x) { |
122 { | 120 return ( |
123 » return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * (aT[9] + | 121 aT[1] + |
124 » » x * (aT[11] + x * (aT[13] + x * (aT[15] + x * (aT[17] + | 122 x * (aT[3] + |
125 » » x * (aT[19] + x * (aT[21] + x * aT[23]))))))))))); | 123 x * (aT[5] + |
| 124 x * (aT[7] + |
| 125 x * (aT[9] + |
| 126 x * (aT[11] + |
| 127 x * (aT[13] + |
| 128 x * (aT[15] + |
| 129 x * (aT[17] + |
| 130 x * (aT[19] + |
| 131 x * (aT[21] + |
| 132 x * aT[23]))))))))))); |
126 } | 133 } |
127 #endif | 134 #endif |
128 | 135 |
129 long double atanl(long double x) | 136 long double atanl(long double x) { |
130 { | 137 union ldshape u = {x}; |
131 » union ldshape u = {x}; | 138 long double w, s1, s2, z; |
132 » long double w, s1, s2, z; | 139 int id; |
133 » int id; | 140 unsigned e = u.i.se & 0x7fff; |
134 » unsigned e = u.i.se & 0x7fff; | 141 unsigned sign = u.i.se >> 15; |
135 » unsigned sign = u.i.se >> 15; | 142 unsigned expman; |
136 » unsigned expman; | |
137 | 143 |
138 » if (e >= 0x3fff + LDBL_MANT_DIG + 1) { /* if |x| is large, atan(x)~=pi/2
*/ | 144 if (e >= 0x3fff + LDBL_MANT_DIG + 1) { /* if |x| is large, atan(x)~=pi/2 */ |
139 » » if (isnan(x)) | 145 if (isnan(x)) |
140 » » » return x; | 146 return x; |
141 » » return sign ? -atanhi[3] : atanhi[3]; | 147 return sign ? -atanhi[3] : atanhi[3]; |
142 » } | 148 } |
143 » /* Extract the exponent and the first few bits of the mantissa. */ | 149 /* Extract the exponent and the first few bits of the mantissa. */ |
144 » expman = EXPMAN(u); | 150 expman = EXPMAN(u); |
145 » if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ | 151 if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ |
146 » » if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) { /* if |x| is small, at
anl(x)~=x */ | 152 if (e < |
147 » » » /* raise underflow if subnormal */ | 153 0x3fff - (LDBL_MANT_DIG + 1) / 2) { /* if |x| is small, atanl(x)~=x */ |
148 » » » if (e == 0) | 154 /* raise underflow if subnormal */ |
149 » » » » FORCE_EVAL((float)x); | 155 if (e == 0) |
150 » » » return x; | 156 FORCE_EVAL((float)x); |
151 » » } | 157 return x; |
152 » » id = -1; | 158 } |
153 » } else { | 159 id = -1; |
154 » » x = fabsl(x); | 160 } else { |
155 » » if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */ | 161 x = fabsl(x); |
156 » » » if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |
x| < 11/16 */ | 162 if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */ |
157 » » » » id = 0; | 163 if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */ |
158 » » » » x = (2.0*x-1.0)/(2.0+x); | 164 id = 0; |
159 » » » } else { /* 11/16 <= |x|
< 19/16 */ | 165 x = (2.0 * x - 1.0) / (2.0 + x); |
160 » » » » id = 1; | 166 } else { /* 11/16 <= |x| < 19/16 */ |
161 » » » » x = (x-1.0)/(x+1.0); | 167 id = 1; |
162 » » » } | 168 x = (x - 1.0) / (x + 1.0); |
163 » » } else { | 169 } |
164 » » » if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.43
75 */ | 170 } else { |
165 » » » » id = 2; | 171 if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.4375 */ |
166 » » » » x = (x-1.5)/(1.0+1.5*x); | 172 id = 2; |
167 » » » } else { /* 2.4375 <= |x
| */ | 173 x = (x - 1.5) / (1.0 + 1.5 * x); |
168 » » » » id = 3; | 174 } else { /* 2.4375 <= |x| */ |
169 » » » » x = -1.0/x; | 175 id = 3; |
170 » » » } | 176 x = -1.0 / x; |
171 » » } | 177 } |
172 » } | 178 } |
173 » /* end of argument reduction */ | 179 } |
174 » z = x*x; | 180 /* end of argument reduction */ |
175 » w = z*z; | 181 z = x * x; |
176 » /* break sum aT[i]z**(i+1) into odd and even poly */ | 182 w = z * z; |
177 » s1 = z*T_even(w); | 183 /* break sum aT[i]z**(i+1) into odd and even poly */ |
178 » s2 = w*T_odd(w); | 184 s1 = z * T_even(w); |
179 » if (id < 0) | 185 s2 = w * T_odd(w); |
180 » » return x - x*(s1+s2); | 186 if (id < 0) |
181 » z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); | 187 return x - x * (s1 + s2); |
182 » return sign ? -z : z; | 188 z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x); |
| 189 return sign ? -z : z; |
183 } | 190 } |
184 #endif | 191 #endif |
OLD | NEW |