1 // boost atanh.hpp header file |
1 // (C) Copyright John Maddock 2005. |
2 |
2 // Use, modification and distribution are subject to the |
3 // (C) Copyright Hubert Holin 2001. |
3 // Boost Software License, Version 1.0. (See accompanying file |
4 // Distributed under the Boost Software License, Version 1.0. (See |
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
5 // accompanying file LICENSE_1_0.txt or copy at |
5 |
6 // http://www.boost.org/LICENSE_1_0.txt) |
6 #ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED |
7 |
7 #define BOOST_MATH_COMPLEX_ATANH_INCLUDED |
8 // See http://www.boost.org for updates, documentation, and revision history. |
8 |
9 |
9 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED |
10 #ifndef BOOST_ATANH_HPP |
10 # include <boost/math/complex/details.hpp> |
11 #define BOOST_ATANH_HPP |
11 #endif |
12 |
12 #ifndef BOOST_MATH_LOG1P_INCLUDED |
13 |
13 # include <boost/math/special_functions/log1p.hpp> |
14 #include <cmath> |
14 #endif |
15 #include <limits> |
15 #include <boost/assert.hpp> |
16 #include <string> |
16 |
17 #include <stdexcept> |
17 #ifdef BOOST_NO_STDC_NAMESPACE |
18 |
18 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; } |
19 |
19 #endif |
20 #include <boost/config.hpp> |
20 |
21 |
21 namespace boost{ namespace math{ |
22 |
22 |
23 // This is the inverse of the hyperbolic tangent function. |
23 template<class T> |
24 |
24 std::complex<T> atanh(const std::complex<T>& z) |
25 namespace boost |
|
26 { |
25 { |
27 namespace math |
26 // |
28 { |
27 // References: |
29 #if defined(__GNUC__) && (__GNUC__ < 3) |
28 // |
30 // gcc 2.x ignores function scope using declarations, |
29 // Eric W. Weisstein. "Inverse Hyperbolic Tangent." |
31 // put them in the scope of the enclosing namespace instead: |
30 // From MathWorld--A Wolfram Web Resource. |
32 |
31 // http://mathworld.wolfram.com/InverseHyperbolicTangent.html |
33 using ::std::abs; |
32 // |
34 using ::std::sqrt; |
33 // Also: The Wolfram Functions Site, |
35 using ::std::log; |
34 // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/ |
36 |
35 // |
37 using ::std::numeric_limits; |
36 // Also "Abramowitz and Stegun. Handbook of Mathematical Functions." |
38 #endif |
37 // at : http://jove.prohosting.com/~skripty/toc.htm |
39 |
38 // |
40 #if defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) |
39 |
41 // This is the main fare |
40 static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L); |
42 |
41 static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L); |
43 template<typename T> |
42 static const T one = static_cast<T>(1.0L); |
44 inline T atanh(const T x) |
43 static const T two = static_cast<T>(2.0L); |
45 { |
44 static const T four = static_cast<T>(4.0L); |
46 using ::std::abs; |
45 static const T zero = static_cast<T>(0); |
47 using ::std::sqrt; |
46 static const T a_crossover = static_cast<T>(0.3L); |
48 using ::std::log; |
47 |
49 |
48 T x = std::fabs(z.real()); |
50 using ::std::numeric_limits; |
49 T y = std::fabs(z.imag()); |
51 |
50 |
52 T const one = static_cast<T>(1); |
51 T real, imag; // our results |
53 T const two = static_cast<T>(2); |
52 |
54 |
53 T safe_upper = detail::safe_max(two); |
55 static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon()); |
54 T safe_lower = detail::safe_min(static_cast<T>(2)); |
56 static T const taylor_n_bound = sqrt(taylor_2_bound); |
55 |
57 |
56 // |
58 if (x < -one) |
57 // Begin by handling the special cases specified in C99: |
59 { |
58 // |
60 if (numeric_limits<T>::has_quiet_NaN) |
59 if(detail::test_is_nan(x)) |
61 { |
60 { |
62 return(numeric_limits<T>::quiet_NaN()); |
61 if(detail::test_is_nan(y)) |
63 } |
62 return std::complex<T>(x, x); |
64 else |
63 else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity())) |
65 { |
64 return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi)); |
66 ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!"); |
65 else |
67 ::std::domain_error bad_argument(error_reporting); |
66 return std::complex<T>(x, x); |
68 |
67 } |
69 throw(bad_argument); |
68 else if(detail::test_is_nan(y)) |
70 } |
69 { |
71 } |
70 if(x == 0) |
72 else if (x < -one+numeric_limits<T>::epsilon()) |
71 return std::complex<T>(x, y); |
73 { |
72 if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity())) |
74 if (numeric_limits<T>::has_infinity) |
73 return std::complex<T>(0, y); |
75 { |
74 else |
76 return(-numeric_limits<T>::infinity()); |
75 return std::complex<T>(y, y); |
77 } |
76 } |
78 else |
77 else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper)) |
79 { |
78 { |
80 ::std::string error_reporting("Argument to atanh is -1 (result: -Infinity)!"); |
79 |
81 ::std::out_of_range bad_argument(error_reporting); |
80 T xx = x*x; |
82 |
81 T yy = y*y; |
83 throw(bad_argument); |
82 T x2 = x * two; |
84 } |
83 |
85 } |
84 /// |
86 else if (x > +one-numeric_limits<T>::epsilon()) |
85 // The real part is given by: |
87 { |
86 // |
88 if (numeric_limits<T>::has_infinity) |
87 // real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x)) |
89 { |
88 // |
90 return(+numeric_limits<T>::infinity()); |
89 // However, when x is either large (x > 1/E) or very small |
91 } |
90 // (x < E) then this effectively simplifies |
92 else |
91 // to log(1), leading to wildly inaccurate results. |
93 { |
92 // By dividing the above (top and bottom) by (1 + x^2 + y^2) we get: |
94 ::std::string error_reporting("Argument to atanh is +1 (result: +Infinity)!"); |
93 // |
95 ::std::out_of_range bad_argument(error_reporting); |
94 // real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2)))) |
96 |
95 // |
97 throw(bad_argument); |
96 // which is much more sensitive to the value of x, when x is not near 1 |
98 } |
97 // (remember we can compute log(1+x) for small x very accurately). |
99 } |
98 // |
100 else if (x > +one) |
99 // The cross-over from one method to the other has to be determined |
101 { |
100 // experimentally, the value used below appears correct to within a |
102 if (numeric_limits<T>::has_quiet_NaN) |
101 // factor of 2 (and there are larger errors from other parts |
103 { |
102 // of the input domain anyway). |
104 return(numeric_limits<T>::quiet_NaN()); |
103 // |
105 } |
104 T alpha = two*x / (one + xx + yy); |
106 else |
105 if(alpha < a_crossover) |
107 { |
106 { |
108 ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!"); |
107 real = boost::math::log1p(alpha) - boost::math::log1p(-alpha); |
109 ::std::domain_error bad_argument(error_reporting); |
108 } |
110 |
109 else |
111 throw(bad_argument); |
110 { |
112 } |
111 T xm1 = x - one; |
113 } |
112 real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy); |
114 else if (abs(x) >= taylor_n_bound) |
113 } |
115 { |
114 real /= four; |
116 return(log( (one + x) / (one - x) ) / two); |
115 if(z.real() < 0) |
117 } |
116 real = -real; |
|
117 |
|
118 imag = std::atan2((y * two), (one - xx - yy)); |
|
119 imag /= two; |
|
120 if(z.imag() < 0) |
|
121 imag = -imag; |
|
122 } |
|
123 else |
|
124 { |
|
125 // |
|
126 // This section handles exception cases that would normally cause |
|
127 // underflow or overflow in the main formulas. |
|
128 // |
|
129 // Begin by working out the real part, we need to approximate |
|
130 // alpha = 2x / (1 + x^2 + y^2) |
|
131 // without either overflow or underflow in the squared terms. |
|
132 // |
|
133 T alpha = 0; |
|
134 if(x >= safe_upper) |
|
135 { |
|
136 // this is really a test for infinity, |
|
137 // but we may not have the necessary numeric_limits support: |
|
138 if((x > (std::numeric_limits<T>::max)()) || (y > (std::numeric_limits<T>::max)())) |
|
139 { |
|
140 alpha = 0; |
|
141 } |
|
142 else if(y >= safe_upper) |
|
143 { |
|
144 // Big x and y: divide alpha through by x*y: |
|
145 alpha = (two/y) / (x/y + y/x); |
|
146 } |
|
147 else if(y > one) |
|
148 { |
|
149 // Big x: divide through by x: |
|
150 alpha = two / (x + y*y/x); |
|
151 } |
|
152 else |
|
153 { |
|
154 // Big x small y, as above but neglect y^2/x: |
|
155 alpha = two/x; |
|
156 } |
|
157 } |
|
158 else if(y >= safe_upper) |
|
159 { |
|
160 if(x > one) |
|
161 { |
|
162 // Big y, medium x, divide through by y: |
|
163 alpha = (two*x/y) / (y + x*x/y); |
|
164 } |
|
165 else |
|
166 { |
|
167 // Small x and y, whatever alpha is, it's too small to calculate: |
|
168 alpha = 0; |
|
169 } |
|
170 } |
|
171 else |
|
172 { |
|
173 // one or both of x and y are small, calculate divisor carefully: |
|
174 T div = one; |
|
175 if(x > safe_lower) |
|
176 div += x*x; |
|
177 if(y > safe_lower) |
|
178 div += y*y; |
|
179 alpha = two*x/div; |
|
180 } |
|
181 if(alpha < a_crossover) |
|
182 { |
|
183 real = boost::math::log1p(alpha) - boost::math::log1p(-alpha); |
|
184 } |
|
185 else |
|
186 { |
|
187 // We can only get here as a result of small y and medium sized x, |
|
188 // we can simply neglect the y^2 terms: |
|
189 BOOST_ASSERT(x >= safe_lower); |
|
190 BOOST_ASSERT(x <= safe_upper); |
|
191 //BOOST_ASSERT(y <= safe_lower); |
|
192 T xm1 = x - one; |
|
193 real = std::log(1 + two*x + x*x) - std::log(xm1*xm1); |
|
194 } |
|
195 |
|
196 real /= four; |
|
197 if(z.real() < 0) |
|
198 real = -real; |
|
199 |
|
200 // |
|
201 // Now handle imaginary part, this is much easier, |
|
202 // if x or y are large, then the formula: |
|
203 // atan2(2y, 1 - x^2 - y^2) |
|
204 // evaluates to +-(PI - theta) where theta is negligible compared to PI. |
|
205 // |
|
206 if((x >= safe_upper) || (y >= safe_upper)) |
|
207 { |
|
208 imag = pi; |
|
209 } |
|
210 else if(x <= safe_lower) |
|
211 { |
|
212 // |
|
213 // If both x and y are small then atan(2y), |
|
214 // otherwise just x^2 is negligible in the divisor: |
|
215 // |
|
216 if(y <= safe_lower) |
|
217 imag = std::atan2(two*y, one); |
|
218 else |
|
219 { |
|
220 if((y == zero) && (x == zero)) |
|
221 imag = 0; |
118 else |
222 else |
119 { |
223 imag = std::atan2(two*y, one - y*y); |
120 // approximation by taylor series in x at 0 up to order 2 |
224 } |
121 T result = x; |
225 } |
122 |
226 else |
123 if (abs(x) >= taylor_2_bound) |
227 { |
124 { |
228 // |
125 T x3 = x*x*x; |
229 // y^2 is negligible: |
126 |
230 // |
127 // approximation by taylor series in x at 0 up to order 4 |
231 if((y == zero) && (x == one)) |
128 result += x3/static_cast<T>(3); |
232 imag = 0; |
129 } |
233 else |
130 |
234 imag = std::atan2(two*y, 1 - x*x); |
131 return(result); |
235 } |
132 } |
236 imag /= two; |
133 } |
237 if(z.imag() < 0) |
134 #else |
238 imag = -imag; |
135 // These are implementation details (for main fare see below) |
239 } |
136 |
240 return std::complex<T>(real, imag); |
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 } |
241 } |
265 |
242 |
266 #endif /* BOOST_ATANH_HPP */ |
243 } } // namespaces |
267 |
244 |
|
245 #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED |