Change cbrt implementation to use a better performing algorithm.

Fix a few warnings along the way.

[SVN r59095]

include/boost/math/distributions/detail/hypergeometric_pdf.hpp

@@ -327,7 +327,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
 template <class T, class Policy>
 inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
 {
-   hypergeometric_pdf_prime_loop_result_entry<T> result = { 1 };
+   hypergeometric_pdf_prime_loop_result_entry<T> result = { 1, 0 };
    hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) };
    return hypergeometric_pdf_prime_loop_imp<T>(data, result);
 }

include/boost/math/special_functions/cbrt.hpp

@@ -10,31 +10,65 @@
 #pragma once
 #endif
 
-#include <boost/math/tools/roots.hpp>
+#include <boost/math/tools/rational.hpp>
+#include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/mpl/divides.hpp>
+#include <boost/mpl/plus.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/type_traits/is_convertible.hpp>
 
 namespace boost{ namespace math{
 
 namespace detail
 {
 
-   template <class T>
-   struct cbrt_functor
-   {
-       cbrt_functor(T const& target) : a(target){}
-       std::tr1::tuple<T, T, T> operator()(T const& z)
-       {
-         T sqr = z * z;
-         return std::tr1::make_tuple(sqr * z - a, 3 * sqr, 6 * z);
-       }
-   private:
-       T a;
-   };
+template <class T>
+struct largest_cbrt_int_type
+{
+   typedef typename mpl::if_<
+      boost::is_convertible<boost::uintmax_t, T>,
+      boost::uintmax_t,
+      unsigned int
+   >::type type;
+};
 
 template <class T, class Policy>
 T cbrt_imp(T z, const Policy& pol)
 {
    BOOST_MATH_STD_USING
+   //
+   // cbrt approximation for z in the range [0.5,1]
+   // It's hard to say what number of terms gives the optimum
+   // trade off between precision and performance, this seems
+   // to be about the best for double precision.
+   //
+   // Maximum Deviation Found:                     1.231e-006
+   // Expected Error Term:                         -1.231e-006
+   // Maximum Relative Change in Control Points:   5.982e-004
+   //
+   static const T P[] = { 
+      static_cast<T>(0.37568269008611818),
+      static_cast<T>(1.3304968705558024),
+      static_cast<T>(-1.4897101632445036),
+      static_cast<T>(1.2875573098219835),
+      static_cast<T>(-0.6398703759826468),
+      static_cast<T>(0.13584489959258635),
+   };
+   static const T correction[] = {
+      static_cast<T>(0.62996052494743658238360530363911),  // 2^-2/3
+      static_cast<T>(0.79370052598409973737585281963615),  // 2^-1/3
+      static_cast<T>(1),
+      static_cast<T>(1.2599210498948731647672106072782),   // 2^1/3
+      static_cast<T>(1.5874010519681994747517056392723),   // 2^2/3
+   };
+
+   if(!boost::math::isfinite(z))
+   {
+      return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
+   }
+
    int i_exp, sign(1);
    if(z < 0)
    {
@@ -44,15 +78,73 @@ T cbrt_imp(T z, const Policy& pol)
    if(z == 0)
       return 0;
 
-   frexp(z, &i_exp);
-   T min = static_cast<T>(ldexp(0.5, i_exp/3));
-   T max = static_cast<T>(ldexp(2.0, i_exp/3));
-   T guess = static_cast<T>(ldexp(1.0, i_exp/3));
-   int digits = (policies::digits<T, Policy>()) / 2;
-   boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
-   guess = sign * tools::halley_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits, max_iter);
-   policies::check_root_iterations("boost::math::cbrt<%1%>", max_iter, pol);
-   return guess;
+   T guess = frexp(z, &i_exp);
+   int original_i_exp = i_exp; // save for later
+   guess = tools::evaluate_polynomial(P, guess);
+   int i_exp3 = i_exp / 3;
+
+   typedef typename largest_cbrt_int_type<T>::type shift_type;
+
+   if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
+   {
+      if(i_exp3 > 0)
+         guess *= shift_type(1u) << i_exp3;
+      else
+         guess /= shift_type(1u) << -i_exp3;
+   }
+   else
+   {
+      guess = ldexp(guess, i_exp3);
+   }
+   i_exp %= 3;
+   guess *= correction[i_exp + 2];
+   //
+   // Now inline Halley iteration.
+   // We do this here rather than calling tools::halley_iterate since we can
+   // simplify the expressions algebraically, and don't need most of the error
+   // checking of the boilerplate version as we know in advance that the function
+   // is well behaved...
+   //
+   typedef typename policies::precision<T, Policy>::type prec;
+   typedef typename mpl::divides<prec, mpl::int_<3> >::type prec3;
+   typedef typename mpl::plus<prec3, mpl::int_<3> >::type new_prec;
+   typedef typename policies::normalise<Policy, policies::digits2<new_prec::value> >::type new_policy;
+   //
+   // Epsilon calculation uses compile time arithmetic when it's available for type T,
+   // otherwise uses ldexp to calculate at runtime:
+   //
+   T eps = (new_prec::value > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
+   T diff;
+
+   if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
+   {
+      //
+      // Safe from overflow, use the fast method:
+      //
+      do
+      {
+         T g3 = guess * guess * guess;
+         diff = (g3 + z + z) / (g3 + g3 + z);
+         guess *= diff;
+      }
+      while(fabs(1 - diff) > eps);
+   }
+   else
+   {
+      //
+      // Either we're ready to overflow, or we can't tell because numeric_limits isn't
+      // available for type T:
+      //
+      do
+      {
+         T g2 = guess * guess;
+         diff = (g2 - z / guess) / (2 * guess + z / g2);
+         guess -= diff;
+      }
+      while((guess * eps) < fabs(diff));
+   }
+
+   return sign * guess;
 }
 
 } // namespace detail

include/boost/math/special_functions/fpclassify.hpp

@@ -106,7 +106,7 @@ inline bool is_nan_helper(T t, const boost::true_type&)
 }
 
 template <class T>
-inline bool is_nan_helper(T t, const boost::false_type&)
+inline bool is_nan_helper(T, const boost::false_type&)
 {
    return false;
 }

minimax/f.cpp

@@ -318,6 +318,11 @@ boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant)
             y = (y == 0) ? 1e-100 : boost::math::sign(y) * 1e-100;
          return (boost::math::log1p(y) - y + y * y / 2) / (y);
       }
+   case 29:
+      {
+         // cbrt over [0.5, 1]
+         return boost::math::cbrt(x);
+      }
    }
    return 0;
 }

performance/Jamfile.v2

@@ -7,13 +7,14 @@
 exe math_performance : 
 distributions.cpp
 main.cpp
+test_cbrt.cpp
 test_erf.cpp
 test_gamma.cpp
 test_ibeta.cpp
 test_igamma.cpp
 test_polynomial.cpp
 test_reference.cpp 
-test_lanczos.cpp  /boost/regex//boost_regex
+/boost/regex//boost_regex
 : 
     <define>BOOST_ALL_NO_LIB
     <link>static

performance/test_cbrt.cpp

@@ -0,0 +1,101 @@
+//  Copyright John Maddock 2007.
+//  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)
+
+#include "required_defines.hpp"
+#include "performance_measure.hpp"
+
+#include <boost/math/special_functions/cbrt.hpp>
+
+double cbrt_test()
+{
+   double result = 0;
+   double val = 1e-100;
+   for(int i = 0; i < 1000; ++i)
+   {
+      val *= 1.5;
+      result += boost::math::cbrt(val);
+   }
+   return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(cbrt_test, "cbrt")
+{
+   double result = cbrt_test();
+
+   consume_result(result);
+   set_call_count(1000);
+}
+
+double cbrt_pow_test()
+{
+   double result = 0;
+   double val = 1e-100;
+   for(int i = 0; i < 1000; ++i)
+   {
+      val *= 1.5;
+      result += std::pow(val, 0.33333333333333333333333333333333);
+   }
+   return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(cbrt_pow_test, "cbrt-pow")
+{
+   double result = cbrt_pow_test();
+
+   consume_result(result);
+   set_call_count(1000);
+}
+
+#ifdef TEST_CEPHES
+
+extern "C" double cbrt(double);
+
+double cbrt_cephes_test()
+{
+   double result = 0;
+   double val = 1e-100;
+   for(int i = 0; i < 1000; ++i)
+   {
+      val *= 1.5;
+      result += ::cbrt(val);
+   }
+   return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(cbrt_test, "cbrt-cephes")
+{
+   double result = cbrt_cephes_test();
+
+   consume_result(result);
+   set_call_count(1000);
+}
+
+#endif
+
+#if defined(__GNUC__) && (__GNUC__ >= 4)
+
+#include <math.h>
+
+double cbrt_c99_test()
+{
+   double result = 0;
+   double val = 1e-100;
+   for(int i = 0; i < 1000; ++i)
+   {
+      val *= 1.5;
+      result += ::cbrt(val);
+   }
+   return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(cbrt_c99_test, "cbrt-c99")
+{
+   double result = cbrt_c99_test();
+
+   consume_result(result);
+   set_call_count(1000);
+}
+
+#endif

src/tr1/llround.cpp

@@ -13,7 +13,7 @@
 
 namespace boost{ namespace math{ namespace tr1{
 
-extern "C" long long BOOST_MATH_TR1_DECL llround BOOST_PREVENT_MACRO_SUBSTITUTION(double x)
+extern "C" boost::long_long_type BOOST_MATH_TR1_DECL llround BOOST_PREVENT_MACRO_SUBSTITUTION(double x)
 {
    return c_policies::llround BOOST_PREVENT_MACRO_SUBSTITUTION(x);
 }

src/tr1/llroundf.cpp

@@ -13,7 +13,7 @@
 
 namespace boost{ namespace math{ namespace tr1{
 
-extern "C" long long BOOST_MATH_TR1_DECL llroundf BOOST_PREVENT_MACRO_SUBSTITUTION(float x)
+extern "C" boost::long_long_type BOOST_MATH_TR1_DECL llroundf BOOST_PREVENT_MACRO_SUBSTITUTION(float x)
 {
    return c_policies::llround BOOST_PREVENT_MACRO_SUBSTITUTION(x);
 }

src/tr1/llroundl.cpp

@@ -13,7 +13,7 @@
 
 namespace boost{ namespace math{ namespace tr1{
 
-extern "C" long long BOOST_MATH_TR1_DECL llroundl BOOST_PREVENT_MACRO_SUBSTITUTION(long double x)
+extern "C" boost::long_long_type BOOST_MATH_TR1_DECL llroundl BOOST_PREVENT_MACRO_SUBSTITUTION(long double x)
 {
    return c_policies::llround BOOST_PREVENT_MACRO_SUBSTITUTION(x);
 }

test/test_cbrt.cpp

@@ -58,6 +58,13 @@ void expected_results()
       "long double",                  // test type(s)
       ".*",                   // test data group
       ".*", 10, 6);                 // test function
+   add_expected_result(
+      ".*",                          // compiler
+      ".*",                          // stdlib
+      ".*",                          // platform
+      "(long )?double|real_concept", // test type(s)
+      ".*",                   // test data group
+      ".*", 2, 2);                 // test function
 }
 
 struct negative_cbrt

test/test_ibeta_inv.cpp

@@ -274,6 +274,7 @@ void test_inverses2(const T& data, const char* type_name, const char* test_name)
 template <class T>
 void test_beta(T, const char* name)
 {
+   (void)name;
    //
    // The actual test data is rather verbose, so it's in a separate file
    //

test/test_sign.cpp

@@ -19,7 +19,7 @@
    using std::setprecision;
 
 template <class RealType>
-void test_spots(RealType /*T*/, const char* type_name)
+void test_spots(RealType /*T*/, const char* /*type_name*/)
 {
    // Basic sanity checks.
    RealType a = 0;