Subject: [Boost-commit] svn:boost r59095 - in trunk: boost/math/distributions/detail boost/math/special_functions libs/math/minimax libs/math/performance libs/math/src/tr1 libs/math/test
From: john_at_[hidden]
Date: 2010-01-17 12:28:36
Author: johnmaddock
Date: 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
New Revision: 59095
URL: http://svn.boost.org/trac/boost/changeset/59095
Log:
Change cbrt implementation to use a better performing algorithm.
Fix a few warnings along the way.
Added:
trunk/libs/math/performance/test_cbrt.cpp (contents, props changed)
Text files modified:
trunk/boost/math/distributions/detail/hypergeometric_pdf.hpp | 2
trunk/boost/math/special_functions/cbrt.hpp | 136 +++++++++++++++++++++++++++++++++------
trunk/boost/math/special_functions/fpclassify.hpp | 2
trunk/libs/math/minimax/f.cpp | 5 +
trunk/libs/math/performance/Jamfile.v2 | 3
trunk/libs/math/src/tr1/llround.cpp | 2
trunk/libs/math/src/tr1/llroundf.cpp | 2
trunk/libs/math/src/tr1/llroundl.cpp | 2
trunk/libs/math/test/test_cbrt.cpp | 7 ++
trunk/libs/math/test/test_ibeta_inv.cpp | 1
trunk/libs/math/test/test_sign.cpp | 2
11 files changed, 135 insertions(+), 29 deletions(-)
Modified: trunk/boost/math/distributions/detail/hypergeometric_pdf.hpp
==============================================================================
--- trunk/boost/math/distributions/detail/hypergeometric_pdf.hpp (original)
+++ trunk/boost/math/distributions/detail/hypergeometric_pdf.hpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -327,7 +327,7 @@
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);
}
Modified: trunk/boost/math/special_functions/cbrt.hpp
==============================================================================
--- trunk/boost/math/special_functions/cbrt.hpp (original)
+++ trunk/boost/math/special_functions/cbrt.hpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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 @@
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
Modified: trunk/boost/math/special_functions/fpclassify.hpp
==============================================================================
--- trunk/boost/math/special_functions/fpclassify.hpp (original)
+++ trunk/boost/math/special_functions/fpclassify.hpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -106,7 +106,7 @@
}
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;
}
Modified: trunk/libs/math/minimax/f.cpp
==============================================================================
--- trunk/libs/math/minimax/f.cpp (original)
+++ trunk/libs/math/minimax/f.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -318,6 +318,11 @@
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;
}
Modified: trunk/libs/math/performance/Jamfile.v2
==============================================================================
--- trunk/libs/math/performance/Jamfile.v2 (original)
+++ trunk/libs/math/performance/Jamfile.v2 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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
Added: trunk/libs/math/performance/test_cbrt.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/performance/test_cbrt.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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
Modified: trunk/libs/math/src/tr1/llround.cpp
==============================================================================
--- trunk/libs/math/src/tr1/llround.cpp (original)
+++ trunk/libs/math/src/tr1/llround.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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);
}
Modified: trunk/libs/math/src/tr1/llroundf.cpp
==============================================================================
--- trunk/libs/math/src/tr1/llroundf.cpp (original)
+++ trunk/libs/math/src/tr1/llroundf.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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);
}
Modified: trunk/libs/math/src/tr1/llroundl.cpp
==============================================================================
--- trunk/libs/math/src/tr1/llroundl.cpp (original)
+++ trunk/libs/math/src/tr1/llroundl.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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);
}
Modified: trunk/libs/math/test/test_cbrt.cpp
==============================================================================
--- trunk/libs/math/test/test_cbrt.cpp (original)
+++ trunk/libs/math/test/test_cbrt.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -58,6 +58,13 @@
"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
Modified: trunk/libs/math/test/test_ibeta_inv.cpp
==============================================================================
--- trunk/libs/math/test/test_ibeta_inv.cpp (original)
+++ trunk/libs/math/test/test_ibeta_inv.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -274,6 +274,7 @@
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
//
Modified: trunk/libs/math/test/test_sign.cpp
==============================================================================
--- trunk/libs/math/test/test_sign.cpp (original)
+++ trunk/libs/math/test/test_sign.cpp 2010-01-17 12:28:34 EST (Sun, 17 Jan 2010)
@@ -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;