215 lines · plain
1// (C) Copyright John Maddock 2005.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED7#define BOOST_MATH_COMPLEX_ATANH_INCLUDED8 9#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED10# include <boost/math/complex/details.hpp>11#endif12#ifndef BOOST_MATH_LOG1P_INCLUDED13# include <boost/math/special_functions/log1p.hpp>14#endif15#include <boost/math/tools/assert.hpp>16 17#ifdef BOOST_NO_STDC_NAMESPACE18namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }19#endif20 21namespace boost{ namespace math{22 23template<class T> 24[[deprecated("Replaced by C++11")]] std::complex<T> atanh(const std::complex<T>& z)25{26 //27 // References:28 //29 // Eric W. Weisstein. "Inverse Hyperbolic Tangent." 30 // From MathWorld--A Wolfram Web Resource. 31 // http://mathworld.wolfram.com/InverseHyperbolicTangent.html32 //33 // Also: The Wolfram Functions Site,34 // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/35 //36 // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."37 // at : http://jove.prohosting.com/~skripty/toc.htm38 //39 // See also: https://svn.boost.org/trac/boost/ticket/729140 //41 42 static const T pi = boost::math::constants::pi<T>();43 static const T half_pi = pi / 2;44 static const T one = static_cast<T>(1.0L);45 static const T two = static_cast<T>(2.0L);46 static const T four = static_cast<T>(4.0L);47 static const T zero = static_cast<T>(0);48 static const T log_two = boost::math::constants::ln_two<T>();49 50#ifdef _MSC_VER51#pragma warning(push)52#pragma warning(disable:4127)53#endif54 55 T x = std::fabs(z.real());56 T y = std::fabs(z.imag());57 58 T real, imag; // our results59 60 T safe_upper = detail::safe_max(two);61 T safe_lower = detail::safe_min(static_cast<T>(2));62 63 //64 // Begin by handling the special cases specified in C99:65 //66 if((boost::math::isnan)(x))67 {68 if((boost::math::isnan)(y))69 return std::complex<T>(x, x);70 else if((boost::math::isinf)(y))71 return std::complex<T>(0, ((boost::math::signbit)(z.imag()) ? -half_pi : half_pi));72 else73 return std::complex<T>(x, x);74 }75 else if((boost::math::isnan)(y))76 {77 if(x == 0)78 return std::complex<T>(x, y);79 if((boost::math::isinf)(x))80 return std::complex<T>(0, y);81 else82 return std::complex<T>(y, y);83 }84 else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))85 {86 87 T yy = y*y;88 T mxm1 = one - x;89 ///90 // The real part is given by:91 // 92 // real(atanh(z)) == log1p(4*x / ((x-1)*(x-1) + y^2))93 // 94 real = boost::math::log1p(four * x / (mxm1*mxm1 + yy));95 real /= four;96 if((boost::math::signbit)(z.real()))97 real = (boost::math::changesign)(real);98 99 imag = std::atan2((y * two), (mxm1*(one+x) - yy));100 imag /= two;101 if(z.imag() < 0)102 imag = (boost::math::changesign)(imag);103 }104 else105 {106 //107 // This section handles exception cases that would normally cause108 // underflow or overflow in the main formulas.109 //110 // Begin by working out the real part, we need to approximate111 // real = boost::math::log1p(4x / ((x-1)^2 + y^2))112 // without either overflow or underflow in the squared terms.113 //114 T mxm1 = one - x;115 if(x >= safe_upper)116 {117 // x-1 = x to machine precision:118 if((boost::math::isinf)(x) || (boost::math::isinf)(y))119 {120 real = 0;121 }122 else if(y >= safe_upper)123 {124 // Big x and y: divide through by x*y:125 real = boost::math::log1p((four/y) / (x/y + y/x));126 }127 else if(y > one)128 {129 // Big x: divide through by x:130 real = boost::math::log1p(four / (x + y*y/x));131 }132 else133 {134 // Big x small y, as above but neglect y^2/x:135 real = boost::math::log1p(four/x);136 }137 }138 else if(y >= safe_upper)139 {140 if(x > one)141 {142 // Big y, medium x, divide through by y:143 real = boost::math::log1p((four*x/y) / (y + mxm1*mxm1/y));144 }145 else146 {147 // Small or medium x, large y:148 real = four*x/y/y;149 }150 }151 else if (x != one)152 {153 // y is small, calculate divisor carefully:154 T div = mxm1*mxm1;155 if(y > safe_lower)156 div += y*y;157 real = boost::math::log1p(four*x/div);158 }159 else160 real = boost::math::changesign(two * (std::log(y) - log_two));161 162 real /= four;163 if((boost::math::signbit)(z.real()))164 real = (boost::math::changesign)(real);165 166 //167 // Now handle imaginary part, this is much easier,168 // if x or y are large, then the formula:169 // atan2(2y, (1-x)*(1+x) - y^2)170 // evaluates to +-(PI - theta) where theta is negligible compared to PI.171 //172 if((x >= safe_upper) || (y >= safe_upper))173 {174 imag = pi;175 }176 else if(x <= safe_lower)177 {178 //179 // If both x and y are small then atan(2y),180 // otherwise just x^2 is negligible in the divisor:181 //182 if(y <= safe_lower)183 imag = std::atan2(two*y, one);184 else185 {186 if((y == zero) && (x == zero))187 imag = 0;188 else189 imag = std::atan2(two*y, one - y*y);190 }191 }192 else193 {194 //195 // y^2 is negligible:196 //197 if((y == zero) && (x == one))198 imag = 0;199 else200 imag = std::atan2(two*y, mxm1*(one+x));201 }202 imag /= two;203 if((boost::math::signbit)(z.imag()))204 imag = (boost::math::changesign)(imag);205 }206 return std::complex<T>(real, imag);207#ifdef _MSC_VER208#pragma warning(pop)209#endif210}211 212} } // namespaces213 214#endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED215