|
1 // boost atanh.hpp header file |
|
2 |
|
3 // (C) Copyright Hubert Holin 2001. |
|
4 // Distributed under the Boost Software License, Version 1.0. (See |
|
5 // accompanying file LICENSE_1_0.txt or copy at |
|
6 // http://www.boost.org/LICENSE_1_0.txt) |
|
7 |
|
8 // See http://www.boost.org for updates, documentation, and revision history. |
|
9 |
|
10 #ifndef BOOST_ATANH_HPP |
|
11 #define BOOST_ATANH_HPP |
|
12 |
|
13 |
|
14 #include <cmath> |
|
15 #include <limits> |
|
16 #include <string> |
|
17 #include <stdexcept> |
|
18 |
|
19 |
|
20 #include <boost/config.hpp> |
|
21 |
|
22 |
|
23 // This is the inverse of the hyperbolic tangent function. |
|
24 |
|
25 namespace boost |
|
26 { |
|
27 namespace math |
|
28 { |
|
29 #if defined(__GNUC__) && (__GNUC__ < 3) |
|
30 // gcc 2.x ignores function scope using declarations, |
|
31 // put them in the scope of the enclosing namespace instead: |
|
32 |
|
33 using ::std::abs; |
|
34 using ::std::sqrt; |
|
35 using ::std::log; |
|
36 |
|
37 using ::std::numeric_limits; |
|
38 #endif |
|
39 |
|
40 #if defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) |
|
41 // This is the main fare |
|
42 |
|
43 template<typename T> |
|
44 inline T atanh(const T x) |
|
45 { |
|
46 using ::std::abs; |
|
47 using ::std::sqrt; |
|
48 using ::std::log; |
|
49 |
|
50 using ::std::numeric_limits; |
|
51 |
|
52 T const one = static_cast<T>(1); |
|
53 T const two = static_cast<T>(2); |
|
54 |
|
55 static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon()); |
|
56 static T const taylor_n_bound = sqrt(taylor_2_bound); |
|
57 |
|
58 if (x < -one) |
|
59 { |
|
60 if (numeric_limits<T>::has_quiet_NaN) |
|
61 { |
|
62 return(numeric_limits<T>::quiet_NaN()); |
|
63 } |
|
64 else |
|
65 { |
|
66 ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!"); |
|
67 ::std::domain_error bad_argument(error_reporting); |
|
68 |
|
69 throw(bad_argument); |
|
70 } |
|
71 } |
|
72 else if (x < -one+numeric_limits<T>::epsilon()) |
|
73 { |
|
74 if (numeric_limits<T>::has_infinity) |
|
75 { |
|
76 return(-numeric_limits<T>::infinity()); |
|
77 } |
|
78 else |
|
79 { |
|
80 ::std::string error_reporting("Argument to atanh is -1 (result: -Infinity)!"); |
|
81 ::std::out_of_range bad_argument(error_reporting); |
|
82 |
|
83 throw(bad_argument); |
|
84 } |
|
85 } |
|
86 else if (x > +one-numeric_limits<T>::epsilon()) |
|
87 { |
|
88 if (numeric_limits<T>::has_infinity) |
|
89 { |
|
90 return(+numeric_limits<T>::infinity()); |
|
91 } |
|
92 else |
|
93 { |
|
94 ::std::string error_reporting("Argument to atanh is +1 (result: +Infinity)!"); |
|
95 ::std::out_of_range bad_argument(error_reporting); |
|
96 |
|
97 throw(bad_argument); |
|
98 } |
|
99 } |
|
100 else if (x > +one) |
|
101 { |
|
102 if (numeric_limits<T>::has_quiet_NaN) |
|
103 { |
|
104 return(numeric_limits<T>::quiet_NaN()); |
|
105 } |
|
106 else |
|
107 { |
|
108 ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!"); |
|
109 ::std::domain_error bad_argument(error_reporting); |
|
110 |
|
111 throw(bad_argument); |
|
112 } |
|
113 } |
|
114 else if (abs(x) >= taylor_n_bound) |
|
115 { |
|
116 return(log( (one + x) / (one - x) ) / two); |
|
117 } |
|
118 else |
|
119 { |
|
120 // approximation by taylor series in x at 0 up to order 2 |
|
121 T result = x; |
|
122 |
|
123 if (abs(x) >= taylor_2_bound) |
|
124 { |
|
125 T x3 = x*x*x; |
|
126 |
|
127 // approximation by taylor series in x at 0 up to order 4 |
|
128 result += x3/static_cast<T>(3); |
|
129 } |
|
130 |
|
131 return(result); |
|
132 } |
|
133 } |
|
134 #else |
|
135 // These are implementation details (for main fare see below) |
|
136 |
|
137 namespace detail |
|
138 { |
|
139 template < |
|
140 typename T, |
|
141 bool InfinitySupported |
|
142 > |
|
143 struct atanh_helper1_t |
|
144 { |
|
145 static T get_pos_infinity() |
|
146 { |
|
147 return(+::std::numeric_limits<T>::infinity()); |
|
148 } |
|
149 |
|
150 static T get_neg_infinity() |
|
151 { |
|
152 return(-::std::numeric_limits<T>::infinity()); |
|
153 } |
|
154 }; // boost::math::detail::atanh_helper1_t |
|
155 |
|
156 |
|
157 template<typename T> |
|
158 struct atanh_helper1_t<T, false> |
|
159 { |
|
160 static T get_pos_infinity() |
|
161 { |
|
162 ::std::string error_reporting("Argument to atanh is +1 (result: +Infinity)!"); |
|
163 ::std::out_of_range bad_argument(error_reporting); |
|
164 |
|
165 throw(bad_argument); |
|
166 } |
|
167 |
|
168 static T get_neg_infinity() |
|
169 { |
|
170 ::std::string error_reporting("Argument to atanh is -1 (result: -Infinity)!"); |
|
171 ::std::out_of_range bad_argument(error_reporting); |
|
172 |
|
173 throw(bad_argument); |
|
174 } |
|
175 }; // boost::math::detail::atanh_helper1_t |
|
176 |
|
177 |
|
178 template < |
|
179 typename T, |
|
180 bool QuietNanSupported |
|
181 > |
|
182 struct atanh_helper2_t |
|
183 { |
|
184 static T get_NaN() |
|
185 { |
|
186 return(::std::numeric_limits<T>::quiet_NaN()); |
|
187 } |
|
188 }; // boost::detail::atanh_helper2_t |
|
189 |
|
190 |
|
191 template<typename T> |
|
192 struct atanh_helper2_t<T, false> |
|
193 { |
|
194 static T get_NaN() |
|
195 { |
|
196 ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!"); |
|
197 ::std::domain_error bad_argument(error_reporting); |
|
198 |
|
199 throw(bad_argument); |
|
200 } |
|
201 }; // boost::detail::atanh_helper2_t |
|
202 } // boost::detail |
|
203 |
|
204 |
|
205 // This is the main fare |
|
206 |
|
207 template<typename T> |
|
208 inline T atanh(const T x) |
|
209 { |
|
210 using ::std::abs; |
|
211 using ::std::sqrt; |
|
212 using ::std::log; |
|
213 |
|
214 using ::std::numeric_limits; |
|
215 |
|
216 typedef detail::atanh_helper1_t<T, ::std::numeric_limits<T>::has_infinity> helper1_type; |
|
217 typedef detail::atanh_helper2_t<T, ::std::numeric_limits<T>::has_quiet_NaN> helper2_type; |
|
218 |
|
219 |
|
220 T const one = static_cast<T>(1); |
|
221 T const two = static_cast<T>(2); |
|
222 |
|
223 static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon()); |
|
224 static T const taylor_n_bound = sqrt(taylor_2_bound); |
|
225 |
|
226 if (x < -one) |
|
227 { |
|
228 return(helper2_type::get_NaN()); |
|
229 } |
|
230 else if (x < -one+numeric_limits<T>::epsilon()) |
|
231 { |
|
232 return(helper1_type::get_neg_infinity()); |
|
233 } |
|
234 else if (x > +one-numeric_limits<T>::epsilon()) |
|
235 { |
|
236 return(helper1_type::get_pos_infinity()); |
|
237 } |
|
238 else if (x > +one) |
|
239 { |
|
240 return(helper2_type::get_NaN()); |
|
241 } |
|
242 else if (abs(x) >= taylor_n_bound) |
|
243 { |
|
244 return(log( (one + x) / (one - x) ) / two); |
|
245 } |
|
246 else |
|
247 { |
|
248 // approximation by taylor series in x at 0 up to order 2 |
|
249 T result = x; |
|
250 |
|
251 if (abs(x) >= taylor_2_bound) |
|
252 { |
|
253 T x3 = x*x*x; |
|
254 |
|
255 // approximation by taylor series in x at 0 up to order 4 |
|
256 result += x3/static_cast<T>(3); |
|
257 } |
|
258 |
|
259 return(result); |
|
260 } |
|
261 } |
|
262 #endif /* defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) */ |
|
263 } |
|
264 } |
|
265 |
|
266 #endif /* BOOST_ATANH_HPP */ |
|
267 |