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

Subject: [Boost-commit] svn:boost r83095 - in trunk/boost/math/special_functions: . detail
From: e_float_at_[hidden]
Date: 2013-02-23 07:50:46

Author: christopher_kormanyos
Date: 2013-02-23 07:50:45 EST (Sat, 23 Feb 2013)
New Revision: 83095
URL: http://svn.boost.org/trac/boost/changeset/83095

Log:
Added support for negative order to cyl_neumann_zero().
Text files modified: 
   trunk/boost/math/special_functions/bessel.hpp                |    62 ++++++++++++++++-----                   
   trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp |   113 ++++++++++++++++++++++++++++++++++++--- 
   2 files changed, 149 insertions(+), 26 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-02-23 07:50:45 EST (Sat, 23 Feb 2013)
@@ -357,12 +357,13 @@
 template <class T, class Policy>
 inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
 {
+   BOOST_MATH_STD_USING // ADL of std names, needed for floor.
+
    static const char* function = "boost::math::cyl_bessel_j_zero<%1%>(%1%, int)";
 
    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
 
-   // Handle the zero'th zero, if requested.
-   // Return NaN if NaN is available or return 0 if NaN is not available.
+   // Handle the non-finite order.
    if (!(boost::math::isfinite)(v) )
    {
      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
@@ -374,12 +375,13 @@
       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol);
    }
 
+   // Get the absolute value of the order.
    const bool order_is_negative = (v < 0);
-
    const T vv((!order_is_negative) ? v : -v);
 
-   const bool order_is_zero     = (vv < half_epsilon);
-   const bool order_is_integer  = ((vv - floor(vv)) < half_epsilon);
+   // Check if the order is very close to zero or very close to an integer.
+   const bool order_is_zero    = (vv < half_epsilon);
+   const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
 
    if(m == 0)
    {
@@ -389,9 +391,10 @@
          return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of J0, but the rank must be > 0 !", m, pol);
       }
 
+      // The zero'th zero of Jv(x) for v < 0 is not defined
+      // unless the order is a negative integer.
       if(order_is_negative && (!order_is_integer))
       {
-         // The zero'th zero of Jv(x) for v < 0 is only defined for negative integer order.
          // For non-integer, negative order, requesting the zero'th zero raises a domain error.
          return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Jv for negative, non-integer order, but the rank must be > 0 !", m, pol);
       }
@@ -401,6 +404,8 @@
    }
 
    // Set up the initial guess for the upcoming root-finding.
+   // If the order is a negative integer, then use the corresponding
+   // positive integer for the order.
    const T guess_root = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess<T, Policy>((order_is_integer ? vv : v), m, pol);
 
    // Select the maximum allowed iterations from the policy.
@@ -415,8 +420,8 @@
       boost::math::tools::newton_raphson_iterate(
          boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv_and_jv_prime<T, Policy>((order_is_integer ? vv : v), order_is_zero, pol),
          guess_root,
-         (std::max)(T(guess_root - 0.3F), T(0)),
-         T(guess_root + 0.3F),
+         (std::max)(T(guess_root - 0.2F), T(0)),
+         T(guess_root + 0.2F),
          my_digits2,
          number_of_iterations);
 
@@ -432,23 +437,48 @@
 template <class T, class Policy>
 inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
 {
+   BOOST_MATH_STD_USING // ADL of std names, needed for floor.
+
    static const char* function = "boost::math::cyl_neumann_zero<%1%>(%1%, int)";
-   // Handle negative order or if the zero'th zero is requested.
+
+   // Handle non-finite order.
    if (!(boost::math::isfinite)(v) )
    {
      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
    }
-   else if(v < 0)
+
+   // Handle if the zero'th root or a negative root is requested.
+   if(m < 0)
    {
-      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be >= 0 !", v, pol);
+      return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol);
    }
-   if(m <= 0)
+
+   const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
+
+   // Get the absolute value of the order.
+   const bool order_is_negative = (v < 0);
+   const T vv((!order_is_negative) ? v : -v);
+
+   // Check if the order is very close to a negative half-integer.
+   const bool order_is_negative_half_integer = (order_is_negative && (((vv * 2U) - floor(vv * 2U)) < boost::math::tools::epsilon<T>()));
+
+   // The zero'th zero of Jv(x) for v < 0 is not defined
+   // unless the order is a negative integer.
+   if((m == 0) && (!order_is_negative_half_integer))
    {
-      return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but must be > 0 !", m, pol);
+      // For non-integer, negative order, requesting the zero'th zero raises a domain error.
+      return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", m, pol);
    }
 
+   // For nrgative half-integers, use the corresponding
+   // spherical Bessel function of positive half-integer order.
+   if(order_is_negative_half_integer)
+      return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol);
+
    // Set up the initial guess for the upcoming root-finding.
-   const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T>(v, m);
+   // If the order is a negative integer, then use the corresponding
+   // positive integer for the order.
+   const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T, Policy>(v, m, pol);
 
    // Select the maximum allowed iterations from the policy.
    boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
@@ -462,8 +492,8 @@
       boost::math::tools::newton_raphson_iterate(
          boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),
          guess_root,
-         (std::max)(T(guess_root - 0.3F), T(0)),
-         T(guess_root + 0.3F),
+         (std::max)(T(guess_root - 0.2F), T(0)),
+         T(guess_root + 0.2F),
          my_digits2,
          number_of_iterations);
 

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-02-23 07:50:45 EST (Sat, 23 Feb 2013)
@@ -252,15 +252,9 @@
           // There is special handling for negative order.
           if(v < 0)
           {
-            // Extract the absolute value of the order.
+            // Create the positive order and extract its positive floor integer part.
             const T vv(-v);
-            const T vv_floor = floor(vv);
-
-            // Check for pure negative integer order.
-            if((vv - vv_floor) < boost::math::tools::epsilon<T>())
-            {
-              return initial_guess(vv, m, pol);
-            }
+            const T vv_floor(floor(vv));
 
             // The to-be-found root is bracketed by the roots of the
             // Bessel function whose reflected, positive integer order
@@ -377,6 +371,25 @@
         }
 
         template<class T, class Policy>
+        class function_object_yv
+        {
+        public:
+          function_object_yv(const T& v,
+                             const Policy& pol) : my_v(v),
+                                                  my_pol(pol) { }
+
+          T operator()(const T& x) const
+          {
+            return boost::math::cyl_neumann(my_v, x, my_pol);
+          }
+
+        private:
+          const T my_v;
+          const Policy& my_pol;
+          const function_object_yv& operator=(const function_object_yv&);
+        };
+
+        template<class T, class Policy>
         class function_object_yv_and_yv_prime
         {
         public:
@@ -419,13 +432,93 @@
           const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&);
         };
 
-        template<class T>
-        T initial_guess(T v, int m)
+        template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
+
+        template<class T, class Policy>
+        T initial_guess(const T& v, const int m, const Policy& pol)
         {
+          BOOST_MATH_STD_USING // ADL of std names, needed for ceil.
+
           // Compute an estimate of the m'th root of cyl_neumann.
 
           T guess;
 
+          // There is special handling for negative order.
+          if(v < 0)
+          {
+            // Create the positive order and extract its positive floor and ceiling integer parts.
+            const T vv(-v);
+            const T vv_ceil (ceil (vv));
+            const T vv_floor(floor(vv));
+
+            // The to-be-found root is bracketed by the roots of the
+            // Bessel function whose reflected, positive integer order
+            // is greater than, but nearest to vv.
+
+            T root_hi;
+            T root_lo;
+
+            if(m == 1)
+            {
+              // The estimate of the first root for negative order is found using
+              // an adaptive range-searching algorithm.
+              if(T(vv - vv_floor) < 0.5F)
+              {
+                root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
+              }
+              else
+              {
+                root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
+              }
+
+              root_lo = T(root_hi - 0.1F);
+
+              const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0);
+
+              while((root_lo > boost::math::tools::epsilon<T>()))
+              {
+                const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0);
+
+                if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
+                {
+                  break;
+                }
+
+                root_hi = root_lo;
+
+                // Decrease the lower end of the bracket using an adaptive algorithm.
+                (root_lo > 0.5F) ? root_lo -= 0.5F : root_lo *= 0.85F;
+              }
+            }
+            else
+            {
+              if(T(vv - vv_floor) < 0.5F)
+              {
+                root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
+                root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
+              }
+              else
+              {
+                root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol);
+                root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
+              }
+            }
+
+            // Perform several steps of bisection iteration to refine the guess.
+            boost::uintmax_t number_of_iterations(12U);
+
+            // Do the bisection iteration.
+            const boost::math::tuple<T, T> guess_pair =
+               boost::math::tools::bisect(
+                  boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv<T, Policy>(v, pol),
+                  root_lo,
+                  root_hi,
+                  boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>,
+                  number_of_iterations);
+
+            return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
+          }
+
           if(m == 1U)
           {
             // Get the initial estimate of the first root.