Subject: [Boost-commit] svn:boost r83167 - trunk/boost/math/special_functions/detail
From: e_float_at_[hidden]
Date: 2013-02-26 15:33:21
Author: christopher_kormanyos
Date: 2013-02-26 15:33:20 EST (Tue, 26 Feb 2013)
New Revision: 83167
URL: http://svn.boost.org/trac/boost/changeset/83167
Log:
Corrected edge cases for cyl_neumann_zero() with negative order.
Text files modified:
trunk/boost/math/special_functions/detail/bessel_jy_zero.hpp | 36 +++++++++++++++++++++++++++++-------
1 files changed, 29 insertions(+), 7 deletions(-)
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-26 15:33:20 EST (Tue, 26 Feb 2013)
@@ -254,6 +254,24 @@
// There is special handling for negative order.
if(v < 0)
{
+ if((m == 1) && (v > -0.5F))
+ {
+ // For small, negative v, use the results of empirical curve fitting.
+ // Mathematica(R) session for the coefficients:
+ // Table[{n, BesselJZero[n, 1]}, {n, -(1/2), 0, 1/10}]
+ // N[%, 20]
+ // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
+ guess = ((((( - T(0.2321156900729)
+ * v - T(0.1493247777488))
+ * v - T(0.15205419167239))
+ * v + T(0.07814930561249))
+ * v - T(0.17757573537688))
+ * v + T(1.542805677045663))
+ * v + T(2.40482555769577277);
+
+ return guess;
+ }
+
// Create the positive order and extract its positive floor integer part.
const T vv(-v);
const T vv_floor(floor(vv));
@@ -291,7 +309,7 @@
}
else
{
- root_lo *= 0.85F;
+ root_lo *= 0.75F;
}
}
}
@@ -319,7 +337,7 @@
{
// Get the initial estimate of the first root.
- if(v < T(2.2))
+ if(v < 2.2F)
{
// For small v, use the results of empirical curve fitting.
// Mathematica(R) session for the coefficients:
@@ -342,7 +360,7 @@
}
else
{
- if(v < T(2.2))
+ if(v < 2.2F)
{
// 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>());
@@ -509,7 +527,7 @@
}
else
{
- root_lo *= 0.85F;
+ root_lo *= 0.75F;
}
}
}
@@ -517,13 +535,17 @@
{
if(T(vv - vv_floor) < 0.5F)
{
- root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
+ root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
+ root_lo += 0.01F;
+ root_hi += 0.01F;
}
else
{
root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol);
root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
+ root_lo += 0.01F;
+ root_hi += 0.01F;
}
}
@@ -546,7 +568,7 @@
{
// Get the initial estimate of the first root.
- if(v < T(2.2))
+ if(v < 2.2F)
{
// For small v, use the results of empirical curve fitting.
// Mathematica(R) session for the coefficients:
@@ -569,7 +591,7 @@
}
else
{
- if(v < T(2.2))
+ if(v < 2.2F)
{
// 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>());