[Boost-commit] svn:boost r82616 - in trunk/boost/math/special_functions: . detail

Subject: [Boost-commit] svn:boost r82616 - in trunk/boost/math/special_functions: . detail
From: e_float_at_[hidden]
Date: 2013-01-25 17:05:50

Author: christopher_kormanyos
Date: 2013-01-25 17:05:49 EST (Fri, 25 Jan 2013)
New Revision: 82616
URL: http://svn.boost.org/trac/boost/changeset/82616

Log:
Re-corrected the behavior of Bessel zeros when used with multiprecision and et_on.
Text files modified: 
   trunk/boost/math/special_functions/bessel.hpp                |    30 ++++++++++++++++++++++--------          
   trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp |    40 +++++++++++++++++++++-------------------
   trunk/boost/math/special_functions/math_fwd.hpp              |    14 --------------                          
   3 files changed, 43 insertions(+), 41 deletions(-)

Modified: trunk/boost/math/special_functions/bessel.hpp
==============================================================================
--- trunk/boost/math/special_functions/bessel.hpp	(original)
+++ trunk/boost/math/special_functions/bessel.hpp	2013-01-25 17:05:49 EST (Fri, 25 Jan 2013)
@@ -303,7 +303,7 @@
 
    if(floor(v) == v)
    {
-      if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
+      if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (4 * fabs(x) * sqrt(sqrt(sqrt(14 * tools::epsilon<T>() * (fabs(x) + fabs(v) / 2) / (5120 * fabs(x))))) > fabs(v)))
       {
          T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
          if((v < 0) && (itrunc(v, pol) & 1))
@@ -332,7 +332,7 @@
    BOOST_MATH_INSTRUMENT_VARIABLE(v);
    BOOST_MATH_INSTRUMENT_VARIABLE(x);
 
-   if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
+   if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (4 * fabs(x) * sqrt(sqrt(sqrt(14 * tools::epsilon<T>() * (fabs(x) + abs(v) / 2) / (5120 * x)))) > abs(v)))
    {
       T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
       if((v < 0) && (v & 1))
@@ -376,7 +376,7 @@
 
    // Handle negative order or if the zero'th zero is requested.
    // Return NaN if NaN is available or return 0 if NaN is not available.
-   if((v < 0) || (m == 0U))
+   if((v < T(0)) || (m == 0U))
       return (std::numeric_limits<T>::has_quiet_NaN ? std::numeric_limits<T>::quiet_NaN() : T(0));
 
    // Set up the initial guess for the upcoming root-finding.
@@ -476,7 +476,7 @@
 
       while(out_it != end_it)
       {
-         *out_it = boost::math::detail::cyl_neumann_zero_imp(v, start_index, pol);
+         *out_it = boost::math::detail::cyl_neumann_zero_imp<T, Policy>(v, start_index, pol);
          ++start_index;
          ++out_it;
       }
@@ -584,13 +584,15 @@
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_zero_imp<value_type>(v, m, pol), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)");
 }
 
 template <class T>
 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, unsigned m)
 {
-   return cyl_bessel_j_zero(v, m, policies::policy<>());
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<T>::is_integer, "Order must be a floating-point type.");
+   return cyl_bessel_j_zero_imp<T, policies::policy<> >(v, m, policies::policy<>());
 }
 
 template <class output_iterator, class T, class Policy>
@@ -603,6 +605,7 @@
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
    detail::cyl_bessel_j_zero_imp<output_iterator, value_type, Policy>(out_it, v, number_of_zeros, start_index, pol);
 }
 
@@ -612,7 +615,11 @@
                               unsigned number_of_zeros,
                               unsigned start_index)
 {
-   detail::cyl_bessel_j_zero_imp(out_it, v, number_of_zeros, start_index, policies::policy<>());
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename detail::bessel_traits<T, T, policies::policy<> >::result_type result_type;
+   typedef typename policies::evaluation<result_type, policies::policy<> >::type value_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
+   detail::cyl_bessel_j_zero_imp<output_iterator, value_type, policies::policy<> >(out_it, v, number_of_zeros, start_index, policies::policy<>());
 }
 
 template <class T, class Policy>
@@ -621,13 +628,15 @@
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_zero_imp<value_type>(v, m, pol), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)");
 }
 
 template <class T>
 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neumann_zero(T v, unsigned m)
 {
-   return cyl_neumann_zero(v, m, policies::policy<>());
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<T>::is_integer, "Order must be a floating-point type.");
+   return cyl_neumann_zero_imp<T, policies::policy<> >(v, m, policies::policy<>());
 }
 
 template <class output_iterator, class T, class Policy>
@@ -640,6 +649,7 @@
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
    detail::cyl_neumann_zero_imp<output_iterator, value_type, Policy>(out_it, v, number_of_zeros, start_index, pol);
 }
 
@@ -649,7 +659,11 @@
                              unsigned number_of_zeros,
                              unsigned start_index)
 {
-   detail::cyl_neumann_zero_imp(out_it, v, number_of_zeros, start_index, policies::policy<>());
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename detail::bessel_traits<T, T, policies::policy<> >::result_type result_type;
+   typedef typename policies::evaluation<result_type, policies::policy<> >::type value_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
+   detail::cyl_neumann_zero_imp<output_iterator, value_type, policies::policy<> >(out_it, v, number_of_zeros, start_index, policies::policy<>());
 }
 
 } // namespace math

Modified: trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp	2013-01-25 17:05:49 EST (Fri, 25 Jan 2013)
@@ -19,6 +19,7 @@
   #include <algorithm>
   #include <boost/math/constants/constants.hpp>
   #include <boost/math/special_functions/math_fwd.hpp>
+  #include <boost/math/special_functions/cbrt.hpp>
   #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
 
   namespace boost { namespace math {
@@ -33,7 +34,7 @@
     namespace bessel_zero
     {
       template<class T>
-      T equation_nist_10_21_19(T v, T a)
+      T equation_nist_10_21_19(const T& v, const T& a)
       {
         // Get the initial estimate of the m'th root of Jv or Yv.
         // This subroutine is used for the order m with m > 1.
@@ -60,9 +61,9 @@
       class equation_as_9_3_39_and_its_derivative
       {
       public:
-        equation_as_9_3_39_and_its_derivative(T zt) : zeta(zt) { }
+        equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }
 
-        boost::math::tuple<T, T> operator()(T z) const
+        boost::math::tuple<T, T> operator()(const T& z) const
         {
           BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt.
 
@@ -86,9 +87,9 @@
       };
 
       template<class T>
-      static T equation_as_9_5_26(T v, T ai_bi_root)
+      static T equation_as_9_5_26(const T& v, const T& ai_bi_root)
       {
-        BOOST_MATH_STD_USING // ADL of std names, needed for pow.
+        BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
 
         // Obtain the estimate of the m'th zero of Jv or Yv.
         // The order m has been used to create the input parameter ai_bi_root.
@@ -107,9 +108,12 @@
         // to refine the value of the estimate of the root of z
         // as a function of zeta.
 
+        const T v_pow_third(boost::math::cbrt(v));
+        const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
+
         // Obtain zeta using the order v combined with the m'th root of
         // an airy function, as shown in  A&S Eq. 9.5.22.
-        const T zeta = pow(v, (T(-2) / 3U)) * (-ai_bi_root);
+        const T zeta = v_pow_minus_two_thirds * (-ai_bi_root);
 
         const T zeta_sqrt = sqrt(zeta);
 
@@ -165,11 +169,10 @@
       namespace cyl_bessel_j_zero_detail
       {
         template<class T>
-        T equation_nist_10_21_40_a(T v)
+        T equation_nist_10_21_40_a(const T& v)
         {
-          BOOST_MATH_STD_USING // ADL of std names, needed for pow.
-
-          const T v_pow_minus_two_thirds = pow(v, T(-2) / 3U);
+          const T v_pow_third(boost::math::cbrt(v));
+          const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
 
           return v * (((((                         + T(0.043)
                           * v_pow_minus_two_thirds - T(0.0908))
@@ -180,7 +183,7 @@
         }
 
         template<class T>
-        T initial_guess(T v, unsigned m)
+        T initial_guess(const T& v, unsigned m)
         {
           // Compute an estimate of the m'th root of cyl_bessel_j.
 
@@ -241,11 +244,11 @@
         class function_object
         {
         public:
-          function_object(T v,
+          function_object(const T& v,
                           const Policy& pol) : my_v(v),
                                                my_pol(pol) { }
 
-          boost::math::tuple<T, T> operator()(T x) const
+          boost::math::tuple<T, T> operator()(const T& x) const
           {
             // Obtain Jv(x) and Jv'(x).
             const T j_v      (boost::math::detail::cyl_bessel_j_imp(  my_v,      x, boost::math::detail::bessel_no_int_tag(), my_pol));
@@ -265,11 +268,10 @@
       namespace cyl_neumann_zero_detail
       {
         template<class T>
-        T equation_nist_10_21_40_b(T v)
+        T equation_nist_10_21_40_b(const T& v)
         {
-          BOOST_MATH_STD_USING // ADL of std names, needed for pow.
-
-          const T v_pow_minus_two_thirds = pow(v, T(-2) / 3U);
+          const T v_pow_third(boost::math::cbrt(v));
+          const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
 
           return v * (((((                         - T(0.001)
                           * v_pow_minus_two_thirds - T(0.0060))
@@ -341,11 +343,11 @@
         class function_object
         {
         public:
-          function_object(T v,
+          function_object(const T& v,
                           const Policy& pol) : my_v(v),
                                                my_pol(pol) { }
 
-          boost::math::tuple<T, T> operator()(T x) const
+          boost::math::tuple<T, T> operator()(const T& x) const
           {
             // Obtain Yv(x) and Yv'(x).
             const T y_v      (boost::math::detail::cyl_neumann_imp(  my_v,      x, boost::math::detail::bessel_no_int_tag(), my_pol));

Modified: trunk/boost/math/special_functions/math_fwd.hpp
==============================================================================
--- trunk/boost/math/special_functions/math_fwd.hpp	(original)
+++ trunk/boost/math/special_functions/math_fwd.hpp	2013-01-25 17:05:49 EST (Fri, 25 Jan 2013)
@@ -1157,24 +1157,10 @@
    inline typename boost::math::detail::bessel_traits<T, T, Policy >::result_type cyl_bessel_j_zero(T v, unsigned m)\
    { return boost::math::cyl_bessel_j_zero(v, m, Policy()); }\
 \
-template <class output_iterator, class T>\
-   inline void cyl_bessel_j_zero_imp(output_iterator out_it,\
-                                     T v,\
-                                     std::size_t number_of_zeros,\
-                                     unsigned start_index)\
-   { boost::math::cyl_bessel_j_zero(out_it, v, number_of_zeros, start_index, Policy(); }\
-\
    template <class T>\
    inline typename boost::math::detail::bessel_traits<T, T, Policy >::result_type cyl_neumann_zero(T v, unsigned m)\
    { return boost::math::cyl_neumann_zero(v, m, Policy()); }\
 \
-template <class output_iterator, class T>\
-   inline void cyl_neumann_zero_imp(output_iterator out_it,\
-                                     T v,\
-                                     std::size_t number_of_zeros,\
-                                     unsigned start_index)\
-   { boost::math::cyl_neumann_zero(out_it, v, number_of_zeros, start_index, Policy(); }\
-\
    template <class T>\
    inline typename boost::math::tools::promote_args<T>::type sin_pi(T x){ return boost::math::sin_pi(x); }\
 \