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

Subject: [Boost-commit] svn:boost r83051 - in trunk/boost/math/special_functions: . detail
From: e_float_at_[hidden]
Date: 2013-02-20 15:40:19

Author: christopher_kormanyos
Date: 2013-02-20 15:40:18 EST (Wed, 20 Feb 2013)
New Revision: 83051
URL: http://svn.boost.org/trac/boost/changeset/83051

Log:
Added support for the zeros of Jv with negative order v.
Added additional checks for the zero'th zero of Jv.
Text files modified: 
   trunk/boost/math/special_functions/bessel.hpp                |    54 +++++--                                 
   trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp |   259 +++++++++++++++++++++++++-------------- 
   2 files changed, 199 insertions(+), 114 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-20 15:40:18 EST (Wed, 20 Feb 2013)
@@ -361,29 +361,47 @@
 
    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
 
-   const bool order_is_zero = ((v > -half_epsilon) && (v < +half_epsilon));
-
-   // Handle negative order or if the zero'th zero is requested.
+   // Handle the zero'th zero, if requested.
    // Return NaN if NaN is available or return 0 if NaN is not available.
-   if(v < 0) 
+   if (!(boost::math::isfinite)(v) )
    {
-      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be >= 0 !", v, pol);
+     return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
    }
-   else if (!(boost::math::isfinite)(v) )
+
+   if(m < 0)
    {
-     return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
+      // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error.
+      return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol);
    }
-   if(m <= 0)
+
+   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);
+
+   if(m == 0)
    {
-      // Special case: The zero'th zero of Jv(x) is zero for v != 0.
-      if((m == 0) && (!order_is_zero))
-         return T(0);
+      if(order_is_zero)
+      {
+         // The zero'th zero of J0(x) is not defined and requesting it raises a domain error.
+         return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of J0, but the rank must be > 0 !", m, pol);
+      }
+
+      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);
+      }
 
-      // Otherwise, the zero'th zero of Jv(x) is not defined and requesting it raises a domain error.
-      return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but must be > 0 !", m, pol);
+      // The zero'th zero does exist and its value is zero.
+      return T(0);
    }
+
    // Set up the initial guess for the upcoming root-finding.
-   const T guess_root = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess<T>(v, m);
+   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.
    boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
@@ -395,9 +413,9 @@
    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
    const T jvm =
       boost::math::tools::newton_raphson_iterate(
-         boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object<T, Policy>(v, pol),
+         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,
-         T(guess_root - 0.3F),
+         (std::max)(T(guess_root - 0.3F), T(0)),
          T(guess_root + 0.3F),
          my_digits2,
          number_of_iterations);
@@ -442,9 +460,9 @@
    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
    const T yvm =
       boost::math::tools::newton_raphson_iterate(
-         boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object<T, Policy>(v, pol),
+         boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),
          guess_root,
-         T(guess_root - 0.3F),
+         (std::max)(T(guess_root - 0.3F), T(0)),
          T(guess_root + 0.3F),
          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-20 15:40:18 EST (Wed, 20 Feb 2013)
@@ -25,12 +25,6 @@
   namespace boost { namespace math {
   namespace detail
   {
-    // Forward declarations of the needed Bessel function implementations.
-    template <class T, class Policy>
-    T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol);
-    template <class T, class Policy>
-    T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol);
-
     namespace bessel_zero
     {
       template<class T>
@@ -83,7 +77,7 @@
         }
 
       private:
-        equation_as_9_3_39_and_its_derivative& operator = (const equation_as_9_3_39_and_its_derivative&);
+        const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&);
         const T zeta;
       };
 
@@ -183,20 +177,142 @@
                           * v_pow_minus_two_thirds + T(1));
         }
 
-        template<class T>
-        T initial_guess(const T& v, int m)
+        template<class T, class Policy>
+        class function_object_jv
+        {
+        public:
+          function_object_jv(const T& v,
+                             const Policy& pol) : my_v(v),
+                                                  my_pol(pol) { }
+
+          T operator()(const T& x) const
+          {
+            return boost::math::cyl_bessel_j(my_v, x, my_pol);
+          }
+
+        private:
+          const T my_v;
+          const Policy& my_pol;
+          const function_object_jv& operator=(const function_object_jv&);
+        };
+
+        template<class T, class Policy>
+        class function_object_jv_and_jv_prime
+        {
+        public:
+          function_object_jv_and_jv_prime(const T& v,
+                                          const bool order_is_zero,
+                                          const Policy& pol) : my_v(v),
+                                                               my_order_is_zero(order_is_zero),
+                                                               my_pol(pol) { }
+
+          boost::math::tuple<T, T> operator()(const T& x) const
+          {
+            // Obtain Jv(x) and Jv'(x).
+            // Chris's original code called the Bessel function implementation layer direct, 
+            // but that circumvented optimizations for integer-orders.  Call the documented
+            // top level functions instead, and let them sort out which implementation to use.
+            T j_v;
+            T j_v_prime;
+
+            if(my_order_is_zero)
+            {
+              j_v       =  boost::math::cyl_bessel_j(0, x, my_pol);
+              j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol);
+            }
+            else
+            {
+                      j_v       = boost::math::cyl_bessel_j(  my_v,      x, my_pol);
+              const T j_v_m1     (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol));
+                      j_v_prime = j_v_m1 - ((my_v * j_v) / x);
+            }
+
+            // Return a tuple containing both Jv(x) and Jv'(x).
+            return boost::math::make_tuple(j_v, j_v_prime);
+          }
+
+        private:
+          const T my_v;
+          const bool my_order_is_zero;
+          const Policy& my_pol;
+          const function_object_jv_and_jv_prime& operator=(const function_object_jv_and_jv_prime&);
+        };
+
+        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 floor.
+
           // Compute an estimate of the m'th root of cyl_bessel_j.
 
           T guess;
 
-          if(m == 0U)
+          // There is special handling for negative order.
+          if(v < 0)
           {
-            // Requesting an estimate of the zero'th root is an error.
-            // Return zero.
-            guess = T(0);
+            // Extract the absolute value of the order.
+            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);
+            }
+
+            // The to-be-found root is bracketed by the roots of the
+            // Bessel function whose reflected, positive integer order
+            // is less than, but nearest to vv.
+
+            T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol);
+            T root_lo;
+
+            if(m == 1)
+            {
+              // The estimate of the first root for negative order is found using
+              // an adaptive range-searching algorithm.
+              root_lo = T(root_hi - 0.1F);
+
+              const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0);
+
+              while((root_lo > boost::math::tools::epsilon<T>()))
+              {
+                const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(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
+            {
+              root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, 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_bessel_j_zero_detail::function_object_jv<T, Policy>(v, pol),
+                  root_lo,
+                  root_hi,
+                  boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>,
+                  number_of_iterations);
+
+            return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
           }
-          else if(m == 1U)
+
+          if(m == 1U)
           {
             // Get the initial estimate of the first root.
 
@@ -242,14 +358,31 @@
 
           return guess;
         }
+      } // namespace cyl_bessel_j_zero_detail
+
+      namespace cyl_neumann_zero_detail
+      {
+        template<class T>
+        T equation_nist_10_21_40_b(const T& v)
+        {
+          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))
+                          * v_pow_minus_two_thirds + T(0.01198))
+                          * v_pow_minus_two_thirds + T(0.260351))
+                          * v_pow_minus_two_thirds + T(0.9315768))
+                          * v_pow_minus_two_thirds + T(1));
+        }
 
         template<class T, class Policy>
-        class function_object
+        class function_object_yv_and_yv_prime
         {
         public:
-          function_object(const T& v,
-                          const Policy& pol) : my_v(v),
-                                               my_pol(pol) { }
+          function_object_yv_and_yv_prime(const T& v,
+                                          const Policy& pol) : my_v(v),
+                                                               my_pol(pol) { }
 
           boost::math::tuple<T, T> operator()(const T& x) const
           {
@@ -257,51 +390,34 @@
 
             const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon));
 
-            // Obtain Jv(x) and Jv'(x).
+            // Obtain Yv(x) and Yv'(x).
             // Chris's original code called the Bessel function implementation layer direct, 
             // but that circumvented optimizations for integer-orders.  Call the documented
             // top level functions instead, and let them sort out which implementation to use.
-            T j_v;
-            T j_v_prime;
+            T y_v;
+            T y_v_prime;
 
             if(order_is_zero)
             {
-              j_v       =  boost::math::cyl_bessel_j(0, x, my_pol);
-              j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol);
+              y_v       =  boost::math::cyl_neumann(0, x, my_pol);
+              y_v_prime = -boost::math::cyl_neumann(1, x, my_pol);
             }
             else
             {
-                      j_v       = boost::math::cyl_bessel_j(  my_v,      x, my_pol);
-              const T j_v_m1     (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol));
-                      j_v_prime = j_v_m1 - ((my_v * j_v) / x);
+                      y_v       = boost::math::cyl_neumann(  my_v,      x, my_pol);
+              const T y_v_m1     (boost::math::cyl_neumann(T(my_v - 1), x, my_pol));
+                      y_v_prime = y_v_m1 - ((my_v * y_v) / x);
             }
 
-            // Return a tuple containing both Jv(x) and Jv'(x).
-            return boost::math::make_tuple(j_v, j_v_prime);
+            // Return a tuple containing both Yv(x) and Yv'(x).
+            return boost::math::make_tuple(y_v, y_v_prime);
           }
 
         private:
           const T my_v;
           const Policy& my_pol;
-          function_object& operator = (const function_object&);
+          const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&);
         };
-      } // namespace cyl_bessel_j_zero_detail
-
-      namespace cyl_neumann_zero_detail
-      {
-        template<class T>
-        T equation_nist_10_21_40_b(const T& v)
-        {
-          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))
-                          * v_pow_minus_two_thirds + T(0.01198))
-                          * v_pow_minus_two_thirds + T(0.260351))
-                          * v_pow_minus_two_thirds + T(0.9315768))
-                          * v_pow_minus_two_thirds + T(1));
-        }
 
         template<class T>
         T initial_guess(T v, int m)
@@ -310,13 +426,7 @@
 
           T guess;
 
-          if(m == 0U)
-          {
-            // Requesting an estimate of the zero'th root is an error.
-            // Return zero.
-            guess = T(0);
-          }
-          else if(m == 1U)
+          if(m == 1U)
           {
             // Get the initial estimate of the first root.
 
@@ -362,49 +472,6 @@
 
           return guess;
         }
-
-        template<class T, class Policy>
-        class function_object
-        {
-        public:
-          function_object(const T& v,
-                          const Policy& pol) : my_v(v),
-                                               my_pol(pol) { }
-
-          boost::math::tuple<T, T> operator()(const T& x) const
-          {
-            const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
-
-            const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon));
-
-            // Obtain Yv(x) and Yv'(x).
-            // Chris's original code called the Bessel function implementation layer direct, 
-            // but that circumvented optimizations for integer-orders.  Call the documented
-            // top level functions instead, and let them sort out which implementation to use.
-            T y_v;
-            T y_v_prime;
-
-            if(order_is_zero)
-            {
-              y_v       =  boost::math::cyl_neumann(0, x, my_pol);
-              y_v_prime = -boost::math::cyl_neumann(1, x, my_pol);
-            }
-            else
-            {
-                      y_v       = boost::math::cyl_neumann(  my_v,      x, my_pol);
-              const T y_v_m1     (boost::math::cyl_neumann(T(my_v - 1), x, my_pol));
-                      y_v_prime = y_v_m1 - ((my_v * y_v) / x);
-            }
-
-            // Return a tuple containing both Yv(x) and Yv'(x).
-            return boost::math::make_tuple(y_v, y_v_prime);
-          }
-
-        private:
-          const T my_v;
-          const Policy& my_pol;
-          function_object* operator = (const function_object&);
-        };
       } // namespace cyl_neumann_zero_detail
     } // namespace bessel_zero
   } } } // namespace boost::math::detail