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

Subject: [Boost-commit] svn:boost r82682 - in trunk/boost/math/special_functions: . detail
From: e_float_at_[hidden]
Date: 2013-02-06 20:24:01

Author: christopher_kormanyos
Date: 2013-02-01 14:39:36 EST (Fri, 01 Feb 2013)
New Revision: 82682
URL: http://svn.boost.org/trac/boost/changeset/82682

Log:
Changed the order of the input parameters for Bessel zeros.
Improved the algorithms for Bessel zeros.
Added and improved the Airy zeros.

Text files modified: 
   trunk/boost/math/special_functions/airy.hpp                   |   171 ++++++++++++++++++++++++++++++++++++    
   trunk/boost/math/special_functions/bessel.hpp                 |    24 ++--                                    
   trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp |   188 +++++++++++++++++++++++++++------------ 
   trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp  |    72 +++++++++-----                          
   trunk/boost/math/special_functions/math_fwd.hpp               |    72 +++++++++-----                          
   5 files changed, 403 insertions(+), 124 deletions(-)

Modified: trunk/boost/math/special_functions/airy.hpp
==============================================================================
--- trunk/boost/math/special_functions/airy.hpp	(original)
+++ trunk/boost/math/special_functions/airy.hpp	2013-02-01 14:39:36 EST (Fri, 01 Feb 2013)
@@ -9,6 +9,8 @@
 
 #include <boost/math/special_functions/bessel.hpp>
 #include <boost/math/special_functions/cbrt.hpp>
+#include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
+#include <boost/math/tools/roots.hpp>
 
 namespace boost{ namespace math{
 
@@ -151,6 +153,131 @@
    }
 }
 
+template <class T>
+T airy_ai_zero_imp(T dummy, unsigned m)
+{
+   BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
+
+   static_cast<void>(dummy);
+
+   // Handle cases when the zero'th zero is requested.
+   // Return NaN if NaN is available or return 0 if NaN is not available.
+   if(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::airy_zero::airy_ai_zero_detail::initial_guess<T>(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)));
+
+   // Use a dynamic tolerance because the roots get closer the higher m gets.
+   T tolerance;
+
+   if     (m <=   10U) { tolerance = T(0.3F); }
+   else if(m <=  100U) { tolerance = T(0.1F); }
+   else if(m <= 1000U) { tolerance = T(0.05F); }
+   else                { tolerance = T(1) / sqrt(T(m)); }
+
+   // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
+   const T am =
+      boost::math::tools::newton_raphson_iterate(
+         boost::math::detail::airy_zero::airy_ai_zero_detail::function_object<T, policies::policy<> >(policies::policy<>()),
+         guess_root,
+         T(guess_root - tolerance),
+         T(guess_root + tolerance),
+         my_digits2,
+         number_of_iterations);
+
+   static_cast<void>(number_of_iterations);
+
+   return am;
+}
+
+template <class output_iterator, class T>
+void airy_ai_zero_imp(T dummy,
+                      unsigned number_of_zeros,
+                      unsigned start_index,
+                      output_iterator out_it)
+{
+   static_cast<void>(dummy);
+
+   for(unsigned i = 0; i < number_of_zeros; ++i)
+   {
+      *out_it = boost::math::detail::airy_ai_zero_imp<T>(T(), start_index + i);
+      ++out_it;
+   }
+}
+
+template <class T>
+T airy_bi_zero_imp(T dummy, unsigned m)
+{
+   BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
+
+   static_cast<void>(dummy);
+
+   // Handle cases when the zero'th zero is requested.
+   // Return NaN if NaN is available or return 0 if NaN is not available.
+   if(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::airy_zero::airy_bi_zero_detail::initial_guess<T>(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)));
+
+   // Use a dynamic tolerance because the roots get closer the higher m gets.
+   // Use a dynamic tolerance because the roots get closer the higher m gets.
+   T tolerance;
+
+   if     (m <=   10U) { tolerance = T(0.3F); }
+   else if(m <=  100U) { tolerance = T(0.1F); }
+   else if(m <= 1000U) { tolerance = T(0.05F); }
+   else                { tolerance = T(1) / sqrt(T(m)); }
+
+   // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
+   const T bm =
+      boost::math::tools::newton_raphson_iterate(
+         boost::math::detail::airy_zero::airy_bi_zero_detail::function_object<T, policies::policy<> >(policies::policy<>()),
+         guess_root,
+         T(guess_root - tolerance),
+         T(guess_root + tolerance),
+         my_digits2,
+         number_of_iterations);
+
+   static_cast<void>(number_of_iterations);
+
+   return bm;
+}
+
+template <class output_iterator, class T>
+void airy_bi_zero_imp(T dummy,
+                      unsigned number_of_zeros,
+                      unsigned start_index,
+                      output_iterator out_it)
+{
+   static_cast<void>(dummy);
+
+   for(unsigned i = 0; i < number_of_zeros; ++i)
+   {
+      *out_it = boost::math::detail::airy_bi_zero_imp<T>(T(), start_index + i);
+      ++out_it;
+   }
+}
+
 } // namespace detail
 
 template <class T, class Policy>
@@ -241,6 +368,50 @@
    return airy_bi_prime(x, policies::policy<>());
 }
 
+template <class T>
+inline typename tools::promote_args<T>::type airy_ai_zero(T dummy, unsigned m)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename tools::promote_args<T>::type result_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<result_type>::is_integer, "Airy dummy parameter must be a floating-point type.");
+   return policies::checked_narrowing_cast<result_type, policies::policy<> >(detail::airy_ai_zero_imp<result_type>(dummy, m), "boost::math::airy_ai_zero<%1%>(%1%)");
+}
+
+template <class output_iterator, class T>
+inline void airy_ai_zero(T dummy,
+                         unsigned number_of_zeros,
+                         unsigned start_index,
+                         output_iterator out_it)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   static_cast<void>(dummy);
+   typedef typename tools::promote_args<T>::type result_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<result_type>::is_integer, "Airy dummy parameter must be a floating-point type.");
+   boost::math::detail::airy_ai_zero_imp<output_iterator, result_type>(result_type(), number_of_zeros, start_index, out_it);
+}
+
+template <class T>
+inline typename tools::promote_args<T>::type airy_bi_zero(T dummy, unsigned m)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename tools::promote_args<T>::type result_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<result_type>::is_integer, "Airy dummy parameter must be a floating-point type.");
+   return policies::checked_narrowing_cast<result_type, policies::policy<> >(detail::airy_bi_zero_imp<result_type>(result_type(), m), "boost::math::airy_bi_zero<%1%>(%1%)");
+}
+
+template <class output_iterator, class T>
+inline void airy_bi_zero(T dummy,
+                         unsigned number_of_zeros,
+                         unsigned start_index,
+                         output_iterator out_it)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   static_cast<void>(dummy);
+   typedef typename tools::promote_args<T>::type result_type;
+   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<result_type>::is_integer, "Airy dummy parameter must be a floating-point type.");
+   boost::math::detail::airy_bi_zero_imp<output_iterator, result_type>(result_type(), number_of_zeros, start_index, out_it);
+}
+
 }} // namespaces
 
 #endif // BOOST_MATH_AIRY_HPP

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-01 14:39:36 EST (Fri, 01 Feb 2013)
@@ -556,10 +556,10 @@
 }
 
 template <class output_iterator, class T, class Policy>
-inline void cyl_bessel_j_zero(output_iterator out_it,
-                              T v,
+inline void cyl_bessel_j_zero(T v,
                               unsigned number_of_zeros,
                               unsigned start_index,
+                              output_iterator out_it,
                               const Policy& pol)
 {
    BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<T>::is_integer, "Order must be a floating-point type.");
@@ -571,12 +571,12 @@
 }
 
 template <class output_iterator, class T>
-inline void cyl_bessel_j_zero(output_iterator out_it,
-                              T v,
+inline void cyl_bessel_j_zero(T v,
                               unsigned number_of_zeros,
-                              unsigned start_index)
+                              unsigned start_index,
+                              output_iterator out_it)
 {
-   return cyl_bessel_j_zero(out_it, v, number_of_zeros, start_index, policies::policy<>());
+   return cyl_bessel_j_zero(v, number_of_zeros, start_index, out_it, policies::policy<>());
 }
 
 template <class T, class Policy>
@@ -597,10 +597,10 @@
 }
 
 template <class output_iterator, class T, class Policy>
-inline void cyl_neumann_zero(output_iterator out_it,
-                             T v,
+inline void cyl_neumann_zero(T v,
                              unsigned number_of_zeros,
                              unsigned start_index,
+                             output_iterator out_it,
                              const Policy& pol)
 {
    BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<T>::is_integer, "Order must be a floating-point type.");
@@ -612,12 +612,12 @@
 }
 
 template <class output_iterator, class T>
-inline void cyl_neumann_zero(output_iterator out_it,
-                             T v,
+inline void cyl_neumann_zero(T v,
                              unsigned number_of_zeros,
-                             unsigned start_index)
+                             unsigned start_index,
+                             output_iterator out_it)
 {
-   return cyl_neumann_zero(out_it, v, number_of_zeros, start_index, policies::policy<>());
+   return cyl_neumann_zero(v, number_of_zeros, start_index, out_it, policies::policy<>());
 }
 
 } // namespace math

Modified: trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp	(original)
+++ trunk/boost/math/special_functions/detail/airy_ai_bi_zero.hpp	2013-02-01 14:39:36 EST (Fri, 01 Feb 2013)
@@ -14,85 +14,155 @@
   #define _AIRY_AI_BI_ZERO_2013_01_20_HPP_
 
   #include <boost/math/constants/constants.hpp>
+  #include <boost/math/special_functions/cbrt.hpp>
 
-  namespace boost { namespace math { namespace detail {
-  namespace airy_zero
+  namespace boost { namespace math {
+  namespace detail
   {
-    template<class T>
-    T equation_as_10_4_105(const T& z)
-    {
-      BOOST_MATH_STD_USING // ADL of std names, needed for pow.
+    // Forward declarations of the needed Airy function implementations.
+    template <class T, class Policy>
+    T airy_ai_imp(T x, const Policy& pol);
+    template <class T, class Policy>
+    T airy_bi_imp(T x, const Policy& pol);
+    template <class T, class Policy>
+    T airy_ai_prime_imp(T x, const Policy& pol);
+    template <class T, class Policy>
+    T airy_bi_prime_imp(T x, const Policy& pol);
 
-      const T one_over_z         = T(1) / z;
-      const T one_over_z_squared = one_over_z * one_over_z;
+    namespace airy_zero
+    {
+      template<class T>
+      T equation_as_10_4_105(const T& z)
+      {
+        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));
+        const T z_pow_third     (boost::math::cbrt(z));
+        const T z_pow_two_thirds(z_pow_third * z_pow_third);
 
-      // Implement the top line of Eq. 10.4.105.
-      const T fz(z_pow_two_thirds * (((((                     + (T(162375596875.0) / 334430208UL)
-                                         * one_over_z_squared - (   T(108056875.0) /   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)));
+        // Implement the top line of Eq. 10.4.105.
+        const T fz(z_pow_two_thirds * (((((                     + (T(162375596875.0) / 334430208UL)
+                                           * one_over_z_squared - (   T(108056875.0) /   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;
-    }
+        return fz;
+      }
 
-    namespace airy_ai_zero_detail
-    {
-      template<class T>
-      T initial_guess(const unsigned m)
+      namespace airy_ai_zero_detail
       {
-        switch(m)
+        template<class T>
+        T initial_guess(const unsigned 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:
+          T guess;
+
+          if(m == 0U)
           {
-            const T t = ((boost::math::constants::pi<T>() * 3U) * T((T(m) * 4U) - 1)) / 8U;
+            // Requesting an estimate of the zero'th root is an error.
+            // Return zero.
+            guess = T(0);
+          }
 
-            return -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+          switch(m)
+          {
+            case 1U: { guess = T(-2.33810741045976703849); break; }
+            case 2U: { guess = T(-4.08794944413097061664); break; }
+            case 3U: { guess = T(-5.52055982809555105913); break; }
+            case 4U: { guess = T(-6.78670809007175899878); break; }
+            case 5U: { guess = T(-7.94413358712085312314); break; }
+            case 6U: { guess = T(-9.02265085334098038016); break; }
+            case 7U: { guess = T(-10.0401743415580859306); break; }
+            case 8U: { guess = T(-11.0085243037332628932); break; }
+            case 9U: { guess = T(-11.9360155632362625170); break; }
+            case 10U:{ guess = T(-12.8287767528657572004); break; }
+            default:
+            {
+              const T t(((boost::math::constants::pi<T>() * 3U) * ((T(m) * 4U) - 1)) / 8U);
+              guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+              break;
+            }
           }
+
+          return guess;
         }
-      }
-    } // namespace airy_ai_zero_detail
 
-    namespace airy_bi_zero_detail
-    {
-      template<class T>
-      T initial_guess(const unsigned m)
+        template<class T, class Policy>
+        class function_object
+        {
+        public:
+          function_object(const Policy pol) : my_pol(pol) { }
+
+          boost::math::tuple<T, T> operator()(const T& x) const
+          {
+            // Return a tuple containing both Ai(x) and Ai'(x).
+            return boost::math::make_tuple(
+              boost::math::detail::airy_ai_imp      (x, my_pol),
+              boost::math::detail::airy_ai_prime_imp(x, my_pol));
+          }
+
+        private:
+          const Policy& my_pol;
+        };
+      } // namespace airy_ai_zero_detail
+
+      namespace airy_bi_zero_detail
       {
-        switch(m)
+        template<class T>
+        T initial_guess(const unsigned 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:
+          T guess;
+
+          if(m == 0U)
           {
-            const T t = ((boost::math::constants::pi<T>() * 3U) * T((T(m) * 4U) - 3)) / 8U;
+            // Requesting an estimate of the zero'th root is an error.
+            // Return zero.
+            guess = T(0);
+          }
 
-            return -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+          switch(m)
+          {
+            case 1U: { guess = T(-1.17371322270912792492); break; }
+            case 2U: { guess = T(-3.27109330283635271568); break; }
+            case 3U: { guess = T(-4.83073784166201593267); break; }
+            case 4U: { guess = T(-6.16985212831025125983); break; }
+            case 5U: { guess = T(-7.37676207936776371360); break; }
+            case 6U: { guess = T(-8.49194884650938801345); break; }
+            case 7U: { guess = T(-9.53819437934623888663); break; }
+            case 8U: { guess = T(-10.5299135067053579244); break; }
+            case 9U: { guess = T(-11.4769535512787794379); break; }
+            case 10U:{ guess = T(-12.3864171385827387456); break; }
+            default:
+            {
+              const T t(((boost::math::constants::pi<T>() * 3U) * ((T(m) * 4U) - 3)) / 8U);
+              guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+              break;
+            }
           }
+
+          return guess;
         }
-      }
-    } // namespace airy_bi_zero_detail
-  } // namespace airy_zero
+
+        template<class T, class Policy>
+        class function_object
+        {
+        public:
+          function_object(const Policy pol) : my_pol(pol) { }
+
+          boost::math::tuple<T, T> operator()(const T& x) const
+          {
+            // Return a tuple containing both Bi(x) and Bi'(x).
+            return boost::math::make_tuple(
+              boost::math::detail::airy_bi_imp      (x, my_pol),
+              boost::math::detail::airy_bi_prime_imp(x, my_pol));
+          }
+
+        private:
+          const Policy& my_pol;
+        };
+      } // namespace airy_bi_zero_detail
+    } // namespace airy_zero
   } // namespace detail
   } // namespace math
   } // namespaces boost

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-01 14:39:36 EST (Fri, 01 Feb 2013)
@@ -94,7 +94,7 @@
 
         // 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
+        // Here, v is larger than about 2.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
@@ -102,9 +102,9 @@
         // 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.
@@ -200,18 +200,20 @@
           {
             // Get the initial estimate of the first root.
 
-            if(v < T(1.2))
+            if(v < T(2.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}]
+              //  Table[{n, BesselJZero[n, 1]}, {n, 0, 22/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);
+              //  Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
+              guess = (((((    - T(0.0008342379046010)
+                           * v + T(0.007590035637410))
+                           * v - T(0.030640914772013))
+                           * v + T(0.078232088020106))
+                           * v - T(0.169668712590620))
+                           * v + T(1.542187960073750))
+                           * v + T(2.4048359915254634);
             }
             else
             {
@@ -221,17 +223,17 @@
           }
           else
           {
-            if(v < T(1.2))
+            if(v < T(2.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>();
+              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);
+              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);
@@ -251,10 +253,24 @@
 
           boost::math::tuple<T, T> operator()(const T& x) const
           {
+            const bool order_is_zero = (   (my_v > -boost::math::tools::epsilon<T>())
+                                        && (my_v < +boost::math::tools::epsilon<T>()));
+
             // 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));
+            T j_v;
+            T j_v_prime;
+
+            if(order_is_zero)
+            {
+              j_v       =  boost::math::detail::cyl_bessel_j_imp(T(0), x, boost::math::detail::bessel_no_int_tag(), my_pol);
+              j_v_prime = -boost::math::detail::cyl_bessel_j_imp(T(1), x, boost::math::detail::bessel_no_int_tag(), my_pol);
+            }
+            else
+            {
+                      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));
+                      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);
@@ -300,18 +316,20 @@
           {
             // Get the initial estimate of the first root.
 
-            if(v < T(1.2))
+            if(v < T(2.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}]
+              //  Table[{n, BesselYZero[n, 1]}, {n, 0, 22/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);
+              //  Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
+              guess = (((((    - T(0.0025095909235652)
+                           * v + T(0.021291887049053))
+                           * v - T(0.076487785486526))
+                           * v + T(0.159110268115362))
+                           * v - T(0.241681668765196))
+                           * v + T(1.4437846310885244))
+                           * v + T(0.89362115190200490);
             }
             else
             {
@@ -321,17 +339,17 @@
           }
           else
           {
-            if(v < T(1.2))
+            if(v < T(2.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>();
+              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);
+              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);

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-02-01 14:39:36 EST (Fri, 01 Feb 2013)
@@ -621,16 +621,17 @@
    typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, unsigned m);
 
    template <class output_iterator, class T>
-   void cyl_bessel_j_zero(output_iterator out_it,
-                                     T v,
-                                     unsigned number_of_zeros,
-                                     unsigned start_index);
+   void cyl_bessel_j_zero(T v,
+                          unsigned number_of_zeros,
+                          unsigned start_index,
+                          output_iterator out_it);
 
    template <class output_iterator, class T, class Policy>
-   void cyl_bessel_j_zero(output_iterator out_it,
-                                     T v,
-                                     unsigned number_of_zeros,
-                                     unsigned start_index, const Policy&);
+   void cyl_bessel_j_zero(T v,
+                          unsigned number_of_zeros,
+                          unsigned start_index,
+                          output_iterator out_it,
+                          const Policy&);
 
    template <class T, class Policy>
    typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zero(T v, unsigned m, const Policy& pol);
@@ -639,16 +640,17 @@
    typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neumann_zero(T v, unsigned m);
 
    template <class output_iterator, class T>
-   void cyl_neumann_zero(output_iterator out_it,
-                                   T v,
-                                   unsigned number_of_zeros,
-                                   unsigned start_index);
+   void cyl_neumann_zero(T v,
+                         unsigned number_of_zeros,
+                         unsigned start_index,
+                         output_iterator out_it);
 
    template <class output_iterator, class T, class Policy>
-   void cyl_neumann_zero(output_iterator out_it,
-                                   T v,
-                                   unsigned number_of_zeros,
-                                   unsigned start_index, const Policy&);
+   void cyl_neumann_zero(T v,
+                         unsigned number_of_zeros,
+                         unsigned start_index,
+                         output_iterator out_it,
+                         const Policy&);
 
    template <class T1, class T2>
    std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x);
@@ -699,6 +701,24 @@
    typename tools::promote_args<T>::type airy_bi_prime(T x);
 
    template <class T, class Policy>
+   typename tools::promote_args<T>::type airy_ai_zero(T dummy, unsigned m);
+
+   template <class output_iterator, class T>
+   void airy_ai_zero(T dummy,
+                     unsigned number_of_zeros,
+                     unsigned start_index,
+                     output_iterator out_it);
+
+   template <class T, class Policy>
+   typename tools::promote_args<T>::type airy_bi_zero(T dummy, unsigned m);
+
+   template <class output_iterator, class T>
+   void airy_bi_zero(T dummy,
+                     unsigned number_of_zeros,
+                     unsigned start_index,
+                     output_iterator out_it);
+
+   template <class T, class Policy>
    typename tools::promote_args<T>::type sin_pi(T x, const Policy&);
 
    template <class T>
@@ -1183,22 +1203,22 @@
    { return boost::math::cyl_bessel_j_zero(v, m, Policy()); }\
 \
 template <class output_iterator, class T>\
-   inline void cyl_bessel_j_zero(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()); }\
+   inline void cyl_bessel_j_zero(T v,\
+                                 unsigned number_of_zeros,\
+                                 unsigned start_index,\
+                                 output_iterator out_it)\
+   { boost::math::cyl_bessel_j_zero(v, number_of_zeros, start_index, out_it, 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(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()); }\
+   inline void cyl_neumann_zero(T v,\
+                                unsigned number_of_zeros,\
+                                unsigned start_index,\
+                                output_iterator out_it)\
+   { boost::math::cyl_neumann_zero(v, number_of_zeros, start_index, out_it, Policy()); }\
 \
    template <class T>\
    inline typename boost::math::tools::promote_args<T>::type sin_pi(T x){ return boost::math::sin_pi(x); }\