boost/math/special_functions/cbrt.hpp

   1 //  (C) Copyright John Maddock 2006.
   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_SF_CBRT_HPP
   7 #define BOOST_MATH_SF_CBRT_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/rational.hpp>
  14 #include <boost/math/policies/error_handling.hpp>
  15 #include <boost/math/special_functions/math_fwd.hpp>
  16 #include <boost/math/special_functions/fpclassify.hpp>
  17 #include <boost/mpl/divides.hpp>
  18 #include <boost/mpl/plus.hpp>
  19 #include <boost/mpl/if.hpp>
  20 #include <boost/type_traits/is_convertible.hpp>
  21
  22 namespace boost{ namespace math{
  23
  24 namespace detail
  25 {
  26
  27 struct big_int_type
  28 {
  29    operator boost::uintmax_t()const;
  30 };
  31
  32 template <class T>
  33 struct largest_cbrt_int_type
  34 {
  35    typedef typename mpl::if_<
  36       boost::is_convertible<big_int_type, T>,
  37       boost::uintmax_t,
  38       unsigned int
  39    >::type type;
  40 };
  41
  42 template <class T, class Policy>
  43 T cbrt_imp(T z, const Policy& pol)
  44 {
  45    BOOST_MATH_STD_USING
  46    //
  47    // cbrt approximation for z in the range [0.5,1]
  48    // It's hard to say what number of terms gives the optimum
  49    // trade off between precision and performance, this seems
  50    // to be about the best for double precision.
  51    //
  52    // Maximum Deviation Found:                     1.231e-006
  53    // Expected Error Term:                         -1.231e-006
  54    // Maximum Relative Change in Control Points:   5.982e-004
  55    //
  56    static const T P[] = {
  57       static_cast<T>(0.37568269008611818),
  58       static_cast<T>(1.3304968705558024),
  59       static_cast<T>(-1.4897101632445036),
  60       static_cast<T>(1.2875573098219835),
  61       static_cast<T>(-0.6398703759826468),
  62       static_cast<T>(0.13584489959258635),
  63    };
  64    static const T correction[] = {
  65       static_cast<T>(0.62996052494743658238360530363911),  // 2^-2/3
  66       static_cast<T>(0.79370052598409973737585281963615),  // 2^-1/3
  67       static_cast<T>(1),
  68       static_cast<T>(1.2599210498948731647672106072782),   // 2^1/3
  69       static_cast<T>(1.5874010519681994747517056392723),   // 2^2/3
  70    };
  71    if((boost::math::isinf)(z) || (z == 0))
  72       return z;
  73    if(!(boost::math::isfinite)(z))
  74    {
  75       return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
  76    }
  77
  78    int i_exp, sign(1);
  79    if(z < 0)
  80    {
  81       z = -z;
  82       sign = -sign;
  83    }
  84
  85    T guess = frexp(z, &i_exp);
  86    int original_i_exp = i_exp; // save for later
  87    guess = tools::evaluate_polynomial(P, guess);
  88    int i_exp3 = i_exp / 3;
  89
  90    typedef typename largest_cbrt_int_type<T>::type shift_type;
  91
  92    BOOST_STATIC_ASSERT( ::std::numeric_limits<shift_type>::radix == 2);
  93
  94    if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
  95    {
  96       if(i_exp3 > 0)
  97          guess *= shift_type(1u) << i_exp3;
  98       else
  99          guess /= shift_type(1u) << -i_exp3;
 100    }
 101    else
 102    {
 103       guess = ldexp(guess, i_exp3);
 104    }
 105    i_exp %= 3;
 106    guess *= correction[i_exp + 2];
 107    //
 108    // Now inline Halley iteration.
 109    // We do this here rather than calling tools::halley_iterate since we can
 110    // simplify the expressions algebraically, and don't need most of the error
 111    // checking of the boilerplate version as we know in advance that the function
 112    // is well behaved...
 113    //
 114    typedef typename policies::precision<T, Policy>::type prec;
 115    typedef typename mpl::divides<prec, mpl::int_<3> >::type prec3;
 116    typedef typename mpl::plus<prec3, mpl::int_<3> >::type new_prec;
 117    typedef typename policies::normalise<Policy, policies::digits2<new_prec::value> >::type new_policy;
 118    //
 119    // Epsilon calculation uses compile time arithmetic when it's available for type T,
 120    // otherwise uses ldexp to calculate at runtime:
 121    //
 122    T eps = (new_prec::value > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
 123    T diff;
 124
 125    if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
 126    {
 127       //
 128       // Safe from overflow, use the fast method:
 129       //
 130       do
 131       {
 132          T g3 = guess * guess * guess;
 133          diff = (g3 + z + z) / (g3 + g3 + z);
 134          guess *= diff;
 135       }
 136       while(fabs(1 - diff) > eps);
 137    }
 138    else
 139    {
 140       //
 141       // Either we're ready to overflow, or we can't tell because numeric_limits isn't
 142       // available for type T:
 143       //
 144       do
 145       {
 146          T g2 = guess * guess;
 147          diff = (g2 - z / guess) / (2 * guess + z / g2);
 148          guess -= diff;
 149       }
 150       while((guess * eps) < fabs(diff));
 151    }
 152
 153    return sign * guess;
 154 }
 155
 156 } // namespace detail
 157
 158 template <class T, class Policy>
 159 inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)
 160 {
 161    typedef typename tools::promote_args<T>::type result_type;
 162    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 163    return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));
 164 }
 165
 166 template <class T>
 167 inline typename tools::promote_args<T>::type cbrt(T z)
 168 {
 169    return cbrt(z, policies::policy<>());
 170 }
 171
 172 } // namespace math
 173 } // namespace boost
 174
 175 #endif // BOOST_MATH_SF_CBRT_HPP
 176
 177
 178
 179