epoc32/include/stdapis/boost/math/complex/asin.hpp
branchSymbian2
changeset 2 2fe1408b6811
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/epoc32/include/stdapis/boost/math/complex/asin.hpp	Tue Mar 16 16:12:26 2010 +0000
     1.3 @@ -0,0 +1,245 @@
     1.4 +//  (C) Copyright John Maddock 2005.
     1.5 +//  Distributed under the Boost Software License, Version 1.0. (See accompanying
     1.6 +//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
     1.7 +
     1.8 +#ifndef BOOST_MATH_COMPLEX_ASIN_INCLUDED
     1.9 +#define BOOST_MATH_COMPLEX_ASIN_INCLUDED
    1.10 +
    1.11 +#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
    1.12 +#  include <boost/math/complex/details.hpp>
    1.13 +#endif
    1.14 +#ifndef BOOST_MATH_LOG1P_INCLUDED
    1.15 +#  include <boost/math/special_functions/log1p.hpp>
    1.16 +#endif
    1.17 +#include <boost/assert.hpp>
    1.18 +
    1.19 +#ifdef BOOST_NO_STDC_NAMESPACE
    1.20 +namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
    1.21 +#endif
    1.22 +
    1.23 +namespace boost{ namespace math{
    1.24 +
    1.25 +template<class T> 
    1.26 +inline std::complex<T> asin(const std::complex<T>& z)
    1.27 +{
    1.28 +   //
    1.29 +   // This implementation is a transcription of the pseudo-code in:
    1.30 +   //
    1.31 +   // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling."
    1.32 +   // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
    1.33 +   // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
    1.34 +   //
    1.35 +
    1.36 +   //
    1.37 +   // These static constants should really be in a maths constants library:
    1.38 +   //
    1.39 +   static const T one = static_cast<T>(1);
    1.40 +   //static const T two = static_cast<T>(2);
    1.41 +   static const T half = static_cast<T>(0.5L);
    1.42 +   static const T a_crossover = static_cast<T>(1.5L);
    1.43 +   static const T b_crossover = static_cast<T>(0.6417L);
    1.44 +   //static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
    1.45 +   static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
    1.46 +   static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);
    1.47 +   static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);
    1.48 +   
    1.49 +   //
    1.50 +   // Get real and imaginary parts, discard the signs as we can 
    1.51 +   // figure out the sign of the result later:
    1.52 +   //
    1.53 +   T x = std::fabs(z.real());
    1.54 +   T y = std::fabs(z.imag());
    1.55 +   T real, imag;  // our results
    1.56 +
    1.57 +   //
    1.58 +   // Begin by handling the special cases for infinities and nan's
    1.59 +   // specified in C99, most of this is handled by the regular logic
    1.60 +   // below, but handling it as a special case prevents overflow/underflow
    1.61 +   // arithmetic which may trip up some machines:
    1.62 +   //
    1.63 +   if(detail::test_is_nan(x))
    1.64 +   {
    1.65 +      if(detail::test_is_nan(y))
    1.66 +         return std::complex<T>(x, x);
    1.67 +      if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
    1.68 +      {
    1.69 +         real = x;
    1.70 +         imag = std::numeric_limits<T>::infinity();
    1.71 +      }
    1.72 +      else
    1.73 +         return std::complex<T>(x, x);
    1.74 +   }
    1.75 +   else if(detail::test_is_nan(y))
    1.76 +   {
    1.77 +      if(x == 0)
    1.78 +      {
    1.79 +         real = 0;
    1.80 +         imag = y;
    1.81 +      }
    1.82 +      else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
    1.83 +      {
    1.84 +         real = y;
    1.85 +         imag = std::numeric_limits<T>::infinity();
    1.86 +      }
    1.87 +      else
    1.88 +         return std::complex<T>(y, y);
    1.89 +   }
    1.90 +   else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
    1.91 +   {
    1.92 +      if(y == std::numeric_limits<T>::infinity())
    1.93 +      {
    1.94 +         real = quarter_pi;
    1.95 +         imag = std::numeric_limits<T>::infinity();
    1.96 +      }
    1.97 +      else
    1.98 +      {
    1.99 +         real = half_pi;
   1.100 +         imag = std::numeric_limits<T>::infinity();
   1.101 +      }
   1.102 +   }
   1.103 +   else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
   1.104 +   {
   1.105 +      real = 0;
   1.106 +      imag = std::numeric_limits<T>::infinity();
   1.107 +   }
   1.108 +   else
   1.109 +   {
   1.110 +      //
   1.111 +      // special case for real numbers:
   1.112 +      //
   1.113 +      if((y == 0) && (x <= one))
   1.114 +         return std::complex<T>(std::asin(z.real()));
   1.115 +      //
   1.116 +      // Figure out if our input is within the "safe area" identified by Hull et al.
   1.117 +      // This would be more efficient with portable floating point exception handling;
   1.118 +      // fortunately the quantities M and u identified by Hull et al (figure 3), 
   1.119 +      // match with the max and min methods of numeric_limits<T>.
   1.120 +      //
   1.121 +      T safe_max = detail::safe_max(static_cast<T>(8));
   1.122 +      T safe_min = detail::safe_min(static_cast<T>(4));
   1.123 +
   1.124 +      T xp1 = one + x;
   1.125 +      T xm1 = x - one;
   1.126 +
   1.127 +      if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
   1.128 +      {
   1.129 +         T yy = y * y;
   1.130 +         T r = std::sqrt(xp1*xp1 + yy);
   1.131 +         T s = std::sqrt(xm1*xm1 + yy);
   1.132 +         T a = half * (r + s);
   1.133 +         T b = x / a;
   1.134 +
   1.135 +         if(b <= b_crossover)
   1.136 +         {
   1.137 +            real = std::asin(b);
   1.138 +         }
   1.139 +         else
   1.140 +         {
   1.141 +            T apx = a + x;
   1.142 +            if(x <= one)
   1.143 +            {
   1.144 +               real = std::atan(x/std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1))));
   1.145 +            }
   1.146 +            else
   1.147 +            {
   1.148 +               real = std::atan(x/(y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1)))));
   1.149 +            }
   1.150 +         }
   1.151 +
   1.152 +         if(a <= a_crossover)
   1.153 +         {
   1.154 +            T am1;
   1.155 +            if(x < one)
   1.156 +            {
   1.157 +               am1 = half * (yy/(r + xp1) + yy/(s - xm1));
   1.158 +            }
   1.159 +            else
   1.160 +            {
   1.161 +               am1 = half * (yy/(r + xp1) + (s + xm1));
   1.162 +            }
   1.163 +            imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
   1.164 +         }
   1.165 +         else
   1.166 +         {
   1.167 +            imag = std::log(a + std::sqrt(a*a - one));
   1.168 +         }
   1.169 +      }
   1.170 +      else
   1.171 +      {
   1.172 +         //
   1.173 +         // This is the Hull et al exception handling code from Fig 3 of their paper:
   1.174 +         //
   1.175 +         if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
   1.176 +         {
   1.177 +            if(x < one)
   1.178 +            {
   1.179 +               real = std::asin(x);
   1.180 +               imag = y / std::sqrt(xp1*xm1);
   1.181 +            }
   1.182 +            else
   1.183 +            {
   1.184 +               real = half_pi;
   1.185 +               if(((std::numeric_limits<T>::max)() / xp1) > xm1)
   1.186 +               {
   1.187 +                  // xp1 * xm1 won't overflow:
   1.188 +                  imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
   1.189 +               }
   1.190 +               else
   1.191 +               {
   1.192 +                  imag = log_two + std::log(x);
   1.193 +               }
   1.194 +            }
   1.195 +         }
   1.196 +         else if(y <= safe_min)
   1.197 +         {
   1.198 +            // There is an assumption in Hull et al's analysis that
   1.199 +            // if we get here then x == 1.  This is true for all "good"
   1.200 +            // machines where :
   1.201 +            // 
   1.202 +            // E^2 > 8*sqrt(u); with:
   1.203 +            //
   1.204 +            // E =  std::numeric_limits<T>::epsilon()
   1.205 +            // u = (std::numeric_limits<T>::min)()
   1.206 +            //
   1.207 +            // Hull et al provide alternative code for "bad" machines
   1.208 +            // but we have no way to test that here, so for now just assert
   1.209 +            // on the assumption:
   1.210 +            //
   1.211 +            BOOST_ASSERT(x == 1);
   1.212 +            real = half_pi - std::sqrt(y);
   1.213 +            imag = std::sqrt(y);
   1.214 +         }
   1.215 +         else if(std::numeric_limits<T>::epsilon() * y - one >= x)
   1.216 +         {
   1.217 +            real = x/y; // This can underflow!
   1.218 +            imag = log_two + std::log(y);
   1.219 +         }
   1.220 +         else if(x > one)
   1.221 +         {
   1.222 +            real = std::atan(x/y);
   1.223 +            T xoy = x/y;
   1.224 +            imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
   1.225 +         }
   1.226 +         else
   1.227 +         {
   1.228 +            T a = std::sqrt(one + y*y);
   1.229 +            real = x/a; // This can underflow!
   1.230 +            imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
   1.231 +         }
   1.232 +      }
   1.233 +   }
   1.234 +
   1.235 +   //
   1.236 +   // Finish off by working out the sign of the result:
   1.237 +   //
   1.238 +   if(z.real() < 0)
   1.239 +      real = -real;
   1.240 +   if(z.imag() < 0)
   1.241 +      imag = -imag;
   1.242 +
   1.243 +   return std::complex<T>(real, imag);
   1.244 +}
   1.245 +
   1.246 +} } // namespaces
   1.247 +
   1.248 +#endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED