123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597 |
- #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/special_functions/detail/bessel_jy_series.hpp>
- #include <boost/math/constants/constants.hpp>
- #include <boost/math/policies/error_handling.hpp>
- #include <complex>
- namespace boost { namespace math {
- namespace detail {
-
-
-
-
-
-
-
-
-
- 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;
- }
-
-
- 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_MATH_ASSERT(fabs(v) <= 0.5f);
- 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;
-
- 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<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
- *Y = -sum;
- *Y1 = -2 * sum1 / x;
- return 0;
- }
-
-
- 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
-
-
-
-
- tolerance = 2 * policies::get_epsilon<T, Policy>();
- tiny = sqrt(tools::min_value<T>());
- C = f = 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<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
- *fv = -f;
- *sign = s;
- return 0;
- }
-
-
-
-
-
-
-
- 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;
-
-
- BOOST_MATH_ASSERT(fabs(x) > 1);
-
-
- T tolerance = 2 * policies::get_epsilon<T, Policy>();
- tiny = sqrt(tools::min_value<T>());
- Cr = fr = -0.5f / x;
- Ci = fi = 1;
-
- T v2 = v * v;
- a = (0.25f - v2) / x;
- br = 2 * x;
- bi = 2;
- temp = Cr * Cr + 1;
- Ci = bi + a * Cr / temp;
- Cr = br + a / temp;
- Dr = br;
- Di = bi;
- 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;
- 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;
- 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;
- }
- policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
- *p = fr;
- *q = fi;
- return 0;
- }
- static const int need_j = 1;
- static const int need_y = 2;
-
-
- template <typename T, typename Policy>
- int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
- {
- BOOST_MATH_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;
- int org_kind = kind;
- T cp = 0;
- T sp = 0;
- 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;
- }
- if (v > static_cast<T>((std::numeric_limits<int>::max)()))
- {
- *J = *Y = policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol);
- return 1;
- }
- n = iround(v, pol);
- u = v - n;
- if(reflect)
- {
- T z = (u + n % 2);
- cp = boost::math::cos_pi(z, pol);
- sp = boost::math::sin_pi(z, pol);
- if(u != 0)
- kind = need_j|need_y;
- }
- if(x == 0)
- {
- if (v == 0)
- *J = 1;
- else if ((u == 0) || !reflect)
- *J = 0;
- else if(kind & need_j)
- *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol);
- else
- *J = std::numeric_limits<T>::quiet_NaN();
- if((kind & need_y) == 0)
- *Y = std::numeric_limits<T>::quiet_NaN();
- else
- {
-
- BOOST_MATH_ASSERT(x != 0);
- }
- return 1;
- }
-
- W = T(2) / (x * pi<T>());
- T Yv_scale = 1;
- if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5)))
- {
-
-
-
-
-
- Jv = bessel_j_small_z_series(v, x, pol);
- Yv = std::numeric_limits<T>::quiet_NaN();
- }
- else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
- {
-
-
-
-
- if(kind&need_j)
- Jv = bessel_j_small_z_series(v, x, pol);
- else
- Jv = std::numeric_limits<T>::quiet_NaN();
- if((org_kind&need_y && (!reflect || (cp != 0)))
- || (org_kind & need_j && (reflect && (sp != 0))))
- {
-
- Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN();
- }
- else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
- {
-
-
-
-
- if(kind&need_j)
- Jv = bessel_j_small_z_series(v, x, pol);
- else
- Jv = std::numeric_limits<T>::quiet_NaN();
- if((org_kind&need_y && (!reflect || (cp != 0)))
- || (org_kind & need_j && (reflect && (sp != 0))))
- {
-
- Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN();
-
- }
- else if(asymptotic_bessel_large_x_limit(v, x))
- {
- if(kind&need_y)
- {
- Yv = asymptotic_bessel_y_large_x_2(v, x, pol);
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN();
- if(kind&need_j)
- {
- Jv = asymptotic_bessel_j_large_x_2(v, x, pol);
- }
- else
- Jv = std::numeric_limits<T>::quiet_NaN();
- }
- else if((x > 8) && hankel_PQ(v, x, &p, &q, pol))
- {
-
-
-
-
-
-
-
-
-
-
-
-
-
- T mod_v = fmod(T(v / 2 + 0.25f), T(2));
- T sx = sin(x);
- T cx = cos(x);
- T sv = boost::math::sin_pi(mod_v, pol);
- T cv = boost::math::cos_pi(mod_v, pol);
- T sc = sx * cv - sv * cx;
- T cc = cx * cv + sx * sv;
- T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x));
- Yv = chi * (p * sc + q * cc);
- Jv = chi * (p * cc - q * sc);
- }
- else if (x <= 2)
- {
- if(temme_jy(u, x, &Yu, &Yu1, pol))
- {
-
- *J = *Y = Yu;
- return 1;
- }
- prev = Yu;
- current = Yu1;
- T scale = 1;
- policies::check_series_iterations<T>(function, n, pol);
- for (k = 1; k <= n; k++)
- {
- T fact = 2 * (u + k) / x;
- if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
- {
- scale /= current;
- prev /= current;
- current = 1;
- }
- next = fact * current - prev;
- prev = current;
- current = next;
- }
- Yv = prev;
- Yv1 = current;
- if(kind&need_j)
- {
- CF1_jy(v, x, &fv, &s, pol);
- Jv = scale * W / (Yv * fv - Yv1);
- }
- else
- Jv = std::numeric_limits<T>::quiet_NaN();
- Yv_scale = scale;
- }
- else
- {
-
- T ratio;
- CF1_jy(v, x, &fv, &s, pol);
-
- T init = sqrt(tools::min_value<T>());
- BOOST_MATH_INSTRUMENT_VARIABLE(init);
- prev = fv * s * init;
- current = s * init;
- if(v < max_factorial<T>::value)
- {
- policies::check_series_iterations<T>(function, n, pol);
- for (k = n; k > 0; k--)
- {
- next = 2 * (u + k) * current / x - prev;
-
-
-
- if (next == 0)
- {
- next = prev * tools::epsilon<T>() / 2;
- }
- prev = current;
- current = next;
- }
- ratio = (s * init) / current;
-
- fu = prev / current;
- }
- else
- {
-
-
-
-
- policies::check_series_iterations<T>(function, n, pol);
- bool over = false;
- for (k = n; k > 0; k--)
- {
- T t = 2 * (u + k) / x;
- if((t > 1) && (tools::max_value<T>() / t < current))
- {
- over = true;
- break;
- }
- next = t * current - prev;
- prev = current;
- current = next;
- }
- if(!over)
- {
- ratio = (s * init) / current;
-
- fu = prev / current;
- }
- else
- {
- ratio = 0;
- fu = 1;
- }
- }
- CF2_jy(u, x, &p, &q, pol);
- T t = u / x - fu;
- gamma = (p - t) / q;
-
-
-
-
-
-
- if(gamma == 0)
- {
- gamma = u * tools::epsilon<T>() / x;
- }
- BOOST_MATH_INSTRUMENT_VARIABLE(current);
- BOOST_MATH_INSTRUMENT_VARIABLE(W);
- BOOST_MATH_INSTRUMENT_VARIABLE(q);
- BOOST_MATH_INSTRUMENT_VARIABLE(gamma);
- BOOST_MATH_INSTRUMENT_VARIABLE(p);
- BOOST_MATH_INSTRUMENT_VARIABLE(t);
- Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
- BOOST_MATH_INSTRUMENT_VARIABLE(Ju);
- Jv = Ju * ratio;
- Yu = gamma * Ju;
- Yu1 = Yu * (u/x - p - q/gamma);
- if(kind&need_y)
- {
-
- prev = Yu;
- current = Yu1;
- policies::check_series_iterations<T>(function, n, pol);
- for (k = 1; k <= n; k++)
- {
- T fact = 2 * (u + k) / x;
- if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
- {
- prev /= current;
- Yv_scale /= current;
- current = 1;
- }
- next = fact * current - prev;
- prev = current;
- current = next;
- }
- Yv = prev;
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN();
- }
- if (reflect)
- {
- if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
- *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
- else
- *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale));
- if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
- *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
- else
- *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
- }
- else
- {
- *J = Jv;
- if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
- *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
- else
- *Y = Yv / Yv_scale;
- }
- return 0;
- }
- }
- }}
- #endif
|