math/complex/atanh.hpp

0001 //  (C) Copyright John Maddock 2005.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
0007 #define BOOST_MATH_COMPLEX_ATANH_INCLUDED
0008
0009 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
0010 #  include <boost/math/complex/details.hpp>
0011 #endif
0012 #ifndef BOOST_MATH_LOG1P_INCLUDED
0013 #  include <boost/math/special_functions/log1p.hpp>
0014 #endif
0015 #include <boost/math/tools/assert.hpp>
0016
0017 #ifdef BOOST_NO_STDC_NAMESPACE
0018 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
0019 #endif
0020
0021 namespace boost{ namespace math{
0022
0023 template<class T>
0024 [[deprecated("Replaced by C++11")]] std::complex<T> atanh(const std::complex<T>& z)
0025 {
0026    //
0027    // References:
0028    //
0029    // Eric W. Weisstein. "Inverse Hyperbolic Tangent."
0030    // From MathWorld--A Wolfram Web Resource.
0031    // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
0032    //
0033    // Also: The Wolfram Functions Site,
0034    // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
0035    //
0036    // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
0037    // at : http://jove.prohosting.com/~skripty/toc.htm
0038    //
0039    // See also: https://svn.boost.org/trac/boost/ticket/7291
0040    //
0041
0042    static const T pi = boost::math::constants::pi<T>();
0043    static const T half_pi = pi / 2;
0044    static const T one = static_cast<T>(1.0L);
0045    static const T two = static_cast<T>(2.0L);
0046    static const T four = static_cast<T>(4.0L);
0047    static const T zero = static_cast<T>(0);
0048    static const T log_two = boost::math::constants::ln_two<T>();
0049
0050 #ifdef _MSC_VER
0051 #pragma warning(push)
0052 #pragma warning(disable:4127)
0053 #endif
0054
0055    T x = std::fabs(z.real());
0056    T y = std::fabs(z.imag());
0057
0058    T real, imag;  // our results
0059
0060    T safe_upper = detail::safe_max(two);
0061    T safe_lower = detail::safe_min(static_cast<T>(2));
0062
0063    //
0064    // Begin by handling the special cases specified in C99:
0065    //
0066    if((boost::math::isnan)(x))
0067    {
0068       if((boost::math::isnan)(y))
0069          return std::complex<T>(x, x);
0070       else if((boost::math::isinf)(y))
0071          return std::complex<T>(0, ((boost::math::signbit)(z.imag()) ? -half_pi : half_pi));
0072       else
0073          return std::complex<T>(x, x);
0074    }
0075    else if((boost::math::isnan)(y))
0076    {
0077       if(x == 0)
0078          return std::complex<T>(x, y);
0079       if((boost::math::isinf)(x))
0080          return std::complex<T>(0, y);
0081       else
0082          return std::complex<T>(y, y);
0083    }
0084    else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
0085    {
0086
0087       T yy = y*y;
0088       T mxm1 = one - x;
0089       ///
0090       // The real part is given by:
0091       //
0092       // real(atanh(z)) == log1p(4*x / ((x-1)*(x-1) + y^2))
0093       //
0094       real = boost::math::log1p(four * x / (mxm1*mxm1 + yy));
0095       real /= four;
0096       if((boost::math::signbit)(z.real()))
0097          real = (boost::math::changesign)(real);
0098
0099       imag = std::atan2((y * two), (mxm1*(one+x) - yy));
0100       imag /= two;
0101       if(z.imag() < 0)
0102          imag = (boost::math::changesign)(imag);
0103    }
0104    else
0105    {
0106       //
0107       // This section handles exception cases that would normally cause
0108       // underflow or overflow in the main formulas.
0109       //
0110       // Begin by working out the real part, we need to approximate
0111       //    real = boost::math::log1p(4x / ((x-1)^2 + y^2))
0112       // without either overflow or underflow in the squared terms.
0113       //
0114       T mxm1 = one - x;
0115       if(x >= safe_upper)
0116       {
0117          // x-1 = x to machine precision:
0118          if((boost::math::isinf)(x) || (boost::math::isinf)(y))
0119          {
0120             real = 0;
0121          }
0122          else if(y >= safe_upper)
0123          {
0124             // Big x and y: divide through by x*y:
0125             real = boost::math::log1p((four/y) / (x/y + y/x));
0126          }
0127          else if(y > one)
0128          {
0129             // Big x: divide through by x:
0130             real = boost::math::log1p(four / (x + y*y/x));
0131          }
0132          else
0133          {
0134             // Big x small y, as above but neglect y^2/x:
0135             real = boost::math::log1p(four/x);
0136          }
0137       }
0138       else if(y >= safe_upper)
0139       {
0140          if(x > one)
0141          {
0142             // Big y, medium x, divide through by y:
0143             real = boost::math::log1p((four*x/y) / (y + mxm1*mxm1/y));
0144          }
0145          else
0146          {
0147             // Small or medium x, large y:
0148             real = four*x/y/y;
0149          }
0150       }
0151       else if (x != one)
0152       {
0153          // y is small, calculate divisor carefully:
0154          T div = mxm1*mxm1;
0155          if(y > safe_lower)
0156             div += y*y;
0157          real = boost::math::log1p(four*x/div);
0158       }
0159       else
0160          real = boost::math::changesign(two * (std::log(y) - log_two));
0161
0162       real /= four;
0163       if((boost::math::signbit)(z.real()))
0164          real = (boost::math::changesign)(real);
0165
0166       //
0167       // Now handle imaginary part, this is much easier,
0168       // if x or y are large, then the formula:
0169       //    atan2(2y, (1-x)*(1+x) - y^2)
0170       // evaluates to +-(PI - theta) where theta is negligible compared to PI.
0171       //
0172       if((x >= safe_upper) || (y >= safe_upper))
0173       {
0174          imag = pi;
0175       }
0176       else if(x <= safe_lower)
0177       {
0178          //
0179          // If both x and y are small then atan(2y),
0180          // otherwise just x^2 is negligible in the divisor:
0181          //
0182          if(y <= safe_lower)
0183             imag = std::atan2(two*y, one);
0184          else
0185          {
0186             if((y == zero) && (x == zero))
0187                imag = 0;
0188             else
0189                imag = std::atan2(two*y, one - y*y);
0190          }
0191       }
0192       else
0193       {
0194          //
0195          // y^2 is negligible:
0196          //
0197          if((y == zero) && (x == one))
0198             imag = 0;
0199          else
0200             imag = std::atan2(two*y, mxm1*(one+x));
0201       }
0202       imag /= two;
0203       if((boost::math::signbit)(z.imag()))
0204          imag = (boost::math::changesign)(imag);
0205    }
0206    return std::complex<T>(real, imag);
0207 #ifdef _MSC_VER
0208 #pragma warning(pop)
0209 #endif
0210 }
0211
0212 } } // namespaces
0213
0214 #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED