diff -r e1b950c65cb4 -r 837f303aceeb epoc32/include/stdapis/boost/math/complex/atanh.hpp
--- a/epoc32/include/stdapis/boost/math/complex/atanh.hpp	Wed Mar 31 12:27:01 2010 +0100
+++ b/epoc32/include/stdapis/boost/math/complex/atanh.hpp	Wed Mar 31 12:33:34 2010 +0100
@@ -1,267 +1,245 @@
-//    boost atanh.hpp header file
+//  (C) Copyright John Maddock 2005.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
-//  (C) Copyright Hubert Holin 2001.
-//  Distributed under the Boost Software License, Version 1.0. (See
-//  accompanying file LICENSE_1_0.txt or copy at
-//  http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
+#define BOOST_MATH_COMPLEX_ATANH_INCLUDED
 
-// See http://www.boost.org for updates, documentation, and revision history.
+#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
+#  include <boost/math/complex/details.hpp>
+#endif
+#ifndef BOOST_MATH_LOG1P_INCLUDED
+#  include <boost/math/special_functions/log1p.hpp>
+#endif
+#include <boost/assert.hpp>
 
-#ifndef BOOST_ATANH_HPP
-#define BOOST_ATANH_HPP
+#ifdef BOOST_NO_STDC_NAMESPACE
+namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
+#endif
 
+namespace boost{ namespace math{
 
-#include <cmath>
-#include <limits>
-#include <string>
-#include <stdexcept>
+template<class T> 
+std::complex<T> atanh(const std::complex<T>& z)
+{
+   //
+   // References:
+   //
+   // Eric W. Weisstein. "Inverse Hyperbolic Tangent." 
+   // From MathWorld--A Wolfram Web Resource. 
+   // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
+   //
+   // Also: The Wolfram Functions Site,
+   // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
+   //
+   // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
+   // at : http://jove.prohosting.com/~skripty/toc.htm
+   //
+   
+   static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
+   static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
+   static const T one = static_cast<T>(1.0L);
+   static const T two = static_cast<T>(2.0L);
+   static const T four = static_cast<T>(4.0L);
+   static const T zero = static_cast<T>(0);
+   static const T a_crossover = static_cast<T>(0.3L);
 
+   T x = std::fabs(z.real());
+   T y = std::fabs(z.imag());
 
-#include <boost/config.hpp>
+   T real, imag;  // our results
 
+   T safe_upper = detail::safe_max(two);
+   T safe_lower = detail::safe_min(static_cast<T>(2));
 
-// This is the inverse of the hyperbolic tangent function.
+   //
+   // Begin by handling the special cases specified in C99:
+   //
+   if(detail::test_is_nan(x))
+   {
+      if(detail::test_is_nan(y))
+         return std::complex<T>(x, x);
+      else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
+         return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi));
+      else
+         return std::complex<T>(x, x);
+   }
+   else if(detail::test_is_nan(y))
+   {
+      if(x == 0)
+         return std::complex<T>(x, y);
+      if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
+         return std::complex<T>(0, y);
+      else
+         return std::complex<T>(y, y);
+   }
+   else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
+   {
 
-namespace boost
-{
-    namespace math
-    {
-#if defined(__GNUC__) && (__GNUC__ < 3)
-        // gcc 2.x ignores function scope using declarations,
-        // put them in the scope of the enclosing namespace instead:
-        
-        using    ::std::abs;
-        using    ::std::sqrt;
-        using    ::std::log;
-        
-        using    ::std::numeric_limits;
-#endif
-        
-#if defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION)
-        // This is the main fare
-        
-        template<typename T>
-        inline T    atanh(const T x)
-        {
-            using    ::std::abs;
-            using    ::std::sqrt;
-            using    ::std::log;
-            
-            using    ::std::numeric_limits;
-            
-            T const            one = static_cast<T>(1);
-            T const            two = static_cast<T>(2);
-            
-            static T const    taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
-            static T const    taylor_n_bound = sqrt(taylor_2_bound);
-            
-            if        (x < -one)
-            {
-                if    (numeric_limits<T>::has_quiet_NaN)
-                {
-                    return(numeric_limits<T>::quiet_NaN());
-                }
-                else
-                {
-                    ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
-                    ::std::domain_error  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-            }
-            else if    (x < -one+numeric_limits<T>::epsilon())
-            {
-                if    (numeric_limits<T>::has_infinity)
-                {
-                    return(-numeric_limits<T>::infinity());
-                }
-                else
-                {
-                    ::std::string        error_reporting("Argument to atanh is -1 (result: -Infinity)!");
-                    ::std::out_of_range  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-            }
-            else if    (x > +one-numeric_limits<T>::epsilon())
-            {
-                if    (numeric_limits<T>::has_infinity)
-                {
-                    return(+numeric_limits<T>::infinity());
-                }
-                else
-                {
-                    ::std::string        error_reporting("Argument to atanh is +1 (result: +Infinity)!");
-                    ::std::out_of_range  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-            }
-            else if    (x > +one)
-            {
-                if    (numeric_limits<T>::has_quiet_NaN)
-                {
-                    return(numeric_limits<T>::quiet_NaN());
-                }
-                else
-                {
-                    ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
-                    ::std::domain_error  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-            }
-            else if    (abs(x) >= taylor_n_bound)
-            {
-                return(log( (one + x) / (one - x) ) / two);
-            }
+      T xx = x*x;
+      T yy = y*y;
+      T x2 = x * two;
+
+      ///
+      // The real part is given by:
+      // 
+      // real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x))
+      // 
+      // However, when x is either large (x > 1/E) or very small
+      // (x < E) then this effectively simplifies
+      // to log(1), leading to wildly inaccurate results.  
+      // By dividing the above (top and bottom) by (1 + x^2 + y^2) we get:
+      //
+      // real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2))))
+      //
+      // which is much more sensitive to the value of x, when x is not near 1
+      // (remember we can compute log(1+x) for small x very accurately).
+      //
+      // The cross-over from one method to the other has to be determined
+      // experimentally, the value used below appears correct to within a 
+      // factor of 2 (and there are larger errors from other parts
+      // of the input domain anyway).
+      //
+      T alpha = two*x / (one + xx + yy);
+      if(alpha < a_crossover)
+      {
+         real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
+      }
+      else
+      {
+         T xm1 = x - one;
+         real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy);
+      }
+      real /= four;
+      if(z.real() < 0)
+         real = -real;
+
+      imag = std::atan2((y * two), (one - xx - yy));
+      imag /= two;
+      if(z.imag() < 0)
+         imag = -imag;
+   }
+   else
+   {
+      //
+      // This section handles exception cases that would normally cause
+      // underflow or overflow in the main formulas.
+      //
+      // Begin by working out the real part, we need to approximate
+      //    alpha = 2x / (1 + x^2 + y^2)
+      // without either overflow or underflow in the squared terms.
+      //
+      T alpha = 0;
+      if(x >= safe_upper)
+      {
+         // this is really a test for infinity, 
+         // but we may not have the necessary numeric_limits support:
+         if((x > (std::numeric_limits<T>::max)()) || (y > (std::numeric_limits<T>::max)()))
+         {
+            alpha = 0;
+         }
+         else if(y >= safe_upper)
+         {
+            // Big x and y: divide alpha through by x*y:
+            alpha = (two/y) / (x/y + y/x);
+         }
+         else if(y > one)
+         {
+            // Big x: divide through by x:
+            alpha = two / (x + y*y/x);
+         }
+         else
+         {
+            // Big x small y, as above but neglect y^2/x:
+            alpha = two/x;
+         }
+      }
+      else if(y >= safe_upper)
+      {
+         if(x > one)
+         {
+            // Big y, medium x, divide through by y:
+            alpha = (two*x/y) / (y + x*x/y);
+         }
+         else
+         {
+            // Small x and y, whatever alpha is, it's too small to calculate:
+            alpha = 0;
+         }
+      }
+      else
+      {
+         // one or both of x and y are small, calculate divisor carefully:
+         T div = one;
+         if(x > safe_lower)
+            div += x*x;
+         if(y > safe_lower)
+            div += y*y;
+         alpha = two*x/div;
+      }
+      if(alpha < a_crossover)
+      {
+         real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
+      }
+      else
+      {
+         // We can only get here as a result of small y and medium sized x,
+         // we can simply neglect the y^2 terms:
+         BOOST_ASSERT(x >= safe_lower);
+         BOOST_ASSERT(x <= safe_upper);
+         //BOOST_ASSERT(y <= safe_lower);
+         T xm1 = x - one;
+         real = std::log(1 + two*x + x*x) - std::log(xm1*xm1);
+      }
+      
+      real /= four;
+      if(z.real() < 0)
+         real = -real;
+
+      //
+      // Now handle imaginary part, this is much easier,
+      // if x or y are large, then the formula:
+      //    atan2(2y, 1 - x^2 - y^2)
+      // evaluates to +-(PI - theta) where theta is negligible compared to PI.
+      //
+      if((x >= safe_upper) || (y >= safe_upper))
+      {
+         imag = pi;
+      }
+      else if(x <= safe_lower)
+      {
+         //
+         // If both x and y are small then atan(2y),
+         // otherwise just x^2 is negligible in the divisor:
+         //
+         if(y <= safe_lower)
+            imag = std::atan2(two*y, one);
+         else
+         {
+            if((y == zero) && (x == zero))
+               imag = 0;
             else
-            {
-                // approximation by taylor series in x at 0 up to order 2
-                T    result = x;
-                
-                if    (abs(x) >= taylor_2_bound)
-                {
-                    T    x3 = x*x*x;
-                    
-                    // approximation by taylor series in x at 0 up to order 4
-                    result += x3/static_cast<T>(3);
-                }
-                
-                return(result);
-            }
-        }
-#else
-        // These are implementation details (for main fare see below)
-        
-        namespace detail
-        {
-            template    <
-                            typename T,
-                            bool InfinitySupported
-                        >
-            struct    atanh_helper1_t
-            {
-                static T    get_pos_infinity()
-                {
-                    return(+::std::numeric_limits<T>::infinity());
-                }
-                
-                static T    get_neg_infinity()
-                {
-                    return(-::std::numeric_limits<T>::infinity());
-                }
-            };    // boost::math::detail::atanh_helper1_t
-            
-            
-            template<typename T>
-            struct    atanh_helper1_t<T, false>
-            {
-                static T    get_pos_infinity()
-                {
-                    ::std::string        error_reporting("Argument to atanh is +1 (result: +Infinity)!");
-                    ::std::out_of_range  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-                
-                static T    get_neg_infinity()
-                {
-                    ::std::string        error_reporting("Argument to atanh is -1 (result: -Infinity)!");
-                    ::std::out_of_range  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-            };    // boost::math::detail::atanh_helper1_t
-            
-            
-            template    <
-                            typename T,
-                            bool QuietNanSupported
-                        >
-            struct    atanh_helper2_t
-            {
-                static T    get_NaN()
-                {
-                    return(::std::numeric_limits<T>::quiet_NaN());
-                }
-            };    // boost::detail::atanh_helper2_t
-            
-            
-            template<typename T>
-            struct    atanh_helper2_t<T, false>
-            {
-                static T    get_NaN()
-                {
-                    ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
-                    ::std::domain_error  bad_argument(error_reporting);
-                    
-                    throw(bad_argument);
-                }
-            };    // boost::detail::atanh_helper2_t
-        }    // boost::detail
-        
-        
-        // This is the main fare
-        
-        template<typename T>
-        inline T    atanh(const T x)
-        {
-            using    ::std::abs;
-            using    ::std::sqrt;
-            using    ::std::log;
-            
-            using    ::std::numeric_limits;
-            
-            typedef  detail::atanh_helper1_t<T, ::std::numeric_limits<T>::has_infinity>    helper1_type;
-            typedef  detail::atanh_helper2_t<T, ::std::numeric_limits<T>::has_quiet_NaN>    helper2_type;
-            
-            
-            T const           one = static_cast<T>(1);
-            T const           two = static_cast<T>(2);
-            
-            static T const    taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
-            static T const    taylor_n_bound = sqrt(taylor_2_bound);
-            
-            if        (x < -one)
-            {
-                return(helper2_type::get_NaN());
-            }
-            else if    (x < -one+numeric_limits<T>::epsilon())
-            {
-                return(helper1_type::get_neg_infinity());
-            }
-            else if    (x > +one-numeric_limits<T>::epsilon())
-            {
-                return(helper1_type::get_pos_infinity());
-            }
-            else if    (x > +one)
-            {
-                return(helper2_type::get_NaN());
-            }
-            else if    (abs(x) >= taylor_n_bound)
-            {
-                return(log( (one + x) / (one - x) ) / two);
-            }
-            else
-            {
-                // approximation by taylor series in x at 0 up to order 2
-                T    result = x;
-                
-                if    (abs(x) >= taylor_2_bound)
-                {
-                    T    x3 = x*x*x;
-                    
-                    // approximation by taylor series in x at 0 up to order 4
-                    result += x3/static_cast<T>(3);
-                }
-                
-                return(result);
-            }
-        }
-#endif /* defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) */
-    }
+               imag = std::atan2(two*y, one - y*y);
+         }
+      }
+      else
+      {
+         //
+         // y^2 is negligible:
+         //
+         if((y == zero) && (x == one))
+            imag = 0;
+         else
+            imag = std::atan2(two*y, 1 - x*x);
+      }
+      imag /= two;
+      if(z.imag() < 0)
+         imag = -imag;
+   }
+   return std::complex<T>(real, imag);
 }
 
-#endif /* BOOST_ATANH_HPP */
+} } // namespaces
 
+#endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED