special_functions/detail/bessel_jy.hpp

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_BESSEL_JY_HPP
0007 #define BOOST_MATH_BESSEL_JY_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/tools/numeric_limits.hpp>
0015 #include <boost/math/tools/type_traits.hpp>
0016 #include <boost/math/special_functions/gamma.hpp>
0017 #include <boost/math/special_functions/sign.hpp>
0018 #include <boost/math/special_functions/hypot.hpp>
0019 #include <boost/math/special_functions/sin_pi.hpp>
0020 #include <boost/math/special_functions/cos_pi.hpp>
0021 #include <boost/math/special_functions/round.hpp>
0022 #include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
0023 #include <boost/math/special_functions/detail/bessel_jy_series.hpp>
0024 #include <boost/math/constants/constants.hpp>
0025 #include <boost/math/policies/error_handling.hpp>
0026
0027 // Bessel functions of the first and second kind of fractional order
0028
0029 namespace boost { namespace math {
0030
0031    namespace detail {
0032
0033       //
0034       // Simultaneous calculation of A&S 9.2.9 and 9.2.10
0035       // for use in A&S 9.2.5 and 9.2.6.
0036       // This series is quick to evaluate, but divergent unless
0037       // x is very large, in fact it's pretty hard to figure out
0038       // with any degree of precision when this series actually
0039       // *will* converge!!  Consequently, we may just have to
0040       // try it and see...
0041       //
0042       template <class T, class Policy>
0043       BOOST_MATH_GPU_ENABLED bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )
0044       {
0045          BOOST_MATH_STD_USING
0046             T tolerance = 2 * policies::get_epsilon<T, Policy>();
0047          *p = 1;
0048          *q = 0;
0049          T k = 1;
0050          T z8 = 8 * x;
0051          T sq = 1;
0052          T mu = 4 * v * v;
0053          T term = 1;
0054          bool ok = true;
0055          do
0056          {
0057             term *= (mu - sq * sq) / (k * z8);
0058             *q += term;
0059             k += 1;
0060             sq += 2;
0061             T mult = (sq * sq - mu) / (k * z8);
0062             ok = fabs(mult) < 0.5f;
0063             term *= mult;
0064             *p += term;
0065             k += 1;
0066             sq += 2;
0067          }
0068          while((fabs(term) > tolerance * *p) && ok);
0069          return ok;
0070       }
0071
0072       // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
0073       // Temme, Journal of Computational Physics, vol 21, 343 (1976)
0074       template <typename T, typename Policy>
0075       BOOST_MATH_GPU_ENABLED int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
0076       {
0077          T g, h, p, q, f, coef, sum, sum1, tolerance;
0078          T a, d, e, sigma;
0079          unsigned long k;
0080
0081          BOOST_MATH_STD_USING
0082             using namespace boost::math::tools;
0083          using namespace boost::math::constants;
0084
0085          BOOST_MATH_ASSERT(fabs(v) <= 0.5f);  // precondition for using this routine
0086
0087          T gp = boost::math::tgamma1pm1(v, pol);
0088          T gm = boost::math::tgamma1pm1(-v, pol);
0089          T spv = boost::math::sin_pi(v, pol);
0090          T spv2 = boost::math::sin_pi(v/2, pol);
0091          T xp = pow(x/2, v);
0092
0093          a = log(x / 2);
0094          sigma = -a * v;
0095          d = abs(sigma) < tools::epsilon<T>() ?
0096             T(1) : sinh(sigma) / sigma;
0097          e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
0098             : T(2 * spv2 * spv2 / v);
0099
0100          T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
0101          T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
0102          T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
0103          f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
0104
0105          p = vspv / (xp * (1 + gm));
0106          q = vspv * xp / (1 + gp);
0107
0108          g = f + e * q;
0109          h = p;
0110          coef = 1;
0111          sum = coef * g;
0112          sum1 = coef * h;
0113
0114          T v2 = v * v;
0115          T coef_mult = -x * x / 4;
0116
0117          // series summation
0118          tolerance = policies::get_epsilon<T, Policy>();
0119          for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
0120          {
0121             f = (k * f + p + q) / (k*k - v2);
0122             p /= k - v;
0123             q /= k + v;
0124             g = f + e * q;
0125             h = p - k * g;
0126             coef *= coef_mult / k;
0127             sum += coef * g;
0128             sum1 += coef * h;
0129             if (abs(coef * g) < abs(sum) * tolerance)
0130             {
0131                break;
0132             }
0133          }
0134          policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
0135          *Y = -sum;
0136          *Y1 = -2 * sum1 / x;
0137
0138          return 0;
0139       }
0140
0141       // Evaluate continued fraction fv = J_(v+1) / J_v, see
0142       // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
0143       template <typename T, typename Policy>
0144       BOOST_MATH_GPU_ENABLED int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
0145       {
0146          T C, D, f, a, b, delta, tiny, tolerance;
0147          unsigned long k;
0148          int s = 1;
0149
0150          BOOST_MATH_STD_USING
0151
0152             // |x| <= |v|, CF1_jy converges rapidly
0153             // |x| > |v|, CF1_jy needs O(|x|) iterations to converge
0154
0155             // modified Lentz's method, see
0156             // Lentz, Applied Optics, vol 15, 668 (1976)
0157             tolerance = 2 * policies::get_epsilon<T, Policy>();
0158          tiny = sqrt(tools::min_value<T>());
0159          C = f = tiny;                           // b0 = 0, replace with tiny
0160          D = 0;
0161          for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
0162          {
0163             a = -1;
0164             b = 2 * (v + k) / x;
0165             C = b + a / C;
0166             D = b + a * D;
0167             if (C == 0) { C = tiny; }
0168             if (D == 0) { D = tiny; }
0169             D = 1 / D;
0170             delta = C * D;
0171             f *= delta;
0172             if (D < 0) { s = -s; }
0173             if (abs(delta - 1) < tolerance)
0174             { break; }
0175          }
0176          policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
0177          *fv = -f;
0178          *sign = s;                              // sign of denominator
0179
0180          return 0;
0181       }
0182       //
0183       // This algorithm was originally written by Xiaogang Zhang
0184       // using std::complex to perform the complex arithmetic.
0185       // However, that turns out to 10x or more slower than using
0186       // all real-valued arithmetic, so it's been rewritten using
0187       // real values only.
0188       //
0189       template <typename T, typename Policy>
0190       BOOST_MATH_GPU_ENABLED int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
0191       {
0192          BOOST_MATH_STD_USING
0193
0194             T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
0195          T tiny;
0196          unsigned long k;
0197
0198          // |x| >= |v|, CF2_jy converges rapidly
0199          // |x| -> 0, CF2_jy fails to converge
0200          BOOST_MATH_ASSERT(fabs(x) > 1);
0201
0202          // modified Lentz's method, complex numbers involved, see
0203          // Lentz, Applied Optics, vol 15, 668 (1976)
0204          T tolerance = 2 * policies::get_epsilon<T, Policy>();
0205          tiny = sqrt(tools::min_value<T>());
0206          Cr = fr = -0.5f / x;
0207          Ci = fi = 1;
0208          //Dr = Di = 0;
0209          T v2 = v * v;
0210          a = (0.25f - v2) / x; // Note complex this one time only!
0211          br = 2 * x;
0212          bi = 2;
0213          temp = Cr * Cr + 1;
0214          Ci = bi + a * Cr / temp;
0215          Cr = br + a / temp;
0216          Dr = br;
0217          Di = bi;
0218          if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
0219          if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
0220          temp = Dr * Dr + Di * Di;
0221          Dr = Dr / temp;
0222          Di = -Di / temp;
0223          delta_r = Cr * Dr - Ci * Di;
0224          delta_i = Ci * Dr + Cr * Di;
0225          temp = fr;
0226          fr = temp * delta_r - fi * delta_i;
0227          fi = temp * delta_i + fi * delta_r;
0228          for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
0229          {
0230             a = k - 0.5f;
0231             a *= a;
0232             a -= v2;
0233             bi += 2;
0234             temp = Cr * Cr + Ci * Ci;
0235             Cr = br + a * Cr / temp;
0236             Ci = bi - a * Ci / temp;
0237             Dr = br + a * Dr;
0238             Di = bi + a * Di;
0239             if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
0240             if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
0241             temp = Dr * Dr + Di * Di;
0242             Dr = Dr / temp;
0243             Di = -Di / temp;
0244             delta_r = Cr * Dr - Ci * Di;
0245             delta_i = Ci * Dr + Cr * Di;
0246             temp = fr;
0247             fr = temp * delta_r - fi * delta_i;
0248             fi = temp * delta_i + fi * delta_r;
0249             if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
0250                break;
0251          }
0252          policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
0253          *p = fr;
0254          *q = fi;
0255
0256          return 0;
0257       }
0258
0259       BOOST_MATH_STATIC const int need_j = 1;
0260       BOOST_MATH_STATIC const int need_y = 2;
0261
0262       // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
0263       // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
0264       template <typename T, typename Policy>
0265       BOOST_MATH_GPU_ENABLED int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
0266       {
0267          BOOST_MATH_ASSERT(x >= 0);
0268
0269          T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
0270          T W, p, q, gamma, current, prev, next;
0271          bool reflect = false;
0272          unsigned n, k;
0273          int s;
0274          int org_kind = kind;
0275          T cp = 0;
0276          T sp = 0;
0277
0278          constexpr auto function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
0279
0280          BOOST_MATH_STD_USING
0281             using namespace boost::math::tools;
0282          using namespace boost::math::constants;
0283
0284          if (v < 0)
0285          {
0286             reflect = true;
0287             v = -v;                             // v is non-negative from here
0288          }
0289          if (v > static_cast<T>((boost::math::numeric_limits<int>::max)()))
0290          {
0291             *J = *Y = policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol);
0292             return 1;  // LCOV_EXCL_LINE previous line will throw.
0293          }
0294          n = iround(v, pol);
0295          u = v - n;                              // -1/2 <= u < 1/2
0296
0297          if(reflect)
0298          {
0299             T z = (u + n % 2);
0300             cp = boost::math::cos_pi(z, pol);
0301             sp = boost::math::sin_pi(z, pol);
0302             if(u != 0)
0303                kind = need_j|need_y;               // need both for reflection formula
0304          }
0305
0306          if(x == 0)
0307          {
0308             if (v == 0)
0309                *J = 1; // LCOV_EXCL_LINE multiprecision case only
0310             else if ((u == 0) || !reflect)
0311                *J = 0;
0312             else if(kind & need_j)
0313                *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity
0314             else
0315                *J = boost::math::numeric_limits<T>::quiet_NaN();  // LCOV_EXCL_LINE, we should never get here, any value will do, not using J.
0316
0317             if((kind & need_y) == 0)
0318                *Y = boost::math::numeric_limits<T>::quiet_NaN();  // any value will do, not using Y.
0319             else
0320             {
0321                // We shoud never get here:
0322                BOOST_MATH_ASSERT(x != 0); // LCOV_EXCL_LINE
0323             }
0324             return 1;
0325          }
0326
0327          // x is positive until reflection
0328          W = T(2) / (x * pi<T>());               // Wronskian
0329          T Yv_scale = 1;
0330          if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5)))
0331          {
0332             //
0333             // This series will actually converge rapidly for all small
0334             // x - say up to x < 20 - but the first few terms are large
0335             // and divergent which leads to large errors :-(
0336             //
0337             Jv = bessel_j_small_z_series(v, x, pol);
0338             Yv = boost::math::numeric_limits<T>::quiet_NaN();
0339          }
0340          else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
0341          {
0342             // Evaluate using series representations.
0343             // This is particularly important for x << v as in this
0344             // area temme_jy may be slow to converge, if it converges at all.
0345             // Requires x is not an integer.
0346             if(kind&need_j)
0347                Jv = bessel_j_small_z_series(v, x, pol);
0348             else
0349                Jv = boost::math::numeric_limits<T>::quiet_NaN();
0350             if((org_kind&need_y && (!reflect || (cp != 0)))
0351                || (org_kind & need_j && (reflect && (sp != 0))))
0352             {
0353                // Only calculate if we need it, and if the reflection formula will actually use it:
0354                Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
0355             }
0356             else
0357                Yv = boost::math::numeric_limits<T>::quiet_NaN();
0358          }
0359          else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
0360          {
0361             // Truncated series evaluation for small x and v an integer,
0362             // much quicker in this area than temme_jy below.
0363             // This code is only used in the multiprecision case, otherwise
0364             // we go via bessel_jn.  LCOV_EXCL_START
0365             if(kind&need_j)
0366                Jv = bessel_j_small_z_series(v, x, pol);
0367             else
0368                Jv = boost::math::numeric_limits<T>::quiet_NaN();
0369             if((org_kind&need_y && (!reflect || (cp != 0)))
0370                || (org_kind & need_j && (reflect && (sp != 0))))
0371             {
0372                // Only calculate if we need it, and if the reflection formula will actually use it:
0373                Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
0374             }
0375             else
0376                Yv = boost::math::numeric_limits<T>::quiet_NaN();
0377             // LCOV_EXCL_STOP
0378          }
0379          else if(asymptotic_bessel_large_x_limit(v, x))
0380          {
0381             if(kind&need_y)
0382             {
0383                Yv = asymptotic_bessel_y_large_x_2(v, x, pol);
0384             }
0385             else
0386                Yv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
0387             if(kind&need_j)
0388             {
0389                Jv = asymptotic_bessel_j_large_x_2(v, x, pol);
0390             }
0391             else
0392                Jv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
0393          }
0394          else if((x > 8) && hankel_PQ(v, x, &p, &q, pol))
0395          {
0396             //
0397             // Hankel approximation: note that this method works best when x
0398             // is large, but in that case we end up calculating sines and cosines
0399             // of large values, with horrendous resulting accuracy.  It is fast though
0400             // when it works....
0401             //
0402             // Normally we calculate sin/cos(chi) where:
0403             //
0404             // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
0405             //
0406             // But this introduces large errors, so use sin/cos addition formulae to
0407             // improve accuracy:
0408             //
0409             T mod_v = fmod(T(v / 2 + 0.25f), T(2));
0410             T sx = sin(x);
0411             T cx = cos(x);
0412             T sv = boost::math::sin_pi(mod_v, pol);
0413             T cv = boost::math::cos_pi(mod_v, pol);
0414
0415             T sc = sx * cv - sv * cx; // == sin(chi);
0416             T cc = cx * cv + sx * sv; // == cos(chi);
0417             T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x));
0418             Yv = chi * (p * sc + q * cc);
0419             Jv = chi * (p * cc - q * sc);
0420          }
0421          else if (x <= 2)                           // x in (0, 2]
0422          {
0423             if(temme_jy(u, x, &Yu, &Yu1, pol))             // Temme series
0424             {
0425                // domain error, this should really have already been handled.
0426                *J = *Y = Yu; // LCOV_EXCL_LINE
0427                return 1;     // LCOV_EXCL_LINE
0428             }
0429             prev = Yu;
0430             current = Yu1;
0431             T scale = 1;
0432             policies::check_series_iterations<T>(function, n, pol);
0433             for (k = 1; k <= n; k++)            // forward recurrence for Y
0434             {
0435                T fact = 2 * (u + k) / x;
0436                if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
0437                {
0438                   scale /= current;
0439                   prev /= current;
0440                   current = 1;
0441                }
0442                next = fact * current - prev;
0443                prev = current;
0444                current = next;
0445             }
0446             Yv = prev;
0447             Yv1 = current;
0448             if(kind&need_j)
0449             {
0450                CF1_jy(v, x, &fv, &s, pol);                 // continued fraction CF1_jy
0451                Jv = scale * W / (Yv * fv - Yv1);           // Wronskian relation
0452             }
0453             else
0454                Jv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
0455             Yv_scale = scale;
0456          }
0457          else                                    // x in (2, \infty)
0458          {
0459             // Get Y(u, x):
0460
0461             T ratio;
0462             CF1_jy(v, x, &fv, &s, pol);
0463             // tiny initial value to prevent overflow
0464             T init = sqrt(tools::min_value<T>());
0465             BOOST_MATH_INSTRUMENT_VARIABLE(init);
0466             prev = fv * s * init;
0467             current = s * init;
0468             if(v < max_factorial<T>::value)
0469             {
0470                policies::check_series_iterations<T>(function, n, pol);
0471                for (k = n; k > 0; k--)             // backward recurrence for J
0472                {
0473                   next = 2 * (u + k) * current / x - prev;
0474                   //
0475                   // We can't allow next to completely cancel out or the subsequent logic breaks.
0476                   // Pretend that one bit did not cancel:
0477                   if (next == 0)
0478                   {
0479                      next = prev * tools::epsilon<T>() / 2;  // LCOV_EXCL_LINE requires specific hardware and rounding to trigger, does get tested on msvc
0480                   }
0481                   prev = current;
0482                   current = next;
0483                }
0484                ratio = (s * init) / current;     // scaling ratio
0485                // can also call CF1_jy() to get fu, not much difference in precision
0486                fu = prev / current;
0487             }
0488             else
0489             {
0490                //
0491                // When v is large we may get overflow in this calculation
0492                // leading to NaN's and other nasty surprises:
0493                //
0494                policies::check_series_iterations<T>(function, n, pol);
0495                bool over = false;
0496                for (k = n; k > 0; k--)             // backward recurrence for J
0497                {
0498                   T t = 2 * (u + k) / x;
0499                   if((t > 1) && (tools::max_value<T>() / t < current))
0500                   {
0501                      over = true;
0502                      break;
0503                   }
0504                   next = t * current - prev;
0505                   prev = current;
0506                   current = next;
0507                }
0508                if(!over)
0509                {
0510                   ratio = (s * init) / current;     // scaling ratio
0511                   // can also call CF1_jy() to get fu, not much difference in precision
0512                   fu = prev / current;
0513                }
0514                else
0515                {
0516                   ratio = 0;
0517                   fu = 1;
0518                }
0519             }
0520             CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
0521             T t = u / x - fu;                   // t = J'/J
0522             gamma = (p - t) / q;
0523             //
0524             // We can't allow gamma to cancel out to zero completely as it messes up
0525             // the subsequent logic.  So pretend that one bit didn't cancel out
0526             // and set to a suitably small value.  The only test case we've been able to
0527             // find for this, is when v = 8.5 and x = 4*PI.
0528             //
0529             if(gamma == 0)
0530             {
0531                gamma = u * tools::epsilon<T>() / x;  // LCOV_EXCL_LINE requires specific hardware and rounding to trigger, does get tested on msvc
0532             }
0533             BOOST_MATH_INSTRUMENT_VARIABLE(current);
0534             BOOST_MATH_INSTRUMENT_VARIABLE(W);
0535             BOOST_MATH_INSTRUMENT_VARIABLE(q);
0536             BOOST_MATH_INSTRUMENT_VARIABLE(gamma);
0537             BOOST_MATH_INSTRUMENT_VARIABLE(p);
0538             BOOST_MATH_INSTRUMENT_VARIABLE(t);
0539             Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
0540             BOOST_MATH_INSTRUMENT_VARIABLE(Ju);
0541
0542             Jv = Ju * ratio;                    // normalization
0543
0544             Yu = gamma * Ju;
0545             Yu1 = Yu * (u/x - p - q/gamma);
0546
0547             if(kind&need_y)
0548             {
0549                // compute Y:
0550                prev = Yu;
0551                current = Yu1;
0552                policies::check_series_iterations<T>(function, n, pol);
0553                for (k = 1; k <= n; k++)            // forward recurrence for Y
0554                {
0555                   T fact = 2 * (u + k) / x;
0556                   if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
0557                   {
0558                      prev /= current;
0559                      Yv_scale /= current;
0560                      current = 1;
0561                   }
0562                   next = fact * current - prev;
0563                   prev = current;
0564                   current = next;
0565                }
0566                Yv = prev;
0567             }
0568             else
0569                Yv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
0570          }
0571
0572          if (reflect)
0573          {
0574             if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
0575                *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);
0576             else
0577                *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale));     // reflection formula
0578             if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
0579                *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);
0580             else
0581                *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
0582          }
0583          else
0584          {
0585             *J = Jv;
0586             if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
0587                *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
0588             else
0589                *Y = Yv / Yv_scale;
0590          }
0591
0592          return 0;
0593       }
0594
0595    } // namespace detail
0596
0597 }} // namespaces
0598
0599 #endif // BOOST_MATH_BESSEL_JY_HPP