[Boost-commit] svn:boost r43800 - in sandbox/math_toolkit: boost/math boost/math/concepts boost/math/distributions boost/math/special_functions libs/math/doc/sf_and_dist libs/math/performance libs/math/test libs/math/test/compile_test

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&#xed;n L&#xf3;pez Mu&#xf1;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;
+}
+