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); }\