ceph/src/boost/boost/math/special_functions/detail/bessel_y1.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_Y1_HPP
   7 #define BOOST_MATH_BESSEL_Y1_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #pragma warning(push)
  12 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
  13 #endif
  14
  15 #include <boost/math/special_functions/detail/bessel_j1.hpp>
  16 #include <boost/math/constants/constants.hpp>
  17 #include <boost/math/tools/rational.hpp>
  18 #include <boost/math/tools/big_constant.hpp>
  19 #include <boost/math/policies/error_handling.hpp>
  20 #include <boost/math/tools/assert.hpp>
  21
  22 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
  23 //
  24 // This is the only way we can avoid
  25 // warning: non-standard suffix on floating constant [-Wpedantic]
  26 // when building with -Wall -pedantic.  Neither __extension__
  27 // nor #pragma diagnostic ignored work :(
  28 //
  29 #pragma GCC system_header
  30 #endif
  31
  32 // Bessel function of the second kind of order one
  33 // x <= 8, minimax rational approximations on root-bracketing intervals
  34 // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968
  35
  36 namespace boost { namespace math { namespace detail{
  37
  38 template <typename T, typename Policy>
  39 T bessel_y1(T x, const Policy&);
  40
  41 template <class T, class Policy>
  42 struct bessel_y1_initializer
  43 {
  44    struct init
  45    {
  46       init()
  47       {
  48          do_init();
  49       }
  50       static void do_init()
  51       {
  52          bessel_y1(T(1), Policy());
  53       }
  54       void force_instantiate()const{}
  55    };
  56    static const init initializer;
  57    static void force_instantiate()
  58    {
  59       initializer.force_instantiate();
  60    }
  61 };
  62
  63 template <class T, class Policy>
  64 const typename bessel_y1_initializer<T, Policy>::init bessel_y1_initializer<T, Policy>::initializer;
  65
  66 template <typename T, typename Policy>
  67 T bessel_y1(T x, const Policy& pol)
  68 {
  69     bessel_y1_initializer<T, Policy>::force_instantiate();
  70
  71     static const T P1[] = {
  72          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),
  73          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),
  74         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7595974497819597599e+11)),
  75          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2144548214502560419e+09)),
  76         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9157479997408395984e+07)),
  77          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2157953222280260820e+05)),
  78         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1714424660046133456e+02)),
  79     };
  80     static const T Q1[] = {
  81          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0737873921079286084e+14)),
  82          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1272286200406461981e+12)),
  83          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7800352738690585613e+10)),
  84          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2250435122182963220e+08)),
  85          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8136470753052572164e+05)),
  86          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.2079908168393867438e+02)),
  87          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
  88     };
  89     static const T P2[] = {
  90          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1514276357909013326e+19)),
  91         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.6808094574724204577e+18)),
  92         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3638408497043134724e+16)),
  93          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0686275289804744814e+15)),
  94         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9530713129741981618e+13)),
  95          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7453673962438488783e+11)),
  96         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1957961912070617006e+09)),
  97          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9153806858264202986e+06)),
  98         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2337180442012953128e+03)),
  99     };
 100     static const T Q2[] = {
 101          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3321844313316185697e+20)),
 102          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.6968198822857178911e+18)),
 103          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0837179548112881950e+16)),
 104          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1187010065856971027e+14)),
 105          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0221766852960403645e+11)),
 106          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.3550318087088919566e+08)),
 107          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0453748201934079734e+06)),
 108          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2855164849321609336e+03)),
 109          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
 110     };
 111     static const T PC[] = {
 112         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
 113         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
 114         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
 115         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
 116         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
 117         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
 118          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
 119     };
 120     static const T QC[] = {
 121         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
 122         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
 123         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
 124         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
 125         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
 126         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
 127          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
 128     };
 129     static const T PS[] = {
 130          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
 131          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
 132          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
 133          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
 134          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
 135          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
 136          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
 137     };
 138     static const T QS[] = {
 139          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
 140          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
 141          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
 142          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
 143          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
 144          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
 145          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
 146     };
 147     static const T x1  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1971413260310170351e+00)),
 148                    x2  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4296810407941351328e+00)),
 149                    x11 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.620e+02)),
 150                    x12 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8288260310170351490e-03)),
 151                    x21 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3900e+03)),
 152                    x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.4592058648672279948e-06))
 153     ;
 154     T value, factor, r, rc, rs;
 155
 156     BOOST_MATH_STD_USING
 157     using namespace boost::math::tools;
 158     using namespace boost::math::constants;
 159
 160     if (x <= 0)
 161     {
 162        return policies::raise_domain_error<T>("boost::math::bessel_y1<%1%>(%1%,%1%)",
 163             "Got x == %1%, but x must be > 0, complex result not supported.", x, pol);
 164     }
 165     if (x <= 4)                       // x in (0, 4]
 166     {
 167         T y = x * x;
 168         T z = 2 * log(x/x1) * bessel_j1(x) / pi<T>();
 169         r = evaluate_rational(P1, Q1, y);
 170         factor = (x + x1) * ((x - x11/256) - x12) / x;
 171         value = z + factor * r;
 172     }
 173     else if (x <= 8)                  // x in (4, 8]
 174     {
 175         T y = x * x;
 176         T z = 2 * log(x/x2) * bessel_j1(x) / pi<T>();
 177         r = evaluate_rational(P2, Q2, y);
 178         factor = (x + x2) * ((x - x21/256) - x22) / x;
 179         value = z + factor * r;
 180     }
 181     else                                // x in (8, \infty)
 182     {
 183         T y = 8 / x;
 184         T y2 = y * y;
 185         rc = evaluate_rational(PC, QC, y2);
 186         rs = evaluate_rational(PS, QS, y2);
 187         factor = 1 / (sqrt(x) * root_pi<T>());
 188         //
 189         // This code is really just:
 190         //
 191         // T z = x - 0.75f * pi<T>();
 192         // value = factor * (rc * sin(z) + y * rs * cos(z));
 193         //
 194         // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4
 195         // which then cancel out with corresponding terms in "factor".
 196         //
 197         T sx = sin(x);
 198         T cx = cos(x);
 199         value = factor * (y * rs * (sx - cx) - rc * (sx + cx));
 200     }
 201
 202     return value;
 203 }
 204
 205 }}} // namespaces
 206
 207 #ifdef _MSC_VER
 208 #pragma warning(pop)
 209 #endif
 210
 211 #endif // BOOST_MATH_BESSEL_Y1_HPP
 212