From: john_at_[hidden]
Date: 2008-03-23 06:32:22
Author: johnmaddock
Date: 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
New Revision: 43800
URL: http://svn.boost.org/trac/boost/changeset/43800
Log:
Added new pow function from Bruno Lalande.
Fixed a few bugs in the non-central distro's that could cause infinite looping.
Added non central distros to the performance test app.
Added:
sandbox/math_toolkit/boost/math/special_functions/pow.hpp (contents, props changed)
sandbox/math_toolkit/libs/math/doc/sf_and_dist/pow.qbk (contents, props changed)
sandbox/math_toolkit/libs/math/test/pow_test.cpp (contents, props changed)
Text files modified:
sandbox/math_toolkit/boost/math/concepts/real_concept.hpp | 4
sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp | 16 ++
sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp | 2
sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp | 211 ++++++++++++++++++++++++++++++++++++++++
sandbox/math_toolkit/boost/math/special_functions.hpp | 1
sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp | 3
sandbox/math_toolkit/libs/math/doc/sf_and_dist/math.qbk | 5
sandbox/math_toolkit/libs/math/doc/sf_and_dist/powers.qbk | 2
sandbox/math_toolkit/libs/math/performance/distributions.cpp | 180 ++++++++++++++++++++++++++++++++++
sandbox/math_toolkit/libs/math/performance/main.cpp | 7
sandbox/math_toolkit/libs/math/test/Jamfile.v2 | 1
sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp | 3
12 files changed, 427 insertions(+), 8 deletions(-)
Modified: sandbox/math_toolkit/boost/math/concepts/real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/real_concept.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -111,6 +111,10 @@
{ return -m_value; }
real_concept const& operator+()const
{ return *this; }
+ real_concept& operator++()
+ { ++m_value; return *this; }
+ real_concept& operator--()
+ { --m_value; return *this; }
private:
real_concept_base_type m_value;
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -45,6 +45,8 @@
// maximum of the poisson weighting term:
//
int k = itrunc(l2);
+ if(k == 0)
+ k = 1;
// Starting Poisson weight:
T pois = gamma_p_derivative(T(k+1), l2, pol);
if(pois == 0)
@@ -74,7 +76,7 @@
{
T term = beta * pois;
sum += term;
- if((fabs(term/sum) < errtol) && (last_term >= term))
+ if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0))
{
count = k - i;
break;
@@ -92,7 +94,7 @@
T term = poisf * betaf;
sum += term;
- if(fabs(term/sum) < errtol)
+ if((fabs(term/sum) < errtol) || (term == 0))
{
break;
}
@@ -122,6 +124,8 @@
// maximum of the poisson weighting term:
//
int k = itrunc(l2);
+ if(k == 0)
+ k = 1;
// Starting Poisson weight:
T pois = gamma_p_derivative(T(k+1), l2, pol);
if(pois == 0)
@@ -154,7 +158,7 @@
T term = poisf * betaf;
sum += term;
- if((fabs(term/sum) < errtol) && (last_term > term))
+ if((fabs(term/sum) < errtol) && (last_term >= term))
{
count = i - k;
break;
@@ -774,6 +778,9 @@
Policy()))
return (RealType)r;
+ if(l == 0)
+ return cdf(beta_distribution<RealType, Policy>(a, b), x);
+
return detail::non_central_beta_cdf(x, 1 - x, a, b, l, false, Policy());
} // cdf
@@ -808,6 +815,9 @@
Policy()))
return (RealType)r;
+ if(l == 0)
+ return cdf(complement(beta_distribution<RealType, Policy>(a, b), x));
+
return detail::non_central_beta_cdf(x, 1 - x, a, b, l, true, Policy());
} // ccdf
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -603,7 +603,7 @@
// Can't do a thing if one of p and q is zero:
//
return policies::raise_evaluation_error<RealType>(function,
- "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
+ "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
}
non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -13,6 +13,7 @@
#include <boost/math/distributions/fwd.hpp>
#include <boost/math/distributions/non_central_beta.hpp> // for nc beta
#include <boost/math/distributions/normal.hpp> // for normal CDF and quantile
+#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/detail/generic_quantile.hpp> // quantile
namespace boost
@@ -476,6 +477,122 @@
return result;
}
+#if 0
+ //
+ // This code is disabled, since there can be multiple answers to the
+ // question, and it's not clear how to find the "right" one.
+ //
+ template <class RealType, class Policy>
+ struct t_degrees_of_freedom_finder
+ {
+ t_degrees_of_freedom_finder(
+ RealType delta_, RealType x_, RealType p_, bool c)
+ : delta(delta_), x(x_), p(p_), comp(c) {}
+
+ RealType operator()(const RealType& v)
+ {
+ non_central_t_distribution<RealType, Policy> d(v, delta);
+ return comp ?
+ p - cdf(complement(d, x))
+ : cdf(d, x) - p;
+ }
+ private:
+ RealType delta;
+ RealType x;
+ RealType p;
+ bool comp;
+ };
+
+ template <class RealType, class Policy>
+ inline RealType find_t_degrees_of_freedom(
+ RealType delta, RealType x, RealType p, RealType q, const Policy& pol)
+ {
+ const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
+ if((p == 0) || (q == 0))
+ {
+ //
+ // Can't a thing if one of p and q is zero:
+ //
+ return policies::raise_evaluation_error<RealType>(function,
+ "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
+ RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
+ }
+ t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true);
+ tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+ //
+ // Pick an initial guess:
+ //
+ RealType guess = 200;
+ std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
+ f, guess, RealType(2), false, tol, max_iter, pol);
+ RealType result = ir.first + (ir.second - ir.first) / 2;
+ if(max_iter >= policies::get_max_root_iterations<Policy>())
+ {
+ policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
+ " or there is no answer to problem. Current best guess is %1%", result, Policy());
+ }
+ return result;
+ }
+
+ template <class RealType, class Policy>
+ struct t_non_centrality_finder
+ {
+ t_non_centrality_finder(
+ RealType v_, RealType x_, RealType p_, bool c)
+ : v(v_), x(x_), p(p_), comp(c) {}
+
+ RealType operator()(const RealType& delta)
+ {
+ non_central_t_distribution<RealType, Policy> d(v, delta);
+ return comp ?
+ p - cdf(complement(d, x))
+ : cdf(d, x) - p;
+ }
+ private:
+ RealType v;
+ RealType x;
+ RealType p;
+ bool comp;
+ };
+
+ template <class RealType, class Policy>
+ inline RealType find_t_non_centrality(
+ RealType v, RealType x, RealType p, RealType q, const Policy& pol)
+ {
+ const char* function = "non_central_t<%1%>::find_t_non_centrality";
+ if((p == 0) || (q == 0))
+ {
+ //
+ // Can't do a thing if one of p and q is zero:
+ //
+ return policies::raise_evaluation_error<RealType>(function,
+ "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
+ RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
+ }
+ t_non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
+ tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+ //
+ // Pick an initial guess that we know is the right side of
+ // zero:
+ //
+ RealType guess;
+ if(f(0) < 0)
+ guess = 1;
+ else
+ guess = -1;
+ std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
+ f, guess, RealType(2), false, tol, max_iter, pol);
+ RealType result = ir.first + (ir.second - ir.first) / 2;
+ if(max_iter >= policies::get_max_root_iterations<Policy>())
+ {
+ policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
+ " or there is no answer to problem. Current best guess is %1%", result, Policy());
+ }
+ return result;
+ }
+#endif
} // namespace detail
template <class RealType = double, class Policy = policies::policy<> >
@@ -507,6 +624,94 @@
{ // Private data getter function.
return ncp;
}
+#if 0
+ //
+ // This code is disabled, since there can be multiple answers to the
+ // question, and it's not clear how to find the "right" one.
+ //
+ static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p)
+ {
+ const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
+ typedef typename policies::evaluation<RealType, Policy>::type value_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float<false>,
+ policies::promote_double<false>,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+ value_type result = detail::find_t_degrees_of_freedom(
+ static_cast<value_type>(delta),
+ static_cast<value_type>(x),
+ static_cast<value_type>(p),
+ static_cast<value_type>(1-p),
+ forwarding_policy());
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ result,
+ function);
+ }
+ template <class A, class B, class C>
+ static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
+ {
+ const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
+ typedef typename policies::evaluation<RealType, Policy>::type value_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float<false>,
+ policies::promote_double<false>,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+ value_type result = detail::find_t_degrees_of_freedom(
+ static_cast<value_type>(c.dist),
+ static_cast<value_type>(c.param1),
+ static_cast<value_type>(1-c.param2),
+ static_cast<value_type>(c.param2),
+ forwarding_policy());
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ result,
+ function);
+ }
+ static RealType find_non_centrality(RealType v, RealType x, RealType p)
+ {
+ const char* function = "non_central_t<%1%>::find_t_non_centrality";
+ typedef typename policies::evaluation<RealType, Policy>::type value_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float<false>,
+ policies::promote_double<false>,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+ value_type result = detail::find_t_non_centrality(
+ static_cast<value_type>(v),
+ static_cast<value_type>(x),
+ static_cast<value_type>(p),
+ static_cast<value_type>(1-p),
+ forwarding_policy());
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ result,
+ function);
+ }
+ template <class A, class B, class C>
+ static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
+ {
+ const char* function = "non_central_t<%1%>::find_t_non_centrality";
+ typedef typename policies::evaluation<RealType, Policy>::type value_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float<false>,
+ policies::promote_double<false>,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+ value_type result = detail::find_t_non_centrality(
+ static_cast<value_type>(c.dist),
+ static_cast<value_type>(c.param1),
+ static_cast<value_type>(1-c.param2),
+ static_cast<value_type>(c.param2),
+ forwarding_policy());
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ result,
+ function);
+ }
+#endif
private:
// Data member, initialized by constructor.
RealType v; // degrees of freedom
@@ -774,6 +979,9 @@
Policy()))
return (RealType)r;
+ if(l == 0)
+ return cdf(students_t_distribution<RealType, Policy>(v), x);
+
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::non_central_t_cdf(
static_cast<value_type>(v),
@@ -817,6 +1025,9 @@
Policy()))
return (RealType)r;
+ if(l == 0)
+ return cdf(complement(students_t_distribution<RealType, Policy>(v), x));
+
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::non_central_t_cdf(
static_cast<value_type>(v),
Modified: sandbox/math_toolkit/boost/math/special_functions.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -52,5 +52,6 @@
#include <boost/math/special_functions/modf.hpp>
#include <boost/math/special_functions/round.hpp>
#include <boost/math/special_functions/trunc.hpp>
+#include <boost/math/special_functions/pow.hpp>
#endif // BOOST_MATH_SPECIAL_FUNCTIONS_HPP
Modified: sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -999,6 +999,9 @@
\
template <class T>\
inline T modf(const T& v, long* ipart){ return boost::math::modf(v, ipart, Policy()); }\
+ \
+ template <int N, class T>\
+ inline typename boost::math::tools::promote_args<T>::type pow(T v){ return boost::math::pow<N>(v, Policy()); }\
#endif // BOOST_MATH_SPECIAL_MATH_FWD_HPP
Added: sandbox/math_toolkit/boost/math/special_functions/pow.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/special_functions/pow.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -0,0 +1,120 @@
+// Boost pow.hpp header file
+// Computes a power with exponent known at compile-time
+
+// (C) Copyright Bruno Lalande 2008.
+// 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)
+
+// See http://www.boost.org for updates, documentation, and revision history.
+
+
+#ifndef BOOST_MATH_POW_HPP
+#define BOOST_MATH_POW_HPP
+
+
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/tools/promotion.hpp>
+#include <boost/mpl/greater_equal.hpp>
+
+
+namespace boost {
+namespace math {
+
+
+namespace detail {
+
+
+template <int N>
+struct positive_power;
+
+template <>
+struct positive_power<0>
+{
+ template <typename T>
+ static typename tools::promote_args<T>::type result(T)
+ { return 1; }
+};
+
+template <>
+struct positive_power<2>
+{
+ template <typename T>
+ static typename tools::promote_args<T>::type result(T base)
+ { return base*base; }
+};
+
+template <int N>
+struct positive_power
+{
+ template <typename T>
+ static typename tools::promote_args<T>::type result(T base)
+ {
+ return (N%2) ? base*positive_power<N-1>::result(base)
+ : positive_power<2>::result(
+ positive_power<N/2>::result(base)
+ );
+ }
+};
+
+
+template <int N, bool>
+struct power_if_positive
+{
+ template <typename T, class Policy>
+ static typename tools::promote_args<T>::type result(T base, const Policy&)
+ { return positive_power<N>::result(base); }
+};
+
+template <int N>
+struct power_if_positive<N, false>
+{
+ template <typename T, class Policy>
+ static typename tools::promote_args<T>::type
+ result(T base, const Policy& policy)
+ {
+ if (base == 0)
+ {
+ return policies::raise_overflow_error<T>(
+ "boost::math::pow(%1%)",
+ "Attempted to compute a negative power of 0",
+ policy
+ );
+ }
+
+ return T(1) / positive_power<-N>::result(base);
+ }
+};
+
+
+template <int N>
+struct select_power_if_positive
+{
+ typedef typename mpl::greater_equal<
+ mpl::int_<N>,
+ mpl::int_<0>
+ >::type is_positive;
+
+ typedef power_if_positive<N, is_positive::value> type;
+};
+
+
+} // namespace detail
+
+
+template <int N, typename T, class Policy>
+inline typename tools::promote_args<T>::type pow(T base, const Policy& policy)
+{ return detail::select_power_if_positive<N>::type::result(base, policy); }
+
+
+template <int N, typename T>
+inline typename tools::promote_args<T>::type pow(T base)
+{ return pow<N>(base, policies::policy<>()); }
+
+
+} // namespace math
+} // namespace boost
+
+
+#endif
Modified: sandbox/math_toolkit/libs/math/doc/sf_and_dist/math.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/sf_and_dist/math.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/sf_and_dist/math.qbk 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -1,13 +1,13 @@
[article Math Toolkit
[quickbook 1.4]
- [copyright 2006, 2007, 2008 John Maddock, Paul A. Bristow, Hubert Holin and Xiaogang Zhang]
+ [copyright 2006, 2007, 2008 John Maddock, Paul A. Bristow, Hubert Holin, Xiaogang Zhang and Bruno Lalande]
[purpose ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22]
[license
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])
]
- [authors [Maddock, John], [Bristow, Paul A.], [Holin, Hubert], [Zhang, Xiaogang]]
+ [authors [Maddock, John], [Bristow, Paul A.], [Holin, Hubert], [Zhang, Xiaogang], [Lalande, Bruno]]
[category math]
[purpose mathematics]
[/last-revision $Date$]
@@ -172,6 +172,7 @@
[def __sqrt1pm1 [link math_toolkit.special.powers.sqrt1pm1 sqrt1pm1]]
[def __powm1 [link math_toolkit.special.powers.powm1 powm1]]
[def __hypot [link math_toolkit.special.powers.hypot hypot]]
+[def __pow [link math_toolkit.special.powers.ct_pow pow]]
[/zeta]
[def __zeta [link math_toolkit.special.zetas.zeta zeta]]
Added: sandbox/math_toolkit/libs/math/doc/sf_and_dist/pow.qbk
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/doc/sf_and_dist/pow.qbk 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -0,0 +1,144 @@
+[section:ct_pow Compile Time Power of a Runtime Base]
+
+The `pow` function effectively computes the compile-time integral
+power of a run-time base.
+
+[h4 Synopsis]
+
+[@../../../../../boost/math/special_functions/pow.hpp `#include <boost/math/special_functions/pow.hpp>`]
+
+ namespace boost { namespace math {
+
+ template <int N, typename T>
+ ``__sf_result`` pow(T base);
+
+ template <int N, typename T, class Policy>
+ ``__sf_result`` pow(T base, const Policy& policy);
+
+ }}
+
+[h4 Rationale and Usage]
+
+Computing the power of a number with an exponent that is known
+at compile time is a common need for programmers. In such cases,
+the usual method is to avoid the overhead implied by
+the `pow`, `powf` and `powl` C functions by hardcoding an expression
+such as:
+
+ // Hand-written 8th power of a 'base' variable
+ double result = base*base*base*base*base*base*base*base;
+
+However, this kind of expression is not really readable (knowing the
+value of the exponent involves counting the number of occurrences of /base/),
+error-prone (it's easy to forget an occurrence), syntactically bulky, and
+non-optimal in terms of performance.
+
+The pow function of Boost.Math helps writing this kind expression along
+with solving all the problems listed above:
+
+ // 8th power of a 'base' variable using math::pow
+ double result = pow<8>(base);
+
+The expression is now shorter, easier to read, safer, and even faster.
+Indeed, `pow` will compute the expression such that only log2(N)
+products are made for a power of N. For instance in the
+example above, the resulting expression will be the same as if we had
+written this, with only one computation of each identical subexpression:
+
+ // Internal effect of pow<8>(base)
+ double result = ((base*base)*(base*base))*((base*base)*(base*base));
+
+Only 3 different products were actually computed.
+
+
+[h4 Return Type]
+
+The return type of these functions is computed using the __arg_pomotion_rules
+when T1 and T2 are different types. For example:
+
+* If T is a `float`, the return type is a `float`.
+* If T is a `long double`, the return type is a `long double`.
+* Otherwise, the return type is a `double`.
+
+[h4 Policies]
+
+[optional_policy]
+
+[h4 Error Handling]
+
+In the case where `pow` is called with a null base and a negative exponent,
+an error occurs since this operation is a division by 0 (it equals to 1/0).
+The error raised is an __overflow_error following the
+[@sf_and_dist/html/math_toolkit/main_overview/error_handling.html
+general policies of error handling in Boost.Math].
+
+The default overflow error policy is `throw_on_error`. A call like `pow<-2>(0)`
+will thus throw a `std::overflow_error` exception. As shown in the
+link given above, other error handling policies can be used:
+
+* `errno_on_error`: Sets `::errno` to `ERANGE` and returns `std::numeric_limits<T>::infinity()`.
+* `ignore_error`: Returns `std::numeric_limits<T>::infinity()`.
+* `user_error`: Returns the result of `boost::math::policies::user_overflow_error`:
+ this function must be defined by the user.
+
+Here is an example of error handling customization where we want to
+specify the result that has to be returned in case of error. We will
+thus use the `user_error` policy, by passing as second argument an
+instance of an overflow_error policy templated with `user_error`:
+
+ // First we open the boost::math::policies namespace and define the `user_overflow_error`
+ // by making it return the value we want in case of error (-1 here)
+
+ namespace boost { namespace math { namespace policies {
+ template <class T>
+ T user_overflow_error(const char*, const char*, const T&)
+ { return -1; }
+ }}}
+
+
+ // Then we invoke pow and indicate that we want to use the user_error policy
+ using boost::math::policies;
+ double result = pow<-5>(base, policy<overflow_error<user_error> >());
+
+ // We can now test the returned value and treat the special case if needed:
+ if (result == -1)
+ {
+ // there was an error, do something...
+ }
+
+Another way is to redefine the default `overflow_error` policy by using the
+BOOST_MATH_OVERFLOW_ERROR_POLICY macro. Once the `user_overflow_error` function
+is defined as above, we can achieve the same result like this:
+
+ // Redefine the default error_overflow policy
+ #define BOOST_MATH_OVERFLOW_ERROR_POLICY user_error
+ #include <boost/math/special_functions/pow.hpp>
+
+ // From this point, passing a policy in argument is no longer needed, a call like this one
+ // will return -1 in case of error:
+
+ double result = pow<-5>(base);
+
+
+[h4 Acknowledgements]
+
+Bruno Lalande submitted this addition to Boost.Math.
+
+'''
+Thanks to Joaquín López Muñoz and Scott McMurray for their help in
+improving the implementation.
+'''
+
+[h4 References]
+
+D.E. Knuth, ['The Art of Computer Programming, Vol. 2: Seminumerical Algorithms], 2nd ed., Addison-Wesley, Reading, MA, 1981
+
+[endsect]
+
+[/
+ Copyright 2008 Bruno Lalande.
+ 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).
+]
+
Modified: sandbox/math_toolkit/libs/math/doc/sf_and_dist/powers.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/sf_and_dist/powers.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/sf_and_dist/powers.qbk 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -268,6 +268,8 @@
[endsect]
+[include pow.qbk]
+
[endsect][/section:powers Logs, Powers, Roots and Exponentials]
Modified: sandbox/math_toolkit/libs/math/performance/distributions.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/distributions.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/distributions.cpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -45,6 +45,20 @@
100000
};
+int small_int_values[] = {
+ 1,
+ 2,
+ 3,
+ 5,
+ 10,
+ 15,
+ 20,
+ 30,
+ 50,
+ 100,
+ 150
+};
+
double real_values[] = {
1e-5,
1e-4,
@@ -58,6 +72,83 @@
100000
};
+#define BOOST_MATH_DISTRIBUTION3_TEST(name, param1_table, param2_table, param3_table, random_variable_table, probability_table) \
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist, name), "dist-" BOOST_STRINGIZE(name) "-cdf")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(param3_table)/sizeof(param3_table[0]);\
+ unsigned d_size = sizeof(random_variable_table)/sizeof(random_variable_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ for(unsigned l = 0; l < d_size; ++l)\
+ {\
+ result += cdf(boost::math:: BOOST_JOIN(name, _distribution) <>(param1_table[i], param2_table[j], param3_table[k]), random_variable_table[l]);\
+ }\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size * d_size);\
+ }\
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist_pdf_, name), "dist-" BOOST_STRINGIZE(name) "-pdf")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(param3_table)/sizeof(param3_table[0]);\
+ unsigned d_size = sizeof(random_variable_table)/sizeof(random_variable_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ for(unsigned l = 0; l < d_size; ++l)\
+ {\
+ result += pdf(boost::math:: BOOST_JOIN(name, _distribution) <>(param1_table[i], param2_table[j], param3_table[k]), random_variable_table[l]);\
+ }\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size * d_size);\
+ }\
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist_quant, name), "dist-" BOOST_STRINGIZE(name) "-quantile")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(param3_table)/sizeof(param3_table[0]);\
+ unsigned d_size = sizeof(probability_table)/sizeof(probability_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ for(unsigned l = 0; l < d_size; ++l)\
+ {\
+ result += quantile(boost::math:: BOOST_JOIN(name, _distribution) <>(param1_table[i], param2_table[j], param3_table[k]), probability_table[l]);\
+ }\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size * d_size);\
+ }
+
#define BOOST_MATH_DISTRIBUTION2_TEST(name, param1_table, param2_table, random_variable_table, probability_table) \
BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist, name), "dist-" BOOST_STRINGIZE(name) "-cdf")\
{\
@@ -190,6 +281,9 @@
BOOST_MATH_DISTRIBUTION1_TEST(students_t, int_values, real_values, probabilities)
BOOST_MATH_DISTRIBUTION2_TEST(weibull, real_values, real_values, real_values, probabilities)
BOOST_MATH_DISTRIBUTION2_TEST(non_central_chi_squared, int_values, int_values, real_values, probabilities)
+BOOST_MATH_DISTRIBUTION3_TEST(non_central_beta, int_values, int_values, real_values, real_values, probabilities)
+BOOST_MATH_DISTRIBUTION3_TEST(non_central_f, int_values, int_values, real_values, real_values, probabilities)
+BOOST_MATH_DISTRIBUTION2_TEST(non_central_t, int_values, small_int_values, real_values, probabilities)
#ifdef TEST_R
@@ -197,7 +291,90 @@
extern "C" {
#include "Rmath.h"
+/*
+double qnchisq(double, double, double, int, int)
+{
+ return 0;
}
+*/
+}
+
+#define BOOST_MATH_R_DISTRIBUTION3_TEST(name, param1_table, param2_table, param3_table, random_variable_table, probability_table) \
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist, name), "dist-" #name "-R-cdf")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(param3_table)/sizeof(param3_table[0]);\
+ unsigned d_size = sizeof(random_variable_table)/sizeof(random_variable_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ for(unsigned l = 0; l < d_size; ++l)\
+ {\
+ result += p##name (random_variable_table[l], param1_table[i], param2_table[j], param3_table[k], 1, 0);\
+ }\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size * d_size);\
+ }\
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist_pdf_, name), "dist-" #name "-R-pdf")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(param3_table)/sizeof(param3_table[0]);\
+ unsigned d_size = sizeof(random_variable_table)/sizeof(random_variable_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ for(unsigned l = 0; l < d_size; ++l)\
+ {\
+ result += d##name (random_variable_table[l], param1_table[i], param2_table[j], param3_table[k], 0);\
+ }\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size * d_size);\
+ }\
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist_quant, name), "dist-" #name "-R-quantile")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(param3_table)/sizeof(param3_table[0]);\
+ unsigned d_size = sizeof(probability_table)/sizeof(probability_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ for(unsigned l = 0; l < d_size; ++l)\
+ {\
+ result += q##name (probability_table[l], param1_table[i], param2_table[j], param3_table[k], 1, 0);\
+ }\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size * d_size);\
+ }
#define BOOST_MATH_R_DISTRIBUTION2_TEST(name, param1_table, param2_table, random_variable_table, probability_table) \
BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist, name), "dist-" #name "-R-cdf")\
@@ -331,6 +508,9 @@
BOOST_MATH_R_DISTRIBUTION1_TEST(t, int_values, real_values, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(weibull, real_values, real_values, real_values, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(nchisq, int_values, int_values, real_values, probabilities)
+BOOST_MATH_R_DISTRIBUTION3_TEST(nf, int_values, int_values, real_values, real_values, probabilities)
+BOOST_MATH_R_DISTRIBUTION3_TEST(nbeta, int_values, int_values, real_values, real_values, probabilities)
+BOOST_MATH_R_DISTRIBUTION2_TEST(nt, int_values, small_int_values, real_values, probabilities)
#endif
Modified: sandbox/math_toolkit/libs/math/performance/main.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/main.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/main.cpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -12,6 +12,7 @@
#include <string>
#include <cstring>
#include <boost/math/tools/config.hpp>
+#include <boost/regex.hpp>
#include <boost/type_traits/is_same.hpp>
#include "performance_measure.hpp"
#include <boost/math/policies/policy.hpp>
@@ -56,6 +57,8 @@
++i;
}
std::cout << "Or use --all to test everything." << std::endl;
+ std::cout << "You can also specify what to test as a regular expression,\n"
+ " for example --dist.* to test all the distributions." << std::endl;
}
void run_tests()
@@ -77,14 +80,14 @@
bool add_named_test(const char* name)
{
bool found = false;
+ boost::regex e(name);
std::set<test_info>::const_iterator a(all_tests().begin()), b(all_tests().end());
while(a != b)
{
- if(std::strcmp(name, (*a).name) == 0)
+ if(regex_match((*a).name, e))
{
found = true;
tests.insert(*a);
- break;
}
++a;
}
Modified: sandbox/math_toolkit/libs/math/test/Jamfile.v2
==============================================================================
--- sandbox/math_toolkit/libs/math/test/Jamfile.v2 (original)
+++ sandbox/math_toolkit/libs/math/test/Jamfile.v2 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -39,6 +39,7 @@
;
run hypot_test.cpp ;
+run pow_test.cpp ;
run log1p_expm1_test.cpp ;
run powm1_sqrtp1m1_test.cpp ;
run special_functions_test.cpp /boost/unit_test//boost_unit_test_framework/<link>static ;
Modified: sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp (original)
+++ sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -245,6 +245,7 @@
long long ll;
boost::math::modf(v1, &ll);
#endif
+ boost::math::pow<2>(v1);
//
// All over again, with a policy this time:
//
@@ -360,6 +361,7 @@
boost::math::llround(v1, pol);
boost::math::modf(v1, &ll, pol);
#endif
+ boost::math::pow<2>(v1, pol);
//
// All over again with the versions in test::
//
@@ -474,6 +476,7 @@
test::llround(v1);
test::modf(v1, &ll);
#endif
+ test::pow<2>(v1);
}
Added: sandbox/math_toolkit/libs/math/test/pow_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/pow_test.cpp 2008-03-23 06:32:21 EDT (Sun, 23 Mar 2008)
@@ -0,0 +1,186 @@
+// Boost pow_test.cpp test file
+// Tests the pow function
+
+// (C) Copyright Bruno Lalande 2008.
+// 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)
+
+
+#include <cmath>
+#include <string>
+#include <iostream>
+
+#include <boost/math/concepts/real_concept.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+#include <boost/typeof/typeof.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/static_assert.hpp>
+
+#include <boost/math/special_functions/pow.hpp>
+
+#include BOOST_TYPEOF_INCREMENT_REGISTRATION_GROUP()
+BOOST_TYPEOF_REGISTER_TYPE(boost::math::concepts::real_concept)
+
+using namespace boost;
+using namespace boost::math;
+
+template <int N, class T>
+void test_pow(T base)
+{
+ typedef typename tools::promote_args<T>::type result_type;
+
+ BOOST_MATH_STD_USING
+
+ if ((base == 0) && N < 0)
+ {
+ BOOST_CHECK_THROW(math::pow<N>(base), std::overflow_error);
+ }
+ else
+ {
+ BOOST_CHECK_CLOSE(math::pow<N>(base),
+ pow(static_cast<result_type>(base), static_cast<result_type>(N)),
+ boost::math::tools::epsilon<result_type>() * 100 * 200); // 200 eps as a %
+ }
+}
+
+template <int N, class T>
+void test_with_big_bases()
+{
+ for (T base = T(); base < T(1000); ++base)
+ test_pow<N>(base);
+}
+
+template <int N, class T>
+void test_with_small_bases()
+{
+ T base = 0.9f;
+ for (int i = 0; i < 10; ++i)
+ {
+ base += base/50;
+ test_pow<N>(base);
+ }
+}
+
+template <class T, int Factor>
+void test_with_small_exponents()
+{
+ test_with_big_bases<0, T>();
+ test_with_big_bases<Factor*1, T>();
+ test_with_big_bases<Factor*2, T>();
+ test_with_big_bases<Factor*3, T>();
+ test_with_big_bases<Factor*5, T>();
+ test_with_big_bases<Factor*6, T>();
+ test_with_big_bases<Factor*7, T>();
+ test_with_big_bases<Factor*8, T>();
+ test_with_big_bases<Factor*9, T>();
+ test_with_big_bases<Factor*10, T>();
+ test_with_big_bases<Factor*11, T>();
+ test_with_big_bases<Factor*12, T>();
+}
+
+template <class T, int Factor>
+void test_with_big_exponents()
+{
+ test_with_small_bases<Factor*50, T>();
+ test_with_small_bases<Factor*100, T>();
+ test_with_small_bases<Factor*150, T>();
+ test_with_small_bases<Factor*200, T>();
+ test_with_small_bases<Factor*250, T>();
+ test_with_small_bases<Factor*300, T>();
+ test_with_small_bases<Factor*350, T>();
+ test_with_small_bases<Factor*400, T>();
+ test_with_small_bases<Factor*450, T>();
+ test_with_small_bases<Factor*500, T>();
+}
+
+
+void test_return_types()
+{
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>('\1')), double>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(L'\2')), double>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(3)), double>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(4u)), double>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(5ul)), double>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(6.0f)), float>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(7.0)), double>::value));
+ BOOST_STATIC_ASSERT((is_same<BOOST_TYPEOF(pow<2>(7.0l)), long double>::value));
+};
+
+
+namespace boost { namespace math { namespace policies {
+template <class T>
+T user_overflow_error(const char*, const char*, const T&)
+{ return 123.456; }
+}}}
+
+
+void test_error_policy()
+{
+ using namespace policies;
+
+ BOOST_CHECK(pow<-2>(
+ 0.0,
+ policy< ::boost::math::policies::overflow_error<user_error> >()
+ )
+ == 123.456);
+}
+
+int test_main(int, char* [])
+{
+ using namespace std;
+
+ cout << "Testing with integral bases and positive small exponents" << endl;
+ test_with_small_exponents<int, 1>();
+ cout << "Testing with integral bases and negative small exponents" << endl;
+ test_with_small_exponents<int, -1>();
+
+ cout << "Testing with float precision bases and positive small exponents" << endl;
+ test_with_small_exponents<float, 1>();
+ cout << "Testing with float precision bases and negative small exponents" << endl;
+ test_with_small_exponents<float, -1>();
+
+ cout << "Testing with float precision bases and positive big exponents" << endl;
+ test_with_big_exponents<float, 1>();
+ cout << "Testing with float precision bases and negative big exponents" << endl;
+ test_with_big_exponents<float, -1>();
+
+ cout << "Testing with double precision bases and positive small exponents" << endl;
+ test_with_small_exponents<double, 1>();
+ cout << "Testing with double precision bases and negative small exponents" << endl;
+ test_with_small_exponents<double, -1>();
+
+ cout << "Testing with double precision bases and positive big exponents" << endl;
+ test_with_big_exponents<double, 1>();
+ cout << "Testing with double precision bases and negative big exponents" << endl;
+ test_with_big_exponents<double, -1>();
+
+ cout << "Testing with long double precision bases and positive small exponents" << endl;
+ test_with_small_exponents<long double, 1>();
+ cout << "Testing with long double precision bases and negative small exponents" << endl;
+ test_with_small_exponents<long double, -1>();
+
+ cout << "Testing with long double precision bases and positive big exponents" << endl;
+ test_with_big_exponents<long double, 1>();
+ cout << "Testing with long double precision bases and negative big exponents" << endl;
+ test_with_big_exponents<long double, -1>();
+
+ cout << "Testing with concepts::real_concept precision bases and positive small exponents" << endl;
+ test_with_small_exponents<concepts::real_concept, 1>();
+ cout << "Testing with concepts::real_concept precision bases and negative small exponents" << endl;
+ test_with_small_exponents<concepts::real_concept, -1>();
+
+ cout << "Testing with concepts::real_concept precision bases and positive big exponents" << endl;
+ test_with_big_exponents<concepts::real_concept, 1>();
+ cout << "Testing with concepts::real_concept precision bases and negative big exponents" << endl;
+ test_with_big_exponents<concepts::real_concept, -1>();
+
+ test_return_types();
+
+ test_error_policy();
+
+ return 0;
+}
+