summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/detail/bessel_jy_asym.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/special_functions/detail/bessel_jy_asym.hpp')
-rw-r--r--boost/math/special_functions/detail/bessel_jy_asym.hpp176
1 files changed, 33 insertions, 143 deletions
diff --git a/boost/math/special_functions/detail/bessel_jy_asym.hpp b/boost/math/special_functions/detail/bessel_jy_asym.hpp
index 0021f8c86a..81f6238e58 100644
--- a/boost/math/special_functions/detail/bessel_jy_asym.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_asym.hpp
@@ -21,61 +21,6 @@
namespace boost{ namespace math{ namespace detail{
template <class T>
-inline T asymptotic_bessel_j_large_x_P(T v, T x)
-{
- // A&S 9.2.9
- T s = 1;
- T mu = 4 * v * v;
- T ez2 = 8 * x;
- ez2 *= ez2;
- s -= (mu-1) * (mu-9) / (2 * ez2);
- s += (mu-1) * (mu-9) * (mu-25) * (mu - 49) / (24 * ez2 * ez2);
- return s;
-}
-
-template <class T>
-inline T asymptotic_bessel_j_large_x_Q(T v, T x)
-{
- // A&S 9.2.10
- T s = 0;
- T mu = 4 * v * v;
- T ez = 8*x;
- s += (mu-1) / ez;
- s -= (mu-1) * (mu-9) * (mu-25) / (6 * ez*ez*ez);
- return s;
-}
-
-template <class T>
-inline T asymptotic_bessel_j_large_x(T v, T x)
-{
- //
- // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
- //
- // Also A&S 9.2.5
- //
- BOOST_MATH_STD_USING // ADL of std names
- T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
- return sqrt(2 / (constants::pi<T>() * x))
- * (asymptotic_bessel_j_large_x_P(v, x) * cos(chi)
- - asymptotic_bessel_j_large_x_Q(v, x) * sin(chi));
-}
-
-template <class T>
-inline T asymptotic_bessel_y_large_x(T v, T x)
-{
- //
- // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
- //
- // Also A&S 9.2.5
- //
- BOOST_MATH_STD_USING // ADL of std names
- T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
- return sqrt(2 / (constants::pi<T>() * x))
- * (asymptotic_bessel_j_large_x_P(v, x) * sin(chi)
- - asymptotic_bessel_j_large_x_Q(v, x) * cos(chi));
-}
-
-template <class T>
inline T asymptotic_bessel_amplitude(T v, T x)
{
// Calculate the amplitude of J(v, x) and Y(v, x) for large
@@ -99,13 +44,14 @@ T asymptotic_bessel_phase_mx(T v, T x)
//
// Calculate the phase of J(v, x) and Y(v, x) for large x.
// See A&S 9.2.29.
- // Note that the result returned is the phase less x.
+ // Note that the result returned is the phase less (x - PI(v/2 + 1/4))
+ // which we'll factor in later when we calculate the sines/cosines of the result:
//
T mu = 4 * v * v;
T denom = 4 * x;
T denom_mult = denom * denom;
- T s = -constants::pi<T>() * (v / 2 + 0.25f);
+ T s = 0;
s += (mu - 1) / (2 * denom);
denom *= denom_mult;
s += (mu - 1) * (mu - 25) / (6 * denom);
@@ -127,10 +73,16 @@ inline T asymptotic_bessel_y_large_x_2(T v, T x)
BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
BOOST_MATH_INSTRUMENT_VARIABLE(phase);
//
- // Calculate the sine of the phase, using:
- // sin(x+p) = sin(x)cos(p) + cos(x)sin(p)
+ // Calculate the sine of the phase, using
+ // sine/cosine addition rules to factor in
+ // the x - PI(v/2 + 1/4) term not added to the
+ // phase when we calculated it.
//
- T sin_phase = sin(phase) * cos(x) + cos(phase) * sin(x);
+ T cx = cos(x);
+ T sx = sin(x);
+ T ci = cos_pi(v / 2 + 0.25f);
+ T si = sin_pi(v / 2 + 0.25f);
+ T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
BOOST_MATH_INSTRUMENT_CODE(sin(phase));
BOOST_MATH_INSTRUMENT_CODE(cos(x));
BOOST_MATH_INSTRUMENT_CODE(cos(phase));
@@ -149,101 +101,39 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x)
BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
BOOST_MATH_INSTRUMENT_VARIABLE(phase);
//
- // Calculate the sine of the phase, using:
- // cos(x+p) = cos(x)cos(p) - sin(x)sin(p)
+ // Calculate the sine of the phase, using
+ // sine/cosine addition rules to factor in
+ // the x - PI(v/2 + 1/4) term not added to the
+ // phase when we calculated it.
//
BOOST_MATH_INSTRUMENT_CODE(cos(phase));
BOOST_MATH_INSTRUMENT_CODE(cos(x));
BOOST_MATH_INSTRUMENT_CODE(sin(phase));
BOOST_MATH_INSTRUMENT_CODE(sin(x));
- T sin_phase = cos(phase) * cos(x) - sin(phase) * sin(x);
+ T cx = cos(x);
+ T sx = sin(x);
+ T ci = cos_pi(v / 2 + 0.25f);
+ T si = sin_pi(v / 2 + 0.25f);
+ T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
return sin_phase * ampl;
}
-//
-// Various limits for the J and Y asymptotics
-// (the asympotic expansions are safe to use if
-// x is less than the limit given).
-// We assume that if we don't use these expansions then the
-// error will likely be >100eps, so the limits given are chosen
-// to lead to < 100eps truncation error.
-//
template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<0>&)
+inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
{
- // default case:
BOOST_MATH_STD_USING
- return 2.25 / pow(100 * tools::epsilon<T>() / T(0.001f), T(0.2f));
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<53>&)
-{
- // double case:
- return 304 /*780*/;
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<64>&)
-{
- // 80-bit extended-double case:
- return 1552 /*3500*/;
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<113>&)
-{
- // 128-bit long double case:
- return 1245243 /*3128000*/;
-}
-
-template <class T, class Policy>
-struct bessel_asymptotic_tag
-{
- typedef typename policies::precision<T, Policy>::type precision_type;
- typedef typename mpl::if_<
- mpl::or_<
- mpl::equal_to<precision_type, mpl::int_<0> >,
- mpl::greater<precision_type, mpl::int_<113> > >,
- mpl::int_<0>,
- typename mpl::if_<
- mpl::greater<precision_type, mpl::int_<64> >,
- mpl::int_<113>,
- typename mpl::if_<
- mpl::greater<precision_type, mpl::int_<53> >,
- mpl::int_<64>,
- mpl::int_<53>
- >::type
- >::type
- >::type type;
-};
-
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&)
-{
- // default case:
- BOOST_MATH_STD_USING
- T v2 = (std::max)(T(3), T(v * v));
- return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&)
-{
- // double case:
- T v2 = (std::max)(T(3), T(v * v));
- return v2 * 33 /*73*/;
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&)
-{
- // 80-bit extended-double case:
- T v2 = (std::max)(T(3), T(v * v));
- return v2 * 121 /*266*/;
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&)
-{
- // 128-bit long double case:
- T v2 = (std::max)(T(3), T(v * v));
- return v2 * 39154 /*85700*/;
+ //
+ // Determines if x is large enough compared to v to take the asymptotic
+ // forms above. From A&S 9.2.28 we require:
+ // v < x * eps^1/8
+ // and from A&S 9.2.29 we require:
+ // v^12/10 < 1.5 * x * eps^1/10
+ // using the former seems to work OK in practice with broadly similar
+ // error rates either side of the divide for v < 10000.
+ // At double precision eps^1/8 ~= 0.01.
+ //
+ return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
}
template <class T, class Policy>
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-11 01:53:20 +0000