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

Subject: [Boost-commit] svn:boost r82599 - in trunk/boost/math/special_functions: . detail
From: e_float_at_[hidden]
Date: 2013-01-24 13:59:15

Author: christopher_kormanyos
Date: 2013-01-24 13:59:14 EST (Thu, 24 Jan 2013)
New Revision: 82599
URL: http://svn.boost.org/trac/boost/changeset/82599

Log:
Added support for the zeros of cylindrical Bessel functions.
Added:
   trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp   (contents, props changed)
   trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp   (contents, props changed)
Text files modified: 
   trunk/boost/math/special_functions/bessel.hpp   |   194 +++++++++++++++++++++++++++++++++++++++ 
   trunk/boost/math/special_functions/math_fwd.hpp |    49 ++++++++++                              
   2 files changed, 241 insertions(+), 2 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-24 13:59:14 EST (Thu, 24 Jan 2013)
@@ -17,6 +17,7 @@
 #include <boost/math/special_functions/detail/bessel_jy.hpp>
 #include <boost/math/special_functions/detail/bessel_jn.hpp>
 #include <boost/math/special_functions/detail/bessel_yn.hpp>
+#include <boost/math/special_functions/detail/bessel_jy_zero.hpp>
 #include <boost/math/special_functions/detail/bessel_ik.hpp>
 #include <boost/math/special_functions/detail/bessel_i0.hpp>
 #include <boost/math/special_functions/detail/bessel_i1.hpp>
@@ -30,6 +31,7 @@
 #include <boost/math/tools/rational.hpp>
 #include <boost/math/tools/promotion.hpp>
 #include <boost/math/tools/series.hpp>
+#include <boost/math/tools/roots.hpp>
 
 namespace boost{ namespace math{
 
@@ -301,7 +303,7 @@
 
    if(floor(v) == 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)))
+      if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
       {
          T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
          if((v < 0) && (itrunc(v, pol) & 1))
@@ -330,7 +332,7 @@
    BOOST_MATH_INSTRUMENT_VARIABLE(v);
    BOOST_MATH_INSTRUMENT_VARIABLE(x);
 
-   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)))
+   if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
    {
       T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
       if((v < 0) && (v & 1))
@@ -367,6 +369,120 @@
    return result * tx;
 }
 
+template <class T, class Policy>
+inline T cyl_bessel_j_zero_imp(T v, unsigned m, const Policy& pol)
+{
+   BOOST_MATH_STD_USING // ADL of std names, needed for log.
+
+   // 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 < 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.
+   const T guess_root = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess<T>(v, m);
+
+   // Select the maximum allowed iterations, being at least 24.
+   boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits<T>::digits10));
+
+   // Select the desired number of binary digits of precision.
+   // Account for the radix of number representations having non-two radix!
+   const int my_digits2 = int(float(std::numeric_limits<T>::digits)
+                              * (  log(float(std::numeric_limits<T>::radix))
+                                 / log(2.0F)));
+
+   // 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),
+         guess_root,
+         T(guess_root - 0.3F),
+         T(guess_root + 0.3F),
+         my_digits2,
+         number_of_iterations);
+
+   static_cast<void>(number_of_iterations);
+
+   return jvm;
+}
+
+template <class output_iterator, class T, class Policy>
+inline void cyl_bessel_j_zero_imp(output_iterator out_it,
+                                  T v,
+                                  unsigned number_of_zeros,
+                                  unsigned start_index,
+                                  const Policy& pol)
+{
+   if(number_of_zeros > 0U)
+   {
+      const output_iterator end_it(out_it + number_of_zeros);
+
+      while(out_it != end_it)
+      {
+         *out_it = boost::math::detail::cyl_bessel_j_zero_imp(v, start_index, pol);
+         ++start_index;
+         ++out_it;
+      }
+   }
+}
+
+template <class T, class Policy>
+inline T cyl_neumann_zero_imp(T v, unsigned m, const Policy& pol)
+{
+   BOOST_MATH_STD_USING // ADL of std names, needed for log.
+
+   // 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 < 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.
+   const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T>(v, m);
+
+   // Select the maximum allowed iterations, being at least 24.
+   boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits<T>::digits10));
+
+   // Select the desired number of binary digits of precision.
+   // Account for the radix of number representations having non-two radix!
+   const int my_digits2 = int(float(std::numeric_limits<T>::digits)
+                              * (  log(float(std::numeric_limits<T>::radix))
+                                 / log(2.0F)));
+
+   // 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),
+         guess_root,
+         T(guess_root - 0.3F),
+         T(guess_root + 0.3F),
+         my_digits2,
+         number_of_iterations);
+
+   static_cast<void>(number_of_iterations);
+
+   return yvm;
+}
+
+template <class output_iterator, class T, class Policy>
+inline void cyl_neumann_zero_imp(output_iterator out_it,
+                                 T v,
+                                 unsigned number_of_zeros,
+                                 unsigned start_index,
+                                 const Policy& pol)
+{
+   if(number_of_zeros > 0U)
+   {
+      const output_iterator end_it(out_it + number_of_zeros);
+
+      while(out_it != end_it)
+      {
+         *out_it = boost::math::detail::cyl_neumann_zero_imp(v, start_index, pol);
+         ++start_index;
+         ++out_it;
+      }
+   }
+}
+
 } // namespace detail
 
 template <class T1, class T2, class Policy>
@@ -462,6 +578,80 @@
    return sph_neumann(v, x, policies::policy<>());
 }
 
+template <class T, class Policy>
+inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_zero(T v, unsigned m, const Policy& pol)
+{
+   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;
+   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<>());
+}
+
+template <class output_iterator, class T, class Policy>
+inline void cyl_bessel_j_zero(output_iterator out_it,
+                              T v,
+                              unsigned number_of_zeros,
+                              unsigned start_index,
+                              const Policy& pol)
+{
+   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;
+   detail::cyl_bessel_j_zero_imp<output_iterator, value_type, Policy>(out_it, v, number_of_zeros, start_index, pol);
+}
+
+template <class output_iterator, class T>
+inline void cyl_bessel_j_zero(output_iterator out_it,
+                              T v,
+                              unsigned number_of_zeros,
+                              unsigned start_index)
+{
+   detail::cyl_bessel_j_zero_imp(out_it, v, number_of_zeros, start_index, policies::policy<>());
+}
+
+template <class T, class Policy>
+inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zero(T v, unsigned m, const Policy& pol)
+{
+   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;
+   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<>());
+}
+
+template <class output_iterator, class T, class Policy>
+inline void cyl_neumann_zero(output_iterator out_it,
+                             T v,
+                             unsigned number_of_zeros,
+                             unsigned start_index,
+                             const Policy& pol)
+{
+   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;
+   detail::cyl_neumann_zero_imp<output_iterator, value_type, Policy>(out_it, v, number_of_zeros, start_index, pol);
+}
+
+template <class output_iterator, class T>
+inline void cyl_neumann_zero(output_iterator out_it,
+                             T v,
+                             unsigned number_of_zeros,
+                             unsigned start_index)
+{
+   detail::cyl_neumann_zero_imp(out_it, v, number_of_zeros, start_index, policies::policy<>());
+}
+
 } // namespace math
 } // namespace boost
 

Added: trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp	2013-01-24 13:59:14 EST (Thu, 24 Jan 2013)
@@ -0,0 +1,100 @@
+//  Copyright (c) 2013 Christopher Kormanyos
+//  Use, modification and distribution are subject to 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)
+//
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+//
+// This header contains implementation details for estimating the zeros
+// of the Airy functions airy_ai and airy_bi on the negative real axis.
+//
+#ifndef _AIRY_AI_BI_ZERO_2013_01_20_HPP_
+  #define _AIRY_AI_BI_ZERO_2013_01_20_HPP_
+
+  #include <boost/math/constants/constants.hpp>
+
+  namespace boost { namespace math { namespace detail {
+  namespace airy_zero
+  {
+    template<class T>
+    T equation_as_10_4_105(const T& z)
+    {
+      BOOST_MATH_STD_USING // ADL of std names, needed for pow.
+
+      const T one_over_z         = T(1) / z;
+      const T one_over_z_squared = one_over_z * one_over_z;
+
+      const T z_pow_two_thirds(pow(z, T(2) / 3U));
+
+      // Implement the top line of Eq. 10.4.105.
+      const T fz(z_pow_two_thirds * (((((                     + (T(162375596875ULL) / 334430208UL)
+                                         * one_over_z_squared - (   T(108056875ULL) /   6967296UL))
+                                         * one_over_z_squared + (       T(77125UL)  /     82944UL))
+                                         * one_over_z_squared - (           T(5U)   /        36U))
+                                         * one_over_z_squared + (           T(5U)   /        48U))
+                                         * one_over_z_squared + (1)));
+
+      return fz;
+    }
+
+    namespace airy_ai_zero_detail
+    {
+      template<class T>
+      T initial_guess(const unsigned m)
+      {
+        switch(m)
+        {
+          case 1U:  return T(-2.33810741045976703849);
+          case 2U:  return T(-4.08794944413097061664);
+          case 3U:  return T(-5.52055982809555105913);
+          case 4U:  return T(-6.78670809007175899878);
+          case 5U:  return T(-7.94413358712085312314);
+          case 6U:  return T(-9.02265085334098038016);
+          case 7U:  return T(-10.0401743415580859306);
+          case 8U:  return T(-11.0085243037332628932);
+          case 9U:  return T(-11.9360155632362625170);
+          case 10U: return T(-12.8287767528657572004);
+          default:
+          {
+            const T t = ((boost::math::constants::pi<T>() * 3U) * T((T(m) * 4U) - 1)) / 8U;
+
+            return -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+          }
+        }
+      }
+    } // namespace airy_ai_zero_detail
+
+    namespace airy_bi_zero_detail
+    {
+      template<class T>
+      T initial_guess(const unsigned m)
+      {
+        switch(m)
+        {
+          case 1U:  return T(-1.17371322270912792492);
+          case 2U:  return T(-3.27109330283635271568);
+          case 3U:  return T(-4.83073784166201593267);
+          case 4U:  return T(-6.16985212831025125983);
+          case 5U:  return T(-7.37676207936776371360);
+          case 6U:  return T(-8.49194884650938801345);
+          case 7U:  return T(-9.53819437934623888663);
+          case 8U:  return T(-10.5299135067053579244);
+          case 9U:  return T(-11.4769535512787794379);
+          case 10U: return T(-12.3864171385827387456);
+          default:
+          {
+            const T t = ((boost::math::constants::pi<T>() * 3U) * T((T(m) * 4U) - 3)) / 8U;
+
+            return -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+          }
+        }
+      }
+    } // namespace airy_bi_zero_detail
+  } // namespace airy_zero
+  } // namespace detail
+  } // namespace math
+  } // namespaces boost
+
+#endif // _AIRY_AI_BI_ZERO_2013_01_20_HPP_

Added: trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp	2013-01-24 13:59:14 EST (Thu, 24 Jan 2013)
@@ -0,0 +1,367 @@
+//  Copyright (c) 2013 Christopher Kormanyos
+//  Use, modification and distribution are subject to 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)
+//
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+//
+// This header contains implementation details for estimating the zeros
+// of cylindrical Bessel and Neumann functions on the positive real axis.
+// Various methods are used to estimate the roots. These include
+// empirical curve fitting and McMahon's asymptotic approximation
+// for small order, and uniform asymptotic expansion for large order.
+//
+#ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_
+  #define _BESSEL_JY_ZERO_2013_01_18_HPP_
+
+  #include <algorithm>
+  #include <boost/math/constants/constants.hpp>
+  #include <boost/math/special_functions/math_fwd.hpp>
+  #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
+
+  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>
+      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.
+        // The order m has been used to create the input parameter a.
+
+        // This is Eq. 10.21.19 in the NIST Handbook.
+        const T mu                  = (v * v) * 4U;
+        const T mu_minus_one        = mu - T(1);
+        const T eight_a_inv         = T(1) / (a * 8U);
+        const T eight_a_inv_squared = eight_a_inv * eight_a_inv;
+
+        const T term3 = ((mu_minus_one *  4U) *     ((mu *    7U) -     T(31U) )) / 3U;
+        const T term5 = ((mu_minus_one * 32U) *   ((((mu *   83U) -    T(982U) ) * mu) +    T(3779U) )) / 15U;
+        const T term7 = ((mu_minus_one * 64U) * ((((((mu * 6949U) - T(153855UL)) * mu) + T(1585743UL)) * mu) - T(6277237UL))) / 105U;
+
+        return a + ((((                      - term7
+                       * eight_a_inv_squared - term5)
+                       * eight_a_inv_squared - term3)
+                       * eight_a_inv_squared - mu_minus_one)
+                       * eight_a_inv);
+      }
+
+      template<typename T>
+      class equation_as_9_3_39_and_its_derivative
+      {
+      public:
+        equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }
+
+        boost::math::tuple<T, T> operator()(const T& z) const
+        {
+          BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt.
+
+          // Return the function of zeta that is implicitly defined
+          // in A&S Eq. 9.3.39 as a function of z. The function is
+          // returned along with its derivative with respect to z.
+
+          const T zsq_minus_one_sqrt = sqrt((z * z) - T(1));
+
+          const T the_function(
+              zsq_minus_one_sqrt
+            - (  acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta)))));
+
+          const T its_derivative(zsq_minus_one_sqrt / z);
+
+          return boost::math::tuple<T, T>(the_function, its_derivative);
+        }
+
+      private:
+        const T zeta;
+      };
+
+      template<class T>
+      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.
+
+        // 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.
+        // Here, v is larger than about 1.2. The estimate is computed
+        // from Abramowitz and Stegun Eqs. 9.5.22 and 9.5.26, page 371.
+        //
+        // The inversion of z as a function of zeta is mentioned in the text
+        // following A&S Eq. 9.5.26. Here, we accomplish the inversion by
+        // performing a Taylor expansion of Eq. 9.3.39 for large z to order 2
+        // and solving the resulting quadratic equation, thereby taking
+        // the positive root of the quadratic.
+        //
+        // In other words: (2/3)(-zeta)^(3/2) approx = z + 1/(2z) - pi/2.
+        // This leads to: z^2 - [(2/3)(-zeta)^(3/2) + pi/2]z + 1/2 = 0.
+        // With this initial estimate, Newton-Raphson iteration is used
+        // to refine the value of the estimate of the root of z
+        // as a function of zeta.
+
+        // 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_sqrt = sqrt(zeta);
+
+        // Set up a quadratic equation based on the Taylor series
+        // expansion mentioned above.
+        const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>());
+
+        // Solve the quadratic equation, taking the positive root.
+        const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U;
+
+        // Establish the range, the digits, and the iteration limit
+        // for the upcoming root-finding.
+        const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1));
+        const T range_zmax = z_estimate + T(1);
+
+        const int digits2_of_t = int(float(std::numeric_limits<T>::digits)
+                                     * (  log(float(std::numeric_limits<T>::radix))
+                                        / log(float(2))));
+
+        const int digits2_for_root = (std::min)(digits2_of_t, std::numeric_limits<double>::digits);
+
+        boost::uintmax_t iteration_count = boost::uintmax_t(std::numeric_limits<T>::digits10 * 2);
+
+        // Calculate the root of z as a function of zeta.
+        const T z = boost::math::tools::newton_raphson_iterate(
+          boost::math::detail::bessel_zero::equation_as_9_3_39_and_its_derivative<T>(zeta),
+          z_estimate,
+          range_zmin,
+          range_zmax,
+          digits2_for_root,
+          iteration_count);
+
+        static_cast<void>(iteration_count);
+
+        // Continue with the implementation of A&S Eq. 9.3.39.
+        const T zsq_minus_one      = (z * z) - T(1);
+        const T zsq_minus_one_sqrt = sqrt(zsq_minus_one);
+
+        // This is A&S Eq. 9.3.42.
+        const T b0_term_5_24 = T(5) / ((zsq_minus_one * zsq_minus_one_sqrt) * 24U);
+        const T b0_term_1_8  = T(1) / ( zsq_minus_one_sqrt * 8U);
+        const T b0_term_5_48 = T(5) / ((zeta * zeta) * 48U);
+
+        const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt);
+
+        // This is the second line of A&S Eq. 9.5.26 for f_k with k = 1.
+        const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt;
+
+        // This is A&S Eq. 9.5.22 expanded to k = 1 (i.e., one term in the series).
+        return (v * z) + (f1 / v);
+      }
+
+      namespace cyl_bessel_j_zero_detail
+      {
+        template<class T>
+        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);
+
+          return v * (((((                         + T(0.043)
+                          * v_pow_minus_two_thirds - T(0.0908))
+                          * v_pow_minus_two_thirds - T(0.00397))
+                          * v_pow_minus_two_thirds + T(1.033150))
+                          * v_pow_minus_two_thirds + T(1.8557571))
+                          * v_pow_minus_two_thirds + T(1));
+        }
+
+        template<class T>
+        T initial_guess(const T& v, unsigned m)
+        {
+          // Compute an estimate of the m'th root of cyl_bessel_j.
+
+          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)
+          {
+            // Get the initial estimate of the first root.
+
+            if(v < T(1.2))
+            {
+              // For small v, use the results of empirical curve fitting.
+              // Mathematica(R) session for the coefficients:
+              //  Table[{n, BesselJZero[n, 1]}, {n, 0, 12/10, 1/10}]
+              //  N[%, 20]
+              //  Fit[%, {n^0, n^1, n^2, n^3, n^4}, n]
+              guess = (((    - T(0.0131733516699981)
+                         * v + T(0.062016703196207))
+                         * v - T(0.163263189752310))
+                         * v + T(1.5412938080110762))
+                         * v + T(2.40485167005651374);
+            }
+            else
+            {
+              // For larger v, use the first line of Eqs. 10.21.40 in the NIST Handbook.
+              guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v);
+            }
+          }
+          else
+          {
+            if(v < T(1.2))
+            {
+              // Use Eq. 10.21.19 in the NIST Handbook.
+              const T a = ((v + T(m * 2U)) - T(0.5)) * boost::math::constants::half_pi<T>();
+
+              guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
+            }
+            else
+            {
+              // Get an estimate of the m'th root of airy_ai.
+              const T airy_ai_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m);
+
+              // Use Eq. 9.5.26 in the A&S Handbook.
+              guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root);
+            }
+          }
+
+          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
+          {
+            // 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));
+            const T j_v_m1   (boost::math::detail::cyl_bessel_j_imp(T(my_v - 1), x, boost::math::detail::bessel_no_int_tag(), my_pol));
+            const T 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 Policy& my_pol;
+        };
+      } // namespace cyl_bessel_j_zero_detail
+
+      namespace cyl_neumann_zero_detail
+      {
+        template<class T>
+        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);
+
+          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, unsigned m)
+        {
+          // Compute an estimate of the m'th root of cyl_neumann.
+
+          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)
+          {
+            // Get the initial estimate of the first root.
+
+            if(v < T(1.2))
+            {
+              // For small v, use the results of empirical curve fitting.
+              // Mathematica(R) session for the coefficients:
+              //  Table[{n, BesselYZero[n, 1]}, {n, 0, 12/10, 1/10}]
+              //  N[%, 20]
+              //  Fit[%, {n^0, n^1, n^2, n^3, n^4}, n]
+              guess = (((    - T(0.0288448856660199)
+                         * v + T(0.1157780677827211))
+                         * v - T(0.2248646596163634))
+                         * v + T(1.4414721779748667))
+                         * v + T(0.89366124647143352);
+            }
+            else
+            {
+              // For larger v, use the second line of Eqs. 10.21.40 in the NIST Handbook.
+              guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v);
+            }
+          }
+          else
+          {
+            if(v < T(1.2))
+            {
+              // Use Eq. 10.21.19 in the NIST Handbook.
+              const T a = ((v + T(m * 2U)) - T(1.5)) * boost::math::constants::half_pi<T>();
+
+              guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
+            }
+            else
+            {
+              // Get an estimate of the m'th root of airy_bi.
+              const T airy_bi_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m);
+
+              // Use Eq. 9.5.26 in the A&S Handbook.
+              guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root);
+            }
+          }
+
+          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
+          {
+            // 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));
+            const T y_v_m1   (boost::math::detail::cyl_neumann_imp(T(my_v - 1), x, boost::math::detail::bessel_no_int_tag(), my_pol));
+            const T 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;
+        };
+      } // namespace cyl_neumann_zero_detail
+    } // namespace bessel_zero
+  } } } // namespace boost::math::detail
+
+#endif // _BESSEL_JY_ZERO_2013_01_18_HPP_

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-24 13:59:14 EST (Thu, 24 Jan 2013)
@@ -614,9 +614,36 @@
    template <class T>
    typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x);
 
+   template <class T, class Policy>
+   typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_zero(T v, unsigned m, const Policy& pol);
+
+   template <class T>
+   typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, unsigned m);
+
+   template <class T1, class T2, class Policy>
+   std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol);
+
+   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);
+
+   template <class T, class Policy>
+   typename detail::bessel_traits<T, T, Policy>::result_type cyl_neuman_zero(T v, unsigned m, const Policy& pol);
+
+   template <class T>
+   typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neuman_zero(T v, unsigned m);
+
    template <class T1, class T2, class Policy>
    std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol);
 
+   template <class output_iterator, class T>
+   inline void cyl_neuman_zero_imp(output_iterator out_it,
+                                   T v,
+                                   std::size_t number_of_zeros,
+                                   unsigned start_index);
+
    template <class T1, class T2>
    std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x);
 
@@ -1139,6 +1166,28 @@
    sph_neumann(unsigned v, T x){ return boost::math::sph_neumann(v, x, Policy()); }\
 \
    template <class T>\
+   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); }\
 \
    template <class T>\