epoc32/include/stdapis/boost/math/complex/atanh.hpp
author William Roberts <williamr@symbian.org>
Wed, 31 Mar 2010 12:33:34 +0100
branchSymbian3
changeset 4 837f303aceeb
parent 3 e1b950c65cb4
permissions -rw-r--r--
Current Symbian^3 public API header files (from PDK 3.0.h)
This is the epoc32/include tree with the "platform" subtrees removed, and
all but a selected few mbg and rsg files removed.
     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    
    40    static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
    41    static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
    42    static const T one = static_cast<T>(1.0L);
    43    static const T two = static_cast<T>(2.0L);
    44    static const T four = static_cast<T>(4.0L);
    45    static const T zero = static_cast<T>(0);
    46    static const T a_crossover = static_cast<T>(0.3L);
    47 
    48    T x = std::fabs(z.real());
    49    T y = std::fabs(z.imag());
    50 
    51    T real, imag;  // our results
    52 
    53    T safe_upper = detail::safe_max(two);
    54    T safe_lower = detail::safe_min(static_cast<T>(2));
    55 
    56    //
    57    // Begin by handling the special cases specified in C99:
    58    //
    59    if(detail::test_is_nan(x))
    60    {
    61       if(detail::test_is_nan(y))
    62          return std::complex<T>(x, x);
    63       else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
    64          return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi));
    65       else
    66          return std::complex<T>(x, x);
    67    }
    68    else if(detail::test_is_nan(y))
    69    {
    70       if(x == 0)
    71          return std::complex<T>(x, y);
    72       if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
    73          return std::complex<T>(0, y);
    74       else
    75          return std::complex<T>(y, y);
    76    }
    77    else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
    78    {
    79 
    80       T xx = x*x;
    81       T yy = y*y;
    82       T x2 = x * two;
    83 
    84       ///
    85       // The real part is given by:
    86       // 
    87       // real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x))
    88       // 
    89       // However, when x is either large (x > 1/E) or very small
    90       // (x < E) then this effectively simplifies
    91       // to log(1), leading to wildly inaccurate results.  
    92       // By dividing the above (top and bottom) by (1 + x^2 + y^2) we get:
    93       //
    94       // real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2))))
    95       //
    96       // which is much more sensitive to the value of x, when x is not near 1
    97       // (remember we can compute log(1+x) for small x very accurately).
    98       //
    99       // The cross-over from one method to the other has to be determined
   100       // experimentally, the value used below appears correct to within a 
   101       // factor of 2 (and there are larger errors from other parts
   102       // of the input domain anyway).
   103       //
   104       T alpha = two*x / (one + xx + yy);
   105       if(alpha < a_crossover)
   106       {
   107          real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
   108       }
   109       else
   110       {
   111          T xm1 = x - one;
   112          real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy);
   113       }
   114       real /= four;
   115       if(z.real() < 0)
   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;
   222             else
   223                imag = std::atan2(two*y, one - y*y);
   224          }
   225       }
   226       else
   227       {
   228          //
   229          // y^2 is negligible:
   230          //
   231          if((y == zero) && (x == one))
   232             imag = 0;
   233          else
   234             imag = std::atan2(two*y, 1 - x*x);
   235       }
   236       imag /= two;
   237       if(z.imag() < 0)
   238          imag = -imag;
   239    }
   240    return std::complex<T>(real, imag);
   241 }
   242 
   243 } } // namespaces
   244 
   245 #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED