ceph/src/boost/libs/math/include/boost/math/complex/atanh.hpp

   1 //  (C) Copyright John Maddock 2005.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
   7 #define BOOST_MATH_COMPLEX_ATANH_INCLUDED
   8
   9 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
  10 #  include <boost/math/complex/details.hpp>
  11 #endif
  12 #ifndef BOOST_MATH_LOG1P_INCLUDED
  13 #  include <boost/math/special_functions/log1p.hpp>
  14 #endif
  15 #include <boost/assert.hpp>
  16
  17 #ifdef BOOST_NO_STDC_NAMESPACE
  18 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
  19 #endif
  20
  21 namespace boost{ namespace math{
  22
  23 template<class T>
  24 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.html
  32    //
  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.htm
  38    //
  39    // See also: https://svn.boost.org/trac/boost/ticket/7291
  40    //
  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 BOOST_MSVC
  51 #pragma warning(push)
  52 #pragma warning(disable:4127)
  53 #endif
  54
  55    T x = std::fabs(z.real());
  56    T y = std::fabs(z.imag());
  57
  58    T real, imag;  // our results
  59
  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       else
  73          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       else
  82          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    else
 105    {
 106       //
 107       // This section handles exception cases that would normally cause
 108       // underflow or overflow in the main formulas.
 109       //
 110       // Begin by working out the real part, we need to approximate
 111       //    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          else
 133          {
 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          else
 146          {
 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       else
 160          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          else
 185          {
 186             if((y == zero) && (x == zero))
 187                imag = 0;
 188             else
 189                imag = std::atan2(two*y, one - y*y);
 190          }
 191       }
 192       else
 193       {
 194          //
 195          // y^2 is negligible:
 196          //
 197          if((y == zero) && (x == one))
 198             imag = 0;
 199          else
 200             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 BOOST_MSVC
 208 #pragma warning(pop)
 209 #endif
 210 }
 211
 212 } } // namespaces
 213
 214 #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED