// Copyright (c) 2006 Xiaogang Zhang | |
// 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) | |
#ifndef BOOST_MATH_BESSEL_JY_HPP | |
#define BOOST_MATH_BESSEL_JY_HPP | |
#ifdef _MSC_VER | |
#pragma once | |
#endif | |
#include <boost/math/tools/config.hpp> | |
#include <boost/math/special_functions/gamma.hpp> | |
#include <boost/math/special_functions/sign.hpp> | |
#include <boost/math/special_functions/hypot.hpp> | |
#include <boost/math/special_functions/sin_pi.hpp> | |
#include <boost/math/special_functions/cos_pi.hpp> | |
#include <boost/math/special_functions/detail/bessel_jy_asym.hpp> | |
#include <boost/math/constants/constants.hpp> | |
#include <boost/math/policies/error_handling.hpp> | |
#include <boost/mpl/if.hpp> | |
#include <boost/type_traits/is_floating_point.hpp> | |
#include <complex> | |
// Bessel functions of the first and second kind of fractional order | |
namespace boost { namespace math { | |
namespace detail { | |
// | |
// Simultaneous calculation of A&S 9.2.9 and 9.2.10 | |
// for use in A&S 9.2.5 and 9.2.6. | |
// This series is quick to evaluate, but divergent unless | |
// x is very large, in fact it's pretty hard to figure out | |
// with any degree of precision when this series actually | |
// *will* converge!! Consequently, we may just have to | |
// try it and see... | |
// | |
template <class T, class Policy> | |
bool hankel_PQ(T v, T x, T* p, T* q, const Policy& ) | |
{ | |
BOOST_MATH_STD_USING | |
T tolerance = 2 * policies::get_epsilon<T, Policy>(); | |
*p = 1; | |
*q = 0; | |
T k = 1; | |
T z8 = 8 * x; | |
T sq = 1; | |
T mu = 4 * v * v; | |
T term = 1; | |
bool ok = true; | |
do | |
{ | |
term *= (mu - sq * sq) / (k * z8); | |
*q += term; | |
k += 1; | |
sq += 2; | |
T mult = (sq * sq - mu) / (k * z8); | |
ok = fabs(mult) < 0.5f; | |
term *= mult; | |
*p += term; | |
k += 1; | |
sq += 2; | |
} | |
while((fabs(term) > tolerance * *p) && ok); | |
return ok; | |
} | |
// Calculate Y(v, x) and Y(v+1, x) by Temme's method, see | |
// Temme, Journal of Computational Physics, vol 21, 343 (1976) | |
template <typename T, typename Policy> | |
int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol) | |
{ | |
T g, h, p, q, f, coef, sum, sum1, tolerance; | |
T a, d, e, sigma; | |
unsigned long k; | |
BOOST_MATH_STD_USING | |
using namespace boost::math::tools; | |
using namespace boost::math::constants; | |
BOOST_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine | |
T gp = boost::math::tgamma1pm1(v, pol); | |
T gm = boost::math::tgamma1pm1(-v, pol); | |
T spv = boost::math::sin_pi(v, pol); | |
T spv2 = boost::math::sin_pi(v/2, pol); | |
T xp = pow(x/2, v); | |
a = log(x / 2); | |
sigma = -a * v; | |
d = abs(sigma) < tools::epsilon<T>() ? | |
T(1) : sinh(sigma) / sigma; | |
e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2) | |
: T(2 * spv2 * spv2 / v); | |
T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v)); | |
T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2); | |
T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv); | |
f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv; | |
p = vspv / (xp * (1 + gm)); | |
q = vspv * xp / (1 + gp); | |
g = f + e * q; | |
h = p; | |
coef = 1; | |
sum = coef * g; | |
sum1 = coef * h; | |
T v2 = v * v; | |
T coef_mult = -x * x / 4; | |
// series summation | |
tolerance = policies::get_epsilon<T, Policy>(); | |
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++) | |
{ | |
f = (k * f + p + q) / (k*k - v2); | |
p /= k - v; | |
q /= k + v; | |
g = f + e * q; | |
h = p - k * g; | |
coef *= coef_mult / k; | |
sum += coef * g; | |
sum1 += coef * h; | |
if (abs(coef * g) < abs(sum) * tolerance) | |
{ | |
break; | |
} | |
} | |
policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol); | |
*Y = -sum; | |
*Y1 = -2 * sum1 / x; | |
return 0; | |
} | |
// Evaluate continued fraction fv = J_(v+1) / J_v, see | |
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 | |
template <typename T, typename Policy> | |
int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol) | |
{ | |
T C, D, f, a, b, delta, tiny, tolerance; | |
unsigned long k; | |
int s = 1; | |
BOOST_MATH_STD_USING | |
// |x| <= |v|, CF1_jy converges rapidly | |
// |x| > |v|, CF1_jy needs O(|x|) iterations to converge | |
// modified Lentz's method, see | |
// Lentz, Applied Optics, vol 15, 668 (1976) | |
tolerance = 2 * policies::get_epsilon<T, Policy>();; | |
tiny = sqrt(tools::min_value<T>()); | |
C = f = tiny; // b0 = 0, replace with tiny | |
D = 0; | |
for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++) | |
{ | |
a = -1; | |
b = 2 * (v + k) / x; | |
C = b + a / C; | |
D = b + a * D; | |
if (C == 0) { C = tiny; } | |
if (D == 0) { D = tiny; } | |
D = 1 / D; | |
delta = C * D; | |
f *= delta; | |
if (D < 0) { s = -s; } | |
if (abs(delta - 1) < tolerance) | |
{ break; } | |
} | |
policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol); | |
*fv = -f; | |
*sign = s; // sign of denominator | |
return 0; | |
} | |
// | |
// This algorithm was originally written by Xiaogang Zhang | |
// using std::complex to perform the complex arithmetic. | |
// However, that turns out to 10x or more slower than using | |
// all real-valued arithmetic, so it's been rewritten using | |
// real values only. | |
// | |
template <typename T, typename Policy> | |
int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) | |
{ | |
BOOST_MATH_STD_USING | |
T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp; | |
T tiny; | |
unsigned long k; | |
// |x| >= |v|, CF2_jy converges rapidly | |
// |x| -> 0, CF2_jy fails to converge | |
BOOST_ASSERT(fabs(x) > 1); | |
// modified Lentz's method, complex numbers involved, see | |
// Lentz, Applied Optics, vol 15, 668 (1976) | |
T tolerance = 2 * policies::get_epsilon<T, Policy>(); | |
tiny = sqrt(tools::min_value<T>()); | |
Cr = fr = -0.5f / x; | |
Ci = fi = 1; | |
//Dr = Di = 0; | |
T v2 = v * v; | |
a = (0.25f - v2) / x; // Note complex this one time only! | |
br = 2 * x; | |
bi = 2; | |
temp = Cr * Cr + 1; | |
Ci = bi + a * Cr / temp; | |
Cr = br + a / temp; | |
Dr = br; | |
Di = bi; | |
//std::cout << "C = " << Cr << " " << Ci << std::endl; | |
//std::cout << "D = " << Dr << " " << Di << std::endl; | |
if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; } | |
if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; } | |
temp = Dr * Dr + Di * Di; | |
Dr = Dr / temp; | |
Di = -Di / temp; | |
delta_r = Cr * Dr - Ci * Di; | |
delta_i = Ci * Dr + Cr * Di; | |
temp = fr; | |
fr = temp * delta_r - fi * delta_i; | |
fi = temp * delta_i + fi * delta_r; | |
//std::cout << fr << " " << fi << std::endl; | |
for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) | |
{ | |
a = k - 0.5f; | |
a *= a; | |
a -= v2; | |
bi += 2; | |
temp = Cr * Cr + Ci * Ci; | |
Cr = br + a * Cr / temp; | |
Ci = bi - a * Ci / temp; | |
Dr = br + a * Dr; | |
Di = bi + a * Di; | |
//std::cout << "C = " << Cr << " " << Ci << std::endl; | |
//std::cout << "D = " << Dr << " " << Di << std::endl; | |
if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; } | |
if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; } | |
temp = Dr * Dr + Di * Di; | |
Dr = Dr / temp; | |
Di = -Di / temp; | |
delta_r = Cr * Dr - Ci * Di; | |
delta_i = Ci * Dr + Cr * Di; | |
temp = fr; | |
fr = temp * delta_r - fi * delta_i; | |
fi = temp * delta_i + fi * delta_r; | |
if (fabs(delta_r - 1) + fabs(delta_i) < tolerance) | |
break; | |
//std::cout << fr << " " << fi << std::endl; | |
} | |
policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol); | |
*p = fr; | |
*q = fi; | |
return 0; | |
} | |
enum | |
{ | |
need_j = 1, need_y = 2 | |
}; | |
// Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see | |
// Barnett et al, Computer Physics Communications, vol 8, 377 (1974) | |
template <typename T, typename Policy> | |
int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol) | |
{ | |
BOOST_ASSERT(x >= 0); | |
T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu; | |
T W, p, q, gamma, current, prev, next; | |
bool reflect = false; | |
unsigned n, k; | |
int s; | |
static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)"; | |
BOOST_MATH_STD_USING | |
using namespace boost::math::tools; | |
using namespace boost::math::constants; | |
if (v < 0) | |
{ | |
reflect = true; | |
v = -v; // v is non-negative from here | |
kind = need_j|need_y; // need both for reflection formula | |
} | |
n = iround(v, pol); | |
u = v - n; // -1/2 <= u < 1/2 | |
if (x == 0) | |
{ | |
*J = *Y = policies::raise_overflow_error<T>( | |
function, 0, pol); | |
return 1; | |
} | |
// x is positive until reflection | |
W = T(2) / (x * pi<T>()); // Wronskian | |
if((x > 8) && (x < 1000) && hankel_PQ(v, x, &p, &q, pol)) | |
{ | |
// | |
// Hankel approximation: note that this method works best when x | |
// is large, but in that case we end up calculating sines and cosines | |
// of large values, with horrendous resulting accuracy. It is fast though | |
// when it works.... | |
// | |
T chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>(); | |
T sc = sin(chi); | |
T cc = cos(chi); | |
chi = sqrt(2 / (boost::math::constants::pi<T>() * x)); | |
Yv = chi * (p * sc + q * cc); | |
Jv = chi * (p * cc - q * sc); | |
} | |
else if (x <= 2) // x in (0, 2] | |
{ | |
if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series | |
{ | |
// domain error: | |
*J = *Y = Yu; | |
return 1; | |
} | |
prev = Yu; | |
current = Yu1; | |
for (k = 1; k <= n; k++) // forward recurrence for Y | |
{ | |
next = 2 * (u + k) * current / x - prev; | |
prev = current; | |
current = next; | |
} | |
Yv = prev; | |
Yv1 = current; | |
if(kind&need_j) | |
{ | |
CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy | |
Jv = W / (Yv * fv - Yv1); // Wronskian relation | |
} | |
else | |
Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. | |
} | |
else // x in (2, \infty) | |
{ | |
// Get Y(u, x): | |
// define tag type that will dispatch to right limits: | |
typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type; | |
T lim; | |
switch(kind) | |
{ | |
case need_j: | |
lim = asymptotic_bessel_j_limit<T>(v, tag_type()); | |
break; | |
case need_y: | |
lim = asymptotic_bessel_y_limit<T>(tag_type()); | |
break; | |
default: | |
lim = (std::max)( | |
asymptotic_bessel_j_limit<T>(v, tag_type()), | |
asymptotic_bessel_y_limit<T>(tag_type())); | |
break; | |
} | |
if(x > lim) | |
{ | |
if(kind&need_y) | |
{ | |
Yu = asymptotic_bessel_y_large_x_2(u, x); | |
Yu1 = asymptotic_bessel_y_large_x_2(T(u + 1), x); | |
} | |
else | |
Yu = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. | |
if(kind&need_j) | |
{ | |
Jv = asymptotic_bessel_j_large_x_2(v, x); | |
} | |
else | |
Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. | |
} | |
else | |
{ | |
CF1_jy(v, x, &fv, &s, pol); | |
// tiny initial value to prevent overflow | |
T init = sqrt(tools::min_value<T>()); | |
prev = fv * s * init; | |
current = s * init; | |
for (k = n; k > 0; k--) // backward recurrence for J | |
{ | |
next = 2 * (u + k) * current / x - prev; | |
prev = current; | |
current = next; | |
} | |
T ratio = (s * init) / current; // scaling ratio | |
// can also call CF1_jy() to get fu, not much difference in precision | |
fu = prev / current; | |
CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy | |
T t = u / x - fu; // t = J'/J | |
gamma = (p - t) / q; | |
Ju = sign(current) * sqrt(W / (q + gamma * (p - t))); | |
Jv = Ju * ratio; // normalization | |
Yu = gamma * Ju; | |
Yu1 = Yu * (u/x - p - q/gamma); | |
} | |
if(kind&need_y) | |
{ | |
// compute Y: | |
prev = Yu; | |
current = Yu1; | |
for (k = 1; k <= n; k++) // forward recurrence for Y | |
{ | |
next = 2 * (u + k) * current / x - prev; | |
prev = current; | |
current = next; | |
} | |
Yv = prev; | |
} | |
else | |
Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. | |
} | |
if (reflect) | |
{ | |
T z = (u + n % 2); | |
*J = boost::math::cos_pi(z, pol) * Jv - boost::math::sin_pi(z, pol) * Yv; // reflection formula | |
*Y = boost::math::sin_pi(z, pol) * Jv + boost::math::cos_pi(z, pol) * Yv; | |
} | |
else | |
{ | |
*J = Jv; | |
*Y = Yv; | |
} | |
return 0; | |
} | |
} // namespace detail | |
}} // namespaces | |
#endif // BOOST_MATH_BESSEL_JY_HPP | |