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/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 dianostic 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 { initboost::math::detail::bessel_y1_initializer::init46 init() 47 { 48 do_init(); 49 } do_initboost::math::detail::bessel_y1_initializer::init50 static void do_init() 51 { 52 bessel_y1(T(1), Policy()); 53 } force_instantiateboost::math::detail::bessel_y1_initializer::init54 void force_instantiate()const{} 55 }; 56 static const init initializer; force_instantiateboost::math::detail::bessel_y1_initializer57 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>("bost::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 213