special_functions/detail/bessel_y1.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_Y1_HPP
0007 #define BOOST_MATH_BESSEL_Y1_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #pragma warning(push)
0012 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
0013 #endif
0014
0015 #include <boost/math/tools/config.hpp>
0016 #include <boost/math/special_functions/detail/bessel_j1.hpp>
0017 #include <boost/math/constants/constants.hpp>
0018 #include <boost/math/tools/rational.hpp>
0019 #include <boost/math/tools/big_constant.hpp>
0020 #include <boost/math/policies/error_handling.hpp>
0021 #include <boost/math/tools/assert.hpp>
0022
0023 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0024 //
0025 // This is the only way we can avoid
0026 // warning: non-standard suffix on floating constant [-Wpedantic]
0027 // when building with -Wall -pedantic.  Neither __extension__
0028 // nor #pragma diagnostic ignored work :(
0029 //
0030 #pragma GCC system_header
0031 #endif
0032
0033 // Bessel function of the second kind of order one
0034 // x <= 8, minimax rational approximations on root-bracketing intervals
0035 // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968
0036
0037 namespace boost { namespace math { namespace detail{
0038
0039 template <typename T, typename Policy>
0040 BOOST_MATH_GPU_ENABLED T bessel_y1(T x, const Policy&);
0041
0042 template <typename T, typename Policy>
0043 BOOST_MATH_GPU_ENABLED T bessel_y1(T x, const Policy&)
0044 {
0045     BOOST_MATH_STATIC const T P1[] = {
0046          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),
0047          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),
0048         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7595974497819597599e+11)),
0049          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2144548214502560419e+09)),
0050         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9157479997408395984e+07)),
0051          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2157953222280260820e+05)),
0052         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1714424660046133456e+02)),
0053     };
0054     BOOST_MATH_STATIC const T Q1[] = {
0055          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0737873921079286084e+14)),
0056          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1272286200406461981e+12)),
0057          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7800352738690585613e+10)),
0058          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2250435122182963220e+08)),
0059          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8136470753052572164e+05)),
0060          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.2079908168393867438e+02)),
0061          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
0062     };
0063     BOOST_MATH_STATIC const T P2[] = {
0064          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1514276357909013326e+19)),
0065         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.6808094574724204577e+18)),
0066         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3638408497043134724e+16)),
0067          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0686275289804744814e+15)),
0068         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9530713129741981618e+13)),
0069          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7453673962438488783e+11)),
0070         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1957961912070617006e+09)),
0071          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9153806858264202986e+06)),
0072         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2337180442012953128e+03)),
0073     };
0074     BOOST_MATH_STATIC const T Q2[] = {
0075          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3321844313316185697e+20)),
0076          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.6968198822857178911e+18)),
0077          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0837179548112881950e+16)),
0078          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1187010065856971027e+14)),
0079          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0221766852960403645e+11)),
0080          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.3550318087088919566e+08)),
0081          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0453748201934079734e+06)),
0082          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2855164849321609336e+03)),
0083          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
0084     };
0085     BOOST_MATH_STATIC const T PC[] = {
0086         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
0087         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
0088         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
0089         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
0090         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
0091         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
0092          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
0093     };
0094     BOOST_MATH_STATIC const T QC[] = {
0095         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
0096         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
0097         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
0098         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
0099         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
0100         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
0101          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
0102     };
0103     BOOST_MATH_STATIC const T PS[] = {
0104          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
0105          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
0106          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
0107          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
0108          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
0109          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
0110          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
0111     };
0112     BOOST_MATH_STATIC const T QS[] = {
0113          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
0114          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
0115          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
0116          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
0117          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
0118          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
0119          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
0120     };
0121     BOOST_MATH_STATIC const T x1  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1971413260310170351e+00)),
0122                    x2  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4296810407941351328e+00)),
0123                    x11 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.620e+02)),
0124                    x12 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8288260310170351490e-03)),
0125                    x21 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3900e+03)),
0126                    x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.4592058648672279948e-06))
0127     ;
0128     T value, factor, r, rc, rs;
0129
0130     BOOST_MATH_STD_USING
0131     using namespace boost::math::tools;
0132     using namespace boost::math::constants;
0133
0134     BOOST_MATH_ASSERT(x > 0);
0135
0136     if (x <= 4)                       // x in (0, 4]
0137     {
0138         T y = x * x;
0139         T z = 2 * log(x/x1) * bessel_j1(x) / pi<T>();
0140         r = evaluate_rational(P1, Q1, y);
0141         factor = (x + x1) * ((x - x11/256) - x12) / x;
0142         value = z + factor * r;
0143     }
0144     else if (x <= 8)                  // x in (4, 8]
0145     {
0146         T y = x * x;
0147         T z = 2 * log(x/x2) * bessel_j1(x) / pi<T>();
0148         r = evaluate_rational(P2, Q2, y);
0149         factor = (x + x2) * ((x - x21/256) - x22) / x;
0150         value = z + factor * r;
0151     }
0152     else                                // x in (8, \infty)
0153     {
0154         T y = 8 / x;
0155         T y2 = y * y;
0156         rc = evaluate_rational(PC, QC, y2);
0157         rs = evaluate_rational(PS, QS, y2);
0158         factor = 1 / (sqrt(x) * root_pi<T>());
0159         //
0160         // This code is really just:
0161         //
0162         // T z = x - 0.75f * pi<T>();
0163         // value = factor * (rc * sin(z) + y * rs * cos(z));
0164         //
0165         // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4
0166         // which then cancel out with corresponding terms in "factor".
0167         //
0168         T sx = sin(x);
0169         T cx = cos(x);
0170         value = factor * (y * rs * (sx - cx) - rc * (sx + cx));
0171     }
0172
0173     return value;
0174 }
0175
0176 }}} // namespaces
0177
0178 #ifdef _MSC_VER
0179 #pragma warning(pop)
0180 #endif
0181
0182 #endif // BOOST_MATH_BESSEL_Y1_HPP
0183