ceph/src/boost/libs/math/include/boost/math/special_functions/detail/bessel_jy.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_BESSEL_JY_HPP
   7 #define BOOST_MATH_BESSEL_JY_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/config.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15 #include <boost/math/special_functions/sign.hpp>
  16 #include <boost/math/special_functions/hypot.hpp>
  17 #include <boost/math/special_functions/sin_pi.hpp>
  18 #include <boost/math/special_functions/cos_pi.hpp>
  19 #include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
  20 #include <boost/math/special_functions/detail/bessel_jy_series.hpp>
  21 #include <boost/math/constants/constants.hpp>
  22 #include <boost/math/policies/error_handling.hpp>
  23 #include <boost/mpl/if.hpp>
  24 #include <boost/type_traits/is_floating_point.hpp>
  25 #include <complex>
  26
  27 // Bessel functions of the first and second kind of fractional order
  28
  29 namespace boost { namespace math {
  30
  31    namespace detail {
  32
  33       //
  34       // Simultaneous calculation of A&S 9.2.9 and 9.2.10
  35       // for use in A&S 9.2.5 and 9.2.6.
  36       // This series is quick to evaluate, but divergent unless
  37       // x is very large, in fact it's pretty hard to figure out
  38       // with any degree of precision when this series actually
  39       // *will* converge!!  Consequently, we may just have to
  40       // try it and see...
  41       //
  42       template <class T, class Policy>
  43       bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )
  44       {
  45          BOOST_MATH_STD_USING
  46             T tolerance = 2 * policies::get_epsilon<T, Policy>();
  47          *p = 1;
  48          *q = 0;
  49          T k = 1;
  50          T z8 = 8 * x;
  51          T sq = 1;
  52          T mu = 4 * v * v;
  53          T term = 1;
  54          bool ok = true;
  55          do
  56          {
  57             term *= (mu - sq * sq) / (k * z8);
  58             *q += term;
  59             k += 1;
  60             sq += 2;
  61             T mult = (sq * sq - mu) / (k * z8);
  62             ok = fabs(mult) < 0.5f;
  63             term *= mult;
  64             *p += term;
  65             k += 1;
  66             sq += 2;
  67          }
  68          while((fabs(term) > tolerance * *p) && ok);
  69          return ok;
  70       }
  71
  72       // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
  73       // Temme, Journal of Computational Physics, vol 21, 343 (1976)
  74       template <typename T, typename Policy>
  75       int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
  76       {
  77          T g, h, p, q, f, coef, sum, sum1, tolerance;
  78          T a, d, e, sigma;
  79          unsigned long k;
  80
  81          BOOST_MATH_STD_USING
  82             using namespace boost::math::tools;
  83          using namespace boost::math::constants;
  84
  85          BOOST_ASSERT(fabs(v) <= 0.5f);  // precondition for using this routine
  86
  87          T gp = boost::math::tgamma1pm1(v, pol);
  88          T gm = boost::math::tgamma1pm1(-v, pol);
  89          T spv = boost::math::sin_pi(v, pol);
  90          T spv2 = boost::math::sin_pi(v/2, pol);
  91          T xp = pow(x/2, v);
  92
  93          a = log(x / 2);
  94          sigma = -a * v;
  95          d = abs(sigma) < tools::epsilon<T>() ?
  96             T(1) : sinh(sigma) / sigma;
  97          e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
  98             : T(2 * spv2 * spv2 / v);
  99
 100          T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
 101          T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
 102          T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
 103          f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
 104
 105          p = vspv / (xp * (1 + gm));
 106          q = vspv * xp / (1 + gp);
 107
 108          g = f + e * q;
 109          h = p;
 110          coef = 1;
 111          sum = coef * g;
 112          sum1 = coef * h;
 113
 114          T v2 = v * v;
 115          T coef_mult = -x * x / 4;
 116
 117          // series summation
 118          tolerance = policies::get_epsilon<T, Policy>();
 119          for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
 120          {
 121             f = (k * f + p + q) / (k*k - v2);
 122             p /= k - v;
 123             q /= k + v;
 124             g = f + e * q;
 125             h = p - k * g;
 126             coef *= coef_mult / k;
 127             sum += coef * g;
 128             sum1 += coef * h;
 129             if (abs(coef * g) < abs(sum) * tolerance)
 130             {
 131                break;
 132             }
 133          }
 134          policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
 135          *Y = -sum;
 136          *Y1 = -2 * sum1 / x;
 137
 138          return 0;
 139       }
 140
 141       // Evaluate continued fraction fv = J_(v+1) / J_v, see
 142       // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
 143       template <typename T, typename Policy>
 144       int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
 145       {
 146          T C, D, f, a, b, delta, tiny, tolerance;
 147          unsigned long k;
 148          int s = 1;
 149
 150          BOOST_MATH_STD_USING
 151
 152             // |x| <= |v|, CF1_jy converges rapidly
 153             // |x| > |v|, CF1_jy needs O(|x|) iterations to converge
 154
 155             // modified Lentz's method, see
 156             // Lentz, Applied Optics, vol 15, 668 (1976)
 157             tolerance = 2 * policies::get_epsilon<T, Policy>();;
 158          tiny = sqrt(tools::min_value<T>());
 159          C = f = tiny;                           // b0 = 0, replace with tiny
 160          D = 0;
 161          for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
 162          {
 163             a = -1;
 164             b = 2 * (v + k) / x;
 165             C = b + a / C;
 166             D = b + a * D;
 167             if (C == 0) { C = tiny; }
 168             if (D == 0) { D = tiny; }
 169             D = 1 / D;
 170             delta = C * D;
 171             f *= delta;
 172             if (D < 0) { s = -s; }
 173             if (abs(delta - 1) < tolerance)
 174             { break; }
 175          }
 176          policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
 177          *fv = -f;
 178          *sign = s;                              // sign of denominator
 179
 180          return 0;
 181       }
 182       //
 183       // This algorithm was originally written by Xiaogang Zhang
 184       // using std::complex to perform the complex arithmetic.
 185       // However, that turns out to 10x or more slower than using
 186       // all real-valued arithmetic, so it's been rewritten using
 187       // real values only.
 188       //
 189       template <typename T, typename Policy>
 190       int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
 191       {
 192          BOOST_MATH_STD_USING
 193
 194             T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
 195          T tiny;
 196          unsigned long k;
 197
 198          // |x| >= |v|, CF2_jy converges rapidly
 199          // |x| -> 0, CF2_jy fails to converge
 200          BOOST_ASSERT(fabs(x) > 1);
 201
 202          // modified Lentz's method, complex numbers involved, see
 203          // Lentz, Applied Optics, vol 15, 668 (1976)
 204          T tolerance = 2 * policies::get_epsilon<T, Policy>();
 205          tiny = sqrt(tools::min_value<T>());
 206          Cr = fr = -0.5f / x;
 207          Ci = fi = 1;
 208          //Dr = Di = 0;
 209          T v2 = v * v;
 210          a = (0.25f - v2) / x; // Note complex this one time only!
 211          br = 2 * x;
 212          bi = 2;
 213          temp = Cr * Cr + 1;
 214          Ci = bi + a * Cr / temp;
 215          Cr = br + a / temp;
 216          Dr = br;
 217          Di = bi;
 218          if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
 219          if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
 220          temp = Dr * Dr + Di * Di;
 221          Dr = Dr / temp;
 222          Di = -Di / temp;
 223          delta_r = Cr * Dr - Ci * Di;
 224          delta_i = Ci * Dr + Cr * Di;
 225          temp = fr;
 226          fr = temp * delta_r - fi * delta_i;
 227          fi = temp * delta_i + fi * delta_r;
 228          for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
 229          {
 230             a = k - 0.5f;
 231             a *= a;
 232             a -= v2;
 233             bi += 2;
 234             temp = Cr * Cr + Ci * Ci;
 235             Cr = br + a * Cr / temp;
 236             Ci = bi - a * Ci / temp;
 237             Dr = br + a * Dr;
 238             Di = bi + a * Di;
 239             if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
 240             if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
 241             temp = Dr * Dr + Di * Di;
 242             Dr = Dr / temp;
 243             Di = -Di / temp;
 244             delta_r = Cr * Dr - Ci * Di;
 245             delta_i = Ci * Dr + Cr * Di;
 246             temp = fr;
 247             fr = temp * delta_r - fi * delta_i;
 248             fi = temp * delta_i + fi * delta_r;
 249             if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
 250                break;
 251          }
 252          policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
 253          *p = fr;
 254          *q = fi;
 255
 256          return 0;
 257       }
 258
 259       static const int need_j = 1;
 260       static const int need_y = 2;
 261
 262       // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
 263       // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
 264       template <typename T, typename Policy>
 265       int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
 266       {
 267          BOOST_ASSERT(x >= 0);
 268
 269          T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
 270          T W, p, q, gamma, current, prev, next;
 271          bool reflect = false;
 272          unsigned n, k;
 273          int s;
 274          int org_kind = kind;
 275          T cp = 0;
 276          T sp = 0;
 277
 278          static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
 279
 280          BOOST_MATH_STD_USING
 281             using namespace boost::math::tools;
 282          using namespace boost::math::constants;
 283
 284          if (v < 0)
 285          {
 286             reflect = true;
 287             v = -v;                             // v is non-negative from here
 288          }
 289          if (v > static_cast<T>((std::numeric_limits<int>::max)()))
 290          {
 291             *J = *Y = policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol);
 292             return 1;
 293          }
 294          n = iround(v, pol);
 295          u = v - n;                              // -1/2 <= u < 1/2
 296
 297          if(reflect)
 298          {
 299             T z = (u + n % 2);
 300             cp = boost::math::cos_pi(z, pol);
 301             sp = boost::math::sin_pi(z, pol);
 302             if(u != 0)
 303                kind = need_j|need_y;               // need both for reflection formula
 304          }
 305
 306          if(x == 0)
 307          {
 308             if(v == 0)
 309                *J = 1;
 310             else if((u == 0) || !reflect)
 311                *J = 0;
 312             else if(kind & need_j)
 313                *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity
 314             else
 315                *J = std::numeric_limits<T>::quiet_NaN();  // any value will do, not using J.
 316
 317             if((kind & need_y) == 0)
 318                *Y = std::numeric_limits<T>::quiet_NaN();  // any value will do, not using Y.
 319             else if(v == 0)
 320                *Y = -policies::raise_overflow_error<T>(function, 0, pol);
 321             else
 322                *Y = policies::raise_domain_error<T>(function, "Value of Bessel Y_v(x) is complex-infinity at %1%", x, pol); // complex infinity
 323             return 1;
 324          }
 325
 326          // x is positive until reflection
 327          W = T(2) / (x * pi<T>());               // Wronskian
 328          T Yv_scale = 1;
 329          if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5)))
 330          {
 331             //
 332             // This series will actually converge rapidly for all small
 333             // x - say up to x < 20 - but the first few terms are large
 334             // and divergent which leads to large errors :-(
 335             //
 336             Jv = bessel_j_small_z_series(v, x, pol);
 337             Yv = std::numeric_limits<T>::quiet_NaN();
 338          }
 339          else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
 340          {
 341             // Evaluate using series representations.
 342             // This is particularly important for x << v as in this
 343             // area temme_jy may be slow to converge, if it converges at all.
 344             // Requires x is not an integer.
 345             if(kind&need_j)
 346                Jv = bessel_j_small_z_series(v, x, pol);
 347             else
 348                Jv = std::numeric_limits<T>::quiet_NaN();
 349             if((org_kind&need_y && (!reflect || (cp != 0)))
 350                || (org_kind & need_j && (reflect && (sp != 0))))
 351             {
 352                // Only calculate if we need it, and if the reflection formula will actually use it:
 353                Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
 354             }
 355             else
 356                Yv = std::numeric_limits<T>::quiet_NaN();
 357          }
 358          else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
 359          {
 360             // Truncated series evaluation for small x and v an integer,
 361             // much quicker in this area than temme_jy below.
 362             if(kind&need_j)
 363                Jv = bessel_j_small_z_series(v, x, pol);
 364             else
 365                Jv = std::numeric_limits<T>::quiet_NaN();
 366             if((org_kind&need_y && (!reflect || (cp != 0)))
 367                || (org_kind & need_j && (reflect && (sp != 0))))
 368             {
 369                // Only calculate if we need it, and if the reflection formula will actually use it:
 370                Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
 371             }
 372             else
 373                Yv = std::numeric_limits<T>::quiet_NaN();
 374          }
 375          else if(asymptotic_bessel_large_x_limit(v, x))
 376          {
 377             if(kind&need_y)
 378             {
 379                Yv = asymptotic_bessel_y_large_x_2(v, x);
 380             }
 381             else
 382                Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
 383             if(kind&need_j)
 384             {
 385                Jv = asymptotic_bessel_j_large_x_2(v, x);
 386             }
 387             else
 388                Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
 389          }
 390          else if((x > 8) && hankel_PQ(v, x, &p, &q, pol))
 391          {
 392             //
 393             // Hankel approximation: note that this method works best when x
 394             // is large, but in that case we end up calculating sines and cosines
 395             // of large values, with horrendous resulting accuracy.  It is fast though
 396             // when it works....
 397             //
 398             // Normally we calculate sin/cos(chi) where:
 399             //
 400             // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
 401             //
 402             // But this introduces large errors, so use sin/cos addition formulae to
 403             // improve accuracy:
 404             //
 405             T mod_v = fmod(T(v / 2 + 0.25f), T(2));
 406             T sx = sin(x);
 407             T cx = cos(x);
 408             T sv = sin_pi(mod_v);
 409             T cv = cos_pi(mod_v);
 410
 411             T sc = sx * cv - sv * cx; // == sin(chi);
 412             T cc = cx * cv + sx * sv; // == cos(chi);
 413             T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x));
 414             Yv = chi * (p * sc + q * cc);
 415             Jv = chi * (p * cc - q * sc);
 416          }
 417          else if (x <= 2)                           // x in (0, 2]
 418          {
 419             if(temme_jy(u, x, &Yu, &Yu1, pol))             // Temme series
 420             {
 421                // domain error:
 422                *J = *Y = Yu;
 423                return 1;
 424             }
 425             prev = Yu;
 426             current = Yu1;
 427             T scale = 1;
 428             policies::check_series_iterations<T>(function, n, pol);
 429             for (k = 1; k <= n; k++)            // forward recurrence for Y
 430             {
 431                T fact = 2 * (u + k) / x;
 432                if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
 433                {
 434                   scale /= current;
 435                   prev /= current;
 436                   current = 1;
 437                }
 438                next = fact * current - prev;
 439                prev = current;
 440                current = next;
 441             }
 442             Yv = prev;
 443             Yv1 = current;
 444             if(kind&need_j)
 445             {
 446                CF1_jy(v, x, &fv, &s, pol);                 // continued fraction CF1_jy
 447                Jv = scale * W / (Yv * fv - Yv1);           // Wronskian relation
 448             }
 449             else
 450                Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
 451             Yv_scale = scale;
 452          }
 453          else                                    // x in (2, \infty)
 454          {
 455             // Get Y(u, x):
 456
 457             T ratio;
 458             CF1_jy(v, x, &fv, &s, pol);
 459             // tiny initial value to prevent overflow
 460             T init = sqrt(tools::min_value<T>());
 461             BOOST_MATH_INSTRUMENT_VARIABLE(init);
 462             prev = fv * s * init;
 463             current = s * init;
 464             if(v < max_factorial<T>::value)
 465             {
 466                policies::check_series_iterations<T>(function, n, pol);
 467                for (k = n; k > 0; k--)             // backward recurrence for J
 468                {
 469                   next = 2 * (u + k) * current / x - prev;
 470                   prev = current;
 471                   current = next;
 472                }
 473                ratio = (s * init) / current;     // scaling ratio
 474                // can also call CF1_jy() to get fu, not much difference in precision
 475                fu = prev / current;
 476             }
 477             else
 478             {
 479                //
 480                // When v is large we may get overflow in this calculation
 481                // leading to NaN's and other nasty surprises:
 482                //
 483                policies::check_series_iterations<T>(function, n, pol);
 484                bool over = false;
 485                for (k = n; k > 0; k--)             // backward recurrence for J
 486                {
 487                   T t = 2 * (u + k) / x;
 488                   if((t > 1) && (tools::max_value<T>() / t < current))
 489                   {
 490                      over = true;
 491                      break;
 492                   }
 493                   next = t * current - prev;
 494                   prev = current;
 495                   current = next;
 496                }
 497                if(!over)
 498                {
 499                   ratio = (s * init) / current;     // scaling ratio
 500                   // can also call CF1_jy() to get fu, not much difference in precision
 501                   fu = prev / current;
 502                }
 503                else
 504                {
 505                   ratio = 0;
 506                   fu = 1;
 507                }
 508             }
 509             CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
 510             T t = u / x - fu;                   // t = J'/J
 511             gamma = (p - t) / q;
 512             //
 513             // We can't allow gamma to cancel out to zero competely as it messes up
 514             // the subsequent logic.  So pretend that one bit didn't cancel out
 515             // and set to a suitably small value.  The only test case we've been able to
 516             // find for this, is when v = 8.5 and x = 4*PI.
 517             //
 518             if(gamma == 0)
 519             {
 520                gamma = u * tools::epsilon<T>() / x;
 521             }
 522             BOOST_MATH_INSTRUMENT_VARIABLE(current);
 523             BOOST_MATH_INSTRUMENT_VARIABLE(W);
 524             BOOST_MATH_INSTRUMENT_VARIABLE(q);
 525             BOOST_MATH_INSTRUMENT_VARIABLE(gamma);
 526             BOOST_MATH_INSTRUMENT_VARIABLE(p);
 527             BOOST_MATH_INSTRUMENT_VARIABLE(t);
 528             Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
 529             BOOST_MATH_INSTRUMENT_VARIABLE(Ju);
 530
 531             Jv = Ju * ratio;                    // normalization
 532
 533             Yu = gamma * Ju;
 534             Yu1 = Yu * (u/x - p - q/gamma);
 535
 536             if(kind&need_y)
 537             {
 538                // compute Y:
 539                prev = Yu;
 540                current = Yu1;
 541                policies::check_series_iterations<T>(function, n, pol);
 542                for (k = 1; k <= n; k++)            // forward recurrence for Y
 543                {
 544                   T fact = 2 * (u + k) / x;
 545                   if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
 546                   {
 547                      prev /= current;
 548                      Yv_scale /= current;
 549                      current = 1;
 550                   }
 551                   next = fact * current - prev;
 552                   prev = current;
 553                   current = next;
 554                }
 555                Yv = prev;
 556             }
 557             else
 558                Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
 559          }
 560
 561          if (reflect)
 562          {
 563             if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
 564                *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
 565             else
 566                *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale));     // reflection formula
 567             if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
 568                *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
 569             else
 570                *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
 571          }
 572          else
 573          {
 574             *J = Jv;
 575             if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
 576                *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
 577             else
 578                *Y = Yv / Yv_scale;
 579          }
 580
 581          return 0;
 582       }
 583
 584    } // namespace detail
 585
 586 }} // namespaces
 587
 588 #endif // BOOST_MATH_BESSEL_JY_HPP
 589